Effects of Freezing Temperature Parameterization on Simulated Sea‐Ice Thickness Validated by MOSAiC Observations

Freezing temperature parameterization significantly impacts the heat balance at sea‐ice bottom and, consequently, the simulated sea‐ice thickness. Here, the single‐column model ICEPACK was used to investigate the impact of the freezing temperature parameterization on the simulated sea‐ice thermodynamic growth during the MOSAiC expedition from October 2019 to September 2020. It is shown that large model errors exist with the standard parameterization and that different formulations for calculating the freezing temperature impact the simulated sea‐ice thickness significantly. Considering the winter mixed layer temperature, a modified parameterization of the freezing point temperature based on Mushy scheme was developed. The mean absolute error (ratio) of simulating sea‐ice thickness for all buoys reduces from 7.4 cm (4.9%) with the “Millero” scheme, which performs the best among the existing schemes in the ICEPACK model, to 4.2 cm (2.9%) with the new developed scheme.


Introduction
The Arctic sea-ice plays a crucial role in the global climate system as it regulates the exchange of energy, water, momentum, and gases between the ocean and the atmosphere.The oceanic heat flux at the ice-ocean interface is a crucial factor in the ice's mass balance and regulates the ice's bottom and lateral growth/melt, which is a crucial variable for understanding the thermodynamics of the sea-ice (Maykut & Untersteiner, 1971;Stanton et al., 2012).
There is evidence that the influence of vertical oceanic heat flux on sea-ice thickness has increased in recent decades (Lin et al., 2022).Observations near the north pole have shown that basal melt doubled between 2008 and 2013 compared to its value from 2000 to 2005 (Perovich et al., 2014).Satellite records show a significant decline in the Arctic sea-ice cover in recent decades, and modeling studies indicate changes in basal growth and melt rates are the main contributor to this decline.Satellite data and models suggest that the increase in oceanic heat flux since the year 2000 reduces winter sea-ice growth through thermodynamic processes (Ricker et al., 2021).With a stand-alone sea-ice model, a sensitivity study reveals that basal melt accounted for more than two-thirds of the mean annual sea ice melting in the Arctic, followed by top melt, while lateral melt accounted for less than 10% (Tsamados et al., 2015).The Coupled Model Intercomparison Project Phase 6 (CMIP6) also confirmed that the two essential processes in modeling the annual cycle of the Arctic sea-ice are basal growth and melting processes.Summer basal melting contributes twice as much to the lost ice mass as surface melting, and the average annual ice thickness increases due to basal growth was four times larger than that due to frazil ice formation (Keen et al., 2021).
Various methods have been developed to calculate the oceanic heat flux into the ice, including direct measurements, the residual method, and the parameterized method (Lin & Zhao, 2019).Direct measurements require frequent high-precision measurements of temperature, salinity, and vertical velocity in the near-surface layer of the ocean beneath the drifting sea-ice, making them costly and difficult to obtain (McPhee, 1992;McPhee & Stanton, 1996;Peterson et al., 2017).The residual method is an indirect method that determines the oceanic heat flux as the residual of the conductive heat flux through the ice and the sea-ice growth/melt rate at the bottom of the ice.Several studies utilize arrays of autonomous ice mass balance buoys deployed in Arctic pack ice to measure sea-ice temperature and ice thickness to calculate oceanic heat flux in this way (Lei et al., 2022;McPhee & Untersteiner, 1982;Perovich et al., 1989;Provost et al., 2017).The parameterized method estimates the oceanic heat flux based on empirical relationships between the ocean and sea-ice variables.This method calculates the difference between the mixed-layer temperature and the freezing temperature, and the friction velocity between the ice and the ocean is statistically determined from the current velocity (McPhee, 1992;McPhee et al., 2008).The parameterization approach is often utilized in numerical models because it is less costly and don't always accompany field campaigns (Peterson et al., 2017).
The parametrizations of the oceanic heat flux, primarily modified by the freezing temperature, can significantly impact the simulation of sea-ice thermodynamic growth.Some freezing temperature parameterizations have been applied in previous studies, such as constant value, a linear function of salinity (Bitz & Lipscomb, 1999), a piecewise linear relation to salinity (Assur, 1960), and determination from smoothed free energy data (Millero, 1978).However, only a few published studies have attempted to evaluate the impact of different freezing temperature parameterizations on sea-ice thermodynamic simulation (Shi et al., 2021;Yu et al., 2022).
The Multidisciplinary Drifting Observatory for the Study of Arctic Climate (MOSAiC) expedition was conducted in the Arctic Ocean from 2019 to 2020 to study the complex interactions between the atmosphere, sea-ice, ocean, and ecosystem for one year.The main goal of the MOSAiC expedition is to improve our understanding of the rapidly changing Arctic climate system and its global implications (Nicolaus et al., 2022;Rabe et al., 2022;Shupe et al., 2022).In this study, MOSAiC observations are used to force and validate the single-column model ICEPACK.We aim to evaluate the impact of the different methods for calculating the freezing temperature on the simulated sea-ice thermodynamic growth.Section 2 describes the data and methods we employ.The simulations of sea-ice thickness at the sites of specific buoys are presented in Section 3.1.Section 3.2 shows comparisons between different freezing-temperature parameterizations in the simulations of all buoys.Section 4 presents conclusions and discussions.

