Modeling of the In ﬂ uence of Sea Ice Cycle and Langmuir Circulation on the Upper Ocean Mixed Layer Depth and Freshwater Distribution at the West Antarctic Peninsula

The Southern Ocean is chronically undersampled due to its remoteness, harsh environment, and sea ice cover. Ocean circulation models yield signi ﬁ cant insight into key processes and to some extent obviate the dearth of data; however, they often underestimate surface mixed layer depth (MLD), with consequences for surface water ‐ column temperature, salinity, and nutrient concentration. In this study, a coupled circulation and sea ice model was implemented for the region adjacent to the West Antarctic Peninsula, a climatically sensitive region which has exhibited decadal trends towards higher ocean temperature, shorter sea ice season, and increasing glacial freshwater input, overlain by strong interannual variability. Hindcast simulations were conducted with different air ‐ ice drag coef ﬁ cients and Langmuir circulation parameterizations to determine the impact of these factors on MLD. Including Langmuir circulation deepened the surface mixed layer, with the deepening being more pronounced in the shelf and slope regions. Optimal selection of an air ‐ ice drag coef ﬁ cient also increased modeled MLD by similar amounts and had a larger impact in improving the reliability of the simulated MLD interannual variability. This study highlights the importance of sea ice volume and redistribution to correctly reproduce the physics of the underlying ocean, and the potential of appropriately parameterizing Langmuir circulation to help correct for biases towards shallow MLD in the Southern Ocean. The model also reproduces Kemper are tested Each of these parameterizations was developed for open ocean and implemented in global models. The of sea ice and a coastline in our regional model an additional challenge given the uncertainties related to the behavior of wind and surface waves in the presence of sea ice, and the different properties of waves in coastal areas. The approach described by Li and Fox ‐ Kemper (2017) parameterizes the in ﬂ uence of LC in the entrainment of dense water from below the mixed layer. of entrainment, entrainment ice drag coef ﬁ cient leads to large differences in the representation of seasonal and interannual sea ice variability. The CTRL simulation, which has a higher drag coef ﬁ cient (2 × 10 − 3 ), shows a bias towards higher SIC and late sea ice retreat when compared to the observations. These biases are ampli ﬁ ed in the simulation with low (5 × 10 − 4 ) drag coef ﬁ cient. The difference between simulations with different drag coef ﬁ cients is much larger during the period of sea ice retreat than during sea ice advance, suggesting that retreat is more in ﬂ uenced by wind action. The smaller differences in SIC between models during sea ice that there are biases in the initial conditions of ocean temperature and salinity used for the model.


