Influence of seasonally varying sea-ice concentration and subsurface ocean heat on sea-ice thickness and sea-ice seasonality for a ‘warm-shelf’ region in Antarctica

Abstract

Processes driving changes in sea-ice seasonality and sea-ice thickness were explored for a ‘warm-shelf’ region along the West Antarctic Peninsula using vertically coupled sea-ice-ocean thermodynamic simulations, with and without assimilated satellite sea-ice observations and moored ocean temperature observations. Simulations with assimilated sea-ice observations permitted investigation of surface [thermodynamic and dynamic (e.g., wind-driven)] processes affecting sea-ice thickness and seasonality. Assimilation of quasi-weekly variability in the depth and temperature of the deep warm pycnocline permitted examination of subsurface processes affecting sea-ice. Simulations using assimilated sea-ice observations (and implied motion) always produced greater surface heat fluxes and overall thinner sea ice. Assimilating seasonal and quasi-weekly variability in the depth and temperature of the pycnocline modified the start of the sea-ice season by −23 to +1 d, and also modified the sea ice thickness/seasonality to be thinner/shorter or thicker/longer at sub-seasonal and seasonal timescales, highlighting a mechanism where a shoaling pycnocline enhanced upward deep-water heat fluxes as transient surface-induced turbulence had a greater effect on a reduced mixed layer volume. The observed interplay of surface, subsurface, and sea-ice modulation of ocean-atmosphere heat transfer underscores the importance of representing the interaction between sea-ice concentration and upper ocean variability in climate projections.

1. Introduction

In polar oceans, heat fluxes from both the atmosphere and subsurface ocean thermodynamically control the growth and melt of sea ice (e.g., Fig. 1). Atmospheric forcing, dominated by shortwave irradiance, grossly follows the seasonal cycle and is the primary driver of sea-ice formation in autumn/winter and melt in spring/summer (Maykut, 1986; Persson and others, 2002). In most of the Antarctic seasonal sea-ice zone seaward of the continental shelf, the potential exists for an upward flux of heat from warm Circumpolar Deep Water (CDW). This deep-ocean heat flux is typically smaller than atmospheric fluxes, yet nevertheless has a large impact on sea-ice production and sea-ice thickness evolution (Martinson, 1990; Martinson and Iannuzzi, 1998; Ackley and others, 2015) and consequently on the seasonality of the sea-ice cover and attendant biogeochemical processes.

Figure 1. Conceptual diagram showing the seasonal progression of heat and salt fluxes and related feedback processes that control sea-ice formation and melt. The vertical axis is not to scale but captures the features and fluxes of a ‘warm’ continental shelf area (e.g., Fig. 2). Red denotes relatively warm waters; blue represents waters at or near the freezing point. For scale, the thermodynamic sea ice thickness typically reaches a maximum of ~70 cm (for a warm shelf region) and snow can accumulate up to ~15 cm. The seasonally warmed and freshened summer mixed layer depth can be 20 to 50 m (dotted black line, far left and far right) while the bottom of the winter water layer can be at 80 to >150 m depth (solid black line; Venables and Meredith, 2014), where the permanent pycnocline marks the transition to warm and salty Circumpolar Deep Water (CDW) further below.

On Antarctic continental shelves, CDW-derived heat can also be an important regulator of sea-ice production, especially for ‘warm-shelf’ regions (e.g., the Bellingshausen and Amundsen shelf seas) where CDW at depth is 1–2°C warmer than for ‘cold-shelf’ regions (e.g., the Weddell and Ross shelf seas) (e.g., Petty and others, 2013; Frew and others, 2019). At the coast on these same warm-shelf regions, CDW-derived heat can also regulate ice shelf thickness (Holland and others, 2020) and the melt of marine glacier termini (Cook and others, 2016) through enhanced melt rates. In contrast, on cool-shelf regions, where CDW is not present or its heat content is greatly attenuated, sea-ice production is typically higher and sea ice thicker (Mahoney and others, 2011). Further, supercooled water emerging from beneath ice shelves and marine glacier termini on cool-shelf regions can contribute to enhanced freezing rates and thicker thermodynamically-grown sea ice (Eicken and Lange, 1989; Hellmer, 2004; Smith and others, 2012; Langhorne and others, 2015).

The seasonal progression of cooling, autumn sea-ice growth, winter sea-ice thickness evolution, and spring sea-ice melt, along with attendant ocean heat fluxes, are depicted conceptually for a warm-shelf region (Fig. 1). In autumn, the cooling atmosphere extracts heat from the seasonally warmed and freshened mixed layer until it reaches the freezing point, at which point sea ice forms. Brine rejection from sea-ice formation mixes surface waters downward, converting the seasonal mixed layer into cold salty Winter Water (WW) that eventually merges with the remnant WW layer below. On warm-shelf regions, warm CDW is present below the WW layer, and continued sea-ice production causes entrainment of CDW heat into the WW layer. This heat delivered to the WW layer and the overlying sea ice regulates sea-ice growth, representing a negative feedback on sea-ice thickness evolution (Martinson, 1990). Some entrained ocean heat might also ventilate to the atmosphere either through open water areas (if sea-ice concentration is less than 100%) and/or through the snow and sea ice cover. In spring, increasing shortwave radiation warms the surface, and if sea-ice concentration is less than 100%, surface warming causes relatively rapid sea-ice melt (relative to areas with 100% sea-ice concentration). This occurs because the relatively dark open ocean has a lower albedo than highly-reflective snow-laden ice. The warmed surface ocean then causes both basal and lateral sea-ice melt, which creates more open water, and in turn more surface warming causing more sea-ice melt (a strong positive feedback). Once sea ice is melted, the resulting fresh water, combined with continued solar heating, creates a warm, fresh and shallow seasonal mixed layer (shallow relative to the remnant WW layer below).

Dynamically, sea-ice growth rates, sea-ice thickness evolution, and sea-ice seasonality are also affected by wind and wave forcing (e.g., Kohout and others, 2014). In the presence of wind and waves, pancake ice forms from the aggregation of frazil ice into disc-like ‘pancakes’. Pancake ice can grow two to three times faster than sea ice under quiescent conditions (Doble and others, 2003). Wind-driven sea-ice movements can also mechanically thicken the sea-ice cover (by rafting and ridging of ice floes), decreasing its sensitivity to further wind and wave forcing and increasing the time needed to subsequently melt-out (e.g., Stammerjohn and others, 2003). In autumn, on-ice (from ocean to ice) winds can lead to delayed ice-edge expansion, or in spring, early ice-edge retreat (and conversely for off-ice winds). Persistent winds can also generate latent heat polynyas, which export large amounts of sea ice in winter while maintaining an area of open water through advection (Smith and others, 1990). When the spring ice-edge retreat is exceptionally early, the amount of ocean heat gained by solar insolation during the expanded ice-free summer can be sufficient to delay sea-ice formation the following autumn (e.g., Stammerjohn and others, 2012). If, on the other hand, the winds are divergent over the pack ice, thus opening the pack ice from within, sea ice can melt faster due to solar-warmed open ocean areas, thus allowing for a faster spring melt-back (Watkins and Simmonds, 1999).

The seasonally-covered ocean just west of the West Antarctic Peninsula (WAP) (Fig. 2) is a warm-shelf region and is one of the fastest warming places on Earth (Meredith and King, 2005; Ducklow and others, 2012; Schmidtko and others, 2014; Stammerjohn and Scambos, 2020) despite also exhibiting strong natural variability (Turner and others, 2016). Coincident changes in winds and sea-ice distribution have led to large ecosystem shifts, including decreased primary productivity and southward shifts in the range of ice-obligate species (Ducklow and others, 2007, 2012; Schofield and others, 2017). Sea ice in particular has changed dramatically along the WAP, with the sea-ice season decreasing from an average of 240 to 180 d over the course of 42 years of satellite observations (1979–2020). To the southwest (in the southern Bellingshausen Sea), the sea-ice season has decreased even more sharply, from an average of 310 to 230 d (Stammerjohn and Maksym, 2017, updated).

Figure 2. Map of the simulation study area showing the Palmer LTER sampling grid, and the mooring location where deep ocean mooring observations and atmospheric forcing reanalyses (ECMWF Interim; Berrisford and others, 2011) were used to perform coupled sea-ice-ocean model experiments. Coastlines and ice contours are derived from the SCAR Antarctic Digital Database (SCAR, 2020) available through the QAntarctica (NPI, 2020) package, and the bathymetry was produced by the Global Multi-Resolution Topography synthesis (Ryan and others, 2009). Currents are reproduced from McKee and Martinson (2020b).