Observation Data
The observations used include atmospheric, oceanic, and ice thickness data from the MOSAiC expedition, all of which are distributed approximately within a range of 15 km (Bliss et al., 2023).The atmospheric data were obtained from the Met City observation site and the three Atmospheric Surface Flux Stations (ASFS).These data are available at one-minute intervals and include variables such as 2 m air temperature, 2 m humidity, 10 m wind speed and direction, atmospheric surface pressure, surface radiation (longwave and shortwave downward radiation), and precipitation.The atmospheric data has undergone careful quality control and processing to ensure accuracy and reliability (Shupe et al., 2022).We employed the average values from Met City and three ASFS as the atmospheric forcing for the ICEPACK.The atmospheric variables included: surface downward shortwave radiation (R sd ), downward longwave radiation (R ld ), 2 m temperature (T 2m ), 2 m specific humidity (Q a ), and 10 m wind speed (U a ).The 2 m wind speed of ASFS is extrapolated to 10 m using a power-law formula with a sea-ice surface roughness length of 1.0 × 10 3 m used in ECMWF-IFS model (Ruti et al., 2008).The oceanic data consist of 10 m sea temperature and salinity obtained from 8 CTD buoys, and 3 m sea temperature, salinity and current velocity obtained from 4 autonomous ocean flux buoys.The average and standard deviation of oceanic data are represented in Figure S1 in Supporting Information S1.All ocean observations are within the mixed layer depth (Rabe et al., 2022).The average of these observed ocean data will be used as ocean forcing to simulate the sea-ice thickness of all buoys.Ice thickness and snow depth are obtained from snow and ice mass balance array (SIMBA).SIMBA buoys are thermistor ice mass balance buoys.These buoys measure a vertical temperature profile through the atmosphere, snow, ice, and ocean from which snow depth and ice thickness and their evolution overtime may be determined.The measurement accuracy of SIMBA is 0.1°C for temperature and 0.02 m for ice thickness and snow depth (Provost et al., 2017).The previous study (Lei et al., 2022) analyzed data from 23 SIMBA buoys during the MOSAiC expedition and investigated the seasonality and timing of the sea-ice mass balance in the Arctic transpolar drift.In our study, 17 buoys were selected based on their operational period, with snow-ice formation excluded.This is because snow depth is prescribed to observed values from SIMBA buoys throughout the simulation in our setup of the sea-ice single-column model ICEPACK, and is uniform for each sub-grid ice class.The ice thickness data obtained from the SIMBA buoy was used to validate the ICEPACK simulation results.SIMBA buoy T63 is used as an example to investigate the impact of freezing temperature parameterization on the sea-ice thickness thermodynamic growth in Section 3.1, because it survived for a longer period of time, including winter, spring, and summer.

ICEPACK
The multi-ice-category single-column model ICEPACK is a state-of-the-art numerical model maintained by the Los Alamos Sea-ice Model (CICE) Consortium.It includes column-based physical processes that affect the seaice thickness within individual grid cell without reference adjacent grid cells.Ice categories are described for subgrid cells by dividing the ICEPACK in a grid cell into discrete thickness categories.The model solves a set of equations representing ice dynamics, thermodynamics, and mechanical redistribution (ridging).The ICEPACK version 1.3.2 was used in this study (Hunke et al., 2022).It incorporates atmospheric and oceanic forcing data as inputs to drive the model simulations.We focus on the influence of parameterization of freezing point temperature on the thermodynamic growth of sea-ice, closing ice dynamics and mechanical redistribution scheme.Model settings, physical parameterizations, and distribution of ice thickness categories for ICEPACK are the same as in a previous study (Gu et al., 2022).In this study, observed snow depth is prescribed to observed values from SIMBA buoys throughout the simulation in ICEPACK to mitigate the impact of simulated snow depth biases on ice thickness, which means that snow depth thermodynamic processes, including snow depth growth and melt, snow ice formation, snow compaction, wind-blown snow, etc., are not considered in this study.

