The evolution of a shallow front in the Arctic marginal ice zone

Arctic sea ice is rapidly changing, with summer sea ice extent and thickness declining more quickly than winter (Perovich et al., 2018). These changes are particularly pronounced in the Chukchi and Beaufort seas of the Western Arctic which have increased open-water in the summer while maintaining winter ice cover, and thus have much higher seasonal ice variability than in the past (Thomson et al., 2016; Stroeve and Notz, 2018; Peng et al., 2018). Accompanying these changes is an increase in width of the summer marginal ice zone (MIZ; Strong and Rigor, 2013): the band of low concentration, highly dynamic ice that separates the open ocean from the ice pack. This increase has led to a renewed interest in understanding the dynamic processes at play in the Beaufort Sea MIZ, with an emphasis on ice-ocean-atmosphere coupling and feedbacks (Lee et al., 2012).

Strong lateral buoyancy gradients in the upper ocean of the MIZ have been identified as having a potentially important role in the heat transport in the ice-ocean-atmosphere system as a result of ice edge upwelling (e.g., Niebauer, 1982; Røed and O’Brien 1983; Häkkinen, 1986) or eddy generation (e.g., Lu et al., 2015). Historically, many of the observational studies of MIZ fronts have taken place in the shallow regions of the Bering Sea (Niebauer, 1982) or the energetic Fram Strait (Johannessen et al., 1983; 1987; Smith and Bird, 1991) where topographic effects may have a significant influence on ocean currents. As the marginal ice zone shifts north into the deeper basins of the Arctic Ocean, these fronts may begin to resemble open-ocean frontal systems. The lateral gradients in upper ocean density observed across the MIZ can have horizontal length scales on the order of 1 km to 10 km (e.g. Buckley et al., 1979). At these scales, submesoscale processes with Rossby and Richardson numbers of order one may dominate the momentum balance (Thomas et al., 2008).

The role of submesoscale dynamics in the MIZ is a subject of recent and active research. Numerical simulations of idealized MIZ have shown spontaneous generation of mixed-layer eddy fields due to submesoscale baroclinic instabilities (SBI) energized by either an existing frontal structure (Manucharyan and Thompson, 2017), or by freshwater input from ice melt (Horvat et al., 2016). Observations and simulations of a submesoscale meltwater front in the Chukchi Sea suggest similar SBI (Lu et al., 2015). These simulations are corroborated by satellite observations showing that MIZ eddy fields are common at scales consistent with SBI (e.g., Bondevik, 2011; Kozlov et al., 2019), and a drifter study has provided in-situ evidence for MIZ submesoscale lateral dispersion (which would be associated with mixed layer eddies (MLE); Mensa et al., 2018). In contrast, high-resolution observations of a submesoscale filament in the Fram Strait instead suggest an alternative instability mechanism as the dominant submesoscale process (von Appen et al., 2018). Timmermans et al. (2012) suggests that the presence of ice-ocean friction may act to modify typical submesoscale processes. Few direct oceanographic measurements of the evolution of MIZ submesoscale fronts have yet been collected, so how the dynamics and forcing are impacted by presence of ice is unclear. Furthermore, while slumping and re-stratification of Arctic submesoscale fronts have been suggested (Timmermans et al., 2012; Timmermans and Winsor, 2013; Wulff et al., 2016), these studies have relied on a limited number of cross-front sections which are unable to resolve any temporal variation.

This paper presents detailed observations of a localized ice-edge front, and addresses its observed evolution in the context of existing submesoscale theory. Following a description of the observations and methods, the evolution of the front is discussed in terms of the dynamical processes that could explain the observed frontogenesis and frontolysis. While these observations with high resolution in both time and space combined with existing submesoscale theory and scaling arguments provide intriguing insight into dynamics of ice-edge fronts, conclusions highlight the need for additional studies.

A ship-based survey took place aboard Norseman II along the edge of the MIZ in the Canada Basin of the Beaufort Sea (Figure 1a). Data were collected over a ~3-day period from 01–04 Oct 2014, so the timing of the survey corresponded to the start of the fall ice-growth season, when the ocean was rapidly cooling and sea ice was beginning to form (the 2014 minimum sea ice extent occurred on 17 Sept 2014; NSIDC Fetterer et al., 2017). In the vicinity of the study site the basin is approximately 3800 m deep and the nearest shelf break is roughly 200 km away, so topographic influence on the local submesoscale dynamics is unlikely. Shipboard measurements were complemented by a number of short deployments (3 to 6 hours each) of two SWIFT drifting buoys (Thomson, 2012). Zippel and Thomson (2016) previously presented data collected during part of this survey in a study about MIZ surface wave dynamics. The experiment was part of the larger Office of Naval Research Marginal Ice Zone program (Lee et al., 2012; Lee and Thomson, 2017).

Oceanographic measurements were collected shipboard using a Teledyne Oceanscience underway profiling system with an underway conductivity-temperature-depth (uCTD) sensor operated on a winch from the stern of the ship (e.g., Rudnick and Klinke, 2007). The uCTD was operated in “yo-yo” mode while the ship was underway (without recovery between casts), thus multiple profiles were executed between instrument recoveries (when probes were swapped and data downloaded). This approach allows for increased temporal and spatial resolution compared to traditional CTD operations; with the ship travelling at 4 kn (~2 ms^–1), profiles collected with a ~4-min period extend to depths of 120–150 m and have ~380 m along-track spacing. The temperature and conductivity sensors both react to variability in the signal on different time frames, so the delay between these becomes a function of fall speed. We corrected for non-uniform fall speed using a technique similar to Ullman and Hebert (2014) and binned the results onto a regular 1-m depth grid. The profile locations were obtained from the navigation files, interpolated to the start time of each cast. Additional post-process quality control was done manually to remove profiles that showed contamination from icing or salinity spiking. The top 4 m of each profile were discarded due to potential ship contamination (the ship has a 4-m draft); however, the strong bulk lateral gradients observed here exist at scales that would be influenced minimally by ship contamination.

