The geographic area and season suggest the potential presence of larger-scale horizontal buoyance gradients. Salinity increases from east to west due to the large river discharge at the eastern end of the gulf (Ylöstalo et al., 2016) and a cross-shore temperature gradient could form due to faster warming in shallow areas (Lips et al., 2014; Lips and Lips, 2017). Furthermore, pycnoclines have varying cross-shore inclination depending on the local wind forcing (Liblik and Lips, 2017).

The data set analyzed in the study was gathered from 9 May to 6 June 2018 using an autonomous underwater vehicle – the Slocum G2 Glider MIA deployed in the Gulf of Finland slope area at 59.57° N–25.05° E, Baltic Sea (Figure 1). The glider collected oceanographic data along the 18 km long transect (S–N direction, sampled 26 times). The Slocum glider can perform YOs without reaching the surface. The mission set-up determined the glider to perform 4 YOs between surfacings (i.e., a cycle). The vehicle started turnaround either 4 m before the surface or 6 m before the seafloor. The glider moved at a horizontal speed of 0.32 ± 0.09 m s^-1. The average distance between the profiles near the surface was 287 ± 39 m (it reduces with depth due to the sawtooth-shaped sampling trajectory), and a profile took 7.8 ± 0.5 min to complete at a characteristic depth of 80–90 m. Therefore, the glider completed a cycle approximately in an hour after covering about a km horizontally.

Temperature and conductivity (salinity) profiles were collected with a pumped Seabird conductivity-temperature-depth recorder (CTD; model G-1451). Technical specifications of the sensors set the accuracy for pressure 0.1% of the full range (0–200 m), for temperature 0.002°C, and for conductivity 0.01 mS cm^-1. The CTD had a sampling rate of 0.5 Hz, and both descending and ascending profiles were recorded. Altogether over 4000 CTD profiles were gathered and analyzed. The raw data was quality controlled, and the coefficients to account for the response time of the sensors and the thermal lag were defined to minimize differences between two consecutive CTD profiles (Appendix). The typical vertical resolution of the raw CTD data was around 0.2 dbar. The half YOs were bin-averaged to a uniform 0.5 dbar vertical grid and arranged as profiles, meaning an average time and coordinates were attributed to a half YO. While the glider data has time-varying orientation in horizontal direction, for this analysis, each transect was interpolated on a regular grid with a latitudinal step of 0.0018° (if not specified otherwise). The horizontal step, equal to 200 m, was slightly larger than the distance between half YOs generally (median 135 ± 56 m) and thus, the interpolated fields were smoother. Note that this median distance between half YOs corresponds to the horizontal distance between the midpoints of the half YOs.

Velocity calculations from the glider observations are limited. Absolute geostrophic velocity can be produced using the hydrographic data and the additional information from the depth-averaged current, which is estimated over a cycle based on dead reckoning and GPS fixes at the surface. The relative geostrophic velocity components were calculated using the dynamic method based on the geostrophic relationships

where the geopotential, ϕ, is proportional to dynamic height, D, as

Coriolis parameter is noted with f and gravitational acceleration with g. Therefore, the horizontal pressure gradient is presented through the geopotential slope of an isobaric surface. The relative geostrophic velocity is evaluated using a dynamic height anomaly relative to a reference pressure. The dynamic height that depends on the density of the seawater was determined from the measured temperature and salinity profiles according to McDougall and Barker (2011) (IOC et al., 2010).

The relative geostrophic velocity component was estimated from the glider CTD data at a lateral scale of 5 km. The choice of reference level (85 dbar) was based on the lowest variability (i.e., standard deviation) of density along isobars, interpreted as the effect of the baroclinic velocity field becoming negligible there. To involve shallower profiles in the geostrophic velocity estimations, the stepped no-motion level method described in Rubio et al. (2009) was used. While the 2D nature of the glider data allows us to estimate only the across-transect component of the geostrophic velocity, the depth-averaged velocity components were rotated to be compatible with cross- and along-transect components. On average, the transect was tilted 14° to the left from the S–N axis. As the glider orientation was relatively stable along the path, the approach of using an average tilt was adopted for simplicity.

