Effects of lateral processes on the seasonal water stratification of the Gulf of Finland: 3-D NEMO-based model study

This paper aims to fill the gaps in knowledge of processes affecting the seasonal water stratification in the Gulf of Finland (GOF). We used a state-of-the-art modelling framework NEMO (Nucleus for European Modelling of the Ocean) designed for oceanographic research, operational oceanography, seasonal forecasting, and climate studies to build an eddy-resolving model of the GOF. To evaluate the model skill and performance, two different solutions were obtained on 0.5 km eddy-resolving and commonly used 2 km grids for a 1-year simulation. We also explore the efficacy of non-hydrostatic effect (convection) parameterizations available in NEMO for coastal application. It is found that the solutions resolving submesoscales have a more complex mixed layer structure in the regions of the GOF directly affected by the upwelling/downwelling and intrusions from the open Baltic Sea. Presented model estimations of the upper mixed layer depth are in good agreement with in situ CTD (BED) data. A number of model sensitivity tests to the vertical mixing parameterization confirm the model’s robustness. Further progress in the submesoscale process simulation and understanding is apparently not connected mainly with the finer resolution of the grids, but with the use of non-hydrostatic models because of the failure of the hydrostatic approach at submesoscale.


Introduction
The Gulf of Finland (GOF) is a 400 km long and 48-135 km wide sub-basin of the Baltic Sea with a mean depth of 37 m and complex bathymetry (see Fig. 1).The large fresh-water input from the Neva River significantly affects the stratification and forms the strong salinity gradient from east to west and from north to south.Sea-surface salinity decreases from 5-6.5 ‰ in the western GOF to about 0-3 ‰ in the easternmost part of the Gulf, where the role of the Neva River is most pronounced (Alenius et al., 1998).In the western GOF, a quasi-permanent halocline is located at a depth of 60-80 m.Salinity in that area can reach values as high as 8-10 ‰ near the sea bed due to the advection of saltier water masses from the Baltic Proper.
The vertical stratification in the GOF as well as in the Baltic Sea is unusual (the thermocline and halocline are usually separated) with a pronounced and relatively stable halocline, whereas the temperature is largely controlled by the seasonal variability of the surface heat fluxes (see, e.g., Hankimo, 1964).During the summer season the water column in the deeper areas of the GOF consists of three layers -the upper mixed layer (UML), the cold intermediate layer, and a saltier and slightly warmer near-bottom layer (see Liblik and Lips, 2012), separated by two pycnoclines -the thermocline at the depths of 10-20 m and the permanent halocline at the depths of 60-70 m.A seasonal thermocline starts to develop in May.The surface mixed layer reaches a maximum depth of 15-20 m by midsummer and an erosion of the thermocline starts in late August due to wind mixing and ther- mal convection.The bottom salinity also shows significant spatio-temporal variability due to irregular saline water intrusions from the Baltic Proper, as well as from changes in river runoff and the precipitation-evaporation balance.There is no permanent halocline in the eastern GOF, where salinity increases approximately linearly with depth (Nekrasov and Lebedeva, 2002;Alenius et al., 2003).
The simulations of the vertical stratification using threedimensional (3-D) numerical models are not so reliable yet (Myrberg et al., 2010).This study shows that the most advanced 3-D circulation models are able to simulate the major features of the hydro-physical fields of the GOF.For example, generally the hind-cast temperatures differ from observations by less than 1-2 • C and the mean error in salinity is less than 1 ‰.Most of the remaining difficulties are connected with problems in adequately representing the dynamics of the mixed layer.The loss of accuracy is most notable in the simulation of the depth and the sharpness of the corresponding thermo-and haloclines.Despite the application of sophisticated turbulent closure schemes and different schemes for vertical mixing, none of the models, analysed in Myrberg et al. (2010), were able to accurately simulate the vertical profiles of temperature and salinity.The latest experiments with turbulence parameterizations of the 3-D hydrodynamic model COHERENS presented in Tuomi et al. (2013) show that the model still underestimates the thermocline depth.Furthermore, the sensitivity of the modelled thermocline depth to the accuracy of the meteorological forcing was studied by increasing the forcing wind speed to better match the measured values of wind speed in the central GOF.The sensitivity test showed that an increase in the wind speed only slightly improved the performance of the turbulence parameterizations in modelling the thermocline depth.
However, a number of studies have reported important effects of the vertical thermohaline structure on the characteristics and processes in the marine ecosystems of the GOF, such as phytoplankton species composition (Rantajarvi et al., 1998) and sub-surface maxima of phytoplankton biomass (Lips et al., 2010), cyanobacteria blooms (Lips et al., 2008), distribution of pelagic fish (Stepputtis et al., 2011), macrozoobenthos abundance (Laine et al., 2007), and oxygen concentrations in the near-bottom layer (Maximov, 2006).
In summary, prediction of the thermohaline structure is a complex problem for the GOF.The spatial variability of the thermohaline structure encompasses a wide range of physical processes at different scales, some of which are still poorly understood (Soomere et al., 2008(Soomere et al., , 2009)).For example, we hypothesize that the local stratification depends very strongly on the across-GOF movements of water masses and that submesoscale eddies generated by baroclinic instability of fronts in the upper layers of the sea not only play an important role in the heterogeneity of spatial distribution of parameters (temperature, nutrients, phytoplankton), but also can contribute to re-stratify the UML, as described in Gent and McWilliams (1990).
In the ocean, submesoscales are scales of motion equal to or less than the Rossby radius of deformation but large enough to be influenced by planetary rotation (Thomas et al., 2007).Recent studies have shown that increasing the horizontal resolution of the model up to 0.5 km (for the GOF Rossby radius aprox.2-4 km) enables models to resolve submesoscale eddies.As a result, surface currents and temperatures show highly detailed patterns that qualitatively match well with the expected features (Zhurbas et al., 2008;Sokolov, 2013).However, there has not yet been any consideration of the influence of eddy motions and across-Gulf movements of water masses on vertical re-stratification of the UML of the GOF.
The motivations behind this study are to provide an insight into the lateral advection processes in the GOF, with particular interest in estimating the contribution of lateral advection processes to the thermocline variations; to assess the impact of horizontal grid resolution on the representation of vertical stratification.