Along the wind-exposed WAP, daily changes in ice-edge location and in sea-ice concentration are largely wind-driven (Stammerjohn and others, 2003; Massom and others, 2006; Massom and others, 2008; Stammerjohn and others, 2008b; Holland and Kwok, 2012). In this scenario, additional open water created by advection of sea ice out of the region can enhance either the ice-albedo feedback in spring by allowing more sunlight to penetrate and warm ice-free waters or enhance ocean heat fluxes in autumn-winter by exposing relatively warm surface waters to rapidly cooling air temperatures and deep wind mixing (Meredith and King, 2005; Stammerjohn and others, 2012; Holland, 2014). The converse of these heat transfer processes is also possible, with a wind-driven influx of sea ice to the region.

Along the WAP, ocean observations also show that localized eddy-driven or regional upwelling of warm CDW on the continental shelf has driven changes in ocean heat content (Martinson and others, 2008; Dinniman and others, 2012; Martinson and McKee, 2012; Graham and others, 2016). The potential influence of subsurface ocean heat on sea ice may be particularly strong on the WAP continental shelf because of vigorous shelf-slope exchange processes (Martinson and McKee, 2012; McKee and Martinson, 2020b). From mooring and other hydrographic observations (e.g., CTD and gliders), it has been shown that CDW is ‘injected’ into mid-depths from the adjacent slope with relatively unattenuated properties compared to other Antarctic shelves, where the thermocline is broader and the temperature maximum is cooler and found at near bottom shelf depths (Martinson and McKee, 2012; McKee and others, 2019). More generally on the WAP, other than by entrainment during brine rejection, cross-pycnocline mixing of CDW heat is believed to occur primarily through shear-driven instability (as opposed to double-diffusive convection; Howard and others, 2004; Brearley and others, 2017). Vertical heat fluxes out of the CDW layer have been estimated by combining observed temperature gradients with diffusivities estimated from parameterizations (using observations of shear and stratification; Howard and others, 2004; Brearley and others, 2017; McKee and others, 2019) or, more recently, from microstructure measurements of dissipation (Scott and others, 2021; Inall and others, 2022), and background values range from less than 1 W m⁻² to about 2 W m⁻² with large spreads, indicating mixing is intermittent. As for the origin of the shear, Howard and others (2004) determined that near-inertial waves were the primary source whereas Brearley and others (2017) showed the origin varies depending on the sea ice cover. With sea ice absent, counterclockwise-polarized near-inertial shear is elevated with the rotary coefficient correlated to local wind stress, together implying episodic wind-events (storminess) drive the vertical heat flux. When sea ice is present, wind input is damped so that inertial waves are suppressed and tides, though weak, are the main source of shear. Other important mixing processes that drive vertical heat movement are tied to bathymetry, including mixing associated with overflows at sills by which colder pycnocline waters are entrained into the CDW (Venables and others, 2017), and hydraulic control at ridges by which persistent, larger heat fluxes (~5 Wm⁻²) can supply CDW heat to the pycnocline after which episodic wind-driven mixing can tap into it (Scott and others, 2021).

The additional heat supplied to surface waters from below is presumed to impact sea-ice presence and thickness in this region of declining sea ice (see Venables and Meredith (2014) for an observational perspective on the depth of CDW and multi-season impacts on subsequent sea-ice formation), however the timing and magnitude of heat delivery from below and above on sea-ice mass balance have not been fully resolved. To investigate the role of subsurface ocean heat on WAP seasonal sea-ice dynamics, we take advantage of year-round mooring measurements of ocean temperature at a 400 m deep mid-shelf location along the WAP (Fig. 2) conducted by the Palmer Long Term Ecological Research (LTER) program. These mooring observations include discrete temperature measurements from near bottom to ~50 m from the surface that show vertical displacements of warm thermocline waters of up to 100 m over the course of a few days (Martinson and McKee, 2012). These ocean temperature measurements are assimilated into a vertically-coupled sea ice and ocean model (described below) to investigate how seasonal and higher-frequency variability of the presence and depth of warm CDW impacts the overlying sea-ice cover. Further, by variably incorporating satellite observations of daily sea-ice concentration, the model is used to explore the relative contributions of surface thermodynamic plus dynamic (e.g. wind-driven) processes and subsurface processes in the seasonal evolution of sea-ice thickness and sea-ice seasonality for this Antarctic warm-shelf region.

2. Methods

2.1. Model

Our focus in this study is to explore heat transfer from the deep ocean and its potential impact on mixed layer dynamics and on sea-ice seasonality and mass balance (e.g., Fig. 1). A regional 3-dimensional (3D) ocean simulation with attendant regional forcing fields could diagnose these coupled processes if the relevant features are accurately resolved (in this case, sub-mesoscale eddies interacting with complex shelf break bathymetry) or observations exist for development of an assimilation system. In the WAP, the only year-round sub-surface observations of CDW are from point-source moorings which can be used to infer horizontal processes, but do not provide a regional-scale view. With these high-resolution vertical observations in mind (and with the added goal of creating a modeling framework that is more accessible than 3D regional models that require significant computing and technical resources), we developed a 1-dimensional (1D) coupled ocean-ice-ecosystem model which explores the thermodynamic coupling between the surface and deep ocean and its role in the seasonal evolution of the upper-ocean marine environment. The KPP-Eco-Ice (KEI) vertical water column model introduced here integrates the KPP mixing model (Large and others, 1994; Doney, 1996), a heat- and mass-conservative sea-ice model (Saenz and Arrigo, 2012, 2014), and a marine ecosystem model (Moore and others, 2002, 2004), and is written in FORTRAN 90. The ecosystem dynamics of the model will be addressed in a separate manuscript. Horizontal dynamical variability is introduced through assimilation of satellite sea-ice concentration and subsurface mooring observation for warm CDW waters, which permits examination of the thermodynamic and dynamic contributions to surface ocean state via perturbation experiments.

The physics of the upper-ocean boundary layer and subsurface are simulated numerically with the K-profile parameterization (KPP) nonlocal, boundary layer scheme developed by Large and others (1994) and further described by Doney (1996), and both of these studies provide the key equations upon which the boundary layer model is built. Here we describe the model physics most relevant for assessing the delivery of subsurface ocean heat to the planetary boundary layer (PBL) and overlying sea-ice cover. The model estimates the time evolution of horizontal velocity components (u and v), potential temperature, and salinity as a function of depth from the conservation equations for momentum, heat, and mass (Doney, 1996, Eqns (10–13)). Vertical turbulent fluxes in the PBL (that decrease to much lower values in the stratified interior) are computed using depth-dependent finite eddy diffusivities times the vertical property gradient and a nonlocal or counter-gradient term (Doney, 1996, Eqn (14)), which is zero except for temperature under unstable conditions (Large and others, 1994). The shape of the boundary layer diffusivity profile is specified as a function of the distance from the surface and scales with the PBL depth and a turbulent velocity scale that increases with surface wind stress and unstable surface buoyancy forcing (i.e., net cooling/evaporation, brine-induced mixing from sea-ice production) (Doney, 1996, Eqn (15)). Vertical diffusion below the PBL in the stratified ocean interior is generated by a local gradient Richardson number instability parameterization computed from the buoyancy frequency and vertical shear (Large and others, Eqns (27 and 28)) and a small background internal wave diffusivity (as opposed to the stability criterion applied during subsurface restoring described below in Section 2.4).

The model does not a priori define a well-mixed layer near the surface, and the PBL depth is taken as the penetration depth of surface generated turbulence. The PBL depth is determined diagnostically based on a stability criterion for the bulk Richardson number and depends upon the stratification, velocity shear, and surface buoyancy forcing (Doney, 1996, Eqn (16)). A bulk Richardson number relative to the surface is calculated from the buoyancy and velocity vector profiles, and the PBL depth is taken as the shallowest depth where the bulk Richardson number equals a specified critical value. An unresolved shear term generated from surface turbulence is added to the resolved vertical shear to give the correct entrainment flux for the pure convective case.

2.2. Model grid, boundaries, and time stepping