Buoyancy is expressed as b y = g ( 1 − ρ ρ 0 ) , where ρ₀ is reference density of 1003 kg m^-3. Along-transect buoyancy gradients, b y = ∂ b ∂ y , with a step of 5 km, were calculated to characterize the observed mesoscale frontal structure. Note that the notation y stands for the along-transect component and notation x for the across-transect component from here on. The stratification of the water column was characterized by the squared Brunt-Väisälä frequency (vertical buoyancy gradient) N 2 = − g ρ ∂ ρ ∂ z , where ρ is density and ∂ ρ ∂ z its vertical gradient.

The analysis of smaller scale variability focused on temperature variations on the isopycnal surfaces. At the same density surfaces, water properties act as passive tracers, and the effect of internal waves is assumed negligible (Rudnick and Cole, 2011). The smaller scale signal was defined as the deviations from the smoothed isopycnal temperature field (moving average of 4 km). The length scale was chosen based on the internal Rossby deformation radius, which is typically 2–4 km in the Gulf of Finland (Alenius et al., 2003). The deviations visualize the patterns created by the smaller scale motions, including the submesoscale range.

The dominating contribution to the energy fluxes at the submesoscale range in the transfer of energy from large to dissipative scales has remained debatable. We calculated the power spectrum of isopycnal temperature in the wavenumber domain, composing the data as a function of cumulative distance. Isopycnal temperature fields on a regular grid were constructed due to irregular sampling in both space and time. First, the temperature profiles were averaged in density bins (step 0.1 kg m^-3), and then interpolated on a grid with a constant horizontal step of 150 m. Averaging in advance allowed us to reduce the effect of spatial and temporal smoothening due to interpolation, which tends to bias spectra towards steepening (e.g., Wortham and Wunsch, 2014). Before taking Fast Fourier Transform, the linear trend was removed and Hanning window applied. The smoothed power spectrum was obtained using 31 elementary frequency bands for averaging. The number of elementary frequency bands can be equalized with the degrees of freedom to find confidence limits (95%).

The wavenumber spectra of isopycnal temperature were calculated in the density range of 3.6–6.0 kg m^-3 with a step of 0.1 kg m^-3 (covered depths from 3–50 m). The spectral slopes were estimated between the spatial scales of 1–10 km. The resolved horizontal length scale was not separated because the short data series limit the resolution of the spectrum and the scales where the transition between distinct dynamics occurs are not known. The glider sampled the 18 km long section repeatedly and thus, the lower wavenumbers than the one corresponding to 10 km were disregarded. The submesoscale in the Baltic Sea can be expected to be active around 1–2 km. The spectrum at <1 km was not analyzed due to greater uncertainties at smaller scales (Nyquist frequency corresponds to the spatial scale of 0.3 km).

Atmospheric forcing is presented through wind, net surface heat flux, and freshwater flux. Parameters defining heat exchange between the atmosphere and the sea surface were extracted on the product grid cells that covered the study area (59.50–59.75° N, 25.00–25.25° E) for May–June 2018 from European Centre for Medium-Range Weather Forecasts atmospheric reanalysis data set (ERA5 hourly data on single levels from 1979 to present, horizontal resolution of 0.25°×0.25°, Hersbach et al., 2020; Hersbach et al., 2022). Hourly u- and v-components of wind at 10 m height, mean surface net short-wave and long-wave radiation flux, mean surface latent and sensible heat flux, mean total precipitation rate, and mean evaporation rate were used.

The net freshwater flux is converted to an equivalent heat flux as also used by Giddy et al. (2021)

Wind components taken from the nearest cell to the glider section (59.75° N, 25.25° E) were smoothed by a Gaussian low pass filter for 6 h. Wind stress is expressed as

where c_d is the drag coefficient (1.2 · 10^-3; Large and Pond, 1981), ρ_a the air density (1.2 kg m^-3), and U_a stands for wind speed. The cumulative wind stress (in N m^-2 h) is the product of wind stress and its duration, demonstrating the impact of the wind during a longer period.

