The sensitivity of landfast sea ice to atmospheric forcing in single-column model simulations: a case study at Zhongshan Station, Antarctica

Single-column sea ice models are used to focus on the thermodynamic evolution of the ice. Generally, these models are forced by atmospheric reanalysis in the absence of atmospheric in situ observations. Here we assess the sea ice thickness simulated by a single-column model (ICEPACK) with in situ observations obtained off Zhongshan Station for the austral winter of 2016. In the reanalysis, the surface air temperature is about 1 ^∘C lower, the total precipitation is about 2 mm d⁻¹ greater, and the surface wind speed is about 2 m s⁻¹ higher compared to the in situ observations. We designed sensitivity experiments to evaluate the simulation bias in sea ice thickness due to the uncertainty in the individual atmospheric forcing variables. Our results show that the unrealistic precipitation in the reanalysis leads to a bias of 14.5 cm in sea ice thickness and 17.3 cm in snow depth. In addition, our data show that increasing snow depth works to gradually inhibit the growth of sea ice associated with thermal blanketing by the snow due to changing the vertical heat flux. Conversely, given suitable conditions, the sea ice thickness may grow suddenly when the snow load gives rise to flooding and leads to snow-ice formation. However, there are still uncertainties related to the model results because superimposed ice and snowdrift are not implemented in the version of the ice model used and because snow-ice formation might be overestimated at locations with landfast sea ice.

1 Introduction

Sea ice plays an essential role in the global climate system by reflecting solar radiation and regulating the heat, moisture, and gas exchanges between the ocean and the atmosphere. In contrast to the rapid decline of sea ice extent and volume in the Arctic (Stroeve et al., 2012; Lindsay and Schweiger, 2015), satellite observations have shown a slight increase in the yearly mean area of Antarctic sea ice since the late 1970s (Parkinson and Cavalieri, 2012) followed by a rapid decline from 2014 (Parkinson, 2019) and a renewed increase in the most recent years (Chemke and Polvani, 2020). Although the sudden decline in Antarctic sea ice is yet to be attributed (Parkinson, 2019), the spatial pattern of Antarctic sea ice changes is suggested to be primarily caused by changes in the atmospheric forcing. For example, the rapid ice retreat in the Weddell Sea from 2015 to 2017 has been associated with the intensification of northerly wind (Turner et al., 2017), while the phase of the southern annular mode (SAM) significantly modulates the sea ice in the Ross Sea and elsewhere, especially in November 2016 (Stuecker et al., 2017; Schlosser et al., 2018; G. Wang et al., 2019).

Landfast sea ice, the immobile fraction of the sea ice, is mainly located near coastal regions of Antarctica, and its change is assumed to be indicative of the evolution of total Antarctic sea ice (Heil et al., 1996; Heil, 2006; Lei et al., 2010; Q. Yang et al., 2016). Unlike drifting sea ice, the change in landfast sea ice is dominated by thermodynamic processes which single-column sea ice models can capture well (Heil et al., 1996; Lei et al., 2010; Y. Yang et al., 2016; Zhao et al., 2017; Liu et al., 2022). Furthermore, a single-column sea ice model is a useful tool to evaluate the impacts of different atmospheric forcings on the sea ice evolution because of the relatively simple structure of the physical processes (Cheng et al., 2013; C. Wang et al., 2019; Merkouriadi et al., 2020). In this study, a state-of-the-art single-column sea ice model, ICEPACK, is chosen to investigate the sensitivity of landfast sea ice to atmospheric forcing for the region off Zhongshan Station in Prydz Bay, East Antarctica (Fig. 1).

Due to the lack of in situ observations, the majority of sea ice studies, especially for the Antarctic, rely on numerical models. Realistic atmospheric forcing is critical for reliable model simulations. Although being criticized for significant deviations from in situ observations (Bromwich et al., 2007; Vancoppenolle et al., 2011; Wang et al., 2016; Barthélemy et al., 2018), atmospheric reanalysis data are assumed to offer reasonable atmospheric forcing for large-scale sea ice models for the Antarctic (Zhang, 2007; Massonnet et al., 2011; Zhang, 2014; Barthélemy et al., 2018). Previous studies reported a large spread between four global atmospheric reanalysis products and in situ observations in the Amundsen Sea Embayment (Jones et al., 2016). Moreover, studies showed that directly using atmospheric reanalysis as forcing for models causes significant biases in the Arctic sea ice simulations (Lindsay et al., 2014; C. Wang et al., 2019). Similar results, accentuated by the sparseness of atmospheric observations entering the reanalysis, can be foreseen for Antarctica. Therefore, the atmospheric forcing needs to be evaluated carefully before simulating Antarctic sea ice. To our knowledge, few studies have given a quantitative evaluation of the effect of different atmospheric forces on sea ice simulations in Antarctica.