Profiles taken with the uCTD showed a strong two-layer structure in the upper ocean (Figure 2). In the open water, a relatively warm and fresh surface layer (T ~ 1.5°C, S_A ~ 27.5 g kg^–1) extended to a depth of 20–25 m, where a sharp pycnocline separated it from colder saltier water below (T ~ –1.5°C, S_A ~ 30 g kg^–1). Transects showed a distinct front of cold and even more fresh water (T ~ –1.5°C, S_A ~ 26.5 g kg^–1) located at the ice edge (Figure 2). Within the front, water temperature approached the in-situ temperature of freezing based on measured salinity, suggesting that the front may be residual meltwater from the sea ice melt season. In the temperature-salinity range observed here, density is driven by changes in salinity (the haline contraction coefficient β ~ 7.8 × 10^–4 kg g^–1, while the thermal expansion coefficient α ~ 3 × 10^–5 K^–1), so the density structure mimics salinity. Adopting a definition of the mixed layer as the density threshold of 0.25 kg m^–3 compared to the shallowest observation (e.g., Cole et al., 2017), the mixed-layer depth was consistent with the surface layer depth (which we define by the 22.15 σ_Θ isopycnal) in open water, while near the ice the mixed layer is shallower due to the presence of the front. While we consider the separate definitions for the surface layer and the mixed layer here, in studies of submesoscale dynamics, descriptions of the “mixed layer” may be based on the locations of the pycnocline and thus would be more consistent with what we refer to as the surface layer (e.g., Fox-Kemper et al., 2008).

The MIZ in the survey area consisted of loose brash ice floes (up to approximately 1 m in diameter) interspersed with small first year floes (up to approximately 5 m in diameter; Figure 3a). Upon arriving in the survey area, a coherent ice edge feature was visually identified on the ship’s navigation radar (Figure 3b, c) and the ship followed that feature as the MIZ deformed and drifted; the ship maintained this feature on the radar during the sampling period. For a roughly 12-hour period on 03 Oct, Norseman II performed short looping patterns while maintaining a relatively constant speed of 4 kn (~2 ms^–1) and following the moving ice (Figure 1b–d). Norseman II is not an ice-capable vessel so for most of the survey the ship only approached the MIZ edge without travelling into the ice; the oceanographic features described here likely extend farther into the MIZ.

The observed MIZ drift pattern was used to identify a local coordinate system with along- and across-drift directions. An origin point was chosen close to the start of the drift. Based on uCTD measurements, this coordinate system is also aligned with the bulk horizontal gradients of temperature and salinity (Figure 1b). The along-drift coordinate axis is approximately parallel with the larger-scale ice edge. Figure 1b shows the path of ship while it was in the survey area, rotated to the local reference frame, and overlaid on a synthetic aperture radar (SAR) image taken by RADARSAT-2 (from MDA Geospatial Services Inc., provided via the U.S. National Ice Center), which shows the ice cover. The angle of the coordinate system was chosen based on approximately matching the direction of the ship track and the bulk gradients of salinity and temperature in the mixed-layer; the results are not sensitive to this angle. The coordinate system is used to define a right-handed x-y axis (Figure 1c), with x aligned in the along-drift (and therefore along-front) direction, and y aligned across-drift/front (so a positive horizontal velocity u is “down-front”, in the direction of a geostrophic frontal jet).

While the movement of regional large-scale ice pack over the survey period was small (~0.05 ms^–1) and generally eastward (NSIDC: Polar Pathfinder; Tschudi et al., 2019), ice in the MIZ was observed to drift to the west and slightly north (indicated by the arrow on Figure 1b). The ice feature being tracked was moving at approximately 0.3 ms^–1 (inferred from the mean along-drift translation of the ship during the looping track). During two deployments bookending this period, both the SWIFT buoys drifted with a speeds ranging from 0.25 ms^–1 to 0.4 ms^–1 towards a mean compass direction of 315°. This drift opposed the local wind direction which originated from 270–300° (Figures 1b and 4b). The observed ice and surface drift direction is consistent with the expected direction of a geostrophic frontal jet formed across the meltwater front. Assuming a level-of-no-motion at the base of the surface layer (z = –h_sl, defined as the 21.15σ_Θ isopycnal, where the strong halocline may decouple surface motions from below), the thermal wind balance $u(0)=f−1∫−hsl0−bydz$⁠, for buoyancy b = –gσ/ρ₀, and f the Coriolis parameter) predicts a velocity of 0.25 ms^–1 within the front for the profile shown in Figure 2. This velocity is roughly consistent with the measured drift velocities from the SWIFTs and inferred velocity from the ship track.

Profiles taken with the uCTD during the looping path were separated into a series of across-drift sections showing temporal evolution of the front in a reference frame following the identified MIZ ice feature (Figure 1c, d). Each section was composed of all casts taken within a 1-hour window, interpolated onto a regular grid in y and z (the across-drift distance and depth), and smoothed with a 2-dimensional Gaussian kernel. The windows overlapped by 30 minutes, so some individual casts may be shared across two different sections; in this way, the sections form a type of moving average and do not map directly to individual “loops” from the path (Figure 1d). Sections containing fewer than 5 casts or spanning a range of less than 1.8 km in the across-drift direction were discarded. We assume that the reference frame following the ice feature provides a good approximation for a Lagrangian frame of the front as a necessary step in describing temporal evolution. Differential ice-ocean velocities and upper-ocean shear may impact this assumption; these effects are discussed later.