Winds and buoyancy fluxes drive mixing in the near-surface layers and affect the development of stratification. In general, wind speed during the study period was <4.5 m s^-1 (Figure 2). At the beginning of the glider mission, E winds prevailed, being relatively strong on 9–12 May and subsiding afterwards (to an average of 2.5 m s^-1). On 18–19 May, stronger northerly winds occurred for a short period, while the last third of May was characterized by alternating E and W winds. The W winds were stronger, on average, 5.4 ± 2.0 m s^-1 compared to E winds of 3.8 ± 0.9 m s^-1. The strong W–NW winds started to prevail on 1 June. The cumulative wind stress perpendicular to the glider transect (approximately positive towards E) picks up the three distinct periods during the mission – from the beginning until 19 May, wind forcing was to W, after that until 31 May, forcing varied between E and W, and from 1 June, it was to E (Figure 2).

The mission was conducted during the formation of the seasonal thermocline. The distributions of temperature, salinity, density, and geostrophic velocity were constructed (Figure 3 shows the composite plots of the upper 25 m from the beginning, middle and end of the mission, 10 days apart). During the mission, the sea surface (0–4 m) temperature increased from 6.7 to 15.2°C, and the salinity decreased from 5.3 to 4.2 g kg^-1. The decrease in salinity indicates the advection of fresher water, which possibly originates from the eastern part of the gulf. The temperatures below the developing seasonal thermocline in the CIL were around 2°C, and the thermocline, characterized by the maximum temperature gradient of about 3°C dbar^-1, was observed at 10–15 m by the beginning of June.

The mesoscale front was defined as an evident inclination in the density distribution, extending over the sampled transect. The horizontal (across-shore) buoyancy gradient formed as a response to the moderate NE–E winds that dominated for ten days prior to 19 May (Figure 2). The NE–E winds force surface waters to move offshore along the southern coast of the Gulf of Finland. Deep water rises in turn, tilting the pycnocline. Figure 4 shows the structure of the water column during the appearance of the front in the study window (~2.5-day period). The variability of buoyancy over the sampled section was small on 20–21 May, apart from the small area at the southern end of the transect, which indicated denser water rising (Figure 4A). However, a westward flow compatible with the geostrophic frontal jet appeared to be developing (Figure 4D). The frontal jet in thermal wind balance flows in the same direction as the wind. The buoyancy increased from south to north at 15 and 25 m depth on 21–23 May, indicating the presence of a horizontal buoyancy gradient (Figures 4B, C). The strong westward flow with velocities over 20 cm s^-1 spanning over the inclined density distribution suggested the presence of the frontal jet (Figures 4E, F).

The appearance of the front in the measurement window could be related to the stronger winds from N–NW on 19–20 May, meaning the generation of the front was not captured. The mesoscale front reaching the surface was not observed. Figure 5C shows smaller scale horizontal buoyancy gradients near the surface. They often had the opposite sign (a decrease of buoyancy from south to north) and were probably related to the fresher water advected into the study area. The mesoscale front persisted until June (Figure 5D), despite rather weak upwelling-favorable conditions alternating with the wind from the opposite direction in the last third of May. It is supported by the geostrophic velocities that exhibit persistence of the negative (westward) flow core. All measured temperature and salinity distributions with density distributions and geostrophic velocities are found in Supplement 1.

Horizontal variations are closely related to the vertical structure. Figures 5A, B illustrates the change in the vertical stratification in the northern and southern part of the transect, compatible with the lighter and denser side of the front, respectively. The beginning of the mission was characterized by weak haline stratification (N² = 0.001 – 0.002 s^-2). The near-surface thermal stratification started to develop on 15–16 May. The appearance of the front on 21–22 May promoted the strengthening of the stratification on the denser side of the front because of the interaction between surface heating and vertical motions associated with the mesoscale front. The stratification was significantly stronger on the denser side of the front than on the lighter side, where actually two pycnoclines were observed after the appearance of the front.

The pycnoclines weakened and strengthened daily on both sides of the front. Pycnoclines on the lighter side appeared to be forced to merge episodically, indicating the alternating of the forcing and/or processes in favor and against the persistence of the front. We suggest that the weakening of the upwelling-favorable forcing arrested the strengthening of the front. Instead of collapsing abruptly, the front hosted ageostrophic frontal processes that equalized the differences in the cross-front direction by slowing the development of the stratification on one side of the front and promoting it on the other. The uniform stratification formed by 3 June over the course of a day, coinciding with the vanishing of the mesoscale front in the study area. Compared to the denser side of the front, the stratification remained slightly weaker on the lighter side.