Ice-Ocean Heat Flux Parameterization
The growth and melting of sea-ice at the ice-ocean interface are governed by an imbalance of the conductive heat flux at the sea-ice bottom and the oceanic heat flux into the ice (Maykut & McPhee, 1995).
Where f bot is the oceanic heat flux into the ice (downward is positive), ρ w is the density of water, c w is the specific heat capacity.T w is the mixed layer temperature under the ice, and T f is the freezing temperature of sea water.c h = 0.006 is the turbulent heat exchange coefficient (Maykut & McPhee, 1995).u * represent the friction velocity and it is limited to be larger than u min = 5 × 10 4 m s 1 in ICEPACK.The ocean current under the ice has a significant impact on friction velocity (McPhee, 1992).However, ICEPACK uses friction velocity as a fixed value to calculate oceanic heat flux.In this study, we implemented a quadratic drag law in ICEPACK to connect the friction velocity and the ocean current (Jenkins et al., 2010).
Where c w = 0.0055 is a dimensionless ice-ocean drag coefficient, and U is the ice velocity relative to the ocean surface velocity.
There are many different approaches to determining the freezing temperature in ICEPACK.The first approach sets T f constant (hereafter referred to as Const): The second approach, following Bitz and Lipscomb (1999), defines T f as a linear function of salinity (hereafter referred to as Linear) Geophysical Research Letters Where S mix is mixed layer salinity.
The third approach determines T f from smoothed free energy data of Doherty and Kester (1974) (Millero (1978), hereafter referred to as Millero): In the fourth approach, following Assur (1960) and Hunke et al. (2022), the freezing temperature is a piecewise linear function of Z (the ratio of the mass of salt (in g) to the mass of pure water (in kg) in brine, hereafter referred to as Mushy): Salinity is the mass of salt (in g) per mass of brine (in kg), so that Combining Equations 6 and 7, we can calculate the freezing temperature by, Where with empirically determined parameters, SSS 0 = 123.667g kg 1 , a 1 = 18.48 g kg 1 K 1 , a 2 = 10.3085g kg 1 K 1 , b 1 = 0 and b 2 = 62.4 g kg 1 .
Mushy is the most comprehensive thermodynamic scheme in ICEPACK (Feltham et al., 2006;Turner et al., 2013).However, the parameterization of Mushy's freezing point temperature underestimates the simulation of sea-ice thickness (refer to Section 3.1 and 3.2 for detailed results).We further evaluated the performance of Mushy parameterization during MOSAiC and found that the difference between mixed layer temperature and freezing point temperature calculated by Mushy was an average of 0.0328°C in winter.However, previous study (Lin & Zhao, 2019;Zhong et al., 2022) found that the ocean freezing temperature of Arctic should be close to mixed layer temperature in winter.To develop a freezing temperature parameterization that is suitable for simulations of Arctic sea-ice, we modified the "Mushy" parameterization (hereafter referred to as "Mushy-m") based on the "Mushy."The "Mushy-m" scheme calculates the mean difference between the freezing temperature and mixed layer temperature within the "Mushy" scheme during the winter and subsequently incorporates this difference into the "Mushy" scheme: The constant offset (0.0328) in "Mushy-m" is the same to calculate the freezing point temperature for all buoys because we use the same mixed layer temperature, salinity, and current velocity as oceanic forcing.

Simulation at the Site SIMBA T63
Mean atmospheric and oceanic observation data are utilized to force ICEPACK.We use the mean absolute error (MAE, MAE = 1 n ∑|h sim h obs |, where n is the number of data point, h sim and h obs are ice thickness simulation and observation, respectively) and MAE ratio (ratio between the MAE and the observation value) as metrics for the evaluation of the simulations.All simulations with the various freezing temperature parameters described above capture the overall evolution of the observation (Figure 1).The sea-ice thickness grows from December to May and starts to decrease in June.When calculating the freezing temperature by using the "Const" parameterization, the simulation underestimates ice thickness in winter (Figure 1a).This is due to a large difference between the mixed layer temperature and the constant freezing temperature of 1.8°C (Figures 1b and 1c), increasing the oceanic heat flux into the ice during winter (Figure 1d), leading to a small sea-ice growth rate (Figure 1f) even though it increases conductive heat flux of ice basal layer (Figure 1e).After the winter, the simulated ice thickness increases rapidly, overestimating the observed ice thickness from mid-April to the end of simulation.In ICE-PACK, the mixed layer temperature should remind higher than the freezing temperature, and the oceanic heat flux remains close to zero from March to June (Figures 1b-1d).The growth rate of sea-ice bottom is related to the difference between conductive heat flux of ice bottom and oceanic heat flux (Figures 1d-1f).The sea-ice melt when the sea-ice growth rate is negative, and vice versa.The MAE for the "Const" simulation is 10.2 cm for the whole simulation period.When using the "Linear" and the "Millero" parameterizations, the simulation exceeds the observations, with MAEs of 17.7 and 12.0 cm, respectively, that is, these parametrizations perform worse than "Const".In the case of the "Mushy" parameterization, the simulation underestimates the ice growth, resulting in an MAE of 11.4 cm.The modified parameterization "Mushy-m" performs well with an MAE of only 4.7 cm."Mushy-m" slightly overestimates in winter due to lower oceanic heat flux but performs very well from mid-June onward.