An acoustic doppler current profiler (ADCP) was mounted in a well on the ship to measure the ocean currents. The ADCP reported measurements in 1-m bins with a sampling period of 30 seconds. The first bin was a depth of 7.3 m below the water surface, so details of velocity within the front were not fully resolved. A VmDas data acquisition system was used to remove ship motion; however, examination of the resulting velocity profiles suggests that some ship contamination was still present in the velocity signal, obscuring small cross-front horizontal velocity gradients. Despite these limitations, we recovered some additional contextual information from the ADCP system during the looping survey (when the ship was moving relatively slowly) through additional processing. We created smoothed, across-drift velocity sections using a method identical to that used for for creating uCTD sections. Because these sections often average over times when the ship was travelling in different directions (Figure 1c), the ship contamination was largely removed by this process.

We characterized temporal evolution of the front by considering a number of section-averaged quantities averaged within the surface layer, which we define by the 22.15σ_Θ isopycnal. We denote these averages by ⟨·⟩_sl, which represent 2-dimensional averages in y-z. This isopycnal choice for the surface layer depth is slightly shallower than the pycnocline separating the surface layer from the water below, which ensures that the sharp gradients within the pycnocline are not skewing averaged quantities. The overall results are not overly sensitive to the choice of isopycnal to describe surface layer depth provided that pycnocline is sufficiently excluded from averages. Additionally, we calculated the bulk frontal slope for each section in two ways. (1) A linear fit was made of some characteristic isopyncals within the front; we contrast results from choosing either the 21.65σ_Θ and 21.85σ_Θ isopyncals. (2) A ratio of the values of the lateral buoyancy gradient (M² = –b_y) to the vertical buoyancy gradient (N² = b_z) averaged within the mixed layer for each section also provides an approximation of isopyncal slope: ⟨M²⟩_sl[⟨N²⟩_sl]^–1 ~ Δz[Δy]^–1.

Atmospheric variables were measured locally by meteorological packages aboard Norseman II and on one of the two SWIFTs. On Norseman II, these variables included air temperature, relative humidity, barometric pressure, and wind speed and direction. Norseman II additionally carried a sonic anemometer mounted on the mast at a height of 10.3 m above the mean water level. The sonic anemometer recorded 3-component relative wind speeds and virtual temperature at 10 Hz sampling rate. These data were processed in 10-minute intervals using the inertial dissipation method described by Yelland et al. (1994) and developed by Large and Pond (1981) in order to calculate the wind stress. Within each 10-minute interval data were de-spiked using the phase-space method of Goring and Nikora (2002), and wind velocity spectra were constructed with 512 point 75% overlapping windows. A –5/3 frequency-power slope was fit within the inertial subrange (identified as 0.6–2 Hz) of each spectra. When wind was from the ship’s stern, erroneous measurements were recorded due to contamination from the ship’s exhaust, so these data were discarded. This matches the procedure used by Zippel and Thomson (2016) for the same measurement data.

We compared measured quantities with atmospheric re-analysis. In an evaluation of different re-analysis products, Lindsay et al. (2014) show that Modern-Era Retrospective analysis for Research and Applications (MERRA) produces the most consistent results with independent observations. Based on those results, we chose to use the updated “version 2” of the MERRA product (MERRA-2; Global Modeling and Assimilation Office (GMAO), 2015; Gelaro et al., 2017), which provides hourly data at 0.5° × 0.625° degree resolution. We compared the on-site measurements with the average of two grid points that fall within the survey area (located at [73°N, 145.62°W] and [73°N, 146.25°W]). Additional comparison was done with two other re-analysis products: the National Center for Climate Prediction and Department of Energy re-analysis version 2 (NCEP-2; Kanamitsu et al., 2002), and the European Centre for Medium-Range Weather Forecasts (ECMWF) interim re-analysis (ERA-Interim; Dee et al., 2011). There was good agreement between the three different re-analysis products during the region and time-frame of record (not shown).

Wind speeds measured aboard Norseman II were generally low, and varied from 0 ms^–1 to a maximum of 10.0 ms^–1 with a mean of 4.3 ms^–1 (Figure 4a) from a compass direction ranging between ~260° to 300° (westerly to northwesterly; Figure 4b). MERRA-2 re-analysis agrees in a bulk sense, showing similar westerly to northwesterly winds with means wind speeds of 3.5 ms^–1, but displays much less variability than locally measured wind. Through most of the survey, the wind direction was approximately aligned with the front, in the “up-front” direction (compare the wind direction in Figure 4 with the approximate frontal orientation indicated by the dashed lines). Shipboard wind measurements show a number of events lasting ~4–6 hours in which the local wind speeds significantly exceeded those reported by re-analysis products, with discrepancies of up to 100% (Figure 4a; peaks at 01 Oct 12:00Z, 02 Oct 13:00Z, 03 Oct 10:00Z, and 03 Oct 20:45Z). The timing of these high-wind events does not correspond to day-night transitions or changes in ship operation (heading or speed, etc.). One of these events occurred while the ship was performing the looping patterns described above and collecting high resolution oceanographic data (denoted by the shaded grey area in Figure 4).