To demonstrate the evolving smaller scale dynamics in the background of the mesoscale front, we analyzed the smaller scale thermohaline variability during the appearance of the horizontal buoyancy gradient. Isopycnal temperature distributions contrast water mass characteristics and/or visualize local diapycnal mixing spots. Figure 4K shows an intrusion manifested as a two-way intersection at 59.675° N – a cold water patch (around 3°C) penetrates through several isopycnal layers from 4.1 to 3.7 kg m^-3 and a warm water patch (around 8°C) penetrates from 3.3 to 3.7 kg m^-3. A similar pattern was observed at 59.65–59.675° N on the following day (Figure 4L), but not the day before (Figure 4J). Therefore, this structure emerged simultaneously with the appearance of the front in the measurement window.

Figures 4G–I demonstrates the smaller scale variations in the temperature by showing the deviations from the surrounding mean (4 km). Two adjacent temperature deviations with opposite signs characterizing the pattern indicated that the observed pattern consisted of two motions in opposite directions. Focusing on the described pattern, the deviations revealed extrema of 1.2°C and -1.4°C at 59.675° N on 21–22 May (Figure 4H). The extrema were larger on the following day (1.3°C and -2.1°C at 59.65–59.675° N, Figure 4I). The sign of the deviation indicated the water to originate from either shallower or deeper layers assuming that higher temperatures had their origin at the sea surface.

The last third of May was characterized by the persistence of the mesoscale front as well as the development of seasonal thermal stratification and the upper mixed layer, where the smaller scale horizontal buoyancy gradients formed. The diurnal cycle of net surface heat flux had a maximum of 400–500 W m^-2 (Figure 6A). The freshwater flux was insignificant compared to it. Previously described tracer patterns, which consisted of warmer water plunging down and colder water penetrating upward, were mapped repeatedly on the background of the mesoscale front. From here on, those intrusions are referred to as a submesoscale pattern.

The structure of the submesoscale patterns suggests that they were traces of cross-front ageostrophic secondary circulation that is known to be generated while a horizontal buoyancy gradient is subject to baroclinic submesoscale, wind-driven and/or baroclinic instability. To identify whether the conditions favored the growth of different instabilities, the potential vorticity, and the equivalent heat fluxes related to restratification processes energized by surface wind forcing or the release of potential energy by baroclinic instability were calculated. Although the heat fluxes demonstrate the impact on the stratification, the elevated values also suggest the potential presence of frontal circulations.

The most pronounced submesoscale patterns captured during 21–23 May could be resulted from different physical processes. The period of the appearance of the mesoscale front in the study window coincides with the most significant episode of the equivalent heat flux related to baroclinic instability (Figure 6B), suggesting the pattern in Figure 4H to result from the frontal circulation generated through the release of potential energy by baroclinic instability. The equivalent heat flux related to wind-induced instability episodically also had comparable magnitudes (Figure 6B). However, this estimate is more likely to be more biased than the estimates suggesting baroclinic instability events. The negative equivalent heat flux from the surface wind forcing suggests destabilization of the water column and possible generation of the overturning circulation, extracting potential vorticity from the pycnocline. The water column structure captured on 23–24 May can be found in Supplement 1.

Downfront winds make the flow susceptible to symmetric instability by reducing potential vorticity. Figure 6C shows the component of wind stress perpendicular to the sampled transect, indicating downfront (upfront) orientation if positive (negative). The unstable flow is characterized by negative potential vorticity. Figure 7A demonstrates the dominance of positive potential vorticity that is in agreement with the development of the seasonal thermocline (relatively strong stratification) near the surface, but a stronger negative signal is found in the second half of the mission. Although the variety of instabilities can develop when q< 0, the ϕ_Rib criteria suggested the presence of only symmetric instability (i.e., 0< Ri_b< 1 and −90°<ϕ_Rib <ϕ_c ). The balanced Richardson number, which was always positive, showed the largest potential for symmetric instability in association with the mesoscale front on 23–24 May (elevated values at 5–10 m, Figures 7B, C). However, the Obukhov length was comparable with the depth of the upper mixed layer and the boundary of the CIL, meaning the criteria –H< z< -h was probably not satisfied. Although potential vorticity was reduced notably in the upper mixed layer in June, the relatively large Obukhov length suggested the convective layer depth to be similar in magnitude to the mixed layer depth, making the generation of symmetric instability still unlikely.

