The Observed Seasonal Cycle of Submesoscale Processes in the Antarctic Marginal Ice Zone

Submesoscale flows are energetic processes in the upper ocean that occur at scales of a few kilometers. These phenomena are important to climate because they can strongly impact vertical motions of heat, carbon, and biogeochemical properties between the surface and deep ocean. Their presence and magnitude are poorly understood in the Antarctic marginal ice zone (MIZ) and under sea ice due to a lack of observations. We present tagged southern elephant seal hydrographic data from the eastern Weddell Sea, over the austral late summer to spring in 2008, to assess the magnitude and seasonality of submesoscale processes in the Antarctic MIZ. We consider the linkages between sea ice cover, ocean stress variability, and mixed layer depth and how these influence the submesoscale fluxes. We find that there is a strong seasonal cycle in the magnitude of the equivalent heat flux from submesoscale flows (240 Wm − 2 ), with peak activity occurring in midwinter. Periods of increased variability in sea ice cover, likely associated with leads, are also linked to restratification through submesoscale fluxes in midwinter. These findings have implications for understanding the magnitude and variability of upper ocean stratification and therefore highlights the importance of developing parameterizations for submesoscale processes in the Antarctic MIZ within global climate models. processes can have a strong effect on the upper ocean and can cause the mixed layer depth to shoal even in the middle of winter. There is also a clear seasonal cycle in the magnitude of the submesoscale activity. This study indicates that submesoscale processes should not be overlooked as a critical component of upper ocean processes in the Antarctic sea ice zone.


Introduction
Globally, the upper ocean plays an important role in water mass transformation, entraining and ventilating deeper waters or downwelling dense surface waters through changes in the mixed layer depth (MLD). The upper ocean acts as a gateway between the atmosphere and ocean interior, so understanding the processes that cause variations in the MLD is of first-order importance for modeling future heat and carbon ocean sequestration (Abernathey et al., 2016;Sallée et al., 2012). Compared to total ocean heat and carbon uptake, the Southern Ocean accounts for respectively 75% and 40% of the global values (DeVries, 2014;Frölicher et al., 2015). However, there is still considerable uncertainty over the temporal and spatial variability in these vertical fluxes, largely due to a lack of adequate observations on relevant scales and seasons, with a bias in current observations to summertime months (Gruber et al., 2019;Sallée et al., 2013). Recent studies using data from Southern Ocean Carbon and Climate Observations and Modeling (SOCCOM) floats revealed that some parts of the Southern Ocean may be a source of atmospheric carbon dioxide during winter months when sea ice is present, although the mechanism for this output of carbon is still unknown due to limited observations in this season (Gray et al., 2018;Gruber et al., 2019).

10.1029/2019JC015587
The growth and melt of sea ice can have a significant impact on the Southern Ocean as a freshwater sink or source and modifies atmosphere-ocean fluxes of momentum and heat. These processes in turn affect both water mass transformation and the transport of heat and carbon through the upper ocean (Abernathey et al., 2016;Asplin et al., 2014;Pellichero et al., 2017Pellichero et al., , 2018. However, the fifth phase of the Coupled Model Intercomparison Project (Taylor et al., 2012) showed global climate models have large variations in trends and distribution of Antarctic sea ice (Roach et al., 2018;Swart & Fyfe, 2013;Turner et al., 2013), with many characterized by insufficient autumn growth (Holmes et al., 2019). The errors in sea ice concentration distribution are in part due to poor representation of thermodynamic sea ice processes in the models, including the effects of lateral melt (Roach et al., 2018), which has been associated with the development of ocean eddies as a result of lateral buoyancy gradients (Horvat et al., 2016). This implies that current coupled climate models misrepresent the importance and magnitude of these lateral buoyancy gradients when sea ice is present in their parameterizations, likely due to the temporal and spatial sparsity of observations of oceanic conditions year-round but especially during wintertime.
Lateral buoyancy gradients, which can be caused by horizontal variations in temperature and salinity, are associated with frontal regions (Boccaletti et al., 2007;Mahadevan et al., 2010;Thomas & Ferrari, 2008). Marginal ice zones (MIZs), the region between the winter sea ice maximum extent and permanent (summer) sea ice, are also frequently characterized by strong lateral gradients in temperature and salinity, caused by sea ice growth or melt, and the interaction between open water and the atmosphere in leads between sea ice floes (Lu et al., 2015;Manucharyan & Thompson, 2017;Timmermans et al., 2012). These lateral buoyancy gradients provide the energy for the development of submesoscale ocean flows (eddies and jets), which typically occur on length scales of kilometers (comparable to the mixed layer deformation radius) and time scales of hours to days (Boccaletti et al., 2007;Fox-Kemper et al., 2008;Horvat et al., 2016). They occur when horizontal buoyancy gradients are equal to or greater than vertical buoyancy gradients and the vertical component of relative vorticity in the water column is comparable to, or greater than, the planetary vorticity (RoRi 2 = 1 Boccaletti et al., 2007). The resulting mixed layer eddies act to redistribute the available potential energy provided by the lateral buoyancy gradients through a vertical buoyancy flux, resulting in a shoaling of the mixed layer (Boccaletti et al., 2007;Klein et al., 2008;Mahadevan & Tandon, 2006). However, these eddies can also be enhanced or destroyed through interaction with submesoscale wind-driven flows (Du Plessis et al., 2019;Mahadevan et al., 2012;Thomas, 2005). The combined effect from submesoscale mixed layer eddies and wind-driven flows can be to either shoal or deepen the mixed layer on short time scales.
Such submesoscale processes thereby change the extent of the mixed layer, which in turn will affect the vertical transport of heat and carbon (McWilliams, 2016). When models are of high-enough resolution to resolve submesoscale processes, vertical velocities and transport increase to up to 5 times larger than when only mesoscale transport is considered Klein & Lapeyre, 2009;Su et al., 2018). In addition to the physical transport of carbon, variations in the mixed layer can impact light and nutrient levels, influencing the rate of biological productivity and the impact of biogeochemical cycles on carbon drawdown into the deep ocean interior Swart et al., 2015;Thomalla et al., 2011). Previous modeling and observational studies at the midlatitudes, with no ice present, indicate a seasonal cycle in submesoscale activity (Hosegood et al., 2013;Thompson et al., 2016), where deeper MLDs in the winter increase the available potential energy, resulting in stronger submesoscale flows Callies et al., 2015). This highlights the importance of year-round observations in order to resolve the seasonal variability.
The Southern Ocean is a region of deep mixed layers, strong fronts, and significant atmospheric forcing through variable heat fluxes and wind stress (Yuan et al., 2009). These characteristics predispose the Southern Ocean to active submesoscale variability, and regional models clearly show an increase in vertical exchange and shoaling that occurs when higher model resolutions are able to resolve submesoscale processes Rosso et al., 2014). Observational studies have also shown evidence of submesoscale activity in the lee of topographic features and associated with eddies or standing meanders linked to the Antarctic Circumpolar Current (ACC) Siegelman, O'Toole, et al., 2019;Viglione et al., 2018). The presence of submesoscale flows, and the balance between the mixed layer eddies and wind-driven flows, plays an important role in springtime stratification and consequently the initiation of springtime blooms in the sub-Antarctic regions (Du Plessis et al., 2017Swart et al., 2015). These studies span from October through to May, providing observations of active submesoscale fields from spring through to autumn. However, they all occur north of the typical extent of winter sea ice, resulting in no observational studies to corroborate the presence, magnitude, and variability of submesoscale flows in the Journal of Geophysical Research: Oceans 10.1029/2019JC015587 Antarctic sea ice covered ocean. Modeling studies and observations from the Arctic suggest that the feedbacks between submesoscale processes and the presence of sea ice are complex and important to understand (Lu et al., 2015;Manucharyan & Thompson, 2017;Timmermans et al., 2012). While the growth or melt of sea ice can influence the development of ocean eddies, the presence of consolidated ice may act to dampen submesoscale processes (Horvat et al., 2016;Manucharyan & Thompson, 2017;Timmermans et al., 2012). Current understanding of submesoscale processes and any seasonal variability is based on the midlatitudes with no sea ice present. It is not known whether the same seasonality will be present when ice is present, and with no models or observational studies focused on the Antarctic ice-covered ocean, this is a critical gap in our knowledge.
In this study, we show that the data from conductivity-temperature-depth (CTD) sensors carried by southern elephant seals is of sufficient temporal and spatial resolution to estimate submesoscale processes and their associated equivalent heat fluxes (section 2). A case study is utilized to assess the ability of this data set to calculate horizontal buoyancy gradients under sea ice and consequent equivalent heat fluxes driven by submesoscale processes (section 3). This analysis is then expanded to all of the available seal data in the Eastern Weddell Sea from 2008 to investigate how the magnitude and relative importance of these processes vary over the seasonal cycle (section 4). Finally, we discuss the linkages between the measured properties, sea ice concentration, and submesoscale fluxes in the Antarctic MIZ (section 5). This represents the first known study showing estimates of wintertime submesoscale processes from a sea ice covered region in the Southern Ocean.