Introduction
The coastal ocean on the western side of the West Antarctic Peninsula (WAP) has undergone dramatic changes over the past several decades. During the second half of the 20th century, the WAP region showed the highest atmospheric warming in the Southern Hemisphere (Vaughan et al., 2003), a decrease in the duration of the sea ice season , and an accelerating retreat of glaciers in the region (Cook et al., 2005Cook & Vaughan, 2010). The surface ocean also exhibited an increase in temperature of more than 1°C, with deeper layers also displaying marked warming trends Meredith & King, 2005;Schmidtko et al., 2014). The increased surface ocean temperature, thought to be caused by increased temperatures in the atmosphere and the decreased sea ice cover, could lead to positive feedback that would further increase the atmospheric warming. Around the turn of the 21st century, however, the atmospheric warming trend leveled off and even reversed direction in some places (Turner et al., 2016). This trend reversal can be attributed to natural climate variability and indicates that the changes observed in the WAP cannot be explained by large-scale global warming alone.
The impact of changes in heat and freshwater input on vertical mixing in the ocean has not been fully quantified. While sea ice melt decreases the density of the surface layer and contributes to stronger near-surface stratification, decreased sea ice cover also leads to increased wind-induced ocean mixing. The effect of sea ice changes on the surface mixed layer depth (MLD) therefore relies on a balance between these two competing effects, which is further complicated by laterally redistributing and mechanically thickening sea ice cover. Persistent lateral motion of the sea ice in one direction can lead to geographic separation of ice production and melt, that is, areas of net sea ice production with increased surface salinity and areas of net sea ice melt with low surface salinity, and thus spatial variability in the stability of the surface mixed layer (Meredith et al., 2008Petty et al., 2014).
Studying the impact of each of these changes on the MLD is important, given that it may affect not only trends in the distribution of water masses in the region but also oceanic primary production, ecosystem dynamics, and the carbon cycle in the WAP (Legge et al., 2017;Vernet et al., 2008). Understanding the impact of increased glacial melt is also important, given that salinity gradients are thought to contribute to the coastal circulation in the WAP (Moffat et al., 2008), and glacial freshwater can be a major source of dissolved iron for primary production (Annett et al., 2017).
One of the challenges of studying the MLD in the WAP is the scarcity of data. Although the Palmer-Long Term Ecological Research (LTER) project has been conducting yearly cruises during January and February since 1993, thus providing an invaluable dataset of over two decades of duration, each cruise constitutes a snapshot relevant to the period of its occupation . Sustained time series with high-frequency sampling exist, including the Rothera Time Series and stations occupied quasi-weekly adjacent to Palmer Station ( Figure 1); however, these point measurements at near-coastal sites do not have full spatial extent across the WAP shelf.
Innovative technologies such as underwater gliders and tagging of marine mammals are helping in the collection of synoptic scale data throughout the shelf that can be used to fill in the gaps in the temporal and spatial coverage of the current datasets (Carvalho et al., 2017;Costa et al., 2008;Couto et al., 2017;Hückstädt et al., 2020;Venables et al., 2017). In addition, there is great scope for high-resolution regional modeling to aid understanding of key processes and the causes of change. Such models can also be used to diagnose the source and amount of freshwater injected into different areas of the WAP shelf and to study the influence of different processes (e.g., wind mixing and freshwater-induced stratification) on MLD and other key ocean properties.
Surface MLD resolved by ocean models are often too shallow, especially in the Southern Ocean (Belcher et al., 2012;Noh & Lee, 2008;Schiller & Ridgway, 2013), and this bias has been attributed partially to the neglect of Langmuir circulation (LC) that results from the interaction between winds and surface waves (Belcher et al., 2012;D'Asaro, 2014). LC comprises a series of shallow vortices that counterrotate on a vertical plane in the surface mixed layer. Strong velocity and vertical shear of the velocity in the Langmuir vortices enhance near-surface turbulence, erode stratification at the top of the thermocline, and thus deepen the surface mixed layer (e.g., Li et al., 2013).
How wind, waves, and Langmuir-induced turbulence might affect MLD and mixing processes in the presence of an ice cover is an active area of research. Although the presence of a consolidated ice cover can greatly attenuate waves, a less consolidated ice cover, which typically characterizes the WAP seasonal sea ice, can permit ocean waves to travel hundreds of kilometers into the ice cover (e.g., Kohout et al., 2014), potentially affecting sea ice growth processes (Doble et al., 2003), sea ice distributions (Kohout et al., 2014), and sea ice melt and near-surface mixing processes (Smith et al., 2018).
As the first part of the Palmer-LTER modeling component to study biophysical and biogeochemical processes on the WAP shelf region, this study aims to understand how different physical factors affect surface MLD in the region, with a view to enable subsequent works on biophysical and biogeochemical processes on the WAP shelf. To test the influence of LC on the MLD in the WAP region, this study tests two different parameterizations of LC, based on  and . While these parameterizations have been tested in global models, this is the first study to implement them in a regional Southern Ocean model. Simulations with different wind-ice drag coefficients are also analyzed to evaluate the influence of sea ice and wind mixing on the MLD at both seasonal and interannual time scales. The simulations are also compared to the available freshwater data (based on oxygen isotope analysis of surface seawater samples collected annually) to determine the relative impact of sea ice production/meltwater and meteoric (glacial and precipitation) freshwater inputs in different areas of the WAP shelf (e.g., Meredith et al., 2017). Results highlight the critical need to correctly incorporate sea ice dynamics to properly simulate the ocean physics in the WAP.
The structure of the paper is as follows. The description of the model used in this study and the parameterizations added to represent LC are described in Section 2. The cruise and satellite data used to validate the model results are described in Section 3. Section 4 contains the results of this study and discusses the observed and modeled advance and retreat of sea ice, and the variability and seasonality of the surface MLD, as well as comparing the vertical profiles of temperature and salinity in the model and observations. The discussion and conclusions of the study are presented in Section 5.

Model and Grid Description
This study uses a hydrostatic version of the Massachusetts Institute of Technology General Circulation Model (MITgcm) and the embedded sea ice and ice shelf modules. MITgcm solves Boussinesq Navier-Stokes equations with the hydrostatic assumption on an Arakawa-C grid and z-level (fixed depth) vertical layers . The sea ice component of the model is based on Hibler-type ice theromodynamics and elastic-viscous-plastic ice rheology (Losch et al., 2010). The ice-shelf component (Losch, 2008) is also used here to include the freshwater input from ice shelf melt. MITgcm has been successfully used to resolve large-scale circulation, sea ice dynamics, and biogeochemistry of the Southern Ocean with horizontal resolutions on the order of a degree of latitude and longitude (Hauck et al., 2015;Holland et al., 2014;Taylor et al., 2013), and physical configurations with higher resolution (of the order of kilometers) have also been used in the Southern Ocean (Regan et al., 2018;Schodlok et al., 2012). (Fretwell et al., 2013) but is complemented with data from the General Bathymetric Chart of the Oceans (GEBCO, www.gebco.net) north of 60°S, which is the northern extent of Bedmap2.

Initial and Boundary Conditions and Surface Forcing
The initial conditions for temperature and salinity were obtained from the World Ocean Atlas (WOA09, Antonov et al., 2010;Boyer et al., 2009;Locarnini et al., 2010). The model has three open boundaries (north, east, and west) with prescribed monthly mean ocean temperature, salinity, and velocity, as well as sea ice concentration (SIC), thickness, and velocity obtained from a large-scale, low-resolution circumpolar MITgcm simulation . On the eastern side of the Antarctic Peninsula, the southern boundary was closed, which could have local effects on the circulation that are beyond the scope of this study but do not affect the currents at the northern tip of the peninsula. Sea ice velocities are not prescribed at the boundaries to avoid spurious ice convergence. The atmospheric forcings were obtained from the ERA-Interim reanalysis (Dee et al., 2011) with a horizontal resolution of 1.5°in latitude and longitude and a 6-hr interval. ERA-Interim has been found to be the most reliable reanalysis for some atmospheric parameters (Bracegirdle & Marshall, 2012). However, it is also known that most reanalyses, including ERA-Interim, underestimate moderate to high wind events in the Southern Ocean (Dong et al., 2020;Yuan, 2004), which could potentially lead to decreased surface mixing (and therefore shallower MLD) and delayed advance and retreat of sea ice, given ice movements are largely driven by wind.
As this model was originally used to study freshwater distribution in the Bellingshausen Sea, the freshwater inputs were carefully chosen to be as realistic as possible, given the constraints imposed by lack of data (Regan et al., 2018). A general surface runoff input was imposed to represent surface melt of land-based ice, and iceberg calving and melting. Values of the total surface runoff, including iceberg decay, are chosen to represent the amount estimated by Van Wessem et al. (2016), and they are distributed uniformly along the coast and linearly decreasing with distance from land out to 100 km.
The estimated runoff to the ocean depends on the climatological and large-scale surface mass balance of freshwater over the continent, and on the estimated iceberg melt, which is highly variable in space and time. In contrast to the steady freshwater input used in Regan et al. (2018), this study introduces a seasonal cycle of the freshwater input with greater values in summer and smaller values in winter. The annual mean glacial input was multiplied by a cosine anomaly in time, which gives a different value for each month and maintains the total input amount over the course of a year. This modification is based on the expectation that most of the freshwater input to the ocean occurs during the summer months. The lack of year-round data makes it impossible to incorporate a more accurate representation of the seasonal cycle of the freshwater sources with presumably high spatial and temporal variability. It is possible that inaccuracies in this imposed freshwater flux could affect the timing of sea ice formation in the coastal areas.

Parameterization of Langmuir Circulation
The model uses the K-profile parameterization (KPP; Large et al., 1994) to parameterize subgrid scale vertical turbulence in the surface boundary layer and ocean interior, and here we consider the influence of LC in the KPP framework as an enhancement factor to the turbulence velocity scale used for the surface boundary layer. The time scale of ocean surface layer responding to winds, forming surface waves, and establishing/ dissipating Langmuir vortices is relatively short (Talley et al., 2011). The impact of LC on mixing and surface MLD could thus be considered by directly coupling surface wave models with atmospheric conditions. However, direct simulation of the ocean surface waves over the spatial and temporal scale of our interests is computationally prohibitive, and there are limited data for validating wave models in the Southern Ocean. For this reason, this study employs empirical wave spectra, following , to estimate the surface wave conditions using winds at 10 m above the sea surface.
Two approaches for including LC, based on the parameterizations of  and , respectively, are tested here. Each of these parameterizations was developed for open ocean and implemented in global models. The presence of sea ice and a coastline in our regional model constitutes an additional challenge given the uncertainties related to the behavior of wind and surface waves in the presence of sea ice, and the different properties of waves in coastal areas. The approach described by  parameterizes the influence of LC in the entrainment of dense water from below the mixed layer. Two main factors affecting the rate of entrainment, that is, entrainment buoyancy flux, w′b′ e , are (i) destabilizing surface buoyancy flux, important in convective turbulent regimes, and (ii) shear instability at the base of the surface mixed layer, dominant under wind-driven regimes. The latter factor is often associated with processes such as inertial oscillations of the surface current. Here w′ is vertical turbulent velocity anomaly and b′ is buoyancy fluctuation. When both factors are present, a proportionality is assumed and a new velocity scale, w x , combining convective velocity scale and water-side friction velocity (Large et al., 1994), is introduced. w x relates to the rate of entrainment so that (1) where h b is the boundary layer depth and is defined as the smallest depth at which the bulk Richardson number reaches a critical value of 0.3. Here u r is a reference speed, u is water speed, b r is a reference buoyancy, b is buoyancy, and U t is a velocity scale that accounts for the influence of subgrid velocity and shear. Several large-scale eddy simulations (e.g., Grant & Belcher, 2009;McWilliams et al., 2014) have suggested that w′b′ e is enhanced in the presence of LC. This would happen because LC enhances turbulent kinetic energy at the base of the mixed layer. Langmuir turbulence could also contribute to the erosion of the thermocline through enhancing shear instability beneath the downwelling regions of the LC cells. However, quantifying the Langmuir enhanced w′b′ e is nontrivial because it is a small portion of the overall turbulent kinetic energy budget in the surface mixed layer ).
The solution proposed by  includes the introduction of a new velocity scale, U t , based on the surface-layer-average Langmuir number: where N is the local buoyancy frequency, w s is the turbulent velocity scale, C υ is a defined parameter based on buoyancy, u* is the velocity scale of pure shear turbulence, w * is the velocity scale of convective turbulence, Ri c is the critical Richardson number, La SL is the surface-layer-averaged Langmuir number: and u s SL is the surface-layer-averaged Stokes drift speed, which is the average speed of a specific fluid parcel over multiple wave periods. While the Stokes drift is not calculated by the model, it is dependent on the high-frequency part of the wave spectrum, which is reasonably well described by the empirical wave spectra. The wave spectra used to calculate the Stokes drift are based on the Phillips spectrum (Phillips, 1985), and the directional spreading of wind waves of Webb and Fox-Kemper (2015) is also considered. In the calculation of wave spectra, the peak frequency is often assumed to be related to the wind speed, so that the wind speed at 10-m height is used for this calculation.
The second approach tested in this study, based on the work of , represents Langmuir turbulence as an enhancement factor applied to the turbulent velocity scale used in the KPP formulation. The idea of an enhancement factor had already been presented by McWilliams and Sullivan (2000), whose enhancement factor is a function of La t . However, the formulation proposed by McWilliams and Sullivan (2000) introduces spurious extra mixing in the extra-tropical regions, possibly because it does not include wind-wave misalignment. The parameterization of  uses a theoretical wave model based on the Stokes drift profile, the boundary layer depth, and the surface friction velocity (see equation 25 of : where u 10 is the wind velocity at 10-m height. The boundary layer depth and surface friction velocity are provided by the model, and the Stokes drift is calculated as previously explained. Both parameterizations of LC depend on wind speed and assume direct work of wind on the ocean surface. However, the presence of sea ice modifies ocean surface waves and the air-ocean drag, which in turn modifies the surface momentum exchange (Andreas et al., 1984;Birnbaum & Lüpkes, 2002). The efficiency of wind in inducing Stoke drift thus varies with SIC, but the relationship between wind forcing and Stokes drift in the presence of sea ice has not been quantified. In this study, therefore, we work with two different assumptions. In the model, ERA-Interim wind forcing is used to calculate the surface Stokes drift and the Stokes transport. To account for the presence of sea ice in the wind work, two different approaches were tested, both consisting of a coefficient to be multiplied by the wind speed at 10 m before it is used in the calculation of surface Stokes drift and transport.
The first approach consists of a modified version of the relationship between SIC and Cd described in Andreas et al. (2010), which does not depend on wind. The study by Andreas et al. (2010) uses data collected in the Arctic over a yearlong experiment and proposes a relationship between SIC and Cd, the drag coefficient of the friction between 10-m wind and ocean: It is seen from this equation ( Figure S1 in the supporting information) that Cd is enhanced at intermediate SIC. This can be explained by the fact that at low to intermediate SIC ice floes will have higher pressure in the upwind and lower pressure in the downwind face, essentially functioning like a sail, which enhances the wind work. At high SIC, ice floes shelter each other and decrease the momentum transfer from air to the ocean. Since frictional velocity is proportional to Cd, we assume that wind effect on Stokes velocity will work in a similar manner.
Therefore, we derived an equation for an enhancement factor (s c ) dependent on SIC to be multiplied by the wind velocity, U 10 , and account for increased momentum transfer at intermediate SIC: so that s c ( Figure S1) equals 1 in the absence of sea ice (no enhancement), peaks at intermediate concentrations and at a rate similar to the Cd enhancement (~35%), and is 0 when SIC approaches 1 (no Langmuir influence in areas completely covered by sea ice due to absence of ocean surface waves).
There are great uncertainties in this relationship between wind influence and SIC, and the uncertainties are, at the moment, impossible to be clarified due to the lack of data. It is possible that this curve would lead to overestimated LC effect on the MLD at mid-SIC concentrations due to the enhanced effect of wind on Stokes velocity.
A second approach assuming that the influence of wind linearly decreases with increasing SIC was implemented to mitigate the potential overestimation of LC, with s c being defined as This approach satisfies the conditions that s c equals 1 in the absence of sea ice and 0 when SIC equals 1. Although none of these approaches are verifiable at this moment due to the lack of data, the comparison among them provides a sense of the locations and times that the MLD is more sensitive to the effect of sea ice on the Stokes drift.

Description of Experiments
Hindcast simulations were integrated between 1984 and 2014, with 1984 to 1990 considered spin-up. Therefore, the model period analyzed covers the period from the beginning of the Palmer-LTER field measurements (1991 onwards, see section 3) until 2014. The model output consists of daily averages of the ocean 10.1029/2020JC016109

Journal of Geophysical Research: Oceans
and sea ice variables. The study examines five simulations with different parameterizations for LC and different values for drag coefficient between air and sea ice, Cd (Table S1 in the supporting information).
The value of Cd is linked to the characteristics of sea ice and shows a wide range of values in both observational and modeling studies, ranging from 5 × 10 −4 to >3 × 10 −3 (e.g., Castellani et al., 2014;Miller et al., 2006;Tsamados et al., 2014;Wamser & Martinson, 1993). The value used in Regan et al. (2018) was 5 × 10 −4 , on the lower end of Cd estimates. The choice of using a low Cd value was intended to adjust the formation of polynyas around Eltanin Bay (south of the Palmer-LTER region) and the freshwater fluxes around George VI Ice Shelf. Lowering the drag coefficient also lowers the impact of wind in transporting sea ice to different regions, hence curbing the formation of polynyas. For the Palmer-LTER region, however, a lower drag coefficient means that more sea ice could be melted locally rather than be transported elsewhere, thus lowering the surface salinity and increasing stability in coastal and shelf regions. Therefore, Cd values of 5 × 10 −4 and 2 × 10 −3 (on the higher end of the observed values) were used in different simulations to assess the impact of sea ice transport in local surface ocean mixing and freshwater budget.
In the control simulation (CTRL), Cd ¼ 2 × 10 −3 , s c is calculated using Equation 7, the parameterization for LC is based on the entrainment buoyancy flux , and only includes the parameterization described in Equation 3. Sensitivity simulation LowCd has a similar configuration to CTRL, but with Cd ¼ 5 × 10 −4 . All of the simulations have an air-ocean drag coefficient of 10 −3 , and the drag coefficient between ice and ocean multiplied by ocean density is 5.5. To assess the effect of different parameterizations for LC, sensitivity simulations without any LC parameterization (NoLC) and with the LC parameterization based on an enhancement factor for wind mixing (EnhLC;  are carried out. All other parameters in NoLC and EnhLC are the same as in CTRL. The difference between the two approaches for s c is investigated by comparing CTRL and simulation Linear_sc, which calculates s c using Equation 8. A table summarizing the different coefficients and parameterizations used in each simulation was added to the supplementary information (Table S1).

Data Used for Skill Assessment
The Palmer-LTER project started in 1991 and has been collecting a range of physical, biological, and chemical data along the WAP since, including semiweekly water-column sampling from Palmer Station from October through March each year and an annual cruise each January-February since 1993 . The data collected during the cruises include CTD (conductivity, temperature and depth) casts, chlorophyll-a concentration, zooplankton abundance, dissolved inorganic carbon and alkalinity, and nutrients. Before 1999, temperature and salinity data were obtained from a SeaBird SEACAT system, and since then the CTD profiles are collected with a dual-pumped Sea-Bird 911+ CTD system. CTD sensors are calibrated by SeaBird before and after each cruise, with standard post-cruise calibration corrections applied to the data. More details about the temperature and salinity data from the Palmer-LTER cruises are described by Martinson et al. (2008). While up to 2008 each cruise sampled the whole Palmer-LTER grid, since 2009 sampling has been limited to one coastal, one shelf, and one slope station in each line, generally complemented with casts performed in process-study cruise locations.
While Palmer Station has increased temporal coverage compared to the cruise data, it is located at the coast and is influenced by the surrounding islands and continental synoptic-scale phenomena that will not be properly simulated by the regional model. Although local phenomena such as the passage of icebergs affect CTD coastal stations sampled during the Palmer-LTER cruises as well, the cruise data are more likely to reflect seasonal and interannual variations that can be captured by the model. Therefore, the simulations are compared to the cruise data but not with the Palmer Station data.
The cruise data used for model validation are from lines 200 to 600 of the Palmer-LTER grid, spanning an area of 500 km along the coast and 250 km across the shelf. The cross-shelf transects are spaced every 100 km in the along-shelf direction, and the stations on each transect are separated by 20 km (Figure 1b). All the data collected are available through the project web page (http://pal.lter.edu/data/). Modeled vertical profiles of salinity and temperature were compared to the cruise CTD. Note, however, that not all CTD stations are sampled every year. As noted above, it is expected that CTD coastal stations are more influenced by continental runoff and phenomena at the synoptic scale that would lead to higher variability that could not be captured by the model, while shelf and slope regions are more representative of oceanic conditions that are less likely to change at short time-scales. The coastal region is less frequently sampled, and interannual variability of MLD in these subregions should be interpreted with caution. The calculated seasonal climatology of MLD, however, is thought to accurately represent the mean state of each region.
In order to decrease the influence of short-term and small-scale processes on the CTD casts during the analyses, we divided the cruise data into north (lines 400, 500, and 600) and south (lines 200 and 300) regions (Figure 1b), and also into coastal, shelf, and slope regions, based on topography. This way, six different subregions are obtained: north coastal (n_cs), north shelf (n_sh), north slope (n_sl), south coastal (s_cs), south shelf (s_sh), and south slope (s_sl).
Onset time of sea ice advance and retreat were calculated using data from the GSFC (Goddard Space Flight Center) Bootstrap algorithm (Comiso, 2017). In this study, following Stammerjohn et al. (2008Stammerjohn et al. ( , 2012 the sea ice year starts in mid-February, from day 46 of the calendar year until day 410. SIC refers to the percentage of each cell area covered by sea ice, ranging from zero (no sea ice cover) to one (totally covered by sea ice). Onset time of sea ice advance was chosen to be the day in which the average SIC of all the stations in each Palmer-LTER grid subregion reaches at least 0.15 for 5 days in a row. The day of advance will be the last on the 5-day series. Onset time of sea ice retreat was calculated as the last day in which SIC is less than 0.15 for the sea ice season. Further details are provided in Stammerjohn et al. (2008Stammerjohn et al. ( , 2012. The sea ice data used to calculate the climatological SIC for each subregion are from the same source, and binned as 8-day means with horizontal resolution of 0.2°in both latitude and longitude. For surface freshwater content, oxygen isotope data, described in Meredith et al. (2008Meredith et al. ( , 2013Meredith et al. ( , 2017, were the primary dataset used for the model skill assessment.

Results
To compare the model results with the Palmer-LTER cruise data, values at model grid points closest to Palmer-LTER stations were analyzed. In the WAP region, temperature is generally very low and changes over a small range, and salinity is the main driver of density changes. Any definition of MLD based on temperature changes is likely to not properly represent stratification, so we chose to use density gradients to compute MLD. We tested different thresholds to define MLD, ranging from a density change of 0.01 to 0.05 kg/m 3 , but due to the low vertical resolution of the model, the differences between the different thresholds were minimal. Following Mitchell and Holm-Hansen (1991), we use the threshold of density change of 0.05 kg/m 3 over a 5-m interval to define surface MLD. Onset times of sea ice advance and retreat in the model data were calculated as per the observational data.
As a first assessment of the skill of the model in reproducing the broad seasonal patterns observed in the surface ocean in the WAP region, we examined summer (December to February) and winter (June to August) climatological surface salinity, temperature, currents, and SIC in simulation CTRL ( Figure S2). The main circulation pattern in the WAP region is a northeastward current of the ACC along the shelf break, which is present in the model results. The horizontal resolution of the model is not sufficient to properly resolve the seasonal southward-flowing Antarctic Peninsula Coastal Current (APCC), which flows in a narrow region along the coast of the WAP (Moffat et al., 2008). The general cross-shelf patterns described by Moffat and Meredith (2018) are represented in the results, such as the summertime sea surface temperature gradient with colder waters at places of recent sea ice melt and the latitudinal gradient with colder waters towards the southern part of the domain. During the winter, sea surface temperature values reflect the sea ice freezing temperature near the coast and shelf region in most of the area.
Simulation CTRL captures an onshore-offshore gradient in summer surface salinity with lower salinities reaching minimums of around 32.5 near the coast, as described by Meredith et al. (2017). While a salinity gradient is present in the winter in the southern continental shelf area of the model domain, it is not as marked as during the summer months, when sea ice melt and increased glacial input contribute to lowering the salinity near the coast, particularly in Marguerite Bay and south of there. Sea ice extends past the shelf break during the winter and is confined to the region south of Marguerite Bay during the summer. Although the model resolution is not high enough to capture the flow of the APCC, the main onshore-offshore and north-south gradients of surface salinity due to sea ice and glacial inputs described by Meredith et al. (2017) are seen in the CTRL results.
The comparison between CTRL and Linear_sc shows how different parameterizations of sea ice influence on wind-driven Stokes drift will affect the sea ice and MLD in the region. The results show that both the climatology and interannual variability of sea ice are similar in the two simulations (Figures 2 and S3). A difference is seen in the simulated MLD, with CTRL showing slightly deeper MLD than Linear_sc between the months of July and November in the shelf and coastal regions. The largest proportional difference is in the northern shelf region during the month of November, when Linear_sc MLD is 6.4% shallower than the CTRL MLD ( Figure 4). As expected, the influence of the parameterization is larger during the period of sea ice retreat, while summer climatology and interannual variability of the MLD are similar between Linear_sc and CTRL. It is therefore important to note that CTRL could overestimate the influence of LC during sea ice retreat and that more data are needed to properly understand how sea ice influences the Stokes drift. Nonetheless, despite the differences between these two simulations, the comparisons with Palmer-LTER data during the summer and analyses of interannual variability are not affected by the parameterization of s c .

Sea ice advance and retreat
To assess how well the model captures sea ice interannual variability, modeled onset times of sea ice advance and retreat over the simulation period are plotted against the observations (Figure 2). Correlation coefficients between modeled and observed onset times of sea ice advance and retreat for each subregion were calculated (Table 1). Sea ice advance is well represented in the model in the southern region, while in the north the model shows earlier advance in the coastal and shelf areas. For sea ice retreat, the CTRL run performs better in the shelf and slope regions than in the coastal regions. In the northern coastal region, the model shows earlier retreat than the observations in most years. In the southern coast the model is able to capture the timing of sea ice retreat observed; however, years of early (2006,2007,2010) and late (2004) retreat are not well represented. Note. All correlations, with the exception of sea ice retreat of LowCd at s_sh, are significant at the 95% confidence level. The location of stations are shown in Figure 1b.
Wind action has a larger influence in the timing of sea ice advance and retreat than LC, evidenced by the larger difference between the CTRL and LowCd results compared to the difference between CTRL and EnhLC/noLC. The comparisons between simulations suggest that wind action has a larger influence on the retreat of sea ice than on the advance. This likely reflects the fact that sea ice advance is driven more by thermodynamics rather than mechanical forcing. The influence of LC on sea ice advance/retreat is much more modest and mainly consists of small differences in the onset time of sea ice advance. Without any parameterization for LC (NoLC), sea ice advances slightly earlier in the coastal areas, and in the slope region in some years.
Interannual variability of the timing of sea ice advance and retreat also seems to be more influenced by wind action than by representation of LC. The correlations between modeled and observed sea ice advance and retreat show that the CTRL simulation has better representation of interannual variability, with correlation coefficients mostly higher than 0.8, with the exception of sea ice advance in the northern coast (0.62) and shelf (0.70) regions. The coefficients calculated for experiments EnhLC and NoLC are similar to those of the CTRL run, and LowCd has the lowest correlations for both sea ice advance and retreat.
In summary, the choice of wind-sea ice drag coefficient has a larger impact on the representation of seasonal and interannual sea ice cycles than the inclusion of LC, with low drag leading to higher SIC and delayed sea ice retreat ( Figure S3). While sea ice retreat seems to be highly affected by wind influence on sea ice, differences in sea ice advance do not seem to be related to wind in all regions except for northern coastal and are likely more affected by changes in surface salinity and temperature between the simulations. Biases in atmospheric forcing (particularly wind forcing) might play a role in the disparities found between the models and observations in some years. Their exact contribution and the associated mechanism are unclear and left for future studies. However, given that the temporal and spatial coverage of in situ wind data is sparse, the magnitude of the wind biases in the reanalysis products remains unknown.

Surface Mixed Layer Depth (MLD) 4.2.1. Variability of Summer MLD
To compare the model results to the Palmer-LTER CTD cruise data, simulated daily MLD in each subregion in the months of January and February were averaged for each year. The results were then compared to the observed mean MLD in each subregion for each cruise (Figure 3). The choice of using the modeled mean MLD in January and February was made because January and February are the months in which the cruises (which last about 6 weeks) take place. Comparing the results on individual dates would not necessarily improve the model-observation comparison, since small-scale phenomena such as passage of submesoscale features and icebergs or a locally strong wind might affect the MLD at each CTD individual station, and these phenomena are not captured by the model. However, with a long enough dataset, and using data from different CTD stations, the Palmer-LTER cruise data should be able to represent the interannual variability of MLD in the WAP to some degree at regional scale.
The CTRL simulation gives mean MLD shallower than observed throughout the domain. The difference between model and observations is particularly large in the southern coastal and northern shelf regions. Less variation is expected in the modeled MLD, given it is represented as a 2-month average, but it also shows lower standard deviation (calculated from daily values). Although the magnitude of the year-to year variations in simulated MLD are much lower than in the observations, the phase of the interannual

Journal of Geophysical Research: Oceans
variability in the CTRL simulation seems to be reasonably well represented in four of the six domains, with the highest correlation with observation for the southern shelf region (0.93; Table S2). Modeled MLD variability is less representative in the northern coastal (correlation ¼ 0.308) and southern slope (correlation ¼ 0.110) regions. The low correlation in the northern coast is expected as MLD in the coastal

Journal of Geophysical Research: Oceans
region is influenced by local phenomena that are not well captured by the model, but it is not clear why the southern slope fails to capture interannual variability. While CTRL simulated and satellite monthly climatological SICs show a similar pattern throughout the model domain (not shown), a possibility is that the southward advection of sea ice during the retreat is less pronounced in the model compared with observations, which would affect the amount of melt in the southern slope. Even if modeled sea ice area in the southern slope is similar to observations, it is possible that sea ice thickness there is overestimated in the model, which would cause the MLD to differ.
The CTRL simulation has the deepest MLD among all the sensitivity simulations. Despite having the same parameterization for LC as CTRL, LowCd MLD is mostly similar or shallower than NoLC and EnhLC, implying that properly representing wind-sea ice drag in the model has as much or more influence on MLD as properly representing LC. Between the two approaches to parameterizing LC, entrainment from the bottom of the mixed layer (CTRL) has a larger impact than applying an enhancement factor (EnhLC).
Overall, the CTRL simulation has better performance in simulating both the depth and interannual variability of the mixed layer, although a bias towards shallow MLD still persists. This indicates that both using a higher drag coefficient and including LC parameterization based on entrainment buoyancy flux are important for resolving the MLD, but properly simulating the sea ice cycle and capturing the redistribution and volume of sea ice melt have a bigger impact in reproducing the interannual variability of the MLD.

MLD Seasonality
Although seasonality of MLD is hard to validate due to lack of data during the cold months of the year, it is known that MLD deepens during the winter due to surface cooling, formation of sea ice, and subsequent brine rejection. Comparisons between the seasonal cycle of MLD in the different model simulations were made to investigate how the changes made to the model parameters affect MLD throughout the year. For that purpose, the monthly climatology of MLD for each subregion was plotted for all model simulations (Figure 4). The CTRL simulation has deeper MLD throughout the year and in all subregions, with maximum depth being reached in August and deeper MLD farther away from the coast.
LowCd has consistently shallower MLD, except for the month of August in the southern slope region. The difference between these two simulations is more accentuated during the spring and early summer, when LowCd overestimates SIC ( Figure S3). Shallower MLD in LowCd occurs even in areas and months in which sea ice is absent and both simulations would have similar physics. Regarding the different parameterizations of LC, applying the enhancement factor (EnhLC) only seems to deepen the mixed layer during the fall and winter in the slope regions, when the difference between EnhLC and noLC is visible but smaller than 5 m. The parameterization used in the CTRL simulation, based on entrainment from below the mixed layer, deepens the MLD throughout the domain with the exception of summer and fall in the northern coast. From the results, we conclude that increased SIC during the spring and summer will shoal the MLD during these seasons, and in this case a shallower MLD can also be seen during fall and winter in most regions. Properly simulating sea ice retreat, therefore, is important for the simulated MLD year-round.

Vertical Structure of the Water Column
The water mass distribution in the surface and subsurface layers of the WAP are marked by the presence of Antarctic Surface Water (AASW), Winter Water (WW), and Circumpolar Deep Water (CDW), respectively, and have been described in several studies (e.g., Klinck et al., 2004;Martinson et al., 2008;Moffat & Meredith, 2018). CDW is a subsurface water mass derived from the ACC, which flows close to the shelf break at the WAP, sporadically intruding onto the shelf. CDW is comparatively warm, with temperatures higher than 1°C, and nutrient rich, playing an important role for primary production. CDW is modified as it flows across the shelf, being generally cooler in the near-shore regions. Above the CDW is the WW, which is formed by deep mixing from brine rejection during winter sea ice production, which then becomes isolated (forming "remnant WW") from the surface by summer stratification. The summer surface stratification is due to increased heat input from solar radiation and freshening during sea ice melt and glacial discharge from land, forming the AASW, which consequently traps the remnant WW between the AASW and the modified CDW, and is represented by a cold layer (temperatures below 1°C) located between 80 and 100 meters.

Journal of Geophysical Research: Oceans
The mean January/February salinity and temperature profiles of each model experiment were plotted for each subregion along with the average Palmer-LTER cruise profile ( Figure 5 for the northern region and Figure 6 for the southern region). The August averages for each experiment were also plotted to represent the winter water column (Figures S4 and S5). The vertical profiles of the CTRL simulation have lower surface salinity and higher surface temperature when compared with the summer observations, which is consistent with late and increased sea ice melt. An increased amount of meltwater would lower surface salinity and increase the stability of the surface layer. With a more stable surface layer, heat from the atmosphere is distributed over a shallower depth, increasing the surface temperature and further increasing the stability of the mixed layer. In the observations, a deeper mixed layer leads to lower surface temperature not only because surface heat is distributed deeper but also because more WW is entrained from below.
Compared to the CTRL simulation, LowCd has lower surface temperature and salinity during the summer, likely because sea ice retreat happens later in LowCd. New meltwater on the surface tends to be cold prior to gaining heat through insolation. Although LowCd brings the surface temperature closer to observations when compared to CTRL, this likely happens for the wrong reason: due to late sea ice retreat, instead of At the surface, the simulation with no LC shows higher temperature and lower salinity in the summer compared to CTRL, increasing the differences from the observations and confirming the positive impact of buoyancy flux entrainment-based LC parameterization on properly simulating the surface ocean layer. Using the enhancement factor to parameterize LC (EnhLC) has similar results to the simulation with noLC. The only difference observed between EnhLC and noLC during the summer is a slightly lower salinity close to the surface in EnhLC. Although EnhLC has larger impact on vertical salinity and temperature during the winter when compared to noLC, the effect of the enhancement factor is still smaller than the parameterization included in the CTRL simulation.
The differences observed in MLD and vertical profiles among the simulations suggest that the influence of late and increased sea ice melt in the coastal and shelf region during the summer persists throughout the year, with lower surface salinity and shallower mixed layer observed during the winter. On the slope, winter MLD is more impacted by the choice of LC parameterization than by the choice of drag coefficient.

Freshwater Seasonality
The freshwater cycle in the WAP has an important role for the coastal circulation (Moffat et al., 2008;Moffat & Meredith, 2018) and for primary production. Glacial meltwater, in particular, is thought to be a major source of dissolved iron (Annett et al., 2017), which is likely the limiting nutrient in the WAP given that nitrogen and silica are rarely depleted (Kim et al., 2016). Studies aimed at quantifying the amount and source of freshwater in the WAP surface ocean have been conducted using oxygen isotope data (Meredith et al., 20082017) to identify whether the freshwater is derived from sea ice melt or from meteoric sources, which is derived from the atmosphere and includes precipitation and glacial melt. Although sea ice formation and melt does not represent a net input of freshwater if integrated over a full season and over a large enough region, an excess of sea ice melt (formation) can be a source (sink) of freshwater at certain locations due to the redistribution of sea ice (i.e., by winds) from the location where in situ growth occurred.
Since the oxygen isotope samples were collected during the Palmer-LTER cruises in January and February, the spatial and temporal patterns of freshwater distribution derived from these data represent the summer distribution only. Given the model used in this study has a reasonable representation of the sea ice cycle ( Figure S3) and includes seasonality in the discharge of glacial freshwater, the results can be used to estimate the influence of different sources of freshwater throughout the year. However, because there is no interannual variability in the glacial runoff in the model, the analysis is restricted to the seasonal climatology of freshwater.
The cruise observations in late summer show that there is a cross-shelf gradient in salinity, with lower salinities around 33.2 near the coast coinciding with higher concentrations of meteoric water . Although meteoric water has a highly variable and heterogeneous spatial distribution, cruise measurements indicate that it is the dominant source of freshwater across the shelf in most years. Sea ice melt concentrations in the near surface are highly variable from year to year . Despite the high variability, the data also suggest that there is more sea ice formation towards the north and more sea ice melt towards the south of the WAP. Sea ice formation, however, leads to deepening of the mixed layer, and the isotopic signature associated with this process is also mixed downward in the water column during this process. Accordingly, there is a seasonal asymmetry in the vertical distribution of isotope-derived freshwater content, which adds complexity to the interpretation of these data .

Journal of Geophysical Research: Oceans
It is seen that although a strong seasonal cycle is included in the glacial runoff, with increased discharge during the summer, when precipitation is added to the glacial runoff, the overall seasonality of meteoric water is not obvious, since the prescribed glacial runoff and precipitation offset each other in the summer and winter, respectively. The CTRL results for sea ice melt indicate that it is the dominant freshwater source between October and December in some subregions. Towards the end of the summer (January and February), as sea ice meltwater approaches zero, meteoric water becomes the dominant source of freshwater. For instance, in regions n_sh and s_sh (Figures 7c and 7d), the meteoric water fluxes in January are 0.09 and 0.065 m/month, respectively, while the sea ice meltwater flux is zero. The prevalence of meteoric water over sea ice meltwater in late summer is in accordance with observations .
CTRL results also show that most of the sea ice formation happens in the coastal region during the winter. Despite melt seasonality being different for each subregion, the overall summer melt, as indicated by the sum of sea ice melt from October to February, is larger in the southern region for coastal (0.6 m north, 0.7 m south), shelf (0.47 m north, 0.53 m south), and slope (0.23 m north, 0.37 m south) stations, consistent with the observational result of more sea ice melt in the southern region.
Sea ice melt happens first in the northern region and in the slope, then progresses towards the south and coast. While this progression is observed in the CTRL simulation, LowCd sea ice melt peaks in December in every subregion with the exception of the northern slope, where melt peaks in November. The largest differences in sea ice melt between CTRL and LowCd are found at the coastal stations in January, where the largest differences in surface temperature and salinity are also observed between these two simulations. The increased melt in the LowCd simulation during the summer confirms the hypothesis that the lower surface temperatures seen in the LowCd experiment (compared to CTRL) are due to recently melted sea ice.
Maps of climatological summer (December to February), winter (June to August), and annual meteoric, net sea ice melt, and total freshwater were also plotted for the CTRL simulation ( Figure 8). They provide a broader assessment of the freshwater distribution throughout the WAP region. The observed onshore-offshore gradient in meteoric water, with higher values towards the coast due to glacial runoff, is present in the summer and also seen in the annual values (Figures 8b and 8e). The summer sea ice melt values are larger towards the coast and southern region, in agreement with the observations. These high-melt areas are also areas of high sea ice formation during the winter. The annual net sea ice melt values, then, suggest net sea ice formation near the coast and in Marguerite Bay, while higher net sea ice melt is seen in the northern shelf. Observations of sea ice melt suggest that meltwater content is higher in the southern region. However, because the sampling takes place in months of high sea ice melt in the southern region, it is possible that the influence of melt is overestimated there. Sea ice formation, in contrast, leads to brine rejection and more vertical mixing, which tends to suppress its signal on the surface. The influence of freezing in the surface could thus be underestimated in the observations.

Discussion and Conclusion
Different ocean-sea ice hindcast simulations were conducted for a region encompassing the WAP with different choices of wind-sea ice drag coefficient and parameterizations for LC to assess the effects of wind mixing and sea ice in the seasonal and interannual representation of MLD and freshwater cycles. Two different approaches to parameterizing LC were tested, based on the studies of  and . These approaches were applied to the KPP mixing parameterization and, respectively, include a LC-based enhancement term for wind mixing throughout the ocean surface mixed layer  and a change in the vertical velocity scale in the critical Richardson number to increase entrainment of subsurface water from below the mixed layer .
We found that the choice of LC parameterization does not significantly change simulated sea ice, but the choice of wind-sea ice drag coefficient leads to large differences in the representation of seasonal and interannual sea ice variability. The CTRL simulation, which has a higher drag coefficient (2 × 10 −3 ), shows a bias towards higher SIC and late sea ice retreat when compared to the observations. These biases are amplified in the simulation with low (5 × 10 −4 ) drag coefficient. The difference between simulations with different drag coefficients is much larger during the period of sea ice retreat than during sea ice advance, suggesting that retreat is more influenced by wind action. The smaller differences in SIC between models during sea ice advance can be at least partially explained by thermodynamic effects, given that lower MLD during the summer leads to higher temperature and low salinity, which could delay the formation of sea ice.
Both sea ice and LC impact the MLD in the WAP. The influence of sea ice on surface MLD is larger in the coastal and shelf regions, while the importance of LC for the MLD is larger than that of sea ice close to the outer sea ice margins (slope) during the winter. Shallower spring and summer MLD in the experiment with low drag coefficient (compared to CTRL) leads to shallower MLD in the fall and winter despite this experiment having similar physics to the CTRL simulation in the ice-free months, suggesting a lasting effect of improperly simulating SIC and seasonality.
A bias towards MLD being shallower than that observed has been seen in several modeling studies. For the LC parameterization, the mixed-layer entrainment approach proposed by  proved to be more effective in deepening the MLD than using an enhancement factor for wind mixing. Although the CTRL simulation includes the parameterization of Li and Fox-Kemper (2017) and a higher wind-sea ice drag coefficient, the CTRL simulation still shows a bias towards shallow MLD compared to observations despite being the best performing simulation and being able to better capture interannual variability. While the interannual variability of MLD is reasonably simulated in most of the domain, the model fails to properly capture it in the northern coast and southern slope subregions. The northern coast is presumably highly influenced by subgrid-scale phenomena that are not fully reproduced by the model, and which can potentially explain the poor model performance there; it is unclear at the moment why the model fails to capture MLD variability in the southern slope.
Interannual variability of the freshwater cycle in the region cannot be studied from the model results given that glacial inputs are based on estimated climatological values. However, because the model reasonably represents sea ice and is thought to have realistic values for meteoric sources, it is possible to do an assessment of the climatological importance of the different freshwater sources in the region in different seasons. Overall, the model is able to reproduce the basic patterns encountered in the observations, such as an onshore-offshore gradient in meteoric freshwater and more sea ice melt in the coastal and southern region (compared to the northern part of the domain) in late summer. The regions of high sea ice melt in late summer are also regions of high sea ice formation during the winter, and annual sea ice melt values indicate net sea ice formation near the coast and in Marguerite Bay, with highest overall net sea ice melt happening in the northern shelf of the domain.
The biases towards shallower MLD could also be related to other limitations in the model configuration. A higher horizontal resolution could improve the model by better resolving eddy activity and coastal circulation; and adding tides could also have an impact in the vertical mixing, although Brearley et al. (2017) find that tidally driven mixing is weak in the WAP shelf. Although ERA-Interim is seen as the most reliable reanalysis over this area (Bracegirdle, 2013), Rodrigo et al. (2013) argue that spatial resolution in the ERA-Interim dataset is not high enough to capture synoptic high-wind events and observe that ERA-Interim underestimates wind speed in the ice shelf areas around Ross Ice Shelf. Lower wind speed could also lead to shallower MLD in the model results, as Brearley et al. (2017) have found that direct impact of storms on ocean surface is responsible for a large portion of the surface ocean mixing in the WAP shelf.
Some of the changes to the model that could potentially improve the results further rely on phenomena for which there is very little data, such as circulation in the coastal ocean and Marguerite Bay, interannual variability of glacial freshwater discharges, and sea ice volume (thickness). Although the largest sea ice thicknesses simulated by the model are of the order of 1 m, which is in accordance with the scarce observations, data describing spatial and temporal changes in sea ice volume are not available. More measurements on these properties would help quantify the freshwater cycle and its influence on surface water stratification. Since most data in the Southern Ocean are collected in areas of open ocean and during the summer, it is possible that there are biases in the initial conditions of ocean temperature and salinity used for the model.
Despite the biases in the model results, this study highlights the influence of sea ice dynamics in the surface ocean mixing and freshwater fluxes, and the importance of a proper choice of parameters when simulating sea ice. Although the model results suggest that using drag coefficients on the higher end of observed values leads to better results in this region of the WAP, it is important to consider that reanalysis products underestimate wind speed and therefore might be inducing the choice of higher drag coefficients to compensate for the model wind speed bias. We also conclude that both mixed layer physics and atmospheric products need to be improved to properly simulate the MLD at seasonal and interannual timescales.