Simulations of All Buoys
Finally, all available SIMBA buoys are utilized to evaluate the impact of freezing temperature parameterization on the simulations.Observed atmosphere and ocean data are employed to force ICEPACK on the individual buoy tracks.The parameterization scheme remains unchanged, but the initial ice thickness and prescribed snow depth are adapted to the condition at each buoy.According to Equation 1, the initial ice thickness and snow depth have no impact on oceanic heat flux and freezing temperature parameterization.Figure 2 shows the monthly mean biases ratio (ratio between the bias and the observation) between the simulations and observations.Initially, the mean biases ratio are small, but they increase toward the end of the simulation period.For all 17 buoys, the "Linear" and "Millero" parameterizations show an almost monotonically increasing positive mean bias over the simulation period, while "Mushy" shows an almost monotonically decreasing negative mean bias.The mean bias of the "Const" parameterization exhibits a negative trend before March 2020, approaching zero in May, followed by a positive trend until the end of the simulation period."Mushy-m" performs best, shows a small positive mean bias in summer, but agrees well with observations in winter and spring.
The MAE and MAE ratio between the simulated sea-ice thickness and the observations from individual buoys is shown in Figure S2 and Table S1 in Supporting Information S1.The largest MAE is 22.3 cm when simulating T70 using "Mushy.""Mushy-m" reduces the MAE significantly to 5.2 cm.For all buoys, "Mushy" and "Linear" exhibit the largest mean MAE of 12.2 and 11.3 cm, respectively."Mushy-m" significantly decreases the mean MAE to 4.2 cm."Millero" also exhibits a small mean MAE (7.4 cm), but it shows a larger simulation bias during spring and summer (Figure 2).The largest MAE of "Millero" is 14.7 cm for the simulation of the buoy T58.Overall, "Mushy-m" exhibits the smallest MAE for most buoys and presents the relative good performance for the buoys T47, T69, T73, T75, and T77.Two potential reasons may affect the performance of "Mushy-m": first, the spatial variability in the observed forcing may contribute to the error (see Figure S1 in Supporting Information S1).We have designed a set of numerical experiments to comprehend the impact of the spatial variability in the ocean forcing on the sea-ice simulation (details refer to Text S1).Second, the ice-ocean drag coefficient is set to a constant value in this study.The drag coefficient has a large uncertainty (Lu et al., 2011), and it significantly affects the friction velocity and the oceanic heat flux, as can be seen in Equations 1 and 2.
Different freezing point temperature parameterization schemes produce different oceanic heat flux and conductive heat flux (Figure 3).Sea-ice bottom growth results from the competition between the conductive heat flux at the ice basal layer and the heat flux from the ocean.Hence, sea-ice thickness simulation performance is jointly determined by the simulated oceanic heat flux and the conductive heat flux.The conductive heat flux changes from upward in winter and spring to downward in summer (Figure 3c), resulting in the difference between the conductive heat flux and the oceanic heat flux from negative to positive, which causes sea-ice to grow in winter and spring and melt in summer.Figure 3e presents the similar length of error bars for different freezing temperature parameterizations in each month.The length of error bars in Figure 3e means the variation in simulated sea-ice growth rate caused by the difference in snow depth and initial ice thickness between different buoys.This suggests that the effects of initial sea-ice thickness on sea-ice growth simulation for each freezing temperature parameterization are consistent.Compared to "Mushy," "Mushy-m" has an increased freezing point temperature throughout the simulation period, decreasing the difference between mixed layer temperature and freezing point temperature and decreasing the upward oceanic heat flux (Figures 3a and 3b).However, the heat loss at sea-ice bottom during growth season for "Mushy-m" is larger than that for "Mushy," while the heat gaining during sea-ice melting season for "Mushy-m" is smaller than that for "Mushy" (Figure 3d).This result indicates the faster sea-ice growth in winter and spring and the slower sea-ice melting in summer for "Mushy-m" compared to "Mushy."The sea-ice growth rate of "Mushy" is underestimated compared to observations in winter and spring, but the ice melting rate performs better in summer.Generally, the greater the difference between the conductive heat flux and oceanic heat flux, the greater the sea-ice growth rate.However, the difference between the conductive heat flux and oceanic heat flux in "Linear" and "Millero" is not the largest, but the rate of sea-ice growth is the largest.The main reason is that the first calculated freezing point temperature by "Linear" and "Millero" in winter is higher than the mixed layer temperature.The ICEPACK compulsively sets the freezing point temperature to be equal to the mixed layer temperature, and the excessive energy (i.e., the first calculated freezing point temperature minus the mixed layer temperature) can lead to the generation of new ice, promoting the growth of sea-ice thickness.Overall, compared to the observed sea-ice growth rate, "Mushy-m" exhibited the smallest average absolute error throughout the simulation period (Figure 3e).