The analysis suggests that the conditions favored the generation of the wind-induced and baroclinic instability, supporting the assumptions for the prevalence of the submesoscale despite the limitations: 1) the calculations relied on the assumption that the glider sampled perpendicular to the horizontal buoyancy gradient; 2) neither the position of the horizontal buoyancy gradients nor the angle at which the glider sampled is known; 3) poorly resolved individual submesoscale features, while sampling at an angle, can go unnoticed due to not satisfying the conditions of the analysis. Further conclusions should be made with caution as several assumptions were made (see section 2.5).

Spatial variability of isopycnal temperature was analyzed in the upper part of the water column. The weakening of spectral energy on deeper isopycnals refers to the decrease in the temperature variability with depth (Figure 8). The slopes were estimated for the length scale from 1–10 km, and the individual σ-layers characterized by similar spectral slopes were averaged. The slopes for averaged spectra in the density ranges of 3.6–3.9 (3–10 m), 4.3–4.6 (10–20 m), 4.9–5.2 (15–25 m), and 5.7–6.0 (35–50 m) kg m^-3 were -2.14, -1.72, -1.84, -1.20, respectively. The slope was the gentlest on the isopycnals positioned in the CIL (5.7–6.0 kg m^-3). The slope near -1 suggests the applicability of the interior-quasigeostrophic theory and the domination of large-scale motions. Upper layers exhibit clearly steeper slopes, suggesting depth-dependence. The atmospheric forcing dominates in the uppermost layers, and energetic smaller scales stir the generated variability relatively efficiently. However, the trend of shallower spectral slopes with increasing depth did not express in the layers associated with the mesoscale front. This indicates the activity of frontal processes and their contribution to the energy flux.

Constantly changing stratification affects the spectra but the extent of it is unknown. Additionally, rather short data series limit the spectral analysis. For that reason, the emphasis is more on the dynamics than the certain values of the slopes, which were near k^-2 in the upper part of the water column though. The observations support the suggestion that the submesoscale flows emerge within the mesoscale front revealed as the enhanced tracer variability in this layer compared to the variability in the layers above or below.

The structure of the frontal submesoscale processes has remained theoretical despite the effort that has been devoted to observing the smaller scale variations. Although recent observational studies well support the presence of the frontal instabilities in the quiescent open-ocean environment, in the Subantarctic region, as well as near boundary currents (e.g., Callies et al., 2015; du Plessis et al., 2017; Yu et al., 2019b), few provide the information on the projection of submesoscale flows in the measurements (Pietri et al., 2013; de Verneil et al., 2019; Pérez et al., 2022). For instance, de Verneil et al. (2019) have demonstrated the orientation of fine-scale phytoplankton patches to be consistent with an ageostrophic secondary circulation generated by frontogenesis, providing, thus, the evidence on submesoscale flows modifying the vertical dynamics of phytoplankton.

According to the dynamical characterization, the coastal sea provides favorable conditions for the submesoscale processes. The transition to the open sea favors the colliding of water masses whether through the meeting of different flow structures or coastal processes like upwelling events, promoting, therefore, the formation of larger and smaller scale horizontal buoyancy gradients. By observing the coastal-offshore area at the end of spring 2018, we captured the appearance of the mesoscale front. The horizontal buoyancy gradient formed due to the wind-induced upwelling forced by the prevailing NE–E winds during the first ten days of the survey.

Simultaneously with the appearance of the mesoscale front, a pronounced submesoscale pattern emerged. The pattern showed upwelling of cold water on the lighter side of the front accompanied by downwelling of warm water on the denser side. The temperatures indicated the origin of the water from deeper or shallower layers, respectively, and, thus, suggested the presence of two adjacent movements in different directions. We suggest the presence of ageostrophic secondary circulation with a width of a couple of km on the profile of the mesoscale front between 5–10 m (Figure 3B). To support this suggestion, we attempt to infer the driving mechanisms that could explain the observed submesoscale activity.