The coastal landfast sea ice in Prydz Bay is generally first-year ice. It usually fractures and is exported or melts out completely between December and the following February, and refreeze occurs from late February onwards (Lei et al., 2010). This seasonal cycle is representative of Antarctic landfast sea ice. This study aims to evaluate the contributions of the various atmospheric forcing variables on landfast sea ice growth. The snow cover exerts influence on the evolution of the vertical sea ice–snow column via a number of mechanisms, including the formation of snow ice added by flooding (Leppäranta, 1983), superimposed ice (Kawamura et al., 1997), and insulating impact (Massom et al., 2001). Understanding the snow depth is a primary concern here.

Two sets of atmospheric forcing have been chosen. The first is spatially interpolated ERA5 onto the location of the observation site, and the second is using in situ atmospheric observations. It is well known that the simulation biases of numerical models are introduced through many shortcomings, including unrealistic surface boundary conditions (here: atmospheric forcing), imperfect physical process formulations, and computational errors. Understanding the uncertainty in sea ice simulations, as well as the sea ice response pattern to atmospheric forcing due to imperfect surface boundaries, is a prerequisite for successful simulations and needs to be assessed first.

This study is arranged as follows: the in situ observations, the numerical model, and the reanalysis are introduced in Sect. 2. The main results are given in Sect. 3, focusing on different kinds of atmospheric forcing on sea ice and snow. Shortcomings, discussions, and conclusions follow in Sects. 4, 5, and 6.

2 Materials and methods

2.1 Meteorological observations

The site of sea ice observations is located in the coastal area off Zhongshan Station (69^∘22^′ S, 76^∘22^′ E; Fig. 1), East Antarctica. The meteorological data were collected at a year-round manned weather observatory run at Zhongshan Station in 2016, which is 1 km inland from the sea ice observation site and 15 m above sea level. Snowfall is measured every 12 h at the Russian Progress II station (located ∼1 km to the southeast of Zhongshan Station). The short- and longwave radiation fluxes were measured every minute with a net radiometer mounted 1.5 m above the surface on a tripod (Q. Yang et al., 2016). Other meteorological variables are available as hourly data, including 2 m air temperature (T_2 m), surface pressure (P_a), specific humidity (calculated from dew-point temperature and P_a), potential temperature (calculated from T_2 m and P_a), air density (calculated by T_2 m and P_a), and 10 m wind speed (U₁₀) (Hao et al., 2019, 2020; Liu et al., 2020).

Figure 1 Location of landfast sea ice surface measurements near Zhongshan Station. The solid triangle denotes the observation site, and the solid circle marks Zhongshan Station. The color on the left represents the terrain.

2.2 Sea ice thickness measurement

A thermistor-chain unit developed by Taiyuan University of Technology (TY) was used to measure sea ice thickness in austral winter 2016. This unit is composed of two parts: the control unit and the thermistor chain. The controller initiates data acquisitions and records and stores the temperature measurements. The thermistor chain is 3 m long with 250 equidistant thermistors. Their sensitivity is 0.063 ^∘C, and the measurement accuracy is ±0.1 ^∘C. The thermistor chain simultaneously records the vertical temperature profile across the near-surface atmosphere, snow cover, sea ice, and surface seawater. The measurement frequency is hourly. Details about the instruments can be found in Hao et al. (2019).

Snow thickness close to the thermistor unit is measured weekly using a ruler with an accuracy of ±0.2 cm. Sea ice thickness is measured with a ruler through a drill hole (5 cm diameter) weekly. The measurement accuracy is ±0.5 cm. The average thickness obtained from three close-by sites is retained. Sea-surface temperature and sea-surface salinity are measured in the drill holes weekly using a Cond 3210 SET1 (Hao et al., 2019).

2.3 Atmospheric reanalysis data

The European Centre for Medium-range Weather Forecasts (ECMWF) released ERA5, the new reanalysis product, in 2017, which is updated in near real time (Hersbach and Dee, 2016; Hersbach et al., 2020). The complete ERA5 data set, extending back to 1950, has been available to the end of 2019 during this study. Compared with the popular ERA-Interim reanalysis, there are several significant improvements in ERA5, including much higher resolutions (both spatially and temporally). ERA5 has global coverage with a horizontal resolution of 31 km by 31 km at the Equator and 10 km by 31 km at the latitude of Zhongshan Station. The ERA5 resolves the vertical atmosphere profile using 137 vertical pressure levels from the surface up to a geopotential height of 0.01 hPa. ERA5 provides hourly analysis and forecast fields and applies a four-dimensional variational data assimilation system (4D-var). ERA5 includes various reprocessed quality-controlled data sets, for example, the reprocessed version of the Ocean and Sea Ice Satellite Application Facilities (OSI SAF) sea ice concentration (Hersbach and Dee, 2016; Hersbach et al., 2020).