Hydrographic Data from Seal Tags
Seal data are acquired from the Marine Mammals Exploring the Oceans Pole to Pole (MEOP) database (Roquet et al., 2013) and has been processed and quality controlled following the protocols described in Siegelman, Roquet, et al. (2019). The data are recorded on CTD sensors attached to the seals as they ascend to ensure near-vertical profiles. The tags were configured to send back 17 data points per dive, using a "broken stick" algorithm, which captures the major inflection points of the profile. The southern elephant seal data set provides a good source of data due to the relatively frequent dive times (multiple dives per day) and dive depths to approximately 500 m, encompassing the deeper wintertime mixed layer, which is typically 80-160 min this region (Pellichero et al., 2017). The data are processed at MEOP, where the measurements are compared to known climatology and ARGO float data colocated with the tag, and vertically interpolated to 1-m bins (Siegelman, Roquet, et al., 2019). The data are not extrapolated to the surface, and all profiles are cut at 6-m depth to ensure a constant depth of the shallowest value. Using this method of validation, the accuracy of the temperature and salinity measurements are estimated as ±0.04 • C and ±0.03 psu, respectively. Nineteen elephant seals were tagged on Bouvet Island in January 2008 (Lowther et al., 2014), and in this study, we use the data sets from 14 of these seals (five seals experienced a loss of the tag or immediately traveled out of the study area). This time series extends between 20 January and 15 November 2008. We focus on this smaller data set of seal tags to allow us to assess their usage for identifying submesoscale processes and to avoid possible variability introduced by including a wider temporal or spatial range. The smaller data set also allows individual assessment of each seal's transit and essential time series variables before they are evaluated statistically as a group. In this study, we will present data using TEOS-10 (McDougall & Barker, 2011) to derive absolute salinity (S A ; g kg −1 ) and conservative temperature (Θ; • C).
The majority of seals transited south toward the continental shelf (Figure 1), before remaining in the sea ice or heading back to Bouvet Island. To remove the effects of transits across the ACC fronts, or the influence of continental shelf water masses and processes (including ice shelf water), we select seal profiles located south of the southern ACC boundary, defined using sea surface height contours (Swart et al., 2010), and north of the 2,000-m bathymetry contour ( Figure 1). This encompasses the seasonally sea ice covered region, which represents the MIZ used in this study. While this still represents a large spatial area (Figure 1), we find that the latitudinal variation has minimal effect on our results, with monthly averages of calculated submesoscale variables only varying by 2% when the latitudinal range is restricted. From the seals within our study area, the modal temporal resolution of profiles (seal dives) was 6 hr 15 min, and horizontal spatial resolution of 9 km, with 50% of dives less than 15.25 km or 8 hr 20 min spacing ( Figure 2) (Boccaletti et al., 2007). We use the full-depth Rossby radius of deformation (Rd) as the upper threshold of dive spacing for this study. For this The dots show each individual profile. Those in light gray are excluded from this study because they lie with the ACC domain or directly next to ice shelves. The seal used as an individual case study (S47) is indicated in blue. Contours show the Southern ACC front (pale pink) and ACC southern boundary (pale green) using SSH contours derived by (Swart et al., 2010). The green diamond indicates Bouvet Island, where the seals were instrumented. data set, the mean Rd is 8.13 ± 2.24 km, and so we use 10.37 km (mean Rd plus one standard deviation) as the cutoff for profiles used for submesoscale estimates. This spatial resolution represents the submesoscale flows directly, in addition to the background buoyancy gradient that is the source of the potential energy to develop the submesoscale flows. By sampling at slightly larger scales, we will be underestimating the magnitude of the submesoscale fronts. This is shown by Swart et al. (2020), where the magnitude of the submesoscale gradients decrease significantly with increasing the sampling length-scale, especially near the mixed layer deformation radius. The mixed layer radius estimated deformation estimated from the seal profiles in this region is 3.7 km (Boccaletti et al., 2007), and we find that 17% of the dives used in this study were separated by this distance or less ( Figure 2c) and so are truly submesoscale resolving. The temporal spacing of the dives resolves the scales required for submesoscales (typically less than 24 hr), with all dives separated by less than 24 hr and only 6% of the profiles separated by time intervals greater than 12 hr ( Figure 2c).
The MLD is defined following the density difference criteria of Δ = 0.03 kg m −3 from a reference depth of 10 m (de Boyer Montégut et al., 2004;Dong et al., 2008).