The COARE bulk algorithms (v3.5; Fairall et al., 2003; Edson et al., 2013) were used to estimate both surface heat fluxes and wind stress using the shipboard measured atmospheric variables. Downward radiative fluxes (long- and short-wave) were not measured, so the COARE algorithms were run using MERRA-2 re-analysis values as inputs for those variables. Despite the differences between re-analysis data and observations present in the wind (Figure 4a), the choice to use re-analysis values for calculating the heat flux is justified by the lower variability and smaller error associated with downwelling radiative fluxes over a variety of re-analysis products compared to other atmospheric variables (Chaudhuri et al., 2014). Sea surface temperature input for the heat flux calculations was taken as the uppermost retained bin of the uCTD casts (at 5-m depth) interpolated onto the time series of atmospheric measurements, so heat fluxes are only available during times when uCTD casts were being taken. This approach provides a “bulk” sea surface temperature. Precipitation was not measured so it was not included in the heat flux estimate. Version 3.5 of the COARE algorithms also includes parametrization to account for sea state using the significant wave height, H_s, and phase speed, c_p, of the waves. We computed heat fluxes using the averages of surface wave properties measured by the two SWIFTs (H_s = 0.38 m, and c_p = 5.8 ms^–1 based on deep-water waves with a measured average peak period of T_p = 3.7s). Given the influence of the MIZ on sea state (e.g., Zippel and Thomson, 2016), whether these are indicative values for the use in the parametrization is unclear; however, in this case the results of application of the COARE routines to the shipboard meteorological data showed limited sensitivity to the inclusion of wave information (heat flux calculated with and without the inclusion of the sea state differed by <1 W² m^–1). The COARE-derived wind stress estimates matched very closely with those measured using the inertial dissipation method.

As expected, there was net ocean cooling throughout the record consistent with the period of ice advance (Figure 4c). Weak incoming short-wave radiation (peaks of ~50 W m^–2) on both 03 and 04 Oct was insufficient to reverse the heat flux direction during daylight hours, although MERRA-2 shows a brief period of heating (~6 h) during daylight on 01 Oct. The net heat flux given by MERRA-2 re-analysis significantly underestimated that calculated by the COARE algorithms during times that the wind speed differed between observations and re-analysis. Differences between the two measures are largely a result of differences in the sensible and latent heat flux components which are both sensitive to wind speed.

Strong gradients across the MIZ suggest that sea-ice-atmosphere interactions drive the processes in these regions. Here, the frontal adjustment appears to be linked to the atmospheric forcing, but dynamics differ from traditional descriptions associated with submesoscale frontal evolution. A number of details of this behaviour are still unresolved.

Across-drift sections show that the front rapidly evolved over an approximately 12- to 18-hour period on 03 Oct (Figure 5a–h). The front steepened over a time period of a few hours with the slope reaching a maximum at approximately 10:00Z on 03 Oct, then subsequently collapsed over a period of roughly 4 h. The qualitative pattern of steepening and collapse is robust across a broad range of isopycnals, as demonstrated by the 21.65σ_Θ and 21.85σ_Θ isopycnal slopes (Figure 5i); however, the magnitude of the change in slope depends on the specific choice of isopycnal. The bulk estimate of slope (based on the ratio of lateral and vertical buoyancy gradients) generally agrees with the slope of 21.65σ isopycnal. The initial steepening and subsequent collapse of isopyncal slope is reflected in the horizontal buoyancy gradient, ⟨M²⟩_sl (Figure 6a), so the evolution can be divided into a frontogenic phase followed by a frontolytic phase. The depth of the mixed layer, as defined by a density difference threshold compared to the shallowest density measurement, does not show similar evolution (Figure 5a–h, grey dashed-dotted line).

While the vertical buoyancy gradient, ⟨N²⟩_sl, evolved over a similar time-frame, its behaviour was more oscillatory with a maximum that appears later and a less pronounced peak compared with ⟨M²⟩_sl (Figure 6b). The overall change in isopycnal slope (Figure 5i) was driven primarily by changes in ⟨M²⟩_sl which has a high variability, while ⟨N²⟩_sl only deviates slightly from its time-mean (⟨M²⟩_sl ranges from ~ –0.3 × 10^–6 s^–2 to1 × 10^–6 s^–2; ⟨N²⟩_sl ranges from ~4 × 10^–4 s^–2 to 5.6 × 10^–4 s^–2; Figure 6a, b).

During the evolution, the along- and across-front velocities in the surface layer (⟨u⟩_sl, ⟨v⟩_sl respectively) also both showed an oscillatory character with a timescale close to the inertial period (Figure 6d, e), and seem to have varied in quadrature. In contrast, the section-averaged vertical shears in the surface layer (⟨u⟩_sl, ⟨v⟩_sl) do not have apparent oscillatory motion, and instead vary together. Through the record, the along-front shear ⟨u_z⟩_sl appears to be in approximate balance with ⟨M²⟩_sl (f^–1), consistent with an along-front thermal wind balance (Figure 6f).

Differential advection by sheared inertial oscillations can lead to oscillations of the stratification ⟨N²⟩_sl (Thomas et al., 2016). While the observed stratification does appear to exhibit such a response (sinusoidal curve fit in Figure 6b), these oscillations would theoretically arise from a stratification budget (⟨N²⟩_sl)_t ~ ⟨v_z⟩_sl⟨M²⟩_sl with constant ⟨M²⟩_sl and oscillatory ⟨v_z⟩_sl (e.g., see Figure 5 in Thomas et al., 2016). We tested if sheared inertial oscillations are setting the stratification by integrating the budget forward with the observed ⟨v_z⟩_sl and ⟨M²⟩_sl. Because ⟨v_z⟩_sl and ⟨M²⟩_sl. typically have the same sign (Figure 6f), the budget leads only to increases in stratification and cannot reproduce the observed oscillatory motion (Figure 6b, blue line). Sheared inertial oscillations are related to the geostrophic adjustment described by Tandon and Garrett (1994), in which an inviscid system that has been disturbed from thermal wind balance (by passing storms, for example) will oscillate around its balanced state. In contrast, our observations show that the front roughly maintains thermal wind balance throughout the adjustment (Figure 6f), suggesting that whatever mechanism is forcing the along-front shear ⟨u_z⟩_sl may also be simultaneous forcing the cross-front buoyancy gradient ⟨M²⟩_sl.

The observed oscillations of the average velocity fields can be characterized with an estimate of the horizontal momentum equations (i.e., a simple slab model):