For comparison and evaluation against the observations in this study, gridded data from ERA5 has been bilinearly interpolated to the observation site (detailed in Sect. 2.1). Directly using atmospheric forcing from coarse grid cells to interpolate to the observation site, although widely accepted in previous studies (e.g., Urraca et al., 2018; C. Wang et al., 2019), may cause errors. We have checked the performance of ERA5 and found that the spatial difference of surface atmospheric variables around the observation site is relatively small, indicating the choice of interpolation techniques will not affect the conclusion of this study.

2.4 ICEPACK

ICEPACK is a column-physics component of the Los Alamos Sea Ice Model (CICE) V6 and is maintained by the CICE Consortium. ICEPACK incorporates column-based physical processes that affect the area and thickness of sea ice. It includes several options for simulating sea ice thermodynamics, mechanical redistribution (ridging), and associated area and thickness changes. In addition, the model supports several tracers, including ice thickness, enthalpy, ice age, first-year ice area, deformed ice area and volume, melt ponds, and biogeochemistry (Hunke et al., 2019). ICEPACK Version 1.1.1 was used in this study, and detailed options of physical parameterizations and model settings for ICEPACK are summarized in Table 1. We employ ICEPACK to distribute the initial ice thickness to each ice thickness category using a distribution function:

$$\begin{array}{}\text{(1)}& p_i={\displaystyle \frac{max\left(\mathrm{2}\times h\times H_{\mathrm{i}}-H_{\mathrm{i}}^{\mathrm{2}},\mathrm{0}\right)}{{\sum }_{i}max\left(\mathrm{2}\times h\times H_{\mathrm{i}}-H_{\mathrm{i}}^{\mathrm{2}},\mathrm{0}\right)}}i=\mathrm{1}\mathrm{\cdots }N,\end{array}$$

where h is the initial ice thickness, H_i is the prescribed ice thickness category (0–0.6, 0.6–1.4, 1.4–2.4, 2.4–3.6, and above 3.6 m; same as for Arctic simulations), and N is the number of ice thickness categories.

Table 1 Detailed options of physical parameterizations and model settings for ICEPACK. SST signifies sea surface temperature.

The atmospheric forcing for the ICEPACK model consists of observations of downward short- and longwave radiation, 2 m air temperature, specific humidity, total precipitation, potential temperature, 2 m air density, and 10 m wind speed. The oceanic forcing includes sea surface temperature, sea surface salinity, and oceanic mixed layer depth. The period concerned in this study is from 22 April, when observed sea ice generally starts to grow, to 22 November 2016. Since there are no observations of the ocean's mixed-layer depth, we set it to 10 m based on a previously published study (Zhao et al., 2019).

3 Results

3.1 Surface atmospheric conditions near the observation site

First, we compare the eight atmospheric variables used to force ICEPACK (surface downward shortwave radiation (R_sd), surface downward longwave radiation (R_ld), surface air temperature (T_a), specific humidity (Q_a), precipitation (P), air potential temperature (Θ_a), air density (ρ_a), and wind speed (U_a) with the respective in situ observations). Table 2 lists the bias (reanalysis minus observation), bias ratio (ratio between the bias and the observation value), the mean value of the in situ observations (Mean_Obs), the correlation coefficient (Corr.), and the root-mean-square deviation (RMSD) between the interpolated ERA5 data and the observations. In general, all eight variables from the two sources closely follow each other (Corr. >0.85), except for P and U_a. In this study, the main attention is on the atmospheric variables T_a, P, and U_a for three reasons. (1) Previous studies (Cheng et al., 2008; Wang et al., 2015) have shown that from all atmospheric forcing variables, uncertainties in T_a, P, and U_a exert a significant impact on the sea ice thickness (Cheng et al., 2008). (2) Surface wind may affect the snow cover in two ways: sublimation due to surface turbulent heat flux (Fairall et al., 2003; Gascoin et al., 2013) and the snowdrift process (Thiery et al., 2012). (3) P and U_a from the reanalysis have the largest bias ratio compared to the in situ observations.