Satellite and Atmospheric Data
Satellite observations of sea ice concentrations and sea ice drift are obtained from the National Snow and Ice Data Center SSMI Climate Record data set (Meier et al., 2017;Tschudi et al., 2016). The data are on a 25 × 25 km grid and smoothed into daily products. Each seal profile is colocated with a sea ice concentration and sea ice drift.
National Center for Environmental Prediction Climate Forecast System Reanalysis (NCEP-CFSR) (Saha et al., 2013) data are used for both meteorological data and for the surface (5 m) ocean currents. The meteorological data include heat flux components (longwave, shortwave, latent, and sensible) and wind data (including 10-m wind speed and wind stress). These data are on hourly temporal resolution and 0.5 • spatial resolution. These are regridded onto the same grid as sea ice concentration and colocated in time and space with the seal profiles.

Calculated Ocean Stress
Wind stress describes the momentum given to the ocean by the shear stress associated with wind blowing across the ocean surface. This is an important term for estimating submesoscale parameters, yet the data presented in this study occur in the sea ice covered ocean where the relationship between wind stress and its input into the ocean is different compared to open waters. We use the term "ocean stress" to describe how the presence of sea ice may affect this momentum exchange, taking into account the momentum transfer from either (or both) air-ocean stress and ice-ocean stress (equation 1).
This results in ocean stress ( ) being calculated as Sea ice concentration is represented as A, ice-ocean stress as io , and atmosphere-ocean stress as ao . The NCEP-CFSR wind-stress product is used for ao , which uses seawater density ( w ), the air-ocean drag coefficient (C d,ao ), and wind speed (U a ). The ice-ocean stress, io , is calculated using seawater density ( w = 1,027 kg m −3 ), the ice-ocean drag coefficient (C d,io ; 1.62 × 10 −3 ; Martinson & Wamser, 1990), sea ice drift (U ice ), and upper ocean velocity (U w ).
We follow the method described in Kim et al. (2017), which uses the assumption that ocean current is negligible compared to sea ice drift (U w = 0), modeling the sea ice as being in free drift. The ice-ocean drag coefficient (C d,io ) used in this study is that of the Martinson and Wamser (1990) study in the Weddell Sea. The drag for Antarctic sea ice, especially in the outer MIZ, will be smaller than that measured in Arctic studies due to the typically thinner, smoother ice characteristics (Tsamados et al., 2014).
Recent studies have also used the Martinson and Wamser (1990) approximation of io as one third of the air-ice stress ( ai ) calculated using wind speed (e.g., Wilson et al., 2019), but while this approximation largely agrees with that in equation 2, we find that this estimates the opposite sign for ocean stress for approximately 4% of dives. The absolute values of ocean stress from ai , io , and io including ocean currents (from NCEP-CFSR reanalysis) differ by less than 7%, and our usage of io with negligible ocean currents represents the smallest of these estimates. We believe that this is the most suitable io approximation to use for this region.

10.1029/2019JC015587
Overall, the use of ocean stress rather than the unmodified wind stress values decreases the stress applied to the ocean when sea ice is present by 58% across all seal profiles. This is also represented as a decrease in the wind-driven submesoscale flow of approximately 50%.