We used a 1 m resolution grid for the water column in the KEI model, extending to 400 m depth, which is well below seasonal mixing (down to ~80–120 m) and the main pycnocline (below ~80–120 m) (e.g., Fig. S1). The bottom ocean boundary was fixed for all tracers, including salinity and temperature, at model initialization. The air-ocean boundary exchanges heat after Doney (1996) and Doney and others (1998), utilizing 2 m air temperature and humidity, wind speed, shortwave and longwave radiation, and precipitation from atmospheric reanalyses (for details, see next section). The surface ocean is assumed to have a constant shortwave albedo of 0.07.

The sea-ice model utilized between 0–26 sea-ice layers, and between 0–14 snow layers, with layer thickness varying in an accordion-fashion after Saenz and Arrigo (2014). The minimum sea-ice or snow layer thickness was 2 cm. Sea ice less than 2 cm thickness, which only occurred during ice melt, was immediately transferred to the liquid state by cooling the mixed layer. The sea-ice surface heat flux follows the method employed in CICE v4 (Hunke and Lipscomb, 2010). In contrast to Saenz and Arrigo (2012), shortwave irradiance is delivered directly to the surface ice or snow layer, and the shortwave albedo of snow and sea ice is fixed at 0.8. This albedo value is at the high-end of measurements and model predictions for snow-covered sea-ice, reflecting the high fraction of snow cover on Antarctic sea ice (Briegleb and Light, 2007).

The KEI model used a base time step of 1 h, which the ocean model and KPP mixing scheme operated at under all circumstances. The sea-ice model advective and fractional coverage calculations were also made hourly, and sea-ice heat and fluid transfer operated with an adaptive time step to maintain numerical stability with a maximum time step of 0.2 h (Saenz and Arrigo, 2014).

2.3. Surface atmospheric and sea-ice forcing

Atmospheric forcing data are derived for each simulation from the ECMWF Interim Reanalysis (ERA Interim), as this dataset reasonably simulates precipitation and thus snow accumulation on Antarctic sea ice (Turner and others, 1999; Nicolaus and others, 2006) but does not account for the redistribution of snow by winds (Leonard and Maksym, 2011). Humidity, temperature, wind speed, and surface pressure were interpolated linearly from 6-hourly values selected from the ERA Interim gridcell nearest to the simulation location (see Section 2.6). Precipitation and short- and long-wave irradiances were calculated from cumulative fluxes at 3-hourly intervals and were assumed constant across each interval. The more recent ERA5 reanalysis was not available when this model was originally constructed and calibrated. Sea-ice concentration at the location of the Palmer LTER mooring was extracted from sea-ice concentration data released by the National Snow Ice Data Center (NSIDC), using the NASA Team algorithm and calculated from passive microwave brightness temperatures from the Special Sensor Microwave Imager/Sounder (SSMI/S) instruments (DiGirolamo and others, 2022).

2.4. Ocean forcing

Moored ocean data come from the moorings maintained by the Palmer LTER project (Martinson and McKee, 2012; Martinson, 2020). For the site and years used in this study, each mooring was equipped with between 11 and 17 fixed-depth thermistors whose depths in the vertical were chosen to yield an integrated ocean heat content that showed an average bias under 0.1% (Martinson and McKee, 2012). Most thermistors were SBE39 sensors from Sea Bird Electronics, with an initial accuracy of 0.002°C and a resolution of 0.0001°C. For the 2011 deployment, 2 of the 17 temperature series were measured by the thermistor on a current meter (Alec Infinity EM). The moorings do not measure salinity, excepting a single sensor during 2007 and 2008. Further mooring sensor details are provided in supplementary Table S1.

Mooring temperature observations were interpolated vertically to 1 m resolution from the thermistor data using a polynomial regression (cubic hermite interpolating polynomial; sensor depths vary in time due to knock-down under strong currents), and then sub-sampled in time (from 10 or 15-minute observations) to match the base-hourly time step of the ocean model to create a continuous subsurface ocean temperature profile (Fig. S1; see Martinson and McKee, 2012). Salinity S _w was restored over a similar depth profile as temperature (see below) using a piecewise linear regression from temperature T _w (Fig. S2) developed from springtime (1993–2012) CTD casts of modified CDW found in the WAP region, and from a mooring-deployed temperature-salinity sensor at ~280 m depth (2007–08) from the Palmer LTER program location 300.100 (Fig. 2):

(1)$$\eqalign{S_w = \;& 0.233\;T_w + 34.338, \;-1.50 < T_w < 0.602, \;\cr S_w = \;& 0.187\;T_w + 34.366, \;\;0.602 \le T_w < 2.00} $$

The linear fit for warmer T_w was generated from the mooring T-S paired data, representative of deeper, less-diluted CDW. The linear fit for cooler T_w is from CTD cast T-S paired data from between 60–200 m depths and captures the greater range of salinities found in WW; T_w = 0.602 is where those fits meet. This is a simplified treatment of salinity, which does not account for known variability in CDW properties, particularly eddy-driven CDW variability discussed further below (McKee and others, 2019). Nonetheless our results show that this T-S parameterization is more than sufficient for distinguishing the relative effects between assimilated vs nonassimilated ocean temperature on sea-ice mass balance, thus serving to demonstrate the importance of subsurface heat on ice-ocean interactions.

During designated simulations, assimilation of the deep ocean profile was accomplished by relaxing temperature values toward the mooring data with a 24-hours restoring period. Longer relaxation periods failed to replicate the full range of observed subsurface temperature variability. Relaxation occurred from either 60 m or the depth of the temperature minimum (whichever was deeper; the upper depth limit of the mooring data depend on ocean current ‘knock down’ of the mooring string [Fig. S1]) downward to 400 m. The upper assimilation depth of 60 m was chosen through model calibration and data availability, as it minimized the addition of spurious temperature stratification from assimilation (this depth lies within WW and temperature change is minimal). Above this depth, which was not subject to assimilation, lies cold ‘winter water’ created during sea ice formation, or in summer, a remnant of this WW layer underlying a solar warmed surface fresh layer (Figs 1 and S1). Using this method, the top 60 m (or more depending on the depth of the temperature minimum) were model-driven in all cases, as opposed to being assimilated, allowing exploration of heat transfer from deeper waters through the surface layer. It is acknowledged, however, that because there is no ocean circulation in the surface layer, the potential effect of laterally-circulated warm and/or fresh water surface anomalies was not simulated.

In contrast, the assimilated subsurface inherently contained laterally-circulated waters as captured by the mooring observations. Ocean eddies for example passed by the mooring as frequently as quasi-weekly and often contained a less attenuated (i.e., warmer and saltier) CDW (Fig. S1) (Martinson and McKee, 2012). Because the salinity parameterization is imperfect, the assimilation process occasionally created an unstable water column with denser water found above lighter restored water. To prevent unrealistic vertical mixing, model temperature and salinity were interpolated across the instability from the bottom-most stable layer to the top-most layer of assimilated mooring data. This assimilation scheme permitted physical simulation of the top 60–100 m of the ocean (depending on the depth of temperature assimilation), allowing continuous estimation of ocean heat fluxes and their effect on sea ice, in addition to applying different scenarios of deep-ocean heat and atmospheric forcing (as described below).

2.5. Model initialization and ice-ocean coupling