The timing of daily variations in T_a is well represented by ERA5, especially for strong cooling events (Fig. 2a). However, ERA5 tends to underestimate warm events by a few degrees, as well as cold events during which differences exceeding 10 ^∘C may occur (Fig. 2d). During the entire observation period in 2016, T_a from ERA5 was 1.2 ^∘C lower than the in situ observations. Also, previous studies reported similar disagreement in T_a between observations and reanalysis in Antarctica (Bracegirdle and Marshall, 2012; Fréville et al., 2014). The cold bias of T_a in the reanalysis was suggested to be caused by the ice surface schemes that cannot accurately describe the ice–atmosphere interactions of strongly stable stratified boundary layers that are frequent in Antarctica.

Table 2 Comparison of atmospheric forcing between ERA5 reanalysis and in situ observations.

The reanalyzed variable with the largest bias ratio from the observations is precipitation (Fig. 2b). Hourly precipitation from ERA5 was accumulated into daily data and compared with the nearest available daily precipitation records from the Progress II station. The maximum daily mean precipitation can reach 19.1 mm d⁻¹ (11 July 2016), with an average of 0.66 mm d⁻¹ from 29 April to 22 November 2016. While ERA5 captures the main precipitation events, it significantly overestimated the magnitude of precipitation events, especially in July. In this month, the mean precipitation rate from ERA5 is 5.83 mm d⁻¹, while the observed is only 1.42 mm d⁻¹. From April to November, the accumulated precipitation from ERA5 is about 300 % larger than that in the in situ observations. Nevertheless, using precipitation from Progress II for Zhongshan Station may be questioned because of the distance of about 1 km to Zhongshan Station. Moreover, the snowdrift due to strong surface wind can affect the precipitation observations and the local accumulated snow mass, which may further cause a significant bias in snow depth between simulation and observations.

Figure 2 Time series of daily (a) surface air temperature, (b) precipitation rate, and (c) wind speed (10 m above the surface). The ERA5 reanalysis data are indicated as red lines. Observations are marked by black lines. Panels (d–f) show the difference (marked by “D”) between ERA5 and the observations (ERA5-observation). The differences are marked by blue lines. The gray lines denote the zero line.

The observed U_a varied from 0.01 to 12.3 m s⁻¹ with an average of 4.2 m s⁻¹ (Fig. 2c). ERA5 captured well the daily and seasonal variation in U_a, but an overestimation of 2.1 m s⁻¹ should be noted, mainly when observed U_a>5 m s⁻¹. One explanation for such an overestimation is that the numerical model underlying ERA5 cannot represent the surface roughness and the katabatic wind in a region with complex orography (Tetzner et al., 2019; Vignon et al., 2019).

3.2 Simulation forced by observed in situ atmospheric variables

The simulation bias of sea ice thickness and snow depth is impacted by many aspects, including unrealistic atmospheric and oceanic forcing and shortcomings in the applied numerical model. In this study, we mainly focus on the influence of imperfect atmospheric forcing.

The sea ice thickness (Obs) measured through a drill hole increases from 29 April (100±2 cm) to 25 October (172±2 cm), remaining level from there on (Fig. 3a). The ice thickness deduced from the TY (Obs_TY) thermistor-chain buoy shows a similar result: sea ice thickness increased from 106 cm on 22 April to 171 cm on 17 November. In November, the sea ice thickness (Obs and Obs_TY) is stationary, indicating a thermodynamic equilibrium between heat loss to the atmosphere and heat gain from the ocean (Q. Yang et al., 2016; Hao et al., 2019).

When forced by atmospheric in situ observations (Sim_Obs), the simulated sea ice thickness agrees well with the observed thickness with a mean bias of less than 1 cm over the growing season. We attribute the excellent simulation result to the fact that the seasonal evolution of landfast sea ice is driven mainly by thermal processes, which ICEPACK captures well.

The average snow depth from observations is 17 cm during the ice-growth season, with low snow depth measured before 11 July (Fig. 3b). After that, the snow depth increases rapidly up to about 37 cm, associated with a precipitation event arising from a single synoptic system. Then it decreases below the seasonal mean (Obs_mean), followed by two secondary maxima (>25 cm) on 8 September and 18 October.

The snow depth in Sim_Obs tracks the observations closely before 2 August (Fig. 3b). Then, the observed snow depth decreased quickly from about 30 cm to about 10 cm, while the Sim_Obs snow depth continued to increase gradually until the onset of surface melting in November. We attribute the observed quick decrease in snow depth to the effect of the snowdrift because the surface wind stayed above 5 m s⁻¹ for most of August (Fig. 2c), giving rise to snowdrift, a process not implemented in the version of ICEPACK used here. The snowdrift might cause a significant spatial difference in accumulated snow patterns (Liston et al., 2018), which may be responsible for the large deviation in snow depth between Sim_Obs and observations. In addition, Sim_Obs underestimated the snow depth on 11 July. As discussed above, using nonlocal observed precipitation from Progress II should be questioned.