Submesoscale Buoyancy Fluxes
Submesocale processes occur when horizontal buoyancy gradients are larger than the vertical buoyancy gradient, resulting in mixed layer eddies, which act to shoal the mixed layer through an overturning streamfunction (Boccaletti et al., 2007;Fox-Kemper et al., 2008;Mahadevan & Tandon, 2006). A greater horizontal buoyancy gradient will provide a larger energy source, leading to a stronger overturning circulation and restratification (Boccaletti et al., 2007;Mahadevan et al., 2010). Wind-driven flows associated with Ekman transport across the lateral buoyancy gradient can result in additional shoaling or mixing (Thomas, 2005). A prevailing down-front wind will act to transport less buoyant over more buoyant waters (resulting in mixing), while an up-front wind will cause a restratification of the water column as more buoyant waters flow over the less buoyant waters.
The combination of these two types of submesoscale flows can result in deepening or shoaling of the MLD, but only when the vertical buoyancy flux they cause is greater than the combined buoyancy flux from the net surface heat flux, freshwater flux (precipitation, evaporation, and brine rejection), and wind forcing. We calculate the vertical buoyancy fluxes associated with submesoscale flows using MLD, along-front ocean stress, vertical stratification (Brunt-Väisälä frequency, N 2 ), and horizontal buoyancy gradient (M 2 ). These vertical buoyancy fluxes can be quantified as an equivalent heat flux in order to compare directly to surface heat fluxes (in W m −2 ) and referred to as a mixed layer eddy flux (MLE) and Ekman buoyancy flux (EBF) throughout.
The vertical buoyancy gradient, or vertical stratification, is used to indicate the stability of the water column; a low N 2 value indicates a susceptibility to mixed layer deepening. As this study is focused on how the water column evolves over a seasonal cycle, we use the average stratification for each individual profile between the surface and the depth of the winter mixed layer isopycnal ( where g is gravitational acceleration, is the mixed layer density, and 0 is the reference density of 1,027 kg m −3 . The horizontal buoyancy gradient (M 2 ) is then calculated as the difference between mixed layer buoyancy on consecutive profiles and normalized by distance between the profiles.
The balanced Richardson number (Ri b ) is a measure of the horizontal buoyancy gradient against the vertical stratification, and as it approaches unity, or is smaller than 1, it can be used as a criterion to indicate that submesoscale processes may occur (Thomas et al., 2013). It is calculated as where f is the coriolis parameter. In this study, we show the percentage of profiles in each season that meet the criterion of Ri b ≤ 1.
The equivalent heat flux associated with MLEs (Q MLE ) is calculated as a function of the MLD and M 2 (Fox-Kemper et al., 2008;Mahadevan et al., 2012): where H is the MLD, C p is the specific heat capacity, is the thermal heat capacity, g is gravitational acceleration, and is the mixed layer density. C E is a coefficient determined through numerical models, and we follow Fox-Kemper et al. (2008) to use 0.06. The vertical structure function, (z), is used to parameterize the variation in the rate of restratification through the mixed layer, where the strongest flux is found in the center of the mixed layer ( (z) = 1), reducing toward the ocean surface and base of the mixed layer ( (z) = 0).

10.1029/2019JC015587
We set (z) to 1 to represent the maximum Q MLE . The values selected for (z) and C E are the most appropriate available but may not optimally represent the environment found under sea ice or in the MIZ. Q MLE will always act as a restratifying flux (always positive), promoting stratification.
The flux caused by the wind-driven flows, Q EBF , is calculated as the balance between the surface ocean stress ( y aligned along-front) with the horizontal buoyancy flux (M 2 ), From this calculation, negative values will indicate a down-front wind (mixing) and positive values an up-front wind (restratification).

Observational Bias
To capture the full magnitude of the horizontal buoyancy gradient, the mixed layer buoyancy must be measured perpendicular to the direction of the front. However, without the ability to direct the seal profiling direction and with no observational knowledge of the mixed layer frontal orientation, we are likely to not be sampling perpendicular to the fronts. Previous studies have shown that straight line transects across features (such as sea ice leads or fronts) are likely to represent a good statistical estimate of the distribution (strength and orientation) when the sample size is large (Key, 1993). Comparable estimates of submesoscale fluxes using ocean glider data have shown that the irregular sampling of buoyancy gradients associated with unknown frontal orientations results in a statistical underestimation of the real horizontal gradient in the range of 51%, 64% or a factor of 1 As we do not have the measured frontal directions, we instead use the assumption that they are oriented perpendicular to a north-south horizontal buoyancy gradient. While the submesoscale fronts are most likely oriented at multiple arbitrary directions (Siegelman, 2020), we make an assumption about orientation in order to provide estimates of the scale of M 2 and associated submesoscale fluxes. The assumption we make is supported by our analysis of the mesoscale sea surface temperature and salinity from the MIMOC climatology (Schmidtko et al., 2013), in addition to documented ocean velocities and transport in this section of the Weddell gyre (Cisewski et al., 2011;Ryan et al., 2016;Schröder & Fahrbach, 1999). Due to this, the calculated EBF may be an overestimate or underestimate of the true vertical buoyancy flux (Du Plessis et al., 2019; Thompson et al., 2016). When all possible orientations of profile direction, frontal direction, and wind stress direction are considered, the ratio of the observed vertical buoyancy gradient to the true vertical buoyancy gradient is typically 0.71 .
This indicates that our results assuming a north-south buoyancy gradient and using zonal ocean stress ( x ) will represent a statistical sampling of horizontal buoyancy gradients and EBF and are most likely to be underestimates of the true submesoscale fluxes.

Regional Setting
The study area in the eastern Weddell Sea-Dronning-Maud Land region (Figure 1) sees one of the largest variations in sea ice edge location around Antarctica (Parkinson & Cavalieri, 2012). The maximum sea ice extent closely correlates with the southern boundary of the ACC in winter (Figure 1; approximately 55 • S) and retreats to approximately 70 • S or the edge of the continent during the melt period in the summer. South of 60 • , the satellite observations of sea ice concentration show that this area is typically sea ice covered between June and November, although in recent years, the Maud Rise polynya has reappeared, thereby significantly impacting the sea ice extent and concentration of this sector of the Southern Ocean (Jena et al., 2019;Swart et al., 2018).
In this region, four water masses are commonly found; Weddell Deep Water (WDW), Circumpolar Deep Water (CDW), Winter Water (WW), and Antarctic Surface Waters (AASW; Figure 3). In the summer, AASW forms through solar radiation and sea ice melt resulting in warming and freshening with typical surface temperatures of around 0 • C. By April, these surface waters are cooling and becoming more saline through heat loss to the atmosphere and brine rejection from sea ice formation ( Figure 3) and eventually lie on the freezing temperature line. This is the formation process of the WW layer, and this water type dominates at the surface for the rest of the data record from this data set. Typically, as the surface waters begin to freshen and warm in the spring, the AASW will be formed at the surface, with WW lying beneath, which is eventually eroded through the summer until the profiles seen for January and February are formed (Figure 3). Below these two water types are the WDW and CDW, warmer and more saline water masses (Figure 3), with a strong halocline that forms between WW and CDW. However, this study focuses on mixed layer processes and stratification, and so the denser water masses of CDW and WDW are treated as bottom boundary water masses to the mixed layer.

Seal S47: A Case Study for Estimating Submesoscale Fluxes
As satellite-retrieved seal tag data have not been used to calculate submesoscale parameters before, we focus on the data from one seal, S47, to test our hypothesis that it is possible to gather information from this source given the available temporal and spatial resolution of sampling achieved by the seals (Figure 2). S47 was tagged on 4 February 2008 and collected measurements until 23 September 2008, resulting in 721 profiles. Of these, 387 profiles lie within our area of study, excluding those where S47 traveled south to the Fimble Ice Shelf between mid-April to end of May (Figure 4, red overhead bar). Further, there are 179 profiles with less than 10.37-km separation that are used in the submesoscale analysis. Of these profiles, 22% are less than the mixed layer radius of deformation. We show the data from outside our area of study for reference.

Seasonal Cycle of Stratification and MLD in the MIZ
The time series from S47 shows the changes in the water column between February and September, with AASW, WW, and CDW all visible throughout the time period (Figure 4). The warm, saline water seen typically between 200 and 400 m lies on the mixing line between WW and CDW, and we only see purer CDW in short intervals in June and July, when S47 is near the continental shelf. The AASW is quickly eroded by the end of March, becoming cooler and more saline as winter approaches. This is associated with a decrease in vertical stratification reflected in the deepening of the MLD (Figures 4a and 4b). Using the colocated sea ice concentration for S47, it shows that the sea ice begins to grow on 28 March and is at a peak concentration between May and July, typically between 85% and 93% (Figure 4a). The sea ice concentration remains above 34% for the rest of the record, except for a short period of 2 days in July when concentrations fell below 20%.
The MLD shows little variability at the beginning of the record (late summer, between February and mid-March). It has a mean depth of 42 m with a standard deviation of only 7.5 m and is associated with strong vertical stratification (4.66 × 10 −5 s −2 ). Following the formation of sea ice, the MLD deepens to an average of 116 ± 18 m between July and September. The MLD is more variable in the winter months, frequently varying up to 40 m within a period of a few days. During this period of sea ice cover, the weakest vertical stratification is observed, with an average of 1.25 × 10 −5 s −2 . These higher frequency changes in the MLD are associated with changes in mixed layer temperature and salinity, particularly visible during August (Figure 4). The slight warming and freshening seen in the mixed layer happens simultaneously with a shoaling of the ML, indicating they could be associated with changes seen in sea ice concentration. The short time scale that they occur on is potentially indicative of submesoscale processes, and it is important to understand what role they may play in upper ocean variability in the MIZ.

Submesoscale Processes Revealed by Seal Tag Data
In order to assess the possible impact of submesoscale processes on this upper ocean variability, we calculate the parameters described in section 2.3 using the seal hydrographic data, together with using the ocean stress calculated from sea ice drift and net heat fluxes from reanalysis data (Figure 5a). The net heat flux is calculated as the sum of longwave and shortwave radiation plus latent and sensible heat (section 2.2), where negative values are out of the ocean. The effect of wind stress on the ocean is greatly dampened by the presence of sea ice, with mean values of ocean stress an order of magnitude smaller once sea ice begins to grow in April (0.18 N m −2 before April to 0.03 N m −2 between April and July; Figure 5a). In the period between July and August, we see peaks in the absolute ocean stress values of up to 0.7 N m −2 , with values then decreasing from mid-August into September to peaks of 0.35 N m −2 . The net heat flux is initially positive but switches to negative (out of ocean) on 11 March and remains negative for the rest of the S47 time series. While there is a clear diurnal oscillation initially, as the polar night approaches, the changes in net heat flux are mostly driven through longwave, sensible, and latent heat fluxes. A significantly colder period occurs between July and August, where values drop below −300 W m −2 , which is associated with the lower sea ice The horizontal buoyancy gradient is fairly constant in magnitude between February and April, with a mean absolute M 2 of 2.62 ± 2.9 × 10 −8 s −2 (Figure 5b). During the period of sea ice formation and establishment, between April and July (excluding the out of region dives), the mean M 2 is 5.29 ± 5.12 × 10 −8 s −2 , indicating increases in both the lateral surface-layer density gradients and the variability of them as the sea ice grows. Over the peak winter period, between August and September, we observe the lowest absolute values of M 2 , with a mean of 1.43 ± 1.51 × 10 −8 s −2 (Figure 5b), although there is a prolonged period between 12 and 18 August where values are consistently elevated above 1 × 10 −8 s −2 . The larger values of M 2 also occur as the MLD is deepening and these properties combine to produce an increase in estimates of MLE fluxes and EBF (Figure 5c). Some of the largest values of calculated MLE and EBF (April-May; Q MLE,EBF > 1,500 W m −2 ) occur when S47 is in front of the ice shelves and likely relate to the input of energetic ice shelf meltwater to the water column (Garabato et al., 2017).
The changes in M 2 between February-April and July-September are not reflected in the magnitude of EBF with mean values of 139 ± 179 W m −2 before April and mean values of 158 ± 460 W m −2 at the end of the record ( Figure 5). The larger standard deviation of the EBF in the July-September period highlights the variability in the magnitude of the EBF and also correlate to stronger ocean stress events, especially seen in mid-August (Figures 5a and 5c). The strong dampening effect of sea ice on ocean stress is reflected Journal of Geophysical Research: Oceans 10.1029 in the low values for EBF (|Q EBF | ≈ 100 W m −2 ) between May and July. In comparison, MLE fluxes are not affected by sea ice cover, instead driven by the increases in horizontal buoyancy gradient and MLD associated with the brine rejection during sea ice formation. Due to this, the MLE fluxes are initially low (O(10)) and begin to increase in April as the MLD starts to deepen and sea ice begins to form (Figure 5c). Although the highest M 2 values occur between April and July, the combination of M 2 and MLD used to calculate MLE fluxes means that the mean values of Q MLE are consistent from April through to the end of the record (approximately 31 W m −2 ). However, both MLD and M 2 vary significantly between August and the end of the data record, and this is reflected in increased standard deviations of M 2 , MLD, and Q MLE in this time period. These variations are reflected in some of the largest values of MLE fluxes from S47, with values reaching over 600 W m −2 .

Event Scale Submesoscale Restratification
Throughout the record from S47, abrupt changes can be seen in the MLD and associated mixed layer properties (Θ, S A , N 2 27.6 , and M 2 ), which can be caused by submesoscale processes causing shoaling or deepening of the mixed layer (Figure 4). In particular, between August and September, the MLD oscillates between about 70 and 130 m (Figures 4a and 4b), which can be associated with variability in the horizontal buoyancy gradient, changes in ocean stress, and resulting variability in submesoscale fluxes.
We focus on an event that occurs between 12 and 18 August (black dots, Figure 5, Figure S1 in the supporting information). During this time, the sea ice concentration is variable with a decreasing trend from 78% to 42% ( Figure 5; Figure S1). The cause of the change in sea ice concentration is not directly observed by this data set, although some positive (into the ocean) heat fluxes are shown during this time from reanalysis data (Figure 5a), which may have induced melting, or it may be due to spatial heterogeneity in sea ice cover and leads opening up between ice floes. Over these 6 days, the ocean stress increases in magnitude (from 0.01 to 0.44 N m −2 ), which could be linked to the decrease in sea ice concentration allowing a greater amount of wind work on the ocean (Figure 5a). Similarly, increases in the magnitude of the horizontal buoyancy gradient may be linked to these changes in sea ice concentration affecting upper ocean salinity and temperature through freshwater fluxes and increased interaction with the atmosphere (Figure 5b).
The changes observed in sea ice concentration, M 2 , and ocean stress will have an effect on the upper ocean submesoscale fluxes and consequently MLD. Initially, there is a deep MLD (130 m), and this, combined with the increases in M 2 , results in the MLE fluxes growing from near 0 to 160 W m −2 . With the decrease in sea ice cover, the EBF fluxes also increase and are consistently positive over the 6-day period with a mean value of 260 W m −2 . When the net heat fluxes, MLE fluxes, and EBF are combined, the equivalent heat flux is net positive over the 6-day period, resulting in a shoaling of the MLD from 130 to 105 m ( Figure S1). Following this restratification event, the sea ice cover increases again, and the net combined equivalent heat flux drops to near 0 W m −2 , with the MLD deepening once more.
This case study highlights how the different upper ocean characteristics can interact to result in shoaling and deepening of the mixed layer over short time scales. In particular, this variability occurred when the sea ice concentration was varying between 40% and 70% coverage, suggesting that the increased sea ice drift associated with lower concentrations can have a significant impact on submesoscale activity.

Seasonality of the Upper Ocean and Submesoscale Processes
Using S47 as a case study demonstrates that it is possible to estimate potential submesoscale fluxes from seal tag data and provides the opportunity to assess the impacts these processes have on the upper ocean. We split the data from S47 into three seasons, summer, autumn, and winter, using the net heat flux (Q net , from reanalysis data) and sea ice concentration for the average location of all seal profiles to be consistent with later analysis (Figure 6; 59.36 • S, 7.35 • E). The seasonal boundaries are defined by the switch from positive to negative net heat flux (summer to autumn, 1 April) and the onset of sea ice growth (autumn to winter, 26 June). The difference in date for these seasonal boundaries varies by less than a week if using the colocated net heat fluxes for S47.

Seasonality using Case Study Seal, S47
To characterize upper ocean and associated submesoscale fluxes in each of the seasons, we group all profiles within each date range described above and assess their distributions using frequency histograms (Figure 7). During the summer, the profiles showed a stratified upper ocean with shallow MLs (mean MLD of 43 m and N 2 27.6 of 4.33 × 10 −5 s −2 ; Figures 7a and 7b), which became deeper through autumn to winter but not necessarily less stratified. The profiles recorded in autumn have the largest range in MLD (38 to 118 m, mean 68 m), but the MLD is on average deeper in winter (117 m). In both autumn and winter, the ocean column is weakly stratified, with a slight increase in the stratification from autumn to winter (1.27 × 10 −5 s −2 to 1.3 × 10 −5 s −2 ; Figure 7a).
The typical values for the horizontal buoyancy gradient increase from summer to autumn, before decreasing in winter, from a mean of 2.39 × 10 −8 s −2 in summer to 4.91 × 10 −8 s −2 in autumn and just 1.36 × 10 −8 s −2 in winter (Figure 7c). However, the ocean stress values are modified, as expected, by the presence of sea ice with lower absolute mean values observed in autumn and winter than in summer (0.08 N m −2 compared with 0.12 N m −2 ; Figure 7d). By combining the seasonality in N 2 and M 2 , we see significant increases in winter of the number of profiles that meet the criterion for submesoscale flows, with 4 times as many profiles in winter having a Richardson number less than or equal to 1 (1.1 % in summer, 4.1 % in winter).
These variations in horizontal buoyancy gradient, MLD, and ocean stress will affect the calculations of the submesoscale parameters. Due to the deepening of the MLD and increase in M 2 , the mean MLE fluxes are larger in winter than in summer, increasing from 9 to 41.2 W m −2 (Figure 7e). The seasonality is not as stark for EBF, due to the significant moderation of the ocean stress due to the presence of sea ice (Figures 5a  and 7d), culminating in summer values of 165.2 W m −2 , which stay almost constant with 164.8 W m −2 in the winter (Figure 7f). EBF values are not further suppressed during periods of sea ice cover and its effect on surface stress due to the maintenance or sometimes enhancement of buoyancy gradients during the winter ( Figure 5).

Seasonality Shown in the Full Seal Data Set
We use the same seasonal boundaries for the full seal data set as those used for the S47 analysis. When all seals are included in the analysis, the time series spans from 19 January to 17 November, resulting in an additional season, spring, being included. The boundary between winter and spring is defined by the switch from negative to positive heat flux, which occurs on 17 October. In total, there are 1,217 profiles in the data set, with summer, autumn, and winter all containing over 340 profiles (Figure 6), providing a robust data set to perform the analysis. Due to the process of tag loss as the seals begin to molt, the spring season only contains 29 profiles ( Figure 6). While the spring profiles are still included in this analysis, the limitation due to the paucity of data during this time period must be considered.
The patterns of change through the seasons observed in the S47 data are evident in the full data set too. From summer through to winter, there is a deepening MLD (from 59 to 117 m) and decrease in the stratification (3.73 × 10 −5 s −2 to 1.04 × 10 −5 s −2 ; Figures 8a and 8b). Similarly, the mean horizontal buoyancy gradient decreases from summer through to winter (2.34 × 10 −8 s −2 to 1.81 × 10 −8 s −2 ), while the ocean stress is modified by the presence of sea ice (Figures 8c and 8d). However, when the M 2 is combined with other factors such as the MLD, the MLE fluxes increase between summer and winter (from a mean of 21.5 ± 93 W m −2 to 83.9 ± 351 W m −2 ). Due to the moderation from sea ice, EBF only has a small increase (|Q EBF | = 197.1 ± 450 W m −2 in summer, |Q EBF | = 230.7 ± 576 W m −2 in autumn and winter). The mean values seen across the full data set are comparable to the seasonal values calculated for S47, except for Q MLE , which is now at peak values during winter (Q MLE > 1,500 W m −2 in September). Also comparable to S47 is the distinct increase in the number of profiles that meet the criterion for submesoscale flows in winter. In summer and autumn, 4.7% of profiles have Ri b ≤ 1, but this doubles in winter to nearly 9% of profiles, linked to both the slight increase in M 2 and the significant decrease in N 2 .
Although there are fewer spring values, they provide an insight into how the upper ocean changes after the deep, weakly stratified winter. In all of the parameters, the spring characteristics follow the opposite trend to the summer-winter trend: a shoaling of the MLD (105 m), decrease in horizontal buoyancy gradient (0.75 × 10 −8 s −2 ), and decreases in the submesoscale fluxes (|Q EBF | = 43.7 W m −2 and Q MLE = 5.1 W m −2 ; Figure 8). The stratification of the upper ocean increases slightly from the winter into the spring (1.04 × 10 −5 to 1.19 × 10 −5 s −2 ). The changes seen in these parameters through the seasons among a total of 1,217 profiles indicate that there is a seasonal cycle in horizontal buoyancy gradient and associated MLE fluxes. Due to the dampening introduced by sea ice on the ocean stress values, the same cycle is not as distinct in EBF values, but interestingly, the main factors determining EBF, M 2 , and wind/ocean stress covary during the seasons and so maintain EBF from summer to autumn to winter.

Feasibility of Seal Tag Data for Resolving Submesoscale Flows
The availability of hydrographic profiles from seal tag data in and at the edge of the sea ice covered domains of the Southern Ocean provide a valuable means by which to assess upper ocean processes, including the impacts of submesoscale flows and variability and their impact on ocean stratification and MLD. The time-space resolution of the seal data set is appropriate to estimate submesoscale gradients and equivalent heat fluxes, especially over seasonal periods where little to no other data exist. The current array of under ice Argo floats, for example, has a repeat cycle rate of 10 days, which does not allow us to calculate submesoscale fluxes, which are known to evolve in the order of hours to a few days and at spatial scales of 10.1029/2019JC015587 O(0.1-10 km). The values for horizontal buoyancy gradient and submesoscale fluxes that are consequently calculated from the seal tags are comparable to contemporaneous studies. A recent study from the Kerguelen Plateau that also used seal tag data reported horizontal buoyancy gradients with peaks up to 2 × 10 −7 s −2 (Siegelman, O'Toole, et al., 2019), and the maximum horizontal buoyancy gradient calculated from this data set is 2.69 × 10 −7 s −2 , despite only 17% of dives fully resolving the submesoscale (Figure 2c). This lower spatial resolution of dives will lead to an underestimation of M 2 (Swart et al., 2020), which also leads to low percentages of profiles meeting the Richardson number criterion for submesoscale flows. Other studies using high-resolution ocean glider data from spring-summer deployments in the Drake Passage and sub-Antarctic regions show similar magnitudes for the contemporaneous time of year (July to March;du Plessis et al., 2017, 2019 andDecember to March;Viglione et al., 2018). The values reported in this study are within the same magnitude as these three concomitant studies, providing confidence that the submesoscale parameters reported using seal data in this study are feasible.
One of the limitations of the satellite-recovered seal tag data compared to other recent submesoscale data sets (such as ocean glider data or recovery of high-resolution seal tag data) is that there are no data on current direction and therefore the orientation of the fronts. This is important for calculating the correct EBF, as the advection of the upper ocean through EBF is driven by the wind stress perpendicular to the front orientation. In this study, we have had to assume an east-west frontal orientation (largely agreeing with the mesoscale buoyancy gradients of this region) and therefore use the zonal wind or ocean stress to calculate the EBF. While these assumptions introduce some uncertainties, they have been shown to provide underestimates of the magnitude of the fluxes (Du Plessis et al., 2019;Thompson et al., 2016). Therefore, for this study, we believe that the EBF calculated are informative to understand the potential impact of such processes to the seasonal cycle of MLD and stratification characteristics in the surface ocean. Future studies should endeavor to obtain data that provide information on the orientation of fine-scale ocean fronts so as to provide improved understanding of the wind-front interaction and its impact on the upper ocean stratification. Gliders, for example, could be a means to provide this observational data given that their estimated surface and depth average currents provide information of the frontal orientation and strength simultaneously (Du Plessis et al., 2019).

Parameterization of the Effect of Sea Ice
In this study, we include data collected by the seals from profiles underneath sea ice. The existing parameterizations for estimating submesoscale fluxes (Fox-Kemper et al., 2008;Mahadevan et al., 2012) are formulated for the open ocean where there is no additional interface between the ocean and atmosphere. We modify the EBF calculation to include ocean stress rather than wind stress, which compensates for the momentum transfer from the atmosphere to the sea ice and then into the ocean by using the sea ice drift. This reduces the EBF by approximately 50%, which indicates that the presence of sea ice significantly dampens these wind-front interactions at the submesoscale. However, it is important to note that there is still considerable uncertainty in the parameterization of ice-ocean stress (Smith et al., 2019;Swart et al., 2019), making this component of the EBF calculations a further source of possible errors. Further observations of small-scale sea ice and upper ocean velocities would be required in the Southern Ocean to improve this parameterization.
However, M 2 was shown to increase from autumn to winter and, together with deep MLDs, results in increases in the MLE fluxes. The parameterization used to estimate MLE fluxes does not contain any adjustments for the presence of sea ice. As the base of the sea ice will present additional drag to the ocean underneath, it is possible that the rate of restratification at the ocean surface will be affected, which will change the shape of (equation 4) through the mixed layer. The magnitude and significance of this effect should be the focus of future research.

Seasonality in Upper Ocean Properties
A recent modeling study suggested that there may not be a large seasonal cycle in the magnitude of submesoscale processes in the Southern Ocean (Su et al., 2018) and instead the submesoscale activity is sustained through deep mixed layers, strong frontal gradients, and wind stress that exist there. While observational studies have shown that large submesoscale fluxes can occur during the summer associated with bathymetry and standing meanders Viglione et al., 2018), no observations exists to determine the presence and magnitude of submesoscale fluxes in the sea ice zone, MIZ, or during midwinter. Using the seal tag data, we can identify seasonal characteristics of MLD, horizontal buoyancy gradients, and submesoscale fluxes between February and November, in a region typically covered by sea ice for at least 4 months of the year.
The seal tag data show a decrease in lateral buoyancy gradients from summer to autumn, followed by a small increase in the winter months. When this is combined with the deepening MLDs, the largest MLE fluxes occur in winter, with the mean values 3 times larger than those in summer. The EBF is near constant, associated with the dampening of wind stress through the presence of sea ice. The changes in the magnitude of horizontal buoyancy gradients and the dampening of EBF both demonstrate how the growth of sea ice plays an important role in the activity of submesoscale processes. While small, there is a clear seasonal signal in the restratifying fluxes from mixed layer eddies.
These seasonal characteristics can be explored further by breaking down the data into monthly means ( Figure 9). Here we focus on the absolute horizontal buoyancy gradients and the magnitude of the submesoscale fluxes. By looking at individual months, it is clear to see the second increase in horizontal buoyancy gradients in midwinter (August-September, Figure 9a) following the initial peak in gradients in February-March. In midwinter, sea ice is present and still variable (e.g., Figure 6), potentially introducing lateral gradients brought on by increased rates of brine rejection. Similarly, leads (gaps in between sea ice floes) may be prolific during this period, creating lateral heterogeneity in ocean-atmosphere interaction. The decrease between April and May in the horizontal buoyancy gradients occurs following the switch to negative heat flux (out of the ocean) in April, resulting in significant, homogeneous heat loss from the surface ocean. In October and November, the magnitude of lateral buoyancy gradients decreases again, possibly due to a more concentrated ice pack minimizing the effect of leads or atmosphere-ocean interaction, reflected in the lack of variability in sea ice concentration in Figure 6.
The standard deviations of these parameters provide an indication of the variability within the data for each month (Figure 9b). When the mean values of lateral buoyancy gradient are high, the standard deviations also increase, indicating that the larger magnitudes of lateral buoyancy gradient occur as short-term events, rather than sustained increased values. However, June shows a different pattern, with comparable mean values of M 2 to May and July but increased standard deviation. The onset of sea ice formation typically occurs in June (Figure 6), and so this increase in standard deviation is likely due to heterogeneous sea ice formation across the study area.

10.1029/2019JC015587
Submesoscale fluxes are calculated using the combination of horizontal buoyancy gradients and MLD or ocean stress, and so the monthly means look similar to those of the horizontal buoyancy gradients (Figure 9a). The effect of deepening MLs and increases in M 2 is clearly visible in the seasonal cycle of Q MLE , with fluxes growing until peak values are reached in August and September (Figure 9a), which are up to 5 times greater than MLE fluxes in the summer and autumn. The dampening effect of sea ice is reflected in the decreasing trend in |Q EBF | values from autumn into winter, with peak values typically constant from February through to August. The largest magnitudes of the combined submesoscale flux (|Q EBF | + Q MLE ) occur in August, driven by the secondary peak in M 2 and significant increase in Q MLE . August and September also have the highest variability in submesoscale fluxes (Figure 9b), indicating frequent short-term events. Including all data between February and November indicates that there may be a seasonal cycle in combined submesoscale fluxes in the sea ice covered ocean of 240 W m −2 .
A remaining question is what happens to the lateral buoyancy gradients and submesoscale fluxes between November and February, where seal tag data is typically missing. There is currently insufficient data to resolve submesoscale processes in this time period, and this should be a focus of future studies in order to capture the sea ice melt season. New field campaigns in the region, using surface and underwater gliders, are beginning to resolve the submesoscale energetics during this period (Swart et al., 2020).

Conclusion
This study presents the first known estimates of horizontal buoyancy gradients and submesoscale fluxes from the ice-covered Southern Ocean through the use of seal tagged data. While the seal-derived hydrography is somewhat limited by the spatial and temporal resolution, it presents the best resolution of wintertime data that currently exists and is of the appropriate resolution to estimate such surface-ocean buoyancy fluxes. The results reveal the dominant seasonal characteristics of submesoscale processes at play in the upper ocean and highlight these important points: 1. Horizontal buoyancy gradients show two peaks throughout the year, at the end of summer and at midwinter. The peaks are associated with the initial onset of cooling and, later, heterogeneity in concentration of the sea ice pack. The decrease in horizontal buoyancy gradients in autumn is likely associated with the dampening of heterogeneity in the upper ocean due to large heat loss between April and July. 2. The second peak in horizontal buoyancy gradients and deep mixed layers result in a seasonal cycle of 240 W m −2 in the combined submesoscale heat fluxes. This is larger than suggested in previous modeling studies (Manucharyan & Thompson, 2017;Su et al., 2018) and highlights that the role of submesoscale flows under sea ice should not be ignored when modeling heat and carbon atmosphere-ocean fluxes. 3. The magnitude of submesoscale fluxes can be large enough to result in restratification of MLD in midwinter.
Due to the limitations and uncertainties associated with the resolution, sea ice parameterizations, and frontal orientations, this study can be considered to represent a statistical representation of the seasonal variability in submesoscale fluxes. With further assessment of the full circumpolar seal tag data set and collection of fine-scale observations (including high-resolution data from retrieved seal tags) between November and February, combined with information of front directions, we will be able to resolve and understand these under sea ice upper ocean processes with the aim of further parameterization into climate models.