where U, V are the surface layer-averaged ageostrophic velocities (U = ⟨u⟩_sl – u_bg with a background [e.g., geostrophic] flow u_bg, and V = ⟨v⟩_sl), and ζ the background vertical vorticity (where we assume no along-front gradients so ζ = –u_y), and h is the depth over which wind stress acts (e.g., the mixed-layer or surface layer depth). Numerically integrating these equations with the observed wind (over a fixed surface layer depth h = 25 m), and using u_bg and ζ as fitting parameters produces good estimates for the observed velocities (Figure 6d, e) with u_bg = 0.02 ms^–1 and ζ = 0.16f if non-zero initial conditions are used for U, V. However, here the non-zero initial conditions necessary to fit the observed response result in the wind forcing terms being of negligible size, so the velocities vary almost exactly as sinusoids with frequency $feff=f(f+ζ)$ (a period of 2π/f_eff ~ 11.1 h). Alternatively, taking the velocity as initially at rest and integrating forward with wind from the beginning of the record (on 01 Oct; Figure 4a) does not capture the phase of the observed oscillations, but can produce oscillations with comparable magnitude (~0.05 ms^–1). Choosing smaller values of h (e.g, the mixed-layer depth as opposed to the surface layer depth) increases the strength of the wind forcing term and produces less desirable fits to the observed velocity.

In describing the frontal evolution, we have assumed that measurements are made in a Lagrangian reference frame following the front. The oscillatory nature of many of the measured fields (Figure 6) suggest the possibility that this is not the case. If instead the front was moving independently from the sampling scheme, then the frontogenic/frontolytic behaviour may be simply a result of advection by an underlying spatially homogeneous inertial velocity field, such as is described by the slab model. Such independent movement would further explain the fact that variations in vertical velocity shear match the variation in lateral buoyancy (Figure 6f). However, because the ice drift opposes the wind direction (Figure 1b), it is likely moving with ocean surface currents. Given the speed of the ice drift, this scenario would require significant decoupling and shear between the near-surface currents and the bulk of the surface layer. Additionally, an explanation would be needed for what was driving those surface currents. Inertial advection of the front past a non-Lagrangian measurement platform does not directly explain the link between the wind stress increase and the increase in lateral buoyancy gradient (see “Frontogenesis” section, below). Furthermore, the development of an eddy field visible in the ice (see “Frontolysis” section, below) suggests that currents driving the ice motion extend beyond a shallow surface layer. While we cannot rule out the possibility that our observations result from simple advection of a steady front, we proceed with the assumption that the primary driver of the observations is unsteady frontal evolution.

Throughout the looping survey, none of the atmospheric variables measured aboard Norseman II (e.g., air temperature, wind stress) showed any measurable cross-front variation, indicating that the front does not drive changes in the atmosphere on $O$(1 km) scales. Nevertheless, we did observe a cross-front variation in the surface heat flux on the order of roughly 10 W m^–2 (Figure 7b) based on composites of the surface heat flux anomaly as a function of across-drift distance (calculated as the heat flux minus the section-averaged heat flux and averaged in 250-m across-drift bins). This variation is driven by the cross-front variation in bulk sea surface temperature (~0.5°C change; Figure 7a). The ocean is losing heat faster away from the ice edge, where the difference between water and air temperature is greater. As the front adjusts, it will mediate the location and magnitude of horizontal gradients in the surface heat flux.

The mean wind speed (~5 ms^–1) and duration (~5 h) of the higher-wind events (when the observed wind deviates from the re-analysis) correspond to an advective length scale of ~90 km (i.e., the length scale if the wind events represent passing storms). Such a scale would only be minimally resolvable by re-analysis products suggesting that these are localized events. Furthermore, this scale is of roughly the same order as the MIZ width, which may indicate that they are laterally constrained by atmospheric heterogeneity across the MIZ. Localized MIZ atmospheric variations associated with the temperature and boundary layer differences across ice-edge have been documented previously (e.g., Grønås and Skeie, 1999; Inoue and Hori, 2011; Guest et al., 2018). Furthermore, Wenta and Herman (2018) found that air circulation within the atmospheric boundary layer is sensitive to ice floe properties on scales typically unresolved by mesoscale weather models.

The submesoscale MLE in the MIZ predicted by models (Lu et al., 2015; Horvat et al., 2016; Manucharyan and Thompson, 2017) have a frontolytic tendency and alone are unable to explain the observed frontal steepening. The frontogenesis is exactly coincident with wind stress increase (compare Figure 6a and 6c), so steepening is driven in some way by surface forcing: either wind stress directly, or associated changes in surface heat flux that accompany the wind. Because the density is set primarily by salinity in this temperature-salinity range, heat flux alone is not expected to have a significant direct affect on the front; however, heat loss from the polar ocean has the secondary impact of causing ice growth which itself is associated with convection.

Some frontal studies evaluate the interaction of surface forcing with MLE by comparing either parametrized magnitudes of the vertical buoyancy flux, ⟨w′b′⟩, associated with each of the drivers of stratification (Mahadevan et al., 2012) or comparing the strength of overturning stream functions (Mahadevan et al., 2010). In this framework, the observed wind stress is unable to lead to isopyncal steepening as the associated buoyancy flux relies on Ekman dynamics (Thomas and Lee, 2005; Mahadevan et al., 2010) which are characterized by subinertial time frames that are longer than the observations. Furthermore, the direction of the wind (“up-front”) would also have a frontolytic sense if it persisted long enough for the Ekman dynamics to manifest.