Using observed meteorological variables as atmospheric forcing in ICEPACK produced unreliable snow depth, while the sea ice thickness was in reasonably good agreement. In other words, the enormous bias in snow depth seems to have little effect on the sea ice thickness in the simulation. This counter-intuitive finding is of great interest to us because it disobeys the general realization that the snow layer significantly modifies the energy exchange on top of the sea ice. Potential causes for this result will be discussed later.

Figure 3 Time series of (a) sea ice thickness and (b) snow depth during the freezing season. Solid black lines with black points show the observations from the drill hole (Obs). Solid green lines show the ice thickness derived from the TY buoy (Obs_TY). Solid red lines show the simulation results under in situ atmospheric forcing (Sim_Obs), and solid blue lines are simulation results under ERA5 forcing (Sim_ERA). In (b), the solid gray line shows the seasonal mean snow depth observations (Obs_mean).

3.3 Simulation forced by ERA5 atmospheric variables

When forced by ERA5 (Sim_ERA), the simulated sea ice thickness shows significant deviations from observations (Fig. 3a). The deviation is only about 1 cm before 11 July, when a heavy precipitation event (∼19 mm d⁻¹) happened. After the precipitation episode, the offset in the sea ice thickness between Sim_ERA and observations was almost constant, about 33 cm.

In contrast to sea ice thickness, the precipitation from ERA5 causes an overestimation in snow depth for the entire simulation period. The snow depth from Sim_ERA is much greater than observations, even before 11 July (Fig. 3b). During the heavy precipitation event (Fig. 2b), the observed snow depth increased from 20 cm to about 40 cm. Although the precipitation rate from ERA5 (∼40 mm d⁻¹) is 2 times larger than the observations, it caused little response in the simulated snow depth. The snow depth increase is near-linear, from about 10 cm to almost 60 cm.

3.4 Sensitivity analysis

To determine which atmospheric variables, including T_a, P, and U_a, are the most crucial in the sea ice simulation, we designed a set of sensitivity simulation experiments named SEN1. The simulation under the forcing from the in situ observed atmospheric variables is the control experiment and is named Sim_Obs. In each experiment of SEN1, one atmospheric variable is replaced by the corresponding variable from ERA5, while all others are identical to those of the control experiment. In Table 3, the averaged bias between the simulation and the observations of the outputs (ice thickness and snow depth) and the bias ratio of forcing atmospheric variables are listed separately.

Table 3 Bias of ice thickness, snow depth, and bias ratio for each forcing variable from Table 2. “All” means using the full set of ERA5 atmospheric forcing.

To determine the sensitivity of sea ice and snow depth near Zhongshan Station to atmospheric forcing, we designed a set of numerical experiments named SEN2. In the control run, the forcing of the simulation directly used the means of observed atmospheric variables (Mean_Obs in Table 4). For a specific atmospheric variable, we build a set of sensitive runs. The focused atmospheric variable changed from its mean (Range in Table 4), and other variables are the same as the control run. Considering the actual range of each observed variable on an interannual scale (M. Van Den Broeke et al., 2004; Jakobs et al., 2020; Roussel et al., 2020), we set the maximum change in T_a, Θ_a, and ρ_a to 2 % and other atmospheric variables to 50 %. Then, we concluded the sensitivity of sea ice and snow to each atmospheric forcing from its corresponding sensitive runs. Because sea ice and snow depth show a quasi-linear response to the change in each specific atmospheric forcing (not shown), the choice of the variable's range will not alter the sensitivity results.

Table 4 The atmospheric forcing (Mean_Obs for the control run and Range for the sensitive run) and sensitivity from SEN2.

Comparing the individual biases in Table 3, it turns out that P and R_ld from ERA5 contribute to the bias in sea ice thickness most strongly. For snow depth, P, U_a, and R_ld contribute the most. In Table 4, the sensitivity of ice thickness and snow depth to each atmospheric variable are listed. Comparing the individual sensitivity, it turns out that the sea ice thickness and snow depth are most sensitive to T_a and Θ_a. However, T_a from ERA5 is close to the in situ observations, so the simulated sea ice thickness and snow depth are hardly impacted (Table 3). The results from SEN1 reveal that the overestimation of P in ERA5 is the primary source of the overestimation of sea ice thickness and snow depth, even with less sensitivity to precipitation (Table 4).