A secondary circulation emerges in the form of an overturning cell as the front is subject to baroclinic, wind-induced, and/or baroclinic submesoscale instability. The submesoscale analysis revealed the large potential role of baroclinic instability on 21–22 May, suggesting that the observed submesoscale features resulted from the secondary circulations generated through the release of potential energy stored in the horizontal buoyancy gradient. The baroclinic instability can act together with surface frontogenesis (Lapeyre et al., 2006) but the surface signature was not observed in the study window. Additionally, the position of the captured feature on top of the mesoscale front is not compatible with frontogenetic secondary circulation that tends to subduct water below the front (Spall, 1995). Baroclinic instability occurs over a period of days (Boccaletti et al., 2007), indicating the possible interaction with the wind-induced instability. The analysis showed the effect of downfront winds that would counteract the restratification on the following days. On 22–23 May, the slope of the isopycnals was relaxed, and a few m deep upper mixed layer had formed (Figure 4F). Destabilizing atmospheric forcing can also favor symmetric instability. However, this is unlikely because the upper mixed layer needs to be deeper than the convective layer for overturning motions to dominate over convective mixing induced by surface forcing (Thomas et al., 2013).

The evolution of smaller scale variability in the background of the mesoscale front suggests the presence of several dynamical processes. Similar tracer patterns with various extent and intensity were mapped repeatedly until the end of the month, but the limitations of the glider observations complicate the interpretations. Fronts are a dynamical phenomenon and the 2D nature of the glider measurements prevents following the evolution of the front. In addition, observational biases should be considered when inferring the dynamical mechanism because of the adoptions for the glider data assuming the sampling perpendicular to the horizontal buoyancy gradient (Thompson et al., 2016). Further, the captured features are subject to soothing due to the effect of time-space aliasing.

Horizontal variations are related to atmospheric forcing and to stirring by larger-scale flow fields. In the Baltic Sea, inclined pycnoclines and coastal upwelling events that appear as a response to winds are frequent (Lehmann et al., 2012; Liblik and Lips, 2017). The Gulf of Finland, a sub-basin of the Baltic Sea, is elongated in the W–E direction. The moderate winds blowing along the gulf are common (Keevallik and Soomere, 2014). The NE–E winds favor upwelling events near the southern coast of the gulf. A corresponding development of hydrographic fields was captured as a mesoscale front in the glider measurement window in May 2018. Although E and W winds alternated after 19 May, the mesoscale front persisted until early June.

The wind has a role in allowing the submesoscale flows to arise due to the disturbance of the balance when the external forces reduce and/or change. We have indicated that on 21–22 May, baroclinic instability appeared to develop during the temporary decrease in wind stress. Interpreting the orientation of the wind with respect to the geostrophic shear is challenging, but a rough estimate would be that the front was oriented approximately in the W–E direction at the sampling site. The sampled section positioned approximately along the S–N axis and the inclination in the density distribution that could be manifested as a front was observable. In accordance with this assumption, the wind changing from E to W direction indicates that mostly the wind had the component either up- or downfront. Thus, this variability of the wind would cause the alignment between the front and the wind direction to be lost or gained episodically (Mahadevan et al., 2010).

The horizontal buoyancy gradients affect the stratification directly through the submesoscale processes. While the stratification is generated and destroyed through diabatic processes and frictional forces, frontal submesoscale can convert gradients from horizontal to vertical (McWilliams, 2016). Lateral variability introduced by the advection of water masses favors the appearance of horizontal buoyancy gradients and, in turn, suggests the emergence of the submesoscale. The glider observations revealed the complex structure of the water column. Atmospheric forcing acting on the sea surface contributed to the formation of the upper mixed layer and seasonal thermal stratification. The appearance of the mesoscale front created favorable conditions for the stratification to strengthen more rapidly in the southern part of the sampled transect as the mesoscale dynamics raised the isopycnals near the surface. We observed the stratification changing daily and developing differently on the denser and lighter side of the front (compatible with the southern and the northern part of sampled transect, respectively; Figures 5A, B). While the denser side was characterized by the strong stratification near the surface, two pycnoclines were observed on the lighter side.