Convection could lead to frontal steepening by mixing waters at the heal of the front. Mahadevan et al. (2010) considers convection only due to surface cooling, which is parametrized by ⟨w′b′⟩_Q = –αQg(ρc_p)^–1 for c_p the heat capacity of water, α the coefficient of thermal expansion, Q and the rate of surface heating/cooling (heat flux out of the ocean is negative). In polar regions, additional convection is possible due to brine rejection during sea ice formation (McPhee, 2008); an associated vertical buoyancy flux, analogous to ⟨w′b′⟩_Q can be constructed for surface salinity input as: ⟨w′b′⟩_S = gβρ_i(ρ₀)^–1ΔS(C_ih_i)_t, where β is the haline contraction coefficient, (C_ih_i) the “effective ice thickness” (C_i is ice concentration and h_i is the ice thickness), ρ_i the ice density, and ΔS the difference in salinity between the ocean and the sea ice. Using reasonable values for sea ice properties (ρ_i = 938 kg m^–3, S_i = 6 g kg^–1, from the Beaufort Sea MIZ in Oct 2015; Smith et al., 2018), and assuming that approximately 5 cm of ice grew during the course of the survey in a 15% ice concentration zone (by visual observation while on site), the total convective buoyancy flux is ⟨w′b′⟩_Q + ⟨w′b′⟩_S ~8.6 × 10^–9 m²s^–3. This value is an order of magnitude below the mean vertical buoyancy gradient associated with MLE: $〈w’b’〉MLE=0.06 hsl2M4f−1~ 9×10−8m2s−3$ (Fox-Kemper et al., 2008), suggesting that convection cannot explain the observed steepening.

Other studies (Dewey et al., 2017; Smith et al., 2018) show that one-dimensional mixing processes often dominate upper ocean physics in the MIZ. We test this with the use of a one-dimensional mixing model: the Price-Weller-Pinkel model (PWP; Price et al., 1986). This model accounts for the combined effects of both wind-driven mixing (due to shear) and convective mixing (due to surface cooling). We ran the model using the surface heat fluxes and wind stresses from the COARE algorithm as forcing for two different vertical profiles. The profiles were separated by a distance of roughly 2 km, so one was located within the front while the other was near the edge of the front. The model applied wind stress over the mixed-layer depth (as defined by a 0.25 kg m^–3 density threshold compared to the shallowest density measurement; as above), so the initial mixed-layer depths were roughly 8 m and 14 m for the two profiles, which is consistent with the depth of the front. In the model, momentum-driven deepening of the mixed layer (due to entrainment) is based on shear between the mixed-layer base and an underlying transition layer (based on Richardson number criteria). The model predicted no appreciable change to the depth of the front for either starting profile (not shown), so the inclusion of shear-driven-mixing (in combination with convection) is still unable to account for the observed frontogenesis. Running the model with different mixed-layer depth definitions (e.g., based on different σ_Θ thresholds, or based on thresholds of Δσ_Θ/Δz) may have resulted in slightly different initial mixed-layer depths, but produced no difference in the qualitative behaviour of the model (i.e., no mixed-layer deepening). To recreate deepening of the front comparable to observations in the model required multiplying the wind stress by a factor of 5.

Frontogenesis is classically explained as a result of horizontal strain fields on a front (Hoskins, 1982, and references therein). In the ocean, mesoscale eddy fields have been shown to provide frontogenic strain fields (e.g., Capet et al., 2008). SAR imagery shows an apparently active eddy field (see Frontolysis section, below), so that there was likely sufficient energy at the mesoscale to provide such a strain field, although our measurements do not allow further resolution of those features. However, if the frontogenic behaviour is a response to a non-local mesoscale strain, then its modulation by the local wind field over short time-scales is unlikely. Instead, we consider processes by which the wind field would induce a strain on the front.

Ignoring mixing and non-conservative terms (justified by the results from the convective/shear-driven deepening experiments), and assuming that along-front fields are uniform (so ∂_x terms are zero), the evolution of lateral buoyancy can be written as

where D_t represents the material derivative (D_t = ∂_t + v∂_y + w∂_z). If advective terms are sufficiently small (which may or may not be a reasonable assumption, as discussed above), changes in ⟨M²⟩_sl are driven by cross-front variations in the ageostrophic velocity fields. Spatially heterogeneous surface stress may drive cross-front variations in the velocity, and thus provide a mechanism by which the front directly responds to the temporally increasing wind stress.

The presence of ice modifies the stress transfer from the atmosphere to the ocean and can lead to heterogeneity in the surface stress input to the ocean. While careful representation of stress transfer across the ice-ocean-atmosphere system requires accounting separately for wind-ice and ice-ocean stresses (e.g., Steele et al., 1989; McPhee, 2008; Cole et al., 2014; 2017), some studies in the MIZ (e.g., Fennel and Johannessen, 1998) consider a single bulk wind stress, with a drag coefficient modulated by the ice. Surface stress measurements complementing that approach show heterogeneity moving from the ocean to the ice, with a local maximum of surface stress in the MIZ even for a spatially homogeneous wind field (this heterogeneity is represented by a spatially varying drag coefficient), where the increase in stress is due to the increased surface roughness elements in the MIZ (Guest et al., 1995; Birnbaum and Lüpkes, 2002; Andreas et al., 2010). Our observations do not show measurable cross-front variations in either wind speed or wind stress over $O$(1 km) scales, but the existence of surface stress variations over the MIZ scale of $O$(10 km) is probable.

Furthermore, the ice may have an impact on the depth to which wind forcing penetrates through modification of the surface turbulent boundary layer. In the slab model presented above we assumed a fixed depth h ~ 25 m over which the wind stress is applied based on the approximate open-ocean surface layer depth. However, studies of under-ice boundary layers have shown that the Ekman layer depth (the steady state depth of wind stress penetration) is less than the mixed-layer depth (McPhee, 2008; Cole et al., 2014; 2017; Randelhoff et al., 2014, and others); for example, measurements from an ice-tethered-profiler in the Beaufort Sea MIZ in Oct 2014 showed an average Ekman layer depth of approximately 10 m in Oct, compared to a mixed-layer depth of 30 m to 40 m (see Figure 10 in Cole et al., 2017).