To clarify the effect of specific forcing further, we replaced the x forcing in Sim_Obs with the corresponding ERA5 variable and named it Sim_ERA_x. Compared with Sim_Obs, Sim_ERA_P overestimates the snow depth from May onwards (Fig. 4b) and shows a significant positive bias in sea ice thickness after 11 July (Fig. 4a). Before 11 July, the sea ice thickness from Sim_ERA_P was even smaller than that from Sim_Obs.

To find out why the snow and sea ice behave differently, we first investigate the net heat flux into the snow surface H_N (positive downward):

$$\begin{array}{}\text{(2)}& H_{\mathrm{N}}=R_{\mathrm{n}}+H_{\mathrm{s}}+H_{\mathrm{l}},\end{array}$$

where R_n, H_s, and H_l are the net surface radiation flux, the sensible heat flux, and the latent heat flux, respectively. All energy fluxes are defined as being positive downward. Because the simulated snow layer in SIM_ERA_P is much deeper than in SIM_Obs, the difference in H_N reflects the modification of the surface energy flux due to the changed snow layer. From Fig. 4d, it can be deduced that the overestimation of snow depth in SIM_ERA_P results in a positive anomaly of H_N before 11 July, which hampers the sea ice growth. Later the difference in H_N becomes relatively small. The dependence of H_N on the snow depth is significant when the snow layer is shallow (<20 cm in this study). If the snow layer is deep enough, its impact on the net surface heat flux ceases.

After 11 July, the difference in sea ice thickness between the two simulations increases quickly from ∼0 to >40 cm (Fig. 4a). We attribute that to flooding with subsequent snow-ice formation (Powell et al., 2005). The continuously deepening snow layer reduces the sea ice freeboard. When heavy snowfall occurs, which frequently happens after 11 July, the snow load pushes the sea ice surface below sea level, and seawater floods onto the sea ice surface, causing the overlying snow to freeze. This snow-ice formation process will form flood ice (snow-ice thickness) at the sea ice surface and rapidly increase the total sea ice thickness (Fig. 4a). The difference (∼100 cm) in accumulated flood ice (Fig. 4c) between Sim_Obs (0.8 cm) and Sim_ERA_P (105.5 cm) is much greater than the difference (∼40 cm) in simulated sea ice thickness (Fig. 4a), while the net surface heat flux compares well after 11 July (Fig. 4d). This difference may be because, as the snow-ice process occurs, the increase in sea ice thickness will reduce the heat loss from the ice cover and inhibit the basal growth of sea ice in winter (Fig. 4e). The flooding-induced snow-ice formation happens at a rate larger than 0.5 cm per hour after 11 July. The snowfall (Fig. 2b) is converted to new snow depth at the top surface (Fig. 4f) using a snow density of 330 kg m⁻³ in ICEPACK (Hunke et al., 2019). Comparing Fig. 4b with Fig. 4f, we find that the change in actual snow depth (11 cm) is much lower than the expected accumulated snowfall (57 cm), indicating that the flooding process reduces about four-fifths of snow depth over sea ice.

Figure 4 Times series of (a) sea ice thickness, (b) snow depth, (c) accumulated flood ice, (d) net surface heat flux, (e) accumulated basal ice growth, and (f) accumulated snowfall. The gray line represents the simulation using precipitation from observations (Sim_Obs). The black line represents the simulation using precipitation from ERA5 (Sim_ERA_P). The color bar represents their difference (Sim_ERA_P − Sim_Obs).

Figure 5 Time series of the simulated (a) sea ice thickness, (b) snow depth, (c) accumulated flood ice, and (d) net surface heat flux in the n experiments of SEN3. The solid black point lines show the in situ observations (Obs). The 11 colored lines denote the 11 sensitivity experiments. When n=0, precipitation is from the in situ observations. When n=10, precipitation is from ERA5.

Figure 6 Box plot of simulation bias (simulation minus observation) of (a) sea ice thickness and (b) snow depth over the daily mean precipitation in the different sensitivity experiments (n increases from left to right). On the x axis, 0.66 mm refers to the experiment with n=0 (in situ precipitation), and 2.66 mm refers to the n=10 experiment (ERA5 precipitation). Two linear regression lines (black and red) are derived for x≤1.06 mm and x>1.06 mm based on the mean of ice thickness and snow depth.

3.5 Additional sensitivity simulations on the precipitation bias