We propose that the mesoscale front hosted ageostrophic frontal processes as the upwelling-favorable forcing weakened. The submesoscale circulations equalized the differences in the cross-front direction by slowing the development of the stratification on one side of the front and promoting it on the other. The overturning circulation induced by baroclinic instability slides dense water under light, acting to slump the front from the horizontal to vertical (Boccaletti et al., 2007). Thus, the baroclinic instability acted toward restratification on 21–22 May. The wind increased on the following days, suggesting the possible interaction with wind-induced instability. The overturning circulations acting toward restratification appeared to be dominating but the temporary return of the NE–E winds supported the persistence of the front. Two pycnoclines were observed on the lighter side of the front again on 27 May. However, the front appeared to shift toward the open sea (northern part), suggesting that partly the frontal circulations managed to restratify the water column. The presence of the submesoscale flows within the horizontal buoyancy gradient could have contributed to the development of the seasonal stratification. The effect of the submesoscale can be especially important during seasonal warming because stratification events can occur more rapidly than solely surface heating could drive (du Plessis et al., 2019). Enhanced stratification favors primary production (Swart et al., 2015).

While measuring submesoscale flow fields is challenging, the underlying dynamics can be deduced by investigating the variability of tracer fields along the same density surfaces (Callies and Ferrari, 2013; Jaeger et al., 2020). Submesoscale flows can cause substantial heterogeneity in tracer distribution because of their dynamics, which by its nature should connect fully two- and three-dimensional motions (McWilliams, 2016). Submesoscale processes have gained more attention in the Baltic Sea only in recent years, though (Lips et al., 2016; Carpenter et al., 2020). Lips et al. (2016) have presented the spectrum of an active tracer in the sub-surface layers (note that the tracer variance along the isobaric surfaces was evaluated). The spectra showed slopes closing to -2 at 1–10 km during periods of high variability of temperature due to upwellings. The discrepancy from -3 that is predicted by quasigeostrophic theory (Charney, 1971) was interpreted as the contribution of the ageostrophic submesoscale processes to the energy cascade.

The studies from the ocean environment are well presented (e.g., Cole and Rudnick, 2012; Callies and Ferrari, 2013; Kunze et al., 2015; Jaeger et al., 2020). However, the processes setting the spectra have remained obscure. In the open ocean, the -2 slope independent of depth has been reported at the submesoscale range (Cole and Rudnick, 2012; Callies and Ferrari, 2013). The main contribution to the energy at the submesoscale comes possibly from unbalanced flows. The characteristics of the Baltic Sea (e.g., the shallow depth, strong stratification, seasonality, negligible tidal forcing) distinguish the sea from the ocean, forming a dynamically complex environment (Tuomi et al., 2012). There are some studies carried out in the stratified environment. Kunze et al. (2015) found tracer variance spectra sloping with a k^-2 at the horizontal length scale 0.03–10 km in the seasonal pycnocline at 20–60 m depth. Evident fronts were not observed at every study site, though. On the other hand, Jaeger et al. (2020) showed steep slopes (k^-3) at 1–10 km in the stratified upper 75 m. Above 10 km scales the spectral slope was closer to -2. They proposed that ageostrophic dynamics could cascade variance in the pycnocline more rapidly than predicted. We, however, found gentler spectral fall-off in the layers associated with the mesoscale front compared to the layers above. This indicates that frontal submesoscale preserves tracer variability longer than the turbulence in the upper layers. The short data series acts as both an advantage and a disadvantage in this study while allowing us to analyze one specific situation but prevents broader conclusions. We will test this suggestion in a future study by incorporating a wide range of glider data sets from different seasons.

Observational spectra are independent of predictions and, thus, extensive data collections contribute to shaping understanding of the energy cascade. The observations have revealed the discrepancies from the quasigeostrophy, but often the limitations prevent identifying the reasons. We conclude that the contribution of the ageostrophic frontal effects is likely in the upper part of the water column (0–25 m in the Baltic Sea). The depth-dependence may be related to the stratified environment as it was noticed here and by Jaeger et al. (2020). Flatter spectral slopes in the deeper part of the water column coincide with stirring by two-dimensional homogeneous turbulence that leads to a tracer spectrum following k^–1 (Callies and Ferrari, 2013).