Conclusions and Discussions
In this study, atmospheric and oceanic observations are used to force the sea-ice single-column model ICEPACK to simulate the sea-ice evolution.We prescribe the snow depth in ICEPACK to observed values from SIMBA buoys to reduce the impact of simulated snow depth biases on ice thickness.Additionally, a simple quadratic drag law is implemented in ICEPACK to account for the influence of oceanic currents on the ice simulation.For the T63 buoy, our results show that the model captures the overall evolution of the sea-ice throughout the entire simulation period.Freezing temperature parameterizations influence the sea-ice thickness simulations significantly.A modified freezing temperature parameterization, called "Mushy-m," is proposed.It reduces the MAE from 11.4 to 4.7 cm for T63 buoy compared to the conventional "Mushy" scheme implemented in ICEPACK.While for all buoys with initial sea ice thickness range of 0.3-2.8m, the "Millero" scheme performs the best among existing schemes with the MAE (ratio) of 7.4 cm (4.9%), and the MAE (ratio) of "Mushy-m" is reduced to 4.2 cm (2.9%).The "Mushy-m" simulations lead to a lower oceanic heat flux in summer, but show much better agreement with the sea-ice thickness observation in winter and spring.However, the constant offset in "Mushym" comes from the difference between "Mushy" freezing point temperature and mixed layer temperature in winter during MOSAIC.The offset need to be changed from site to site or simulation to simulation to make sure the freezing point temperature in "Mushy" is close to mixed layer temperature in winter.
There is still room for improvement.The ice-ocean drag coefficient plays a crucial role in determining the rate of heat transfer between the ocean and the sea-ice cover (Lu et al., 2011).In this study, the drag coefficient is set to a constant value of 0.0055 according to previous studies (McPhee, 2008;McPhee & Stanton, 1996).For our future studies, this value can potentially be replaced with measurements of turbulent stress.Furthermore, the snow accumulates on the surface of the sea-ice during winter, forming an insulating layer that reduces heat transfer from the ocean to the atmosphere.It also impacts the ice albedo (reflectivity), influencing the amount of solar radiation absorbed or reflected (Callaghan et al., 2011).In this study, observed snow depth is prescribed in ICEPACK to mitigate the impact of simulation bias of snow depth.In future studies, the snow depth should be actively simulated but for that better atmospheric precipitation data are needed.

Figure 1 .
Figure 1.Time series of (a) sea-ice thickness, (b) freezing temperature, (c) the different between mixed layer temperature and freezing temperature, (d) oceanic heat flux, (e) conductive heat flux of ice basal layer and (f) sea-ice growth rate.The black line in (a) and (b) represent ice thickness from SIMBA T63 and mixed layer temperature from CTD and autonomous ocean flux buoys buoys, respectively.Other color lines represent the simulation using different parameterization of freezing point temperature.

Figure 2 .
Figure 2. Mean bias ratio between the sea-ice thickness simulations and observations for each month.Each color represents the mean bias ratio using different parameterizations.The 17 buoys' simulated minimum or maximum bias is shown by error bars, and their mean is represented by points.The black dash line represents the zero line.

Figure 3 .
Figure 3.Time series of monthly average (a) freezing point temperature and mixed layer temperature, (b) mixed layer temperature minus freezing point temperature, (c) oceanic heat flux and conductive heat flux of ice basal layer, (d) conductive heat flux minus oceanic heat flux, and (e) sea-ice growth rate bias.Black color represents the observed value, while other colors represent parameterization of different freezing temperature.The 17 buoys' simulated minimum or maximum conductive heat flux and sea-ice growth rate bias for each month is shown by error bars, and their mean is represented by points.The black dash line represents the zero line.