The precipitation from ERA5 shows the most significant deviation compared to the in situ observations and contributes the largest sea ice and snow simulation bias. To determine the cause of differences in the sea ice and snow response to precipitation, we set up 10 sensitivity experiments named SEN3 (Fig. 5). In the nth experiment, n×10 % of the daily difference between P from ERA5 and the in situ observations is added to the observed P on that day. This procedure gradually increases the magnitude of the precipitation in the experiments, while the timing of the daily precipitation events remains almost unchanged.

We define the bias as the difference between simulations and observations from 27 July to the end of November. Different start or end dates of this period do not change this result. The bias of both sea ice thickness and snow depth linearly grows with increasing precipitation (Fig. 6). The simulation bias of the sea ice thickness is relatively small before the precipitation increases by about 1 mm per day. We suggested that the snow-ice formation is small (Fig. 5c), and the insulation of the snow layer (Fig. 5d) hampers the sea ice growth. In fact, the simulated sea ice thickness even decreases (at a rate of −3.4 cm (mm d⁻¹)⁻¹) when the added precipitation is <1 mm d⁻¹. When the added precipitation is >1 mm d⁻¹, the simulated sea ice thickness quickly increases at a rate of 22 cm (mm d⁻¹)⁻¹.

In contrast, the simulated snow depth increases rapidly at 23.9 cm (mm d⁻¹)⁻¹ when the enforced precipitation remains small but at a rate of 6.5 cm when the added precipitation is large. This is because more snow is converted into flood ice, and the snow-ice formation process strongly overrules the larger insulation effect from the snow layer, promoting sea ice growth.

The snow-ice process is based on Archimedes' principle. Therefore, the threshold value (1 mm d⁻¹) is related to the density value of ice, snow, and water in the model parameterization, as well as the sea ice thickness and snow depth. If sea ice, snow density, and initial snow depth decrease, or seawater density and initial ice thickness increase, the threshold will increase, and vice versa. These different effects of increases in precipitation on the snow and sea ice growth are illustrated in Fig. 7, emphasizing the role of flooding via snow-ice formation. When the snow layer is shallow, increases in precipitation will quickly deepen the snow layer and inhibit the growth of sea ice thickness due to the insulation of snow. The decrease in the surface net heat flux is the dominant factor. While the snow layer is deep and a lot of precipitation is present, the flooding process induces snow-ice formation, and the sea ice grows quickly, while the snow depth increases only slowly.

Figure 7 Schematic diagram for (a) low precipitation and (b) large precipitation events illustrating the precipitation effect on sea ice growth. The orange arrows represent surface net heat flux, and different colored boxes indicate the layer of snow, flood ice, and sea ice.

4 Shortcomings

The simulated ice thickness and snow depth deviate from the observations in this study (Fig. 3). We list the shortcomings that could affect the simulation: (1) superimposed ice is not considered in this study; (2) the snow-ice formation might be overestimated on the landfast sea ice in ICEPACK; and (3) the snowdrift process has not been involved in the version of ICEPACK used here.

Superimposed ice is present in early autumn when the snow starts to melt (Kawamura et al., 1997) and contributes significantly to sea ice growth (up to 20 % of mass) (Granskog et al., 2004). Superimposed ice usually corresponds to liquid precipitation or melted snow that permeates downward to form a fresh slush layer and refreezes. The superimposed ice is implemented in ICEPACK via the melt pond parametrization but has not been considered in this study. Therefore, the simulation may underestimate sea ice thickness and overestimate snow depth compared to the observations in November (Fig. 3a). We will apply the melt pond scheme in follow-up research work.

Flooding-induced snow-ice formation is common in the Antarctic Ocean because of the thin ice and heavy snowfall (Kawamura et al., 1997). It can contribute to considerable ice mass (12 %–36 %) and reduce the snow depth by up to 42 %–70 %, depending on the season and location (Jeffries et al., 2001). The parameterization of the flooding process in ICEPACK is based on Archimedes' principle for the pack ice, which might be problematic for the coastal landfast sea ice. With a much larger volume and shallower seawater around than the pack sea ice, part of the coastal landfast sea ice might contact the sea bed rather than float in the sea. Thus, the flooding should be much weaker even with weighted snow cover. Moreover, the change in density of ice due to the flooding process is significant (Saloranta, 2000) but not well considered in ICEPACK. For example, a slushy layer of 10 cm depth would refreeze within 3 d from observations (Provost et al., 2017), while the process only needs 1 d in ICEPACK. Hence, the landfast sea ice growth due to snow-ice formation needs improvement in ICEPACK, especially when the input precipitation is significantly exaggerated, e.g., the ERA5 forcing.