Initial salinity and temperature vertical profiles were derived from a combination of mooring observations (Fig. S1) and CTD casts (Fig. S3) from the Palmer LTER program acquired during the annual austral summer research cruise along the WAP. Mooring observations were averaged in time over the two weeks prior to the annual CTD cast. The reason for this initialization procedure was to start all simulations for a given year (described below) with an equivalent (in terms of each year's annual cycle) vertical structure and heat content, allowing layer variability in the deep pycnocline (defined here as the minimum depth of the warm CDW) to drive model changes. Initialization values above the topmost depth of valid mooring observations to the surface (typically −60 m) were taken directly from CTD measurements. These CTD casts are acquired sometime in late January to early February, at a time when no sea ice was present. Model runs were started within 2 hours of the original CTD cast time and continued forward for a period of 550 d to capture the subsequent winter and summer evolution.

In the KEI model, approximating changes to the ice-covered fraction necessarily includes approximations and assumptions of the heat fluxes between the ice and ocean. The 1-dimensional nature of our model excludes the important processes of mechanical ridging and rafting. The single estimate of ice thickness (i.e., at each time step) cannot effectively simulate changes in ice amount and heat fluxes that would result from a distribution of ice thicknesses. For example, since typical Antarctic sea ice is weighted towards thinner ice types, thinner ice would grow and melt more quickly than the mean ice thickness, the latter being effectively simulated here (Bitz and others, 2001). For the purpose of examining the influence of ocean heat on sea ice, we have chosen computational methods that allow the sea ice to mimic observed ice cover changes, rather than parameterize the ice physics at the floe-level. A more realistic approach would be to simulate the ice-covered fraction by using a Lagrangian framework with multiple ice classes or categories that would allow a full ice thickness distribution to evolve (e.g., Schramm and others, 1997); however, without flow thickness or ridging/rafting observations to gauge model results, there is little data to understand the contribution of these processes to overall heat fluxes. Nevertheless, our more conceptual approach allows us to test the variable influence of deep ocean heat over time (in relation to heat delivered from above via solar warming of the mixed layer) and its effect on the seasonal evolution of sea-ice mass balance.

New ice formation is treated as follows. Once the mixed layer cools to near the freezing point, any additional atmospheric heat extraction results in supercooled water. The heat deficit of the supercooled water is then transformed into new ice, which we assume consists of a constant thickness of consolidated frazil ice. In reality, during the initial autumn expansion of Antarctic sea ice, large storm-generated swells can cause mechanical thickening of these frazil layers, allowing them to reach up to 40–50 cm in thickness before the pack ice solidifies (Doble and others, 2003; Kohout and others, 2014). Additional ridging and rafting of thin ice may increase ice thickness, while the opening of leads facilitates new ice growth and larger overall ice volumes. As an approximation between these various new ice types, we chose an initial ice thickness of 30 cm. This initial ice thickness also produces seasonal mean ice thicknesses between 0.6–0.8 m in the model, which appears reasonable given in situ observations of ice thickness observed along the WAP (0.37–1.02 m) (Perovich and others, 2004; Fritsen and others, 2008).

Considerations for sea-ice melt include the following. Steele and others (1992) demonstrated computationally that lateral melting at ice floe edges must be a tiny fraction of overall sea-ice melt due to the small surface area of ice floe edges compared to bottom surface area. This may be technically true, however the rapid retreat of the sea-ice edge in springtime involves additional mechanisms, primarily wind/wave-induced divergence (Kohout and others, 2014). This allows a feedback where increased shortwave radiance absorption is increased by a greater amount of open-water that in turn increases lateral melt, a process that also includes the fracturing of large floes into smaller ones, resulting in significantly greater edge exposure and lateral melt (Roach and others, 2018; Smith and others, 2022).

To simulate ice melt and ice-edge retreat in KEI, we apply a boundary layer melting potential F _bl (heat content of surface ocean boundary layer above freezing temperature, delivered over a time step) in a 2-step process. First, the basal heat flux F _w is calculated across a thin boundary layer at the ice-ocean interface (McPhee, 1992; Maykut and McPhee, 1995) as:

(2)$$F_w = \rho _wc_wc_{hw}u^\ast ( {T_w-T_f} ) $$

where ρ _w is the seawater density, c _w = 3.96 J g⁻¹ K⁻¹ is the specific heat of seawater, c _hw = 0.006 is a heat transfer coefficient (McPhee, 1992; Maykut and McPhee, 1995), u* is the friction velocity at the ice-water interface, and T _w and T _f are the temperature and freezing temperature of the seawater boundary layer. While the KPP mixing model does estimate local wind-driven lateral acceleration and inertial rotation, without momentum and other forcings in the 1-dimensional model domain, a mechanistic determination of u* is not possible and its value is fixed at 0.05 m s⁻¹ (Hunke and Lipscomb, 2010). The ocean to ice flux F _w is added to the conductive heat flux F _c through the sea ice (and overlaying snow, if any) to the basal ice surface (Saenz and Arrigo, 2012, 2014) to find the total basal ice heat flux, and the thickness of basal ice melt or growth dh/dt is determined as:

(3)$$\displaystyle{{{\rm d}h} \over {{\rm d}t}} = \displaystyle{{F_w + F_c} \over {-q_b}}\;$$

where q _b is the heat of fusion at the bottom layer of sea ice (Bitz and Lipscomb, 1999), and Fc results from the sea-ice heat flux integration (Saenz and Arrigo, 2012).

Second, because Eqn (2) may not allow transfer of all the boundary layer melting potential F _bl to the sea ice, any remaining heat potential is used to melt the ice laterally. The fraction of ice removed during melt df _i/dt is determined as:

(4)$$\displaystyle{{{\rm d}f_i} \over {{\rm d}t}} = \displaystyle{{F_{bl}-F_w} \over {-q_{tot}}}$$

where q _tot, calculated similarly to q _b but using all vertical sea ice and snow layers and their thickness and densities, represents the heat of fusion of the full sea-ice pack with units of J m⁻².

This second process is somewhat at odds with observed local-scale sea-ice melt processes where water temperatures can rise above freezing since the approach just described requires that surface ocean temperature remain near the freezing point until all ice has melted. On the one hand, keeping the surface ocean temperature near freezing slows melt-back and ice-edge retreat, since higher temperatures would accelerate basal ice melt. On the other hand, lateral melt increases the open water fraction, increasing surface-ocean heating through greater exposure to shortwave irradiance, and thus increasing basal melt rate (but at a slower rate relative to not keeping the surface ocean near freezing). Tests showed that our approach to the melting process was sensitive to observed surface forcing, and adequately simulated the timing of ice melt in accordance with satellite sea-ice observations, such that relatively long and short ice seasons were replicated faithfully. More complicated schemes of including or parameterizing ice edge melt processes in a 1-dimensional model did not seem warranted.

2.6. Simulation location and simulation types

The simulation location (300.100; 66.50°S, 69.87°W) corresponds to one of the Palmer LTER mooring locations (Fig. 2), which is located at grid station 100 along the 300 line of the Palmer LTER sampling grid (station numbers increase in 20-km increments from on-to-offshore, and grid lines increase in 100-km increments from southwest to northeast, i.e., along but not exactly perpendicular to the peninsula). The 300.100 station is located mid-shelf with bottom depth ~487 m, ~80 km west of Adelaide Island and ~60 km east of the shelf break (Fig. 2; Martinson and McKee, 2012). This site was selected for monitoring because it lies along the bank of a submarine canyon (Marguerite Trough) that bisects the continental shelf. The canyon helps to steer intrusions of CDW that appear on the shelf (McKee and Martinson, 2020b) and that are detected as transient upward and downward displacements of warm water that occur typically below the surface boundary layer. The 300 gridline, of which this station is a part, experiences a relatively cooler climate with a longer duration of ice cover than higher-numbered grid lines (e.g., 400) to the northeast, including locations at the coast (e.g., Palmer Station and Vernadsky/Faraday Station).

For a select number of years (2007, 2008, 2011, determined by the availability of mooring observations), we ran a suite of four simulations. The four simulations are characterized as follows: ‘FREE’ (no assimilation of satellite sea ice or ocean mooring data), ‘IC’ (ice conforming based on satellite sea-ice observations), ‘TC’ (thermocline conforming based on ocean mooring temperature observations), and ‘ITC’ (ice and thermocline conforming). The four different simulations over the 3 different years (with good yearly contrasts in both surface and subsurface forcing) allow us to isolate the effects of local air-sea fluxes, seasonally varying sea-ice concentration, and eddy-induced vertical exchanges of ocean heat, and to evaluate their respective controls on upper-ocean vertical structure and sea-ice seasonality and the evolution of sea-ice thickness.

The FREE simulations allow the surface forcing alone to drive sea-ice formation and melt and to determine the mixed layer evolution. In the IC simulations, the fractional sea-ice cover conforms to SSM/I satellite estimations, which introduces ice advection as a modifier of ice/ocean dynamics. Ice entering the model domain is assumed to be identical to existing ice (potential advective convergence effects on ice thickness were ignored). TC simulations incorporate ocean mooring observations as described in the ocean forcing section. ITC simulations used both ice and thermocline assimilation methods.

3. Results

3.1. Yearly and seasonal contrasts

The three simulation years (2007, 2008 and 2011) are characterized by different surface forcing conditions as reflected by monthly mean surface air temperature and winds (Figs. 3a, b), and by the observed daily sea-ice concentration and sea-ice seasonality (Fig. 3c). Surface air temperatures were generally warmer during 2008, while the winter minimum was differently timed and of varying magnitude in all three years (earliest and lowest in 2007, latest and weakest in 2011). During autumn, winds were relatively strong at the mooring location (>~9 m s⁻¹) in May-June 2008 and May 2011, while during spring, the strongest winds were in October 2007 (Fig. 3). Region-wide, winds were strong and mostly northwesterly throughout most of 2008 (Fig. S4b), whereas in 2007 strong northwesterly winds occurred primarily in August-September (Fig. S4a), and in March, May, October-December in 2011 (Fig. S4c).

Figure 3. Direct comparison of surface conditions between the three modeled years (2007, 2008, 2011). (a) Mean monthly 2-m air temperature and (b) wind speed from ECWMF Interim reanalyses in the vicinity of the Palmer LTER mooring location (Fig. 2). (c) NSIDC Daily Sea-Ice concentration (NASA Team algorithm; DiGirolamo and others, 2022) for the same three years from the nearest gridpoint to the mooring location. (d) Simulated sea-ice and snow thickness evolution using both assimilated sea-ice and ocean temperature observations at the mooring location (i.e., the ITC simulation; see also Fig. 4, bottom row).

Figure 4. Time series of simulated sea-ice (downward from 0 on y-axis, and shaded) and snow (upward on y-axis) thicknesses (meters) for the four simulation types, across three different years. Sea ice salinity is indicated by color.

Over the satellite era (1979–2020) both sea-ice concentration and the sea-ice season along the WAP show strong decreasing trends that are superimposed on pronounced year-to-year variability. The three years studied here (2007, 2008, 2011), though close together in time, strongly reflect this year-to-year variability (e.g., Fig. 3c, S4). Sea-ice variability is primarily driven by winds, and both wind speed and direction (whether on-ice or off-ice) drive changes in sea-ice seasonality and sea-ice concentration as illustrated in Figure S4. Of the 3 years, 2011 is most representative of near-mean sea-ice conditions (as observed over the satellite era), while 2008 is representative of the shortest sea-ice season (and lowest concentrations) yet observed.

We define ice-on as the first day of the ice season when sea-ice concentration reaches 15% or above and persists (above 15%) for 5 consecutive days thereafter (to eliminate spurious weather-related false-ice appearances, as well as the ephemeral advection of ice in/out of a given location on daily timescales). Conversely, ice-off is the day at the end of the ice season when sea-ice concentration drops below 15% and persists (below 15%) for 5 consecutive days. On average, the northward progression of the sea-ice edge reaches the mooring location by June 4th (based on ice-on [dates at the mooring location over 1979–2020), but in these three years, the ice-on date arrived substantially later by 3–5 weeks. In spring, the southward retreating ice edge on average passes the mooring location by November 30th (based on ice-off dates at the mooring location over 1979–2018), but in 2007 and 2008, the spring ice-off date was substantially earlier by 4 and 6 weeks, respectively, while in 2011 the timing coincided with the long-term average (November 30).

These three years are also distinctly different in terms of the simulated snow and sea-ice thickness. Across all four simulation types, 2008 has the thinnest sea-ice, while 2011 has the thickest (Fig. 4). The FREE simulations show the thickest sea ice relative to the other simulations (Fig. 4, top row), except notably the TC simulation in 2011. Although sea-ice thickness starts to thin by summer in the FREE simulations, the sea-ice fails to completely melt by end of summer. The thickest ice corresponds to the weakest ocean-ice heat fluxes (FREE simulations in all 3 years and the TC simulation in 2011) (Table 1; Fig. S5). Compared to the FREE simulations, the 2007 and 2008 TC simulations show reduced sea-ice thickness, more variable bottom melt (Fig. 4, second row) and generally stronger ocean-ice heat fluxes (Fig. S5), but with 2011 showing the opposite. Also, for both the FREE and TC simulations (Fig. 4, top two rows), high snow accumulation leads to seawater flooding at the snow-ice interface, resulting in high-salinity near the top of the sea ice (most notable in 2008 but present in all 3 years).

Table 1. Winter ice season duration and sea ice mass balance and relevant forcing variables

To allow comparison, means are calculated during the observed sea-ice period for all simulations. SAT is surface air temperature; THF is thermal heat flux.

The depth of the warm pycnocline varied between 80–130 m in both the IC and ITC simulations, but with differences in mean seasonal depths (Figs 5, 6; Table 1). The assimilated (ITC) ocean observations demonstrate pronounced high frequency variations in both the depth and temperature of pycnocline waters (Fig. S1). The heat content of the CDW layer on the WAP is a function of the CDW properties and the layer thickness, which do not necessarily vary in synchrony. A long-term warming trend of the CDW water mass is well-documented (e.g., Martinson and others, 2008; Schmidtko and others, 2014), whereas a positive trend in the heat content (Martinson and others, 2008) has leveled off in the early 2000's, since then being governed more by interannual variability (McKee, 2019, Fig. 1.6). The interannual variability is likely driven by local (Martinson and others, 2008; Spence and others, 2014) and remote (Spence and others, 2017; McKee and Martinson, 2020a) changes in the long-shore wind and its impact on isopycnal slope at the shelf edge or on-shore transport, which influence layer thickness. For context, 2011 is characterized by an anomalously deep pycnocline and thin CDW layer while 2008 is characterized by a relatively warm and thick layer. For all 3 years, the high frequency vertical excursions of subsurface water in the ITC simulations are revealed as distinct warm anomalies in difference plots between the ITC and IC ocean temperature simulations, particularly in 2008 (Fig. S6), which is consistent with the mean winter ITC pycnocline depth being shallower in 2008 than in 2007 or 2011 (88 m vs 99 and 97 m, respectively, Table 1, Fig. 6).

Figure 5. Vertical water column temperature profiles through time, comparing FREE (a, c, e subplots) and ITC (b, d, f subplots) simulations for 3 different years. Vertical dashed lines in June and July indicate the sea-ice onset period (black = FREE, red = TC, green = IC/ITC [observed]). Purple vertical dashed lines indicate the end of the satellite-observed (IC/ITC [observed]) sea-ice period. Sea ice does not melt in the FREE simulations. The density-derived mixed layer depth (in cyan) and KPP-defined turbulent mixing depth (gray line) are also indicated (see Section 3.5).

Figure 6. Comparison of (KPP-defined) turbulent mixing depth and deep pycnocline depth (density-defined depth of the warm CDW) between the IC and ITC simulations for 2007 (subplot a), 2008 (subplot b), and 2011 (subplot c). The vertical black dashed lines enclose the satellite-observed sea-ice period.

Directional (positive) component net heat fluxes (W m⁻²) at key physical boundaries illustrate seasonal differences in the net heat flux driving sea-ice evolution (Fig. 7). In all three years the autumn period, prior to the ice-on date, is dominated by net heat loss to the atmosphere from the seasonally warmed PBL (Figs 7a–c, pink stacks). Variability in the autumn PBL-atmosphere net heat flux is driven mostly by variability in atmospheric temperature and winds, which in turn drive cooling rates. With continued autumn cooling the PBL deepens, eroding the seasonal mixed layer, which now also releases some seasonally-gained heat to the pycnocline below (in orange), i.e., to the remnant WW layer. The relatively large heat injection into the remnant WW in the beginning of April coincides with a distinct warming of the remnant WW layer (Fig. 5). After early April, a small upward heat flux from the pycnocline to the PBL (in blue) occurs (in the ITC simulation) due to the seasonal deepening of the PBL (Fig. 5). Just prior to the ice-on date there is a larger amount of heat injected from the pycnocline to the PBL in 2007 and 2008 (Figs 7a, b in blue), due to an initial, but ephemeral, sea-ice growth and brine rejection event (e.g. free-melt cycling).

Figure 7. Stacked directional net heat fluxes for different simulation years at physical interfaces (subplots a, b, c) and difference between simulations (subplots d, e, f, g) in W m⁻². Vertical heights represent weekly-averaged flux magnitudes (or differences of flux magnitudes). The vertical black dashed lines enclose the satellite-observed sea-ice period. Subplots a, b, and c represent simulation ITC, and are comparable between years. The ITC fluxes shown represent the best approximation by the 1D KEI model of the dynamic and thermodynamic process that drive sea-ice evolution. Subplots d, e, and f represent the difference of simulations ITC-IC, and are comparable between years. The differences between ITC and IC simulations represent estimates of the flux contributions due to dynamic subsurface forcing. Subplot g is distinct from the others; it shows the difference between simulation ITC-TC for 2007, and is indicative of the flux contributions due to dynamic surface forcing.

During the period following the ice-on date, heat fluxes become more complex, as cooling induced by the atmosphere upon the PBL and sea-ice (Fig. 7 pink, purple) is countered by warming of the PBL and then sea-ice from below, i.e., across the deep warm pycnocline (Fig. 7 blue, green). In 2008, weakly average fluxes, including fluxes that freeze and melt sea-ice, are smaller overall than in the other two years. There are instances in the ITC simulation when there is a transfer of heat from the PBL downward as well (in orange) that are associated with downward excursions of the pycnocline (Figs 5, 6) that are also relatively cooler (Fig. S6) and less saline. The change in stratification, although relatively episodic (in relation to the vertical excursion), is associated with a redistribution of heat downward (fluxes in orange) via the dilation/contraction of pycnocline waters.

Throughout the winter ice-covered period, the magnitude and ratio of these heat fluxes vary in concert with the vertical excursions of the pycnocline (see Fig. 6) and the mixing of heat upward or in some cases downward, together with the variable sea-ice concentrations (Fig. 3c). By early September, sea-ice concentration starts to rapidly decrease, though somewhat intermittently (Fig. 3c), and by mid-September solar radiation begins to gradually increase, thus initiating some small downward fluxes from the atmosphere to the PBL (in brown) or sea ice (in red) and from the PBL to the pycnocline (in orange). Simultaneously, upward heat fluxes (in blue, green, pink) associated with intermittent sea-ice production but at decreasing sea-ice concentration (Figs 3 c, d) also continue until the ice-off date in early November. As expected, after the ice-off date, the heat fluxes are dominated by the downward heat flux from the atmosphere (in brown). In the following sections, differences between the simulations in regard to the seasonal timing, magnitude and source of heat will be discussed for each of the years.

Differences in directional fluxes between the ITC and IC simulations serve as an estimate of the contribution of subsurface ocean forcing upon the coupled atmosphere, ocean and sea-ice column (Figs 7–f). In 2008, there is very little difference between the ITC and IC simulations among outgoing atmospheric fluxes (Fig. 7e pink, purple), indicating the apparent changes in subsurface ocean forcing in the ITC simulation (Fig. 7e blue, green, orange) are not ‘feeding back’ and causing changes to surface fluxes. In 2007 and 2011, the subsurface heat assimilation drives fluxes of distinctly different magnitudes during the sea-ice period. In 2007 after ice-on, reduced heat fluxes across the deep warm pycnocline (Fig. 7d blue) are associated with sea-ice cooling (Fig. 7d purple). This pattern switches 4 weeks after ice-on, with much greater heat flux upward across the deep warm pycnocline (Fig. 7d blue) in the ITC simulation, which is primarily transferred to sea-ice (Fig. 7d green). In 2011, assimilated subsurface temperatures drive reduced heat transfer to the PBL and bottom of the sea ice (Fig. 7f blue, green); heat flux to the atmosphere through the sea ice is similarly reduced (Fig. 7f purple).

3.2. Timing of autumn ice growth onset and ice edge advance

The three years are distinguished by very different summer open water seasons preceding the timing of ice growth onset (Table 2). The shortest open water period was prior to the 2007 sea-ice season (221 d) and the longest was prior to the 2011 sea-ice season (264 d). A longer open water summer season permits greater solar warming of the upper ocean. The atmosphere-only forced FREE simulations reflect the open water season differences, with ice growth onset occurring earliest in 2007 (June 28), and latest in 2011 (July 20). Additionally a longer open water season may lead to deeper ocean mixed layers (due to longer exposure to wind mixing), and this too is reflected in the mixed layer depths (MLDs) of the FREE simulations. MLDs are shallower in 2007 and deeper in 2011 prior to the ice-on date (Table 2).

Table 2. Autumn sea ice onset dates and relevant forcing variables. SAT is surface air temperature

The ice-on date and MLD for 2007 and 2008 are similar for the FREE and IC simulations (Table 2). In 2011 the ice-on date is distinctly earlier by 25 d (June 25 vs July 20) and the MLD shallower (56 m vs 72 m) in the IC simulation due to the reduced exposure to winds beneath the earlier ice cover. Given that both the FREE and IC simulations were initiated with the same ocean profile and ran with the same atmospheric surface forcing, this suggests the ice-on date for 2011 was forced by subsurface ocean dynamics. Assimilated subsurface observations (TC simulation) drive the timing of the ice-on date to within 2 d to the observed in 2011 (6/27 vs 6/25). Subsurface temperature variability indicates that waters with decreased upper ocean heat content moved into the study area compared to what is predicted by the FREE or IC simulations (Figs. S5, S6c). In contrast, the FREE, IC and TC simulations in 2008 had ice-on dates all within 2 d of each other, similar MLDs (Table 2), and small inter-simulation heat flux differences (Fig. 7) compared to the other years, suggesting that the subsurface ocean was responding primarily to local atmospheric forcing in 2008.

In 2007, the ice-on date in the TC simulation is 11–12 d earlier than in the FREE and IC simulations (6/17 vs 6/28 and 6/29, respectively; Table 2). The ocean mixed layer is also shallower in the TC vs FREE and IC simulations (46 m vs 55 and 54 m) indicating that the subsurface ocean appears to have affected ocean mixed layer dynamics in fall 2007 (specifically on MLD). The shallower mixed layer in the TC simulation is coincident with observed (assimilated) shoaling of the pycnocline, and is consistent with previous reports that variability in pycnocline depth can affect mixed layer dynamics in this area, at both seasonal and quasi-weekly timescales (Figs 5, 6) (Martinson and McKee, 2012; McKee and others, 2019). Generally, a shallower mixed layer cools faster thus helping to explain the earlier ice-on date for the TC vs FREE and IC simulations in 2007 (Table 2; see Venables and Meredith, 2014).

It is also worth noting that all three years experienced anomalously strong regional, peninsula-ward winds and southward displaced ice-edges during May-June (Fig. S4) such that the ice-on dates for all three years (29 June, 6 July and 25 June, respectively; Table 2) were late compared to the long-term (1979–2020) mean (4 June). However, in 2007, the ocean had cooled by 17 June and was already growing sea-ice, but winds kept the location clear until 12 d later (29 June). By mechanically delaying an ice-on date, our simulations indicate a considerable amount of ocean heat can be additionally extracted from the ocean to the atmosphere (up to 170 W m⁻²; Fig. 7a, the peak prior to ice-on date). The amount of ocean heat ventilated just prior to the ice-on date is also greater in the IC vs ITC simulations (Fig. 7).

3.3. Seasonal evolution of sea-ice thickness during winter

In the FREE simulations, yearly differences in winter sea-ice thickness are forced only by surface cooling and mean depth of warm pycnocline waters. Although the sea-ice season is unrealistically long in the FREE simulations, the overall thinner winter sea ice in 2008 vs 2007 vs 2011 (Fig. 4 top row, Table 1) is consistent with relatively warmer mean winter temperature in 2008 (coolest in 2011) and shallower mean pycnocline in 2008 (deepest in 2011; Table 1).

When the observed subsurface ocean temperatures are assimilated (TC simulations), the episodic intrusions of warmer (less modified) CDW higher in the water column (Figs 5, 6) permit greater-intensity bottom melt/freeze episodes via a linked mechanism. Shallower warm water in the TC simulation permits mixing processes to deliver more heat to the basal sea ice, which melts and thins the sea ice considerably (Table 1; Fig. 4). In turn, thinner sea ice allows greater subsequent sea ice growth during further freeze cycles due to the reduced insulation of thinner sea-ice (and snow) cover over a given period. This was most apparent in 2007, where 1.12 m of sea ice grew in the FREE simulation (0.58 m mean thickness) while 1.48 m of sea ice grew in the TC simulation (0.37 m mean thickness) over the same period (Table 1). A very different picture is presented in 2011 with sea-ice thickening beginning more slowly (than in the FREE simulation) but reaching thicker sea ice in late winter (Fig. 4), consistent with overall cooler assimilated subsurface ocean temperatures in 2011 (e.g., Fig. S6).

The sea-ice thicknesses in IC simulations are strikingly different from simulations without sea-ice concentration assimilation (FREE, TC), largely due to the abruptly shorter ice season (Fig. 4, third row from top; Table 1) and more subtly, the greater variability in daily sea ice concentration (Fig. S5, top panels). Despite the assimilated differences in sea-ice concentration and seasonality, sea-ice thickness evolved similarly in the FREE and IC simulations in 2007 and 2008 over the period of assimilated satellite-observations. In year 2011, the IC simulation was forced to have an earlier ice onset compared with the FREE simulation (by 25 d, Table 2). The upper ocean still had considerable residual heat at the start of the ice-on date in the IC simulation (as discussed in the previous section), and the mean winter depth of the pycnocline was shallower by 21 m (Table 1), such that when the ice-on date occurred (based on satellite observations), it went through several basal melt cycles before it slowly thickened.

The evolution of sea-ice thickness in the 2007 and 2008 ITC vs TC simulations (Fig. 4 bottom row) both show the effects of the assimilated vertical excursions of warm pycnocline waters on melt/freeze cycles in sea-ice thickness. In 2011, sea-ice thickness evolution differed from the other two years where by mid-October, sea-ice thickening ceases in the ITC simulation, whereas sea ice thickness continued to increase in the TC simulation. This difference is due to the more variable sea ice concentrations in the ITC simulation that result in a thinner snow cover. With less insulative snow cover and more variable sea-ice concentrations, the bottom melt/freeze cycles imparted by subsurface ocean variability become more accentuated. There is also a stronger contrast between initial sea-ice thickening in the ITC vs IC simulation, the former showing thicker sea ice (by 7 cm, Table 1, Fig. 4), especially early in the ice season. The thicker sea ice in the 2011 ITC simulation is due to the assimilated subsurface ocean being overall cooler than the simulated subsurface ocean (Fig. S6), permitting more rapid sea-ice growth at the start of sea-ice growth season. The 2011 ITC simulation also shows a large amount of total sea-ice growth compared to its mean thickness; ice production more than doubled compared to the IC simulation (Table 1).

3.4. Seasonal evolution of spring sea-ice break-up and melt

In the FREE simulations, sea-ice and snow thicknesses are too thick to melt-out completely the following year, except for 2008 which does finally melt-out, albeit 5 months later than observed (Table 3). In the TC simulations, adding deep ocean forcing helps in thinning the sea-ice cover, particularly in 2007 and 2008 (Fig. 4), but again, only in 2008 does the sea-ice completely melt the following summer (though still 2 months later than observed; Table 3), indicating that subsurface ocean heat variability is not sufficient to drive large/first-order changes in the sea-ice duration. The lack of complete sea-ice melt in majority of the FREE and TC simulations indicates the climate of the WAP is cool enough to require mechanical (advection, leads) and ice-edge effects to permit sea-ice retreat in some years.

Table 3. Spring sea-ice melt/break-up dates and relevant forcing variables. SAT is surface air temperature

In simulations with assimilated sea-ice concentration, despite being much thinner, sea ice does not completely melt by the end of specified ice-off date. In 2007 and 2008 IC and ITC simulations, the sea-ice thickness at ice-off is comparable to midwinter (Fig. 4 third, fourth rows), indicating the majority of sea-ice melt occurs elsewhere following advection of the ice pack from the study area. In 2011, IC and ITC simulations show rapid sea-ice thinning before ice-off, indicating much more of the sea ice melts locally.

3.5. Mechanisms of assimilated subsurface heat transfer

To explore these effects mechanistically, we first note that the KPP model distinguishes between a boundary layer depth (based on instantaneous depth of surface driven turbulence; solid gray line in Fig. 5) and a mixed layer depth (determined by density stratification; cyan line in Fig. 5). The depth of these two metrics shows the mixed layer to be deeper than the boundary layer over much of the simulations (Fig. 5), indicating a relatively thick zone of weak, but nonzero stratification. The effect of the upward vertical excursions of the pycnocline in the TC/ITC vs FREE/IC simulations equates to an additional 20–60 W m⁻² of net upward ocean heat flux in 2007 and 2008, and 5–20 W m⁻² upward ocean heat flux in 2011 (following an enhanced freezing period; Fig. 7). Higher heat fluxes from the PBL to the pycnocline are made possible by a shoaling pycnocline, as transient atmospheric heating has a greater effect on a reduced mixed layer volume.

We also note that the KPP model has high vertical diffusivities in the boundary layer, with peak values as high as 0.1 m² s⁻¹ and mixed layer mean values of ~0.03 m² s⁻¹ (during winter), but typically the temperature gradient in the boundary layer is very small. Higher diffusivities (higher than the background value of 10⁻⁵ m² s⁻¹) can penetrate the upper pycnocline, an important aspect of the KPP model dynamics, due either to deep turbulent mixing by strong winds or brine-induced vertical mixing by sea-ice production. To quantify the effect of these deep mixing events, we computed mean diffusivities and turbulent heat fluxes in the upper pycnocline in a band ± 20 m about the deep pycnocline and then assessed the combined effect on heat transfer from the ocean to the sea ice (Table 1). The heat transfer from the ocean to the sea ice is higher in TC/ITC vs the FREE simulations, particularly in 2007 and 2008 (Table 1), but also in 2011 for the ITC simulation. Interestingly, the mean diffusivities calculated from the same depth interval as the turbulent heat fluxes are ~50% lower in the TC simulations, due to the generally sharper density gradient present in the assimilated temperature profiles. Heat transfer to the sea ice is driven by the difference in temperature between the surface ocean and the sea-ice temperature, therefore episodic vertical intrusions of typically warmer pycnocline waters assimilated in the TC/ITC (vs FREE) simulations allowed greater heat transfer to the sea ice by causing greater episodic increases in temperature of the (often thinner/smaller volume) PBL, despite both lower turbulent potential and deeper mean winter pycnocline depths (Table 1).

Under the same dynamical surface conditions (ITC vs IC, Fig. 7g, which repeats Fig. 7d), the addition of subsurface dynamical variability (the ‘TC’ in the ITC simulation) leads to a seasonal shift from lesser heat flux upward from the deep pycnocline in autumn to greater heat flux in winter. As described above, the lesser heat flux in autumn (relative to the IC simulation) is due to a shallowing of the ocean mixed layer in response to a shoaling, and sharpening, of the pycnocline (Figs 5, 6) that in turn dampens deep mixing and the entrainment of heat from below. However, by early winter, greater heat flux across the deep warm pycnocline is facilitated by the episodic vertical excursions in the pycnocline and brine rejection from sea ice (Fig. 6).

Under the same dynamical subsurface forcing (ITC vs TC), the addition of surface dynamical variability (the ‘IC’ in the ITC simulation) leads to greater heat flux to the atmosphere in autumn, in association with later ice-on dates and because of generally lower sea-ice concentration (Fig. 7g). Greater heat flux to the atmosphere continues into the winter as a result of more variable and lower sea-ice concentration (via advective assimilation), but greater deep-ocean heat flux to the PBL and sea ice also become evident, resulting from episodic upticks and downticks in sea-ice concentration that in turn release subsurface ocean heat upward. The upticks in sea-ice concentration are associated with increases in sea ice production, brine rejection and deep mixing, while the downticks in sea-ice concentration facilitate wind-driven deep mixing.

4. Discussion

4.1. Effects of subsurface ocean forcing on sea-ice seasonality and sea-ice thickness

By comparing three different sea-ice seasons and four differently prescribed modeled ice-ocean simulations, we identified two effects of subsurface ocean forcing on sea-ice seasonality and sea-ice thickness evolution not previously well-described: (1) periodic shoaling of the pycnocline in maintaining a (comparatively) shallow boundary layer, particularly in autumn (Figs 5, 6), and its subsequent role in modulating ice-on dates (Table 2), and (2) quasi-weekly vertical excursions of pycnocline waters that modify the turbulent heat exchange particularly during the sea-ice covered period, leading to both positive (2011) and negative (2007, 2008) perturbations to sea-ice thickness (Figs 5, 6, and S6).

The implications of subsurface forcing of sea-ice phenology are well illustrated by the specific combination of surface and ocean forcing in autumn of 2007, which led to extensive connection between sub-pycnocline waters and the atmosphere. During autumn, cooling of the surface waters leads to sea-ice formation, but wind-driven advection of newly formed sea ice not only affects the timing of ice-edge expansion but also air-ocean heat fluxes (Massom and others, 2006; Stammerjohn and others, 2008b; Stammerjohn and others, 2011; Turner and others, 2013; Venables and Meredith, 2014). Subsurface ocean heat can then be transferred to the surface through turbulent diffusion, or, more actively, through convective mixing from brine rejection during sea-ice growth (e.g., Martinson and Iannuzzi, 1998). Based on our model simulations with assimilated satellite sea-ice concentration and subsurface ocean observations (ITC runs), ~50 W m⁻² of deep ocean heat was ventilated over a 12-d delay in the ice-on date in 2007 in the TC vs ITC simulation, demonstrating how the particulars of seasonality in the WAP can variably connect the deep ocean to atmosphere. The overall potential for ocean heat ventilation during anomalously late ice-on dates is consistent with other studies, and for the WAP area in particular, the trend towards a later ice-on date and increased ocean heat flux in late autumn has been linked to amplified atmospheric warming in autumn-winter (Meredith and King, 2005; Stammerjohn and others, 2008a, 2008b; Turner and others, 2013).

Venables and Meredith (2014) outline a related physical feedback mechanism where deep mixing and subsequent summer surface water heat entrainment due to low winter sea-ice cover further contribute to a reduction in the sea-ice cover the following autumn. We illustrate the mechanism by which deep mixing during winter also can thin and reduce (current) sea-ice cover by tapping a different heat source: shoaling CDW waters. Venables and Meredith (2014) note that mCDW water in Ryder Bay (coast-ward of the 300.100 mooring site) is decreasing in temperature, indicating these mechanisms for modulating sea ice may interact. The CDW observed further out on the shelf appears to show less signs of cooling; however, given the pace of change of climate in the WAP, further study of the interplay of surface and subsurface heat delivery to the shrinking sea-ice cover in the WAP is warranted.

Sea-ice thickness and salinity were most highly modified by including sea-ice concentration assimilations (IC, ITC simulations), affirming that surface dynamical processes implicitly included in satellite-observed sea-ice concentrations (such as wind advection, lead formation, and other rheological properties) imparted a larger effect on sea ice than subsurface dynamical forcing (e.g., Bitz and others, 2005). Nonetheless, subsurface dynamical forcing also modulated sea-ice thickness and salinity, and while surface forcing alone always resulted in decreased sea-ice thickness, subsurface forcing imparted either sea-ice thickness increases or decreases. Sea-ice thickness exerts a first-order control on the sea-ice rheology (Hibler, 1979; Feltham, 2008); changes in mean sea-ice thickness due to subsurface forcing could therefore influence processes such as lead formation/divergence and ridging that feedback into surface dynamical forcing in this region. The additional brine convection associated with the substantial increase in total ice growth observed in some simulations alters the redistribution of salt in the water column, where it has attendant effects in mixed-layer dynamics and circulation (Matsumura and Hasumi, 2008; Lui and others, 2019).

Additionally, subsurface ocean forcing on sea ice may impact sea ice ecology. The amount of snow-ice formation determines the amount of habitat for the infiltration of algal communities, the most highly productive sea-ice communities (Arrigo, 2003, 2014; Saenz and Arrigo, 2014). Rapid melt-freeze cycling induced by the periodicity of subsurface heat flux to the sea-ice bottom could either enhance algal productivity via brine convection (Mundy and others, 2007) or remove lower-ice-associated algal communities through melt.

4.2. Implications of our modeling approach and study area

The tool developed and used here to investigate the potential contributions of different surface and especially subsurface ocean forcing on sea-ice seasonality and thickness is a coupled 1-dimensional column of snow, sea ice and ocean, and by design lacks 3-dimensional ocean and sea-ice circulation and attendant effects. Dynamical variability is therefore introduced to the model by assimilating satellite observations of sea-ice concentration and mooring observations of subsurface ocean temperature. We acknowledge that our modeling framework is simplified, and that we only partially account for horizontal variability in our simulations with assimilated observations. However, using a 1-dimensional approach (vs a more complex 3-dimensional approach) allowed separation and assessment of the relative contributions of strictly thermodynamically-imposed vs dynamically-imposed variability on a system that is strongly coupled in the vertical. The 1-D model also complements full 3-D ocean-sea ice simulations where errors relative to local mooring and CTD observations reflect a complicated balance of uncertainties in surface forcing, lateral boundary conditions, and model sea-ice and ocean physical dynamics (e.g., Schultz and others, 2020).

We also acknowledge that our results involving the assimilated subsurface ocean temperatures are dependent on the scheme we used to interpolate salinity. The parameterized temperature-salinity (T-S) relationship used in this study (Eqn (1)) is based on the range in T-S observed for modified CDW (mCDW; CDW water that has mixed with another water mass) as sampled at a mid-shelf location (i.e., near the mooring site) by the Palmer LTER program every January, as well as from a single mid-depth moored salinity sensor. These data likely do not include the full range in mCDW T-S that results from the episodic and seasonal variability in mCDW at this location. Episodic variability is due in part to intermittent delivery of upper CDW from the continental slope and its subsequent vertical and lateral mixing at peripheries of mesoscale structures to form mCDW (McKee and others, 2019). Seasonal variability in mooring temperatures show a warmer WW, cooler CDW, and a weakened vertical temperature gradient in April-May (Fig. S1). Potentially strong winds in these months (e.g., Fig. S4) would enhance vertical mixing of heat, and their directionality (southwestward) could reduce shelf-slope exchange (McKee and Martinson, 2020a) and hence they provide a seasonal signal to T-S properties of mCDW not captured in our parameterization.

Ocean temperature observations from this area of the WAP region show that the remnant WW layer is often eroded by end of summer/early autumn due in part by the mixing-down of surface heat and fresh water and/or the mixing up of subsurface ocean heat and saltier water. Because we do not have coincident salinity observations, the slightly warmed WW could either be fresher or saltier depending on the degree of deep mixing. During the autumn transition, what is shown in the TC and ITC simulations is a consistently shallower PBL as compared to the IC and FREE simulations (Figs 6, S5). How much of this is due to the seasonal shoaling of the pycnocline observed for this location, during autumn in particular (Figs 5. 6), vs imposed by our interpolated salinity is unknown. However, once sea-ice growth is well underway, the mixing depth becomes far more variable in all simulations, thus erasing any marked differences in mean mixing depths between the IC and ITC simulations (Fig. 6).

5. Conclusions

Quasi-weekly and seasonal variability in the depth and temperature of the pycnocline in WAP waters demonstrated potentiality to drive both generally positive (2007, 2008) and negative (2011) deep-ocean heat transfer to surface processes compared to a passive subsurface ocean. Sea ice exerts a strong positive feedback through its high reflectance of incoming solar radiation thus enabling it to seasonally persist longer than it would otherwise (Curry and others, 1995; Perovich and others, 2007). Martinson (1990) showed that another feedback exists below Antarctic sea ice, where a limit to sea-ice growth exists depending on the depth of warm pycnocline waters. In the WAP, mooring observations reveal quasi-weekly variability in both the depth and temperature of pycnocline waters, with vertical excursions of 20 m or more on daily time scales, superimposed on more smoothly varying low-frequency modes of variability (Figs 5. S1). Here we show this variability in the depth and temperature of pycnocline waters, particularly on ‘warm shelf’ regions, are complexed with the above processes to impart significant changes in sea-ice seasonality and the winter evolution of sea-ice thickness.

Given the scarcity of ocean observations for continental shelf seas ~ Antarctica, it has been difficult to assess the role of subsurface ocean forcing on seasonal sea-ice concentration variability. Regional and global coupled models often perform poorly when simulating sea-ice seasonality, extent, concentration, and thickness (Blanchard-Wrigglesworth and others, 2015; Guemas and others, 2014; Wayand and others, 2019). Sun and Eisenman (2021) demonstrate that improving advection of sea ice though assimilation leads to more accurate simulations of sea-ice phenology; analogously these model results demonstrate that high-frequency subsurface forcing has effects that will influence sea-ice characteristics and attendant climate-regulating processes, air-sea fluxes and deep ocean ventilation. Combined representation and/or parameterization of the coupled vertical surface and subsurface processes explored here will improve our ability to model and predict Antarctic sea-ice seasonality and thickness, and its global impact on heat budgets, ocean circulation and marine ecosystems.