Approach
The traditional point of view is that the eddy diffusion dominates in the horizontal direction and in the vertical direction mixing due to eddies is limited, and small-scale processes such as turbulence provide the majority of mixing.Based on this idea, most commonly a 1-D approach is used to set up vertical mixing by tuning a turbulent scheme.The GOF, an enclosed basin with complex bathymetry and strong stratification mixed layer dynamics, can be strongly affected by lateral advective processes.To investigate this phenomenon, we present a state-of-the art 3-D model of the GOF with high vertical and two different horizontal resolutions.Shelf sea modelling is characterized by a demand for many different configurations to meet multiple science and user needs.NEMO (Nucleus for European Modelling of the Ocean) gives the capability to rapidly configure shelf sea models using appropriate high resolutions and parameterizations for the representation of coastal dynamics.

General model set-up
Our study is based on a 3-D thermo-hydrodynamic model build on the NEMO code initially designed for the open ocean and adopted by our team for the GOF (NEMO GOF).
The NEMO is a 3-D hydrostatic, baroclinic primitive equation model toolkit laid out horizontally on the Arakawa Cgrid (Madec et al., 1998;Madec, 2012).The NEMO is developing in a framework of a community of European institutes and benefits from the recent scientific and technical developments implemented in most ocean modelling platforms.The NEMO implementation for the GOF uses the TVD (Madec et al., 1998) advection scheme in the horizontal direction, the piecewise parabolic method (PPM) in the vertical direction (Liu and Holt, 2010), the non-linear variable volume (VVL) scheme for the free surface.In the horizontal plane, the model uses the standard Jacobean formulation for the pressure gradient, the viscosity and diffusivity formulation with a constant coefficient for momentum and tracer diffusion.The horizontal viscosity and diffusivity operators are rotated to be aligned with the density iso-surfaces to accurately reproduce density flows.
There are NEMO set-ups for the Baltic Sea that have been recently published by Hordoir et al. (2013Hordoir et al. ( , 2015)).The GOF set-up was developed in parallel to the Baltic Sea model and aimed to introduce resolution able to resolve the submesoscale processes in a horizontal direction and insure accurate representation of the vertical structure by increasing the vertical resolution to 1 m.General model set-up for the GOF shares most of the parameterization and schemes with Baltic Sea model.
In this paper, we used a gridded bathymetric data set with a resolution of 0.25 nm for the GOF (Andrejev, 2010).Choosing different grid resolutions of the model is formally equivalent to the choice of an appropriate averaging operator (lowpass filtering at the grid step) and an approach to estimate the contribution of smaller scales to the general motion.To assess the impact of submesoscale motion on the vertical stratification, two configurations of NEMO GOF were generated by utilizing different horizontal resolutions and the same vertical resolution of 1 m.Both configurations have 94 verti-cal levels, but 1 min zonal and 2 min meridional resolution (∼ 2 km) in a standard configuration and 0.25 min zonal and 0.5 min meridional resolution (∼ 0.5 km) in a finer resolution configuration.The parameters of configurations were kept as identical as possible.The main exception is the coefficients of horizontal diffusivity and viscosity, which were set to the minimum values guaranteeing numerical stability.
Numerical experiments were started from rest and initialized with temperature and salinity fields from the operational model of the Baltic Sea HIROMB (Funkquist, 2001).The computational domain covers the entire GOF with the open boundary set at 23 • 30 E longitude (see Fig. 1), boundary conditions also being taken from HIROMB.According to the inter-comparison of several model results for the GOF (Myrberg et al., 2010), HIROMB was rated as the best model for the western part of the GOF.The operational status of the model gave us additional benefits.The model was forced by the surface-forcing data set HIRLAM (http://hirlam.org)(using the CORE bulk-forcing algorithm) and climatic rivers runoff (Stalnacke et al., 1999).We used SMHI version of HI-ROMB with HIRLAM atmospheric fields included in output files as a part of a standard operational product of SMHI.Temporal resolution for the atmospheric forcing and boundary conditions is 1 h.

Parameterization of convective flows
One of the possible mechanisms by which the lateral motion affects the stratification is a shear-induced convection: situation in which heavy water may be advected on top of lighter water.This mechanism has been observed, e.g. in the bottom boundary layer of lakes (Lorke et al., 2005) and on the continental shelf (Rippeth et al., 2001).Evidently, the shearinduced convection can take place throughout the water column, for example, during upwelling.In nature, convective processes quickly re-establish the static stability of the water column (Umlauf, 2005).These processes have been removed from the model via the hydrostatic assumption so they must be parameterized.
Convective mixing can be parameterized in NEMO by (1) a computationally efficient solution "TKE (turbulent kinetic energy) scheme" in combination with convective adjustment procedures (a non-penetrative convective adjustment or an enhanced vertical diffusion) and (2) a physically more accurate GLS (generic length scale) scheme.
The TKE scheme is a turbulence closure scheme proposed by Bougeault and Lacarrére (1989) originally developed in a model for the atmospheric boundary layer.In the Mellor and Yamada (1974) hierarchy it is a 1.5-level closure and consists of a prognostic closure for the TKE and an algebraic formulation for the mixing length scale.The time evolution of TKE is the result of the production of TKE through vertical shear, its suppression through stratification, its vertical diffusion, and dissipation of the Kolmogorov (1942) type: where N is the local buoyancy frequency, l ε and l k are the dissipation and mixing length scales, u and v are the horizontal velocity components, k is the layer number, e 3 = 1 m is the vertical scale factor, P rt is the Prandtl number and K m and K ρ are the vertical eddy viscosity and diffusivity coefficients.The parameter C k is known as a stability function and is defined as a constant in the TKE scheme.The constants C k = 0.1 and C ε = 0.7 are specified to deal with vertical mixing at any depth (Gaspar et al., 1990).K e is the eddy diffusivity coefficient for the TKE.In NEMO K e = K m .For computational efficiency, the original formulation of the turbulent length scales proposed by Gaspar et al. (1990) has been simplified to the following first-order approximation This simplification valid in a stable stratified region with constant values of the buoyancy frequency has two major drawbacks: it makes no sense for locally unstable stratification and the computation no longer uses all the information contained in the vertical density profile.To overcome these drawbacks, NEMO TKE scheme implementation adds an extra assumption concerning the vertical gradient of the computed length scale.Therefore, the length scales are first evaluated as in Eq. ( 4) and then bounded such that In order to impose the constraint of Eq. ( 5), NEMO introduces two additional length scales: l up and l dwn .The length scales l up and l dwn are respectively the upward and downward distances to which a fluid parcel is able to travel from current z-level k, converting its TKE into the potential energy by working against the stratification, and they can be evaluated as where nk is the number of levels in vertical and l (k) is computed using Eq.(4); i.e. Finally, The GLS scheme is formally equivalent to the TKE scheme, except it uses (1) a prognostic equation for the generic length scale φ and (2) expressions for the complex stability functions instead constants.We used k − ε turbulent closure scheme (Rodi, 1987) with φ = C 3 0µ e 3/2 l −1 , where C 0µ is a constant depending on the choice of the stability function (Galperin et al., 1988;Kantha and Clayson, 1994).
This prognostic length scale is valid for convective situations and arbitrarily increases diffusivity to represent convection (Umlauf andBurchard, 2003, 2005): Here C 1 , C 2 , C 3 , σ φ are constants for the k − ε turbulent closure scheme.They are equal 1.44, 1.92, 1.0, 1.3 respectively.C µ and C µ are calculated from the stability function.
As known, the equation fails in stably stratified flows, and for this reason almost all authors apply a clipping of the length scale as an ad hoc remedy.With this clipping, the maximum permissible length scale is determined by A value of C lim = 0.53 is often used (Galperin et al., 1988).Umlauf and Burchard (2005) showed that the value of the clipping factor is of crucial importance for the entrainment depth predicted in stably stratified flows.Another value is 0.26, several authors have suggested limiting the dissipative length-scale in the presence of stable stratification even down to 0.07 (Holt and Umlauf, 2008).
In addition, convective mixing can be parameterized in NEMO by an enhancement to the eddy viscosity and diffusivity (ED), if for N 2 <0, K m and K ρ are locally set to the value of 100 m 2 s −1 .
We performed comparative tests of above-listed convection parameterizations to investigate their principal applicability for shear-induced convective situations.

Numerical experiments
The modelling period was chosen from 1 April to 31 August 2011 when pronounced thermocline occurs.The thermocline starts its formation in early May when the surface Ocean Sci., 12, 987-1001, 2016 www.ocean-sci.net/12/987/2016/heating and turbulent mixing are dominant processes.Note that year 2011 was characterized by strong upwelling events in the beginning and in the end of the modelling period.In Sect.3.1 the GLS, TKE, and ED mixing parameterizations are compared in a series of sensitivity experiments.The choice of closure scheme and the effects of a varying Galperin limit were investigated against MODIS sea-surface temperature (SST) to get the best reproduction of SST pattern.
In Sect.3.2 we present results of the model runs compared with available CTD (BED) data to study the performance of the chosen parameterizations to represent the UML evolution.Also, the ability of the model to correctly capture such features as fronts was tested against SST images for different resolutions in beginning of August 2011 when there were cloud-free images.

Sensitivity to vertical mixing parameterizations
In this section we study closure schemes and enhanced diffusion parameterization performance for convective situations caused by upwelling near the Estonian coast starting on 12 May.Figure 2 shows a cross section of the GOF for the density field (black isolines) overlaid by the vertical eddy diffusivity coefficient (colour filled).
Fragment A of Fig. 2 illustrates the mechanism instability formation.It is a hypothetical solution obtained with constant eddy diffusivity coefficients set to the minimum possible for this case, i.e. values of 10 −4 -10 −5 m 2 s −1 , and ED switched off.All south-north cross sections present the situation mainly formed by an upwelling event near the Estonian coast (left side of the cross section).Due to the presence of permanent density gradient from Estonian to the Finish coast and a strong offshore current caused by upwelling, dense waters originated from the Estonian side overlay fresher lighter water in the downwelling area near the Finish coast.
Fragment B illustrates the performance of the ED procedure setting the eddy viscosity and diffusivity coefficients equal to 100 m 2 s −1 in the areas of unstable stratification.According to this experiment, the maximum depth of convection penetration is equal to 10 m in the centre of GOF and reaches up to 25 m near the Finish coast.
Fragment C illustrates the performance of solution with the TKE closure scheme including previously described modifications introduced in NEMO.As seen, the solution demonstrates high values of eddy diffusion coefficients in the areas of unstable stratification.The depth of the mixed layer is not limited by the convection penetration depth (see Fig. 2b) and formed as a result of a joint action of current velocity shear, buoyancy, and TKE diffusion and dissipation (see Eq. 1).
Fragment D shows the combined effect of cases B and C. As seen from comparison of Fig. 2d and c, the solution with the modified TKE scheme captures most of the existing instabilities.ED (Fig. 2b) triggered only in some small areas in the centre of the mixed layer and did not affect the actual mixing depth.
Fragments E and F present the performance of the solution with the GLS closure scheme with Galperin limits of 0.53 and 0.26 respectively.A solution with GLS parameterization with switched-off length-scale limitation was also obtained but turned out to be practically equal to case E. UML depth in these solutions is comparable to that in the cases C and D, confirming the success of TKE modifications in NEMO.
The above tests confirm that both TKE and GLS closure schemes used in NEMO are able to catch the convection induced by upwelling.As seen in Fig. 2, an instability of the vertical column initiates dramatic increase in vertical diffusivity coefficients up to 0.04 m 2 s −1 TKE (Fig. 2c and d) or 0.036 m 2 s −1 GLS (Fig. 2e and f) from the background value set to 10 −6 m 2 s −1 .The TKE scheme forms a core with stronger mixing in the area of downwelling but at the same time the UML depth is comparable in both cases.Switched on ED does not modify the UML depth predicted by turbulent closure schemes.
Evaluation of the actual performance of presented alternative parameterizations of convective processes is a complex task requiring high spatial and temporal resolution of in situ data that is not available at the moment.The SST derived from the satellite thermal infrared imagery during cloud-free conditions provides significant information for monitoring of the relevant key ocean structures, such as fronts, eddies, and upwelling.At the same time, the SST fields can be used as an indicator of vertical mixing processes.SST fields can be considered as integral of subsurface dynamic, but, for example, we cannot directly estimate a depth of the thermocline from them.Alternatively the comparison of the modelled frontal structure at the sea surface and MODIS data during an upwelling event (lifting water from under the UML) could indicate how well the model reproduces stratification.As soon as we would get a realistic stratification, the surface pattern of simulated SST will also be in agreement with remotely observed SST.
Results of the comparison of modelled (various mixing parameterizations and resolutions) and MODIS-derived SST are presented at Fig. 3.The model shows that maximum upwelling development occurs on 14 May when the upwelling front reaches the centre of the GOF and is characterized by a maximum temperature difference across the front of up to 5 • C. Unfortunately, due to heavy cloudiness, the satellite images captured only the relaxation phase of the upwelling dated on 20 May.
As seen, the model performs better if the GLS scheme is used and the value of C lim is 0.53 (Galperin's value).The stronger length-scale limitation leads to underestimation of mixing and increased SST values compared to MODIS data.On the other hand, the solution obtained with TKE scheme underestimates mixing; nevertheless, it is not too far from the observations.The best performance takes place at the higher resolution and GLS scheme used when the solution is in a good agreement with the MODIS SST (Fig. 3b).Based on presented sensitivity tests, the GLS mixing scheme was chosen and the length scale limiting was fixed as C lim = 0.53.

General model performance
To evaluate the general model performance, we used in situ data for temperature and salinity obtained during the Russian State Hydrometeorological University expedition dated from 20 July 2011 to 5 August 2011.The comparison of model and data has been performed for the last decade in July just before the UML starts to degrade due to heating and wind conditions (Fig. 4).CTD data were grouped into three sets of profiles representing western (23.5-26 longitude, 10 profiles), central (26-28.2longitude, 12 profiles) and eastern (28.2-30 longitude, 12 profiles) parts of the GOF.According to the presented (Fig. 4) averaged CTD profiles (black curves), the UML is much deeper in the western part of the GOF and considerably shallower and sharper in the central and eastern parts.This UML behaviour, typical for the GOF, was captured quite well by all the model realizations (coloured curves).Standard deviation of CTD data given as error bars presents the variability range of in situ data.All presented solutions with different parameterizations are in good agreement with the data in terms of the UML depth, whereas the fine spatial resolution slightly better represents the nature in the western part of GOF.In the eastern part of the GOF strongly influenced by the Neva outflow the modelled thermocline is about 5 m deeper than observed.This is mainly due to prescribing climatic boundary conditions at the river mouth not allowing for the differences in individual years and complicated hydrodynamics of the estuary.
One more comparison between model and data is presented in Fig. 5 where the modelled SST for the two resolutions is given vs. MODIS SST on 2 August 2011.At this time it was possible to fix the upwelling again near the southern coast of the GOF.In the high-resolution model solution the temperature of cold water rising to the surface drops down to 6 • C, which is consistent with the satellite SST.In the case of coarse resolution the upwelling effect is less pronounced: the lowest temperature in the core region is about 10 • C. Solutions with both resolutions reproduce spatial patterns of upwelling.Although the coarse-resolution solution gives a more flattened upwelling front (shown by the isotherm of 19.5 • C), high-resolution solution is more rugged due to reproduced submesoscale features that corresponds well with observed SST.
Results of model comparison with SST and in situ data confirm the robustness of the developed model, which allows us to use it in a more detailed evaluation of the vertical structure formation mechanisms of the sea and its temporal evolution.

Results
During the upwelling/downwelling event in May modelled on both grids simulates a substantial re-stratification of the UML.The re-stratification is characterized by a sharpening and at the same time a deepening of the thermocline down to 40 m near the Finish coast and export of cold water to the surface near the Estonian coast (Fig. 6). Figure 6a and b show maps of the turbocline depth on the 16 May 2011.The turbocline depth is defined as the depth at which the vertical eddy diffusivity coefficient falls below a given value (here taken equal to background value of 5 cm 2 s −1 ) and can be interpreted as a maximum penetration depth of the turbulent motion in the surface layer.
According to Fig. 6a and b presenting solutions on 2 and 0.5 km grids respectively, the turbocline depth reaches the maximum in the areas near the Finnish coast where the convection is a dominant factor in vertical mixing.We can note the significant differences in the spatial patterns of the turbocline for fine and rough resolutions.Solution on the 0.5 km grid shows a deeper and more complex thermocline pattern.It can be explained by the fact that small-scale frontal structures induced by strong horizontal gradients and captured by the fine-resolution model lead to convective instabilities (Boccaletti et al., 2007) acting to locally re-stratify UML.The model with 2 km resolution cannot resolve submesoscale frontal features and high values (compare to fine resolution) of lateral diffusion coefficients act to smooth the front in other words decreasing potential energy of the front.Unfortunately, few data are available for validation of these differences.Locations of CTD profiles on 16 May are marked as points I, II, III in Fig. 6a and c. Figure 6 (I, II, III) shows the vertical profiles of temperature at locations near the Finish coast.At panel (I) the UML depth for the 2 km resolution model (dashed black line) is shallower than the observed UML depth (solid black line) by 13 m.At the same time, observations and the 0.5 km resolution model (grey line) temperature are almost collocated, and UML depth reaches 40 m.At panel (II) the modelled UML depth is overestimated, but the misfit reaches 7 m for the 2 km resolution model and only 3 m for the 0.5 km resolution model.
We cannot compare the UML depth from the results presented at panel III since none of the models were able to reproduce lateral intrusions observed.The low model performance at this point can be explained by the proximity of the frontal zone between coastal and deep water masses due to the upwelling.We assume that a small error in the predicted location of the front can lead to a serious misfits in the vertical profile.Note also that point (III) is located in a zone of rapid turbocline depth variations (see Fig. 6a and b).This fact confirms a complex front structure, which is formed by the set of randomly spaced small-scale features.The deterministic model can only predict their appearance but not the exact location.
Figure 7 presents evolution of the thermocline through the season.Left panels present the maximum depth of the turbocline and thermocline for the May when the thermocline was formed.Right panels present the same but for the period from 1 June to 28 July.This period ends just before the upwelling in July-August from which the UML erosion begins.Thermocline depth was defined as the depth of 3.5 • C isotherm (see Fig. 4).As seen in the presented data, turbulent mixing during the upwelling in May was the strongest throughout the season (see Fig. 7b).At the same time, the increase in the 3.5 • C isotherm depth up to 45 m during June-July is not accomplished by any considerable turbulent activity (maximum turbocline depth during June-July do not exceed 20 m for the most of the area of the GOF).Taking into consideration the low value of the background vertical diffusivity coefficient (10 −6 m 2 s −1 ), this fact highlights the importance of the advective processes for the formation of the shape and depth of the thermocline.Advective processes resulting in deepening of the isotherm are initiated by intrusion of warm The sensitivity of the model solution to increased horizontal resolution is manifested in the different intrusion propagations to the east (compare the right plots on Fig. 7d and  f).Density fronts associated with the intrusion are a source of baroclinic instability, which are resolved differently by the 0.5 km eddy permitting configuration (Fig. 7c ) compared to the 2 km configuration (Fig. 7e).

Discussion and conclusions
We used the state-of-the-art modelling framework NEMO, initially developed for the open ocean, to build an eddy-resolving model of the GOF.To evaluate the model skill and performance two different solutions where obtained: a commonly used 2 km grid and 0.5 km eddy-resolving fine grid.
With the resolution of 0.5 km the model starts to resolve submesoscale eddies.In the ocean, submesoscales are scales of motion equal to or less than the baroclinic Rossby radius of deformation.For the GOF the baroclinic Rossby radius varies between 2 and 4 km and we need at least 4 points to resolve the eddy.According to Gent and McWilliams (1990), the eddies can act to re-stratify the UML of the ocean, causing the vertical transport through the thermocline.
By moving from 2 to 0.5 km it is logical to expect an intensification of vertical movements induced by smaller vortices resolution.Figure 8 presents the comparison of vertical velocity absolute values for 2 km and 500 m resolutions.The fields are averaged for the depth of 5 m and a 5-day period in May characterized by high-intensity wind-induced dynamics.The main features of the horizontal distribution of the vertical velocity, including the regions of extreme values are similar in both cases.However, on a finer grid structures resembling meanders, currents, and filaments appeared in the middle of the bay at the Estonian coast as well as near the Finland coast where there is a set of point maxima.Both of these small-scale features are absent at coarse grid.It is important to note that the difference in the vertical velocity field appear mainly in the upper mixed layer of the sea.Below the pycnocline the vertical velocity patterns in both cases are very similar.Thus, marked differences could be attributed to the vortex centres of submesoscale eddies, but this assumption is not confirmed by visual horizontal velocity field analysis: explicit vortices are absent in ultraviolet (UV) horizontal field.An alternative hypothesis links these features with local elevations of the bottom topography.
An additional effect of resolved lateral submesoscale processes was investigated in Sect. 4. It was shown that submesoscale motion affects the plume propagation caused by salty water intrusion to the GOF from the Baltic Sea.Generally speaking this process has found to be dominant in the formation of the shape of termocline through the summer season, while the depth of UML was formed by an intensive mixing during spring upwelling.In both cases advective processes act as the main "driving force".
The presented model demonstrates a substantial improvement in the basin stratification compared to previous numerical studies.The traditional point of view is that the smallscale processes, such as turbulence, provide the majority of mixing in the vertical direction.Most commonly the 1-D approach is used to set up vertical mixing by tuning a turbulent scheme.The GOF, an enclosed basin with complex bathymetry and strong stratification mixed layer dynamics can be strongly affected by lateral advective processes.Adequate representation of lateral processes by the model let us decrease the role of background constants in a turbulent mixing scheme (we set them to minimum possible values).This simplifies the traditional trade-off between the depth and sharpness of the thermocline.Setting the background values of vertical eddy viscosity and diffusivity to 10 −5 and 10 −7 m 2 s −1 respectively, lets us keep the sharp form of the thermocline and halocline while the UML depth corresponds to observations.Since the time period of the runs was rather short (less than 1 year) and the model had not been used before, it is obvious that the values of some parameters might have been somewhat improperly chosen for use in this study.Through fine tuning of the model better results could probably be obtained.However, the focus in this study was to examine the differences arising from different horizontal resolutions, the fact that model parameters were similar in each case should be considered to be far more important than the quantitative agreement between observations and model results.Actually, it was shown that the model results for both resolutions are in reasonable agreement with available observations.In some cases the 0.5 km model performs better and at the same time there are areas not covered by observations where we can note more substantial difference between models.It is found that simulations resolving submesoscale are characterized by the deeper UML with more complex structure in the regions of the GOF directly affected by the upwelling/downwelling.
The GOF is a highly dynamic region with lateral currents causing temperature contrasts and/or rapid temporal variations on the surface.From the satellite picture we can identify whether the model properly reproduces the frontal structure at the surface.For example, the temperature drop during an upwelling event and resulting temperature contrast at the surface reached 2.5 • C. We assume it to be a considerably more substantial signal compared to known uncertainties of satellite SST measurements (0.4 • C; https: //podaac.jpl.nasa.gov).The usage of results of hydrodynamic modelling together with SST information can provide an extended analysis and deeper understanding of the upwelling process.Re-stratification of the UML caused by upwelling results in changes of the SST pattern that can be observed from satellites.From the comparison of modelled and observed satellite SST we can identify whether the model reproduces the stratification itself and as a result properly reproduces the frontal structure at the surface.
Refinement of the model resolution below the level of 0.5 km would be of limited benefit in a hydrostatic model.For the purpose of deep investigation of submesoscale processes in GOF such as transport across the UML and on/offshore the non-hydrostatic formulation is needed.It lets us avoid "artificial smoothing" of the velocity field.Other possible improvements of the model performance, which we are planning for the next steps, will include sensitivity tests for the different boundary conditions with higher spatial resolution at the open boundary and surface, and utilization of recently available data with high spatial coverage from the expeditions during the Gulf of Finland Year 2014.
The Baltic Environmental Database -BED -was initiated in 1990 as part of the research project "Large-scale Environmental Effects and Ecological Processes in the Baltic Sea".The basic idea behind the Baltic Environmental Database has been to make available data sets on the conditions in the Baltic Sea and on forcing functions.Presently, the database is placed at the Baltic Nest Institute (BNI) at the Stockholm University http://www.balticnest.org/bed.MODIS Sea Surface Temperature product generated from the sensors installed on two satellites Terra and Aqua by NASA Goddard Space Flight Center.The gridded global product with spatial resolution 1 km available from http://dx.doi.org/10.5067/TERRA/MODIS_OC.2014.0 and http://dx.doi.org/10.5067/AQUA/MODIS_OC.2014.0.

Figure 1 .
Figure 1.The bathymetry of the Baltic Sea.Red line -open boundary of the model domain, yellow line -location of the meridional cross section for Fig. 2.

Figure 3 .
Figure 3. SST on 20 May 2011: (a) MODIS SST; (b) GLS with a Galperin limit 0.53, and horizontal resolution 0.5 km; (c) GLS with a Galperin limit 0.53 and horizontal resolution 2 km; (d) GLS with a Galperin limit 0.26 and horizontal resolution 2 km; (e) TKE with convective adjustment and horizontal resolution 2 km; (f) GLS with a Galperin limit 0.07 and horizontal resolution 2 km.

Figure 4 .
Figure 4. Averaged vertical profiles of temperature and salinity in western (a, d), central (b, e), and eastern (c, f) parts of the GOF for the period 20 July-5 August 2011.Grey lines -CTD data with standard deviation corridors, solid and dashed black lines -are modelled on grids 0.5 and 2 km.

Figure 5 .
Figure 5. SST maps of the GOF on 2 August 2011: (a) MODIS data, (b) and (c) modelled SST on grids 0.5 and 2 km respectively.

Figure 6 .
Figure 6.Modelled turbocline depth (m) in the GOF on 20 May 2011: (a) and (b) horizontal distributions on grids 0.5 and 2 km respectively; (I), (II), and (III) -vertical profiles of temperature at the locations marked on maps (a) and (b).

Figure 8 .
Figure 8. Vertical velocity absolute values (log scale) averaged for the depth of 5 m and a 5-day interval: (a) 2 km model grid, (b) 500 km model grid.