Surface drifting snow particles play an essential role in the surface mass balance (M. R. Van den Broeke et al., 2004). Figure 3b shows that the observed snow depth has quickly decreased from 32 cm on 2 August to 15.5 cm on 10 August, which should be attributed to the snowdrift because the surface wind is >8 m s⁻¹ in most of this period (Fig. 2c). Friction velocity becomes sufficiently high to overcome the gravity and bonds between snow particles in this strong wind and raise the snow particles from the surface (van den Broeke et al., 2006; Thiery et al., 2012; Tanji et al., 2021). However, the mean surface wind in ERA5 is convergent around the observation site during the intense wind period (Fig. 8), which might not be expected for snow depth to decrease due to snowdrift. The coarse resolution of the atmospheric reanalysis might not produce a realistic surface wind field, which is primarily determined by the local topography (Van Den Broeke et al., 1999; Frezzotti et al., 2005). In addition, surface sublimation of drifting snow particles, which is most significant in warm, dry, and windy weather (Thiery et al., 2012), plays an important role in surface mass balance (M. R. Van den Broeke et al., 2004) but has not been involved in ICEPACK yet.

Figure 8 The mean ERA5 surface wind and divergence from 2 to 10 August. The black line represents the coastline, and the red point represents the observation site.

5 Discussions

The surface wind can affect the snow depth by changing the surface heat fluxes (Fairall et al., 2003). Compared with Sim_Obs, Sim_ERA_W gives a $-\mathrm{2.5}\times {\mathrm{10}}^{\mathrm{4}}$ W m⁻² lower latent heat flux (positive downward) on average (Fig. 9b), i.e., a larger sublimation (Fig. 9c), and a reduction of about −3.4 cm of the snow depth (Fig. 9a). Therefore, the overestimation in the surface wind from ERA5 partly neutralizes the effect of overestimated precipitation.

Figure 9 Times series of (a) snow depth, (b) accumulated latent heat flux, and (c) accumulated snow sublimation. The gray line represents the simulation using wind from the observations (Sim_Obs). The black line represents the simulation using wind from ERA5 (Sim_ERA_W). The color bar represents their difference (Sim_ERA_W − Sim_Obs).

The oceanic forcing also plays an essential role in sea ice evolution (Uotila et al., 2019). Heat flux from the ocean boundary layer changes the sea ice energy balance (Maykut and Untersteiner, 1971). The ocean heat flux is mainly impacted by summer insolation through open leads, thin ice, melt ponds (Perovich and Maykut, 1990), and upward heat transfer through vertical turbulent mixing (McPhee et al., 1999). Because the oceanic observations under sea ice are challenging, most sea ice models directly use some empirical values, like the default value in CCSM3, to build the ocean boundary condition (e.g., Y. Yang et al., 2016; Turner and Hunke, 2015). Although some oceanic variables, like the water temperature and salinity, are from observations, others refer to previous studies, like the mixed layer depth. The uncertainty in oceanic forcing might be as important as the atmospheric ones, which will be focused on in our coming work.

6 Conclusions

This work uses the single-column sea ice model ICEPACK forced by the ERA5 atmospheric reanalysis and atmospheric in situ observations to simulate snow depth and sea ice thickness at Zhongshan Station, Antarctic. We find that forced by atmospheric variables from in situ observations, ICEPACK can reasonably simulate the sea ice thickness evolution, but it significantly overestimates the snow depth after the heavy snowfall on 11 July. When using atmospheric forcing from ERA5, sea ice thickness simulation is close to observations before 11 July but suddenly increases after the snowfall event.

From the sensitivity experiments, we find that the significant deviation in the precipitation of ERA5 contributes to the largest bias in both sea ice thickness and snow depth even though the precipitation is moderately sensitive to sea ice thickness (−0.032 cm %⁻¹) and snow depth (0.135 cm %⁻¹). On average, about 2 mm d⁻¹ more precipitation in ERA5 is found during the observation period, which produces about 14.5 cm excess sea ice thickness and 17.3 cm more snow depth.

We further explore the physical mechanism of the effect of precipitation on ice thickness. Snow-ice formation can be triggered by a heavy snowfall episode, like on 11 July. It efficiently produces ice at the sea ice surface, decelerates the snow accumulation, and inhibits sea ice's basal growth. When the snowfall is weak, the snow layer thickens quickly and hampers the sea ice growth through its insulation effect. When the snowfall increases to a certain degree (∼1 mm d⁻¹), it will trigger a continuous flooding process, accelerating the sea ice growth and slowing down the snow layer thickening.