In previous studies, non-uniform surface stress across the MIZ is associated with ice-edge upwelling/downwelling (e.g., Niebauer, 1982; Fennel and Johannessen, 1998), and eddy generation and shedding (Häkkinen, 1986). However, these studies consider longer subinertial timescales, so that Ekman dynamics can be set up (upwelling/downwelling at the ice edge is associated with convergence/divergence in the Ekman layer, and eddy shedding is a non-linear interaction with Ekman velocity). If the along-front wind stress was maintained, the system may evolve towards one of those states; however, we see highly dynamic transient behaviour that was not discussed in those studies.

Following relaxation of the wind, the front rapidly collapsed (Figure 5). Submesoscale frontolysis is driven by instabilities such as SBI, which lead to the development of eddies in the mixed layer and re-stratification. A comparison of SAR images prior to the survey and at the end of the survey show a change in character of the ice edge (Figure 8). Prior to the survey (29 Sept; Figure 8a), the ice edge was generally compact, while over the course of the survey there was considerable deformation of the ice edge. A SAR image from late on 03 Oct (Figure 8b) shows that this deformation resulted in a number of coherent filaments and vortical structures over a wide range of scales. Within the MIZ, sea ice is particularly mobile (Lund et al., 2018) and can act as tracer for surface ocean currents (e.g., Johannessen et al., 1987; Manucharyan and Thompson 2017), so this deformation is indicative of an underlying ocean eddy field.

The appearance of these eddies may be due to the development of an eddy field over this time frame (i.e., the growth of MLE; Manucharyan and Thompson 2017). Linear stability analysis (Stone, 1966) gives that the fastest growing SBI mode has a wavelength of ℓ_MLE = 2π$U$(f^–1)[(2/5)(Ri + 1)]^1/2, where $U$ is a characteristic velocity, and $Ri=N2(uz2+vz2)−1$ is the Richardson number (which is $O$(1) in submesoscale flows). The linear model is valid during the initial stages of instability, before the inverse cascade transfers energy to longer wavelengths (Fox-Kemper et al., 2008). Median values for Ri in the mixed layer for the different sections were $O$(1) (section medians varied approximately from 0.3 to 1.2, although spatially varying values varied over several orders of magnitude). Using Ri = 1, and the surface velocity as a characteristic velocity scale, $U$ ~ 0.25 ms^–1, gives the wavelength of the most unstable SBI mode as ℓ_MLE ~ 10 km. This scale is consistent with many of the observed eddies within the SAR image (Figure 8b) suggesting the possibility that MLE and associated processes are present here. This presents evidence for MLE-driven frontolysis; however, such frontolysis is also associated with restratification (Fox-Kemper et al., 2008). The apparent oscillatory behaviour of ⟨N²⟩_sl (Figure 6b) does not present strong evidence for MLE restratification. Furthermore, the speed of frontolysis may be faster than would be associated with MLE re-stratification (which might evolve over several days (Fox-Kemper et al., 2008).

Alternatively, frontolysis may be separate from the eddy activity. Instead of developing over the course of the survey, an ocean eddy field may have already been active but decoupled from the sea ice (for example, due to the influence of wind). Figure 8b shows the MIZ at a time when the wind had been effectively stopped for a period of hours (Figure 4a; at 16:20Z on 03 Oct), so it is possible that the ocean currents arranged the ice into the visible eddy field over that time. For example, Kozlov et al. (2019) observed eddies in the satellite record within the MIZ that may be associated with gyre processes. From the onset of the survey, however, the observed ice drift opposed the direction of local winds, suggesting that it was already moving along with the upper ocean. Cole et al. (2017) have suggested at these low ice concentrations, the transfer of momentum from wind to ice is relatively inefficient.

The primary contribution of this study is a set of high resolution observations of an ice-edge located front taken over a period of days, which allowed us to investigate its temporal evolution. These observations show a dynamic adjustment of the front over fast timescales that has not been described previously in the literature. The front appears to respond to a wind event that drives a frontogenic behaviour; following relaxation of the wind, the front rapidly relaxes. While important metrics are derived (e.g., Figure 6), our measurements are insufficient to fully explain what is driving this behaviour. We hypothesize that the presence of sea in the MIZ may be important factor by driving spatial heterogeneity in the surface stress and momentum transfer from the atmosphere to the ocean. Furthermore, the wind event driving the frontal adjustment may be a localized event that itself is impacted by the MIZ.

Drag from differential ice-ocean velocity and variations in surface roughness in the presence of ice have both been shown to modify wind stress transfer to the ocean (McPhee, 2008). In considering wind-ice and ice-ocean stresses separately, Cole et al. (2017) indicated that at low ice concentrations the transfer of momentum from the wind to the ice is inefficient, but also suggested that ice-ocean momentum transfer might be highest in low ice concentrations. Atmospheric measurements show an increase in surface stress moving from the open ocean to the MIZ (Guest et al., 1995; Birnbaum and Lüpkes, 2002; Andreas et al., 2010). In mixed ice and open water conditions (e.g., in the MIZ), estimates of bulk drag coefficients have been constructed based on open-ocean and ice-based coefficients weighted by ice concentration (Zippel and Thomson, 2016). However, in low-concentration MIZs, observations of the total atmosphere-ice-ocean momentum transfer are sparse, and do not yet give a full picture of all processes. For example, recent work suggests that in an energetic sea state the damping of short waves by ice in the MIZ may decrease the effective transfer velocity from wind to ocean (Smith and Thomson, 2019). Overlaid on the complicated picture of wind stress transfer to the ocean in MIZs are localized weather patterns, such as low-level jets, which may result from the transitions of atmospheric properties across the MIZ (e.g., Guest et al., 2018; Liu and Schweiger, 2019).

The use of the COARE algorithms has been shown previously as insufficient for predicting surface fluxes in the MIZ (Yu et al., 2017); however, in this case the strong match between the wind stress estimated with COARE versus that calculated using the inertial dissipation method (Figure 6c) suggests that the COARE results here are valid. The validity of the COARE algorithms here is likely because the ship track was generally travelling only on the open-water side of the ice edge. Travel into the ice may have necessitated a different approach for calculating surface fluxes (e.g., with alternative parameterizations for turbulent transfer coefficients from Andreas et al., 2010). We used the COARE-derived wind stress and heat fluxes as forcing for the PWP model. Previous application of PWP to the MIZ has used a hybrid approach where surface fluxes were derived from a three-dimensional Arctic model in order to account for ice growth/melt, and impacts from ice on surface forcing (Dewey et al., 2017). However, that study examined much longer trends (annual variations) where melt/growth processes may have been significant. For the short model runs we performed (evolving over hours), we do not expect the sea ice to modify our application of the PWP model significantly, with the exception of its role in modifying surface stress. We can test the effect of sea ice on the PWP model results through application of some scaling factor on the surface stress. For example, Guest et al. (2018) predicted a factor 2.5 times greater drag coefficient in the MIZ compared to open water. However, in our simulations a factor of 5 increase was required for any appreciable deepening of the front that would be comparable to the observations.

Lateral buoyancy fluxes due to MLE have been theorized to be a potential driving process for setting the upper ocean structure in the MIZ (Lu et al., 2015; Horvat et al., 2016; Manucharyan and Thompson, 2017). SAR imagery shows an active eddy field at the study location (Figure 8b) with a number of features with scales consistent with the fastest growing mode from linear instability theory (Stone, 1966). However, whether MLE processes are dominant here is not yet clear. The interaction of surface forcing (and in particular wind stress) with MLE restratification has not yet been investigated for the MIZ. In the open ocean, wind aligned in an “up-front” direction leads to frontolytic behaviour (Thomas and Lee, 2005), as opposed to the frontogenesis observed here. However, that effect is driven by Ekman buoyancy flux, which is driven necessarily by Ekman dynamics that need a longer time-scale to develop. For open-ocean fronts the competing effects of re-stratification from MLE and de-stratification from convection has been investigated previously based on the relative strength of parametrized vertical buoyancy flux ⟨w′b′⟩ associated with MLE and surface cooling. We have extended that approach by developing an equivalent convective ⟨w′b′⟩ to account for brine rejection during surface cooling.

The apparent response of an ice-edge meltwater front to variations in localized wind forcing emphasizes the need to study the combined ice-ocean-atmosphere system rather than considering these components individually. This study highlights the challenges inherent to studying these processes in the MIZ. Shipboard sampling of these meltwater fronts is particularly difficult due to the relatively small horizontal and vertical spatial scales. Sampling near-surface stratification and velocity from ships with large drafts relative to the depth of the upper ocean mixed layer is particularly challenging. Wulff et al. (2016) have demonstrated that autonomous underwater vehicle platforms are capable of capturing near-surface profiles within a meltwater front with a high enough horizontal resolution to properly characterize a number of frontal properties and without the inherent problems associated with ship wake contamination. However, whether such vehicles can sample with an appropriate temporal resolution to capture the type of dynamic adjustment we have observed here is not clear. Furthermore, these platforms are unable to measure local wind patterns and thus must rely on re-analysis products, which we have shown to miss local ice-edge variations.

Episodic ice melt during the fall ice-growth season due to strong winds has been observed previously in the MIZ due to one-dimensional mixing processes (Smith et al., 2018). However, the impact of wind on mixing is strongly mediated by underlying oceanographic conditions; as we have seen here, one-dimensional mixing may not be valid in the presence of a strong frontal structure in the MIZ. Modelling studies indicate that sea ice may be particularly sensitive to vertical and lateral heat fluxes associated with MLE (Fox-Kemper et al., 2011; Horvat et al., 2016), though the present study is unable to answer questions about the impact of MIZ frontal processes on ice growth/melt.

While the specifics of the frontal steepening are not fully resolved, the described contrast with open-ocean conditions indicates that different frontogenic mechanisms may be of importance in the MIZ and provides a first step in understanding the forced dynamics of submesoscale fronts in polar regions. As the sea ice edge retreats farther north over the deep basin of the Beaufort Sea during the summer ice minimum, these processes may become increasingly important in controlling the heat fluxes during the fall freeze-up period.

Data from the Ocean Science underway profiling system (uCTD) along with shipboard meteorological and ADCP data are available at http://hdl.handle.net/1773/45389. Data from SWIFTs are available at apl.uw.edu/swift.

We would like to thank the captain and crew of Norseman II and our colleagues in the Marginal Ice Zone program. Additionally, we would like to acknowledge support from the Office of Naval Research (ONR) for funding and the National Ice Centre (NIC) for providing satellite imagery. RADARSAT-2 Data and Products ©MDA Geospatial Services Inc. (2014) – All Rights Reserved. RADARSAT is an official mark of the Canadian Space Agency.

Funding for the MIZ program was from the Office of Naval Research (ONR) grants N000141210180 and N000141210113. Additional funding is through ONR Early Student Support (322).

The authors have no competing interests to declare.

Samuel Brenner prepared the manuscript and performed the primary analysis. Luc Rainville and Jim Thomson collected the data and performed initial processing. They also contributed tools and ideas to the analysis, assisted with manuscript preparation, and provided feedback and supervision throughout. Craig Lee provided equipment for data collection and assisted with experiment design. He also contributed to data analysis and interpretation.