Modeling present and future ice covers in two Antarctic lakes

Abstract

Antarctic lakes with perennial ice covers provide the opportunity to investigate in-lake processes without direct atmospheric interaction, and to study their ice-cover sensitivity to climate conditions. In this study, a numerical model – driven by radiative, atmospheric and turbulent heat fluxes from the water body beneath the ice cover – was implemented to investigate the impact of climate change on the ice covers from two Antarctic lakes: west lobe of Lake Bonney (WLB) and Crooked Lake. Model results agreed well with measured ice thicknesses of both lakes (WLB – RMSE= 0.11 m over 16 years of data; Crooked Lake – RMSE= 0.07 m over 1 year of data), and had acceptable results with measured ablation data at WLB (RMSE= 0.28 m over 6 years). The differences between measured and modeled ablation occurred because the model does not consider interannual variability of the ice optical properties and seasonal changes of the lake's thermal structure. Results indicate that projected summer air temperatures will increase the ice-cover annual melting in WLB by 2050, but that the ice cover will remain perennial through the end of this century. Contrarily, at Crooked Lake the ice cover becomes ephemeral most likely due to the increase in air temperatures.

1. Introduction

Lakes are well known as sentinels of climate change as many of their physical, chemical and biological properties are sensitive to variations in meteorological and hydrological conditions (Laird and others, 1996; Adrian and others, 2009; Schindler, 2009). In endorheic basins, lakes are the lowest points of watersheds; linking the climate processes effects across landscapes (Castendyk and others, 2016), being this integration characteristic for this type of lakes (LaBaugh and others, 1997). Polar and sub-polar lakes are considered some of the most sensitive lakes to climate change (Vincent and others, 1998; Lyons and others, 2006; Vincent and others, 2008; Williamson and others, 2009; Paquette and others, 2015). In particular, perennially ice-covered, meromictic Antarctic lakes are exceptional candidates to be indicators of climate variability and change as a result of their isolation from direct human and animal impact (Wharton and others, 1992; Smith, 1999; Reid and Crout, 2008).

The McMurdo Dry Valleys (MDVs) of Victoria Land, Antarctica (Fig. 1a), contain several lakes with these features, and are distinctly sensitive to climate variations as small perturbations can lead to extreme variations in the hydrologic regime; also known as polar amplification phenomena (Dana and others, 1998; Fountain and others, 1998, 1999). This site was discovered by Scott (1905), who made observations of lake levels and their ice covers. Nowadays, large meteorological and physical data of different lakes exist thanks to the expeditions made by the New Zealand Antarctic Program, and the United States, McMurdo Dry Valleys, Long Term Ecological Research (MCM-LTER) Program (Castendyk and others, 2016).

Fig. 1. (a) Map from a portion of the MDVs, showing WLB (Taylor Valley); the McMurdo Sound; and the Taylor Glacier, which is an outer glacier from the Antarctic Plateau. (b) Map of Vestfold Hills indicating the position of Crooked Lake. Squares represent the location of the meteorological stations, and triangles define the centroids of the global cells, from the climatological models, where the blue ones were used in the simulations of each lake.

The lakes of the MDVs have been studied over 80 years. Specifically, the features of Lakes Bonney, Fryxell, Hoare and Vanda have been widely recorded (Spigel and Priscu, 1998; Gibson, 1999). Ragotzkie and Likens (1964) studied Lakes Bonney and Vanda, noting their different albedo and thermal structures, establishing the basis of how distinctive the lakes within MDVs can be. With meteorological information collected between 1986 and 2000, Doran and others (2002a) made a climate analysis of the MDVs. They investigated the behavior of Taylor (Fig. 1a), Wright and Victoria Valleys. They found that Lake Vanda, located in Wright Valley, was the warmest place in summer with 74.7 degree days above freezing (DDAF) in average, whereas Lake Brownworth, also in Wright Valley, was the coldest site in summer with only 6.2 DDAF in average. Castendyk and others (2016) discussed the changes to the physical characteristics of Lake Vanda, mainly the loss of ice cover and rise of the lake level. They showed that the ice cover thinning at this lake was due to the warming of the upper layers of the water column; whereas the lake level rise was closely related to the regional climate change, which increases air temperatures and the frequency of foehn winds. Obryk and others (2016) studied the evolution of the ice covers of Lakes Bonney and Fryxell. They associated the long-term trends of the ice thickness not only to climate change, but also to the thermal structure of each lake. They found that the lake deep heat storage speeds up or retards the response of ice covers to climate changes.

Throughout the 20th century there has been a general trend of lake levels rise at the MDVs (Chinn, 1993), except for the 1986–2000 period, in which the lake levels dropped. The latter period coincided with a decrease in summer temperatures (Doran and others, 2002b). This cooling period was abruptly terminated with the anomalously warm summer of 2001/02, referred to as the ‘flood year’ (Doran and others, 2008; Gooseff and others, 2017). Doran and others (2008) showed that the ‘flood year’ was associated with an increase of 2.4^°C in mean summer temperature and an increase of 51.4 in mean DDAF, relative to the previous year (2000/01). They also found that the increase of DDAF was closely related to an increase in the frequency of down valley winds, particularly in inland sites of the MDVs, such as Lake Bonney. Therefore, the summer wind regime also had a relevant role in the flood year; melting ice and raising lake levels. On the other hand, Gooseff and others (2017) investigated the impact of the ‘flood year’ on physical and biological properties of perennial ice-covered lakes from the MDVs. Between 2001 and 2002, the maximum mean air temperature (~3^°C) coincided with one of the highest mean summer solar fluxes recorded (~242 W m⁻²). Thus, a decadal summer cooling phase (− 1.2^°C per decade) ended suddenly in 2002 provoking the greatest amount of glacial meltwater since 1969. As a consequence, the perennial ice covers began to thin significantly until today. For example, the ice cover of the west lobe of Lake Bonney (WLB), located in Taylor Valley (Fig. 1a), thickened from 1987 to 2001 at a rate of 0.85 m per decade, and since 2002 it began to thin at a rate of 0.61 m per decade (Gooseff and others, 2017). These results show that it is likely that long-term changes of ice thickness are linked to abrupt events in response to climate change.

Another Antarctic place with lakes of such characteristics is Vestfold Hills of Princess Elizabeth Land, which has the largest concentration of meromictic lakes in Antarctica and possibly in the whole world (Burke and Burton, 1988; Gallagher and others, 1989; Gibson, 1999). Lakes at Vestfold Hills are generally saline or hypersaline near the coast, and fresh, as Crooked Lake (Fig. 1b), near the Antarctic Plateau (Burke and Burton, 1988). In 2003, insights of the ice cover seasonal behavior from Crooked Lake were found by Reid and Crout (2008). They established that in such lake, the water temperature was affected mainly by air temperature changes.

Recently, Obryk and others (2019) studied the ice thickness evolution of WLB using climate scenarios based on historical extremes of locally observed solar radiation, surface air temperature, relative humidity and wind speed froma 16-year meteorological record (1996–2012). Whereas we relied on General Circulation Model (GCM) projections for our future climate scenarios, Obryk and others (2019) selected the years that exhibited the maximum and minimum values of each meteorological variable (based on the annual mean), and then built future scenarios with different combinations of those years. Obryk and others (2019) found that the ice cover of WLB vanished between ~2025 and ~2055 when scenarios with maximum values of solar radiation, i.e. those observed during the ‘flood year’, were simulated. On the contrary, scenarios with both solar radiation and surface air temperature at minimum produced a persistent perennial ice cover with thicknesses between ~4 and ~6 m at 2060.

In this work, we analyze if changes in climate can lead to substantial responses of two Antarctic ice-covered lakes: WLB and Crooked Lake. To accomplish this objective, an ice-lake model was developed and validated. Then, air temperature projections from the Community Earth System Model-Large Ensemble (CESM-LE) under a business-as-usual scenario (RCP8.5) were used to simulate the evolution of the ice covers of WLB and Crooked Lake under climate change, as opposed to the work of Obryk and others (2019) that used extreme trends of meteorological information from a 16-year record. The aim of choosing these particular lakes was to show insights on how climate change will affect Antarctic ice covers with distinct meteorological forcing, physical properties and lake thermal structure. Additionally, the CESM-LE (Kay and others, 2015) was chosen above other general circulation models because it has a good spatial resolution; it is an earth system model that couples all the earth system components, including the cryosphere, which allows to have a better representation of temperatures over the Antarctic region.

2. Sites description

This study focuses on two Antarctic lakes: WLB and Crooked Lake. WLB is located in Taylor Valley of the MDVs, as shown in Figure 1a. Crooked Lake is situated in Vestfold Hills (Fig. 1b).

The MDVs of Victoria Land, Antarctica contain ice-free terrains with perennially ice-covered lakes supplied by glacial melt (Ragotzkie and Likens, 1964; Castendyk and others, 2016; Obryk and others, 2016). Additionally, the MDVs are within a polar desert, where mean annual temperatures range from −14.8 to −30.0^°C, and the precipitation is less than 50 mm snow water equivalent (SWE) per year (Doran and others, 2002a; Fountain and others, 2010; Castendyk and others, 2016). The interannual variability of their ice covers ranged between 3 and 6 m over the past two decades (Obryk and others, 2016). Lake Bonney (77^° 44′ S, 161^° 10′ E) lies at the lowest part of the MDVs (Hoare and others, 1964; Doran and others, 2002a). It is influenced by strong, warm and dry katabatic winds from the Polar Plateau (Obryk and others, 2016). Doran and others (2002a) determined that Lake Bonney is the second warmest site (during summer) in the MDVs, having on average 34.3 DDAF between 1993 and 2000.

Like the MDVs, the Vestfold Hills of Princess Elizabeth Land, Antarctica, contain ice-free terrains. An important difference with the MDVs is that some lakes of this place may lose their ice cover for a brief period during the austral summer (Laybourn-Parry and Pearce, 2007). Gibson (1999) observed that the ice covers from Vestfold Hills had ice thickness between 0.5 and 2 m. Vestfold Hills is also a polar desert, where mean annual temperatures range from −8.9 to −12.0^°C, and the precipitation is about 64.2 mm SWE per year (Barnes-Keoghan, 2016). Crooked Lake (68^° 37′ S, 78^° 22.3′ E) is one of the largest (~12.6 km²) and deepest (140 m) freshwater lakes in Antarctica (Gibson, 1999; Reid and Crout, 2008). It is affected by strong katabatic winds from the Antarctic Plateau (Reid and Crout, 2008). Crooked Lake had on average 43.3 DDAF on the period 1996–2003 (Barnes-Keoghan, 2016).

3. Materials and methods

3.1. Model description

Several models have been used to solve the thermodynamics of ice shelves. Since the classical analytical model proposed by Stefan (1891), many models of different complexity have been developed (e.g. Maykut and Untersteiner, 1971; Semtner, 1976; Launiainen and Cheng, 1998; Bitz and Lipscomb, 1999; Reid and Crout, 2008; Hunke and others, 2015; Obryk and others, 2016). Hunke and others (2015) developed the Los Alamos Sea Ice Model (CICE), which solves the thermodynamics of ice shelves considering radiative, conductive and turbulent heat fluxes, among other physical processes. The model implemented in the present study is based on the work of Reid and Crout (2008), and on CICE with the following assumptions: (i) only the thermodynamic aspects of ice are taken into account since the ice covers studied here are stable (Hibler, 1979; Lepparanta, 2015); (ii) the model does not account for the effects of salinity – at the bottom of the ice cover of WLB the salinity is nearly zero, whereas Crooked Lake is a freshwater lake (Ragotzkie and Likens, 1964; McKay and others, 1985; Reid and Crout, 2008); (iii) throughout the year there is no snow cover over the ice due to the low precipitation rates (Castendyk and others, 2016); (iv) there is no lateral melting, even when melt ponds are created around the edges in the MDVs lakes (Fountain and others, 1999). Nonetheless, they represent a small percentage of the total lake surface (Spigel and Priscu, 1998; Dugan and others, 2013); and (v) the incoming radiative fluxes are calculated based on equations of an Antarctic site (Reid and Crout, 2008; Obryk and others, 2016).

The model solves the 1-D heat transfer equation along the vertical direction considering heat conduction within the ice:

(1)$$ {\partial \over \partial t}\left(\rho_{{\rm i}} c_{{\rm i}}T_{{\rm i}}\left(z,t\right)\right)=-{\partial \over \partial z}\left(-k_{{\rm i}}{\partial T_{{\rm i}}\left(z,t \right) \over \partial z}+I(z,t)\right) $$

where z is the depth (positive downwards), t is the time, T _i is the ice temperature, ρ_i is the density of the ice, c _i is the specific heat capacity of ice, k _i is the thermal conductivity of ice and I(z, t) is a heat source (i.e. positive when heat is gained by the ice). The lateral heat conduction is neglected since the diffusion length scale in fresh ice is about 6 m in a year (Lepparanta, 2015). Also, Eqn (1) neglects the advection due to the upwards motion of the ice. Advection is typically ignored in these types of models, although sometimes it can result in a discrepancy of ~10% (McKay, 2004). The term I(z, t) is an internal heat source that represents the attenuation of the shortwave radiation that penetrates into the ice (Maykut and Untersteiner, 1971; Semtner, 1976), and it is modeled with the Beer–Lambert law (McKay and others, 1985; Obryk and others, 2016):

(2)$$ I(z,t)=\beta (1-\alpha) F_{{\rm sw}} \exp^{-\kappa z}$$

where β is the fraction of absorbed shortwave radiation that penetrates the ice (following Lepparanta, 2015, this parameter was set to 0.45), α is the albedo, F _sw is the shortwave radiation, and κ is the extinction coefficient. Note that the Beer–Lambert law neglects scattering and ignores heterogeneities in the ice cover (McKay and others, 1994), but is a commonly used approximation in these systems (Obryk and others, 2016).

The surface temperature (top boundary condition) is calculated from the ice-cover surface energy balance, which considers the radiative fluxes, the atmospheric turbulent fluxes and the conductive heat flux at the top of the ice (Fig. 2):

(3)$$ (1-\beta)(1-\alpha)F_{{\rm sw}}+F_{{\rm lw}\downarrow}-F_{{\rm lw}\uparrow}+F_{{\rm s}}+F_{{\rm l}}=F_{{\rm ct}}$$

where F _lw↓ is the incoming longwave radiation, F _lw↑ is the outgoing longwave radiation, F _s is the sensible heat flux, F _l is the latent heat flux and F _ct is the conductive heat flux at the top of the ice. In this study, to determine the turbulent sensible and latent fluxes, the atmospheric stability is taken into account (Monin and Obukhov, 1954; Lumley and Panofsky, 1964; Launiainen and Cheng, 1998). The nomenclature of all the parameters is presented in Appendix A, and the parameterization of all these fluxes is described in the Supplementary Material.

Fig. 2. Schematic of heat fluxes on an ice-covered lake, and the temperature profile in the interior of the ice (dashed line). Fluxes toward the ice cover are positive.

The melting rate at the top of the ice cover is modeled using Eqn (4):

(4)$$ -\rho_{\rm i}L_{\rm f}\bigg[{{\rm d} h \over {\rm d} t}\bigg]_{z=0} = \Bigg\{ \matrix{F_{\rm net}-F_{\rm c} & if \;\; F_{\rm net}\gt F_{\rm ct} \cr \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!0 & {\rm otherwise} \Bigg.} $$

where h is the ice thickness of the top layer, L _f is the latent heat of freezing and F _net is the net heat flux from the atmosphere (left side of Eqn (3)). When melting is occurring, the temperature of the ice layer is kept constant at 0^°C. If the net heat flux from the atmosphere is greater than the conductive heat flux at the top of the ice, there is an excess of energy that will melt the top layer of the ice; otherwise, the ice thickness of the top layer remains constant. Any additional energy that remains will melt the other ice layers from the top to the bottom, and it is assumed that the melted ice (i.e. liquid water) reaches the upper layers of the lake water column. The model also takes into account the sublimation of the ice when latent heat is released to the atmosphere, according to what is proposed by Hunke and others (2015):

(5)$$ -\rho_{{\rm i}}\left(L_{{\rm s}}-c_{{\rm i}}T\right)\left[{{\rm d} h \over {\rm d} t}\right]_{z=0}=F_{{\rm l}} $$

where L _s is the latent heat of sublimation. Adding the rates of ice loss in Eqns (4) and (5), the ablation rate at the top of the ice can be found. Note that at each layer of our model, sublimation and melting processes can occur independently, in the same way as in the CICE (Hunke and others, 2015), where the excess energy is used from top to bottom, and Eqns (4) and (5) allow the partitioning between melting and sublimation. At the bottom of the ice cover, ice can grow or melt depending on the energy balance between the conductive heat flux that occurs at the bottom of the ice cover (F _cb) and the sensible heat flux from the lake (F _w) (Hunke and others, 2015), as shown in Figure 2. As the melted ice reaches the upper layers of the water column, in our model F _w also accounts for the liquid that is formed by melting and for any other process that could enhance sensible heat flux from the lake, such as an additional solar radiation penetration that reaches the upper layers of the water column. The rate of ice melt or growth at the bottom of the ice cover is described by:

(6)$$ -\rho_{{\rm i}}L_{{\rm f}}\left[{{\rm d} h \over {\rm d} t}\right]_{z=H}=F_{{\rm cb}}+F_{{\rm w}}$$

where h is the ice thickness of the bottom layer. The sensible heat at the ice–water interface, F _w, can be considered as a constant if its value is small (Maykut and Untersteiner, 1971), otherwise it should be modeled with the turbulent boundary layer theory as a bulk flux (Maykut and McPhee, 1995; Launiainen and Cheng, 1998; Reid and Crout, 2008)

Equations (1)–(6) were solved using the Backward Time-Centered Space (BTCS) implicit method, as it is a well-known formulation to solve the 1-D heat equation (Kereyu and Gofe, 2016), coupled with the moving boundaries of the ice cover. In this way, the ice cover was discretized into 50 equally spaced layers, where each one of them had a temperature that evolved in time. A 3 h time step was used in all the simulations. The initial temperature profile of the ice was assumed to be linear between the surface temperature (being the initial one equal to the air temperature at the beginning of the simulation) and the fixed fresh water freezing temperature at the bottom (Reid and Crout, 2008; Obryk and others, 2016). Since some fluxes of Eqn (3) depend on the surface temperature, this equation was solved using the iterative method described in Hunke and others (2015). After the ice cover had melt or grown, the layers inside it did not have the same thicknesses anymore. To keep a uniform spatial discretization, the layers were updated based on the new value of the total ice thickness, H. Finally, to ensure energy conservation, the enthalpies of the new layers were recalculated (Hunke and others, 2015). A detailed description of the developed algorithm is given in the Supplementary Material.

3.2. Meteorological data

The information needed to run the model are meteorological data: relative humidity, air temperature, wind speed, shortwave radiation, cloud cover; optical properties of ice: albedo, extinction coefficient; and water temperature beneath the ice (or a known value of F _w).

For WLB, the model was run with 16 years of information. The climate data were obtained from Lake Bonney's weather station (Fig. 1a) (Doran and others, 1995, 2002a; Fountain and Doran, 2014).There is no information of cloud cover at Lake Bonney. Thus, in a similar way to that presented by Obryk and others (2016), the cloud cover was defined as a random number with a uniform distribution generated at daily intervals. The albedo was obtained from an interpolation based on 5 months (September–February) of the albedo measurements made by Fritsen and Priscu (1999). The generated albedo function was extrapolated to the rest of the year (spring and fall), to have an annual distribution. Because of the minor importance of shortwave radiation in spring and fall, and its absence in winter; the error associated to the extrapolation is considered to be small (Obryk and others, 2016). The extinction coefficient of ice was set to 0.85, based on the mean of a range of values estimated at Lake Bonney in mid summer (Priscu, 1991; Howard-Williams and others, 1998). The sensible heat flux from the lake was considered constant, and was adjusted minimizing the RMSE between the modeled and measured values of ice thickness.

For Crooked Lake, the model was run with meteorological information of 2003 (Reid and Crout, 2008); and the water temperature of the lake and the optical properties of the ice were obtained by Palethorpe and others (2004). Daily air temperature data were retrieved from the Australian Antarctic Data Centre (Barnes-Keoghan, 2016).

3.3. Validation data

Ice thickness of WLB were measured by Priscu (2014). The validation data contain ice thickness measurements over 16 years, only in the austral summers (except in 2008 when data were recorded also during austral autumn) because of logistical restraints (Obryk and others, 2016). Continuous ablation data over 6 years were used to validate the model (Dugan and others, 2013; Doran, 2014).

During 2003 in Crooked Lake, ice thickness data were obtained by Palethorpe and others (2004). However, between 1 February and 18 May, a data gap exists because the ice cover was too thin to support the automatic sensing probe (Reid and Crout, 2008).

3.4. Climate projections

Future simulations were forced by climatological projections from the CESM-LE Project (Kay and others, 2015). CESM-LE includes 40 ensemble members, each simulating from 1920 to 2100 (except ensemble member 001, which starts from 1850) (Hurrel and others, 2013; Kay and others, 2015). CESM couples an atmospheric model (CAM5; Gettelman and others, 2010), an ocean model (POP2; Smith and others, 2010), a land model (CLM4; Lawrence and others, 2011) and a sea ice model (CICE; Hunke and others, 2015).The recent implementation of CESM-LE showed improvements compared to earlier versions, influencing the simulated climate feedbacks at high latitudes (Holland and others, 2012; Hurrel and others, 2013).

The CESM-LE output used was the daily air temperature from all the ensemble members forced with a business-as-usual scenario of climate change; RCP8.5 (Meinshausen and others, 2009). As WLB and Crooked Lake are located in polar desert regions, precipitation is negligible (Gibson, 1999; Castendyk and others, 2016; Obryk and others, 2016), and consequently changes in future precipitation rates were not analyzed. A hybrid delta approach, similar to that used by Hausner and others (2014), was carried out considering the output from the 40 ensemble members to generate the future air temperatures at both WLB and Crooked Lake. A difference was established between the observed and the CESM-LE output daily air temperatures of the baseline period (1996–2011 for WLB, and 1996–2003 for Crooked Lake). For WLB, a moving average with a 16-year time frame was used, whereas at Crooked Lake the time frame was of 8 years. Then, the 10th, 50th and 90th percentiles for the delta values were found, and they were applied to the daily averaged air temperatures of the baseline period, from the meteorological data of WLB and Crooked Lake, respectively. The outliers from the resulting time series were removed using a 3D phase space method (Mori and others, 2007). These data were used to run the model toward the future. To evaluate the impact of climate change on the ice cover (thickness evolution, surface temperature, ablation, bottom melt and ice growth rates), three time periods were considered: historical (1996–2011 for WLB and only 2003 for Crooked Lake); mid-term (2040–2055); and long-term (2084–2099). In the projections at WLB, relative humidity, wind speed and radiation were modeled using a representative climate year that was calculated from the historical period; whereas at Crooked Lake, future simulations were run with the same meteorological data from 2003, except for the water temperature, which was determined as a function of the air temperature through the linear regression proposed by Reid and Crout (2008). The optical properties of ice, in both lakes, were held constants.

3.5. Statistics for model assessment

3.5.1. Nash–Sutcliffe coefficient of efficiency

To assess the performance of the proposed ice model, the Nash–Sutcliffe coefficient of efficiency (E) was used (Bennett and others, 2013). The E compares the efficacy of the model to a model that only uses the mean of the observed data. It ranges from − ∞ to 1, in which the latter indicates the best agreement. A value of zero indicates that the performance is equally good if the mean observed data are used, while a value less than zero shows that the model prediction is worst than using the mean observed data. E is calculated as:

(7)$$ E=1-{\sum_{i=1}^{n}{\left(O_{i}-P_{i}\right)^{2}} \over \sum_{i=1}^{n}{\left(O_{i} -\bar{O}_{i}\right)}}$$

where O _i are the observed (measured) data, $\bar {O}_{i}$ is the mean of the observed data, P _i are the predicted values from the model and n is the number of data points used for the calculation.

3.5.2. Root mean square error

To evaluate the difference between the observed and modeled data, the root mean square error (RMSE) was used (Bennett and others, 2013):

(8)$$ {\rm RMSE}=\sqrt{{1 \over n}\sum_{i=1}^{n}{\left(O_{i}-P_{i}\right)^{2}}}$$

3.5.3. Linear regressions

Linear regressions were found using the MATLAB function regress, which also report the coefficient of determination (r ²) and the p-value. The future trends mentioned below are calculated based on the median projected climate delta.

3.6. Sensitivity analysis

A model sensitivity was performed by changing in ±10% the original values of the input parameters, except for air temperature, which was increased or decreased by 1 K. To evaluate the sensitivity, a simple sensitivity index (Si) was adopted from Hoffmann and Gardner (1983):

(9)$$ {\rm Si}=\left\vert 1-{D_{\min} \over D_{\max}}\right\vert$$

where D _min and D _max are the model output (H) when a parameter was decreased or increased, respectively. Values closer to 1 indicates high sensitivity to changes, whereas Si < 0.01 means no sensitivity to changes (Obryk and others, 2016).

4. Results

4.1. Model validation

The ice model was validated with 16 years of measured ice thickness from WLB. For this lake, there are no measurements of water temperature below the ice cover, so the sensible heat flux at the ice–water interface was adjusted by minimizing the RMSE between the modeled and measured ice thickness. The statistical metrics of eight simulations with different values of F _w are summarized in Table 1. The F _w fluxes were held constant in the model because their values are small (Maykut and Untersteiner, 1971). The results of three of these simulations (F _w = 3.5, 5.5 and 7.0 W m⁻²) are shown in Figure 3a. The best fit was obtained when F _w = 5.5 W m⁻² with an RMSE of 0.21 m (red solid line in Fig. 3a), and an E of 0.64.

Fig. 3. (a) Three simulations of the WLB ice cover evolution for 3.5 W m⁻² < F _w < 7.0 W m⁻². The circles correspond to the measured ice thicknesses and the red solid line is the best modeled ice thickness with F _w = 5.5 W m⁻². The error bars represent the standard deviation of the ice thickness measurements. (b) Modeled ice thickness (solid line) forcing a change of the sensible heat flux at the ice bottom on 26 October 2001. The drastic ice thinning in March 2008 could not be modeled. In that season, the ice was observed to be flooded with lake waters.

Table 1. RMSE and E for all the modeled F _w at the ice–water interface

Although the best fit was found with a constant value of F _w through the whole simulation (F _w = 5.5 W m⁻²), the trend of increasing ice thickness between 1996 and 2002 could not be predicted using a constant sensible heat flux. Therefore, we performed a second adjustment minimizing the RMSE of the sensible heat flux at the ice–water interface considering two periods divided by the day with the maximum ice thickness measured, i.e. from 1 January 1996 to 25 October 2001; and from 26 October 2001 to 1 January 2012. F _w at the first period was forced to be smaller to be consistent with findings reported elsewhere (Doran and others, 2008; Gooseff and others, 2017). Gooseff and others (2017) pointed out that in the austral summer of 2001/02, a great amount of glacial meltwater was produced, which increased the stream flow into WLB. Doran and others (2008) found that this episode was closely related to an increase in the frequency of down valley winds. This meltwater episode thinned the perennial ice cover of WLB, most likely because a combination of the meteorological conditions and the changes in the sensible heat flux from the shallow water beneath the ice cover. The best adjustment was with F _w = 3.8 W m⁻² and F _w = 5.9 W m⁻² for the first and second time periods, respectively. The latter value remained constant for the climate projections. The overall modeled ice thickness agreed the best with the observed data (Fig. 3b) with an RMSE of 0.11 m, and an E of 0.87. Nonetheless, the drastic ice thinning that occurred in March 2008 could not be represented correctly most likely because the ice was observed to be flooded with lake water, as described by Obryk and others (2016).

To complete the validation in WLB, the modeled ablation was compared with 6years of measured ablation (Fig. 4). A fit of r ² = 0.96 (p<0.001) and an RMSE of 0.28 m were obtained. To understand how the different parameters affected the modeled ice thickness, a sensitivity analysis was performed to the climate variables, optical properties of ice, F _sw and F _w (Table 2). Results suggest that ice thickness evolution is very sensitive to changes in solar radiation, and how it is absorbed and transmitted through the ice. Also, the air temperature plays an important role in ice thickness changes. Therefore, it is expected that future temperature increases due to climate change affect the variations of ice covers. Ice thickness evolution is equally sensitive to the turbulent heat flux from the water body as to the air temperature. Thus, long-term trends of ice thickness are affected by the heat below the ice cover (Obryk and others, 2016).

Fig. 4. Comparison of the modeled and measured ice ablation from WLB over 6 years.

Table 2. Si values for the validation of 16 years in WLB. The original values of the parameters were changed in ±10%, except the air temperature, which was increased or decreased by 1 K. The F _w was of 3.8 W m⁻² from 1 January 1996 to 25 October 2001; and of 5.9 W m⁻² from 26 October 2001 to 1 January 2012. The units of each variable are presented in the nomenclature (Appendix)

Also, the model was validated with the ice thickness measurements from Crooked Lake. Since this lake had data of water temperature below the ice cover, the correlation proposed by Maykut and McPhee (1995) was used to model the sensible heat flux at the ice–water interface (see Supplementary Material). The modeled and averaged measured ice thickness for the year 2003 is shown in Figure 5. Results of ice thickness had an r ² of 0.94 (p<0.001) and an RMSE of 0.07 m, which are within the range of results found by Reid and Crout (2008): r ² = 0.89 and RMSE = 0.09 m. The model represented well the seasonal behavior of the ice cover within a year, although there is a data gap that occurred during the winter period. Model results predicted that the ice grew slower than the measured ice thickness, reaching a maximum modeled thickness that was slightly smaller compared to observations. This behavior most likely occurred because the model did not consider seasonal variations on the surface of the ice; in winter the surface is smoother and flatter, and in summer it is broken and ridged (Reid and Crout, 2008).

Fig. 5. Modeled and measured ice thickness in Crooked Lake.

4.2. Future projections of ice cover thicknesses

4.2.1. West Lake Bonney

The evolution of the mean annual ice thickness, the mean annual air temperature and the mean summer air temperature between 1996 and 2099 at WLB are presented in Figure 6. Under future climate conditions, the mean summer air temperatures at WLB rose at a rate of 0.52^°C per decade (p<0.001). Maximum mean air summer temperatures occurred in different years, and were different for each climate delta. The 2087–2088 austral summer, under the 90th percentile climate delta, was the warmest with a mean air temperature of −2.11^°C. Note that Doran and others (2008) highlighted that the summer wind regime had a relevant role in lake level increase and ice cover melting. They found that the warm summer of 2001/02 was strongly associated with an increase in DDAF between the ‘non-flood’ and ‘flood’ seasons. Doran and others (2008) also demonstrated a strong correlation between maximum summer air temperature and DDAF. Unpublished analysis of the data from all the MCM-LTER lake stations show that there is no significant difference in the correlation between DDAF and mean summer air temperature, and the correlation between DDAF and maximum summer air temperature (McKay, Personal Communication). In our work, mean projected air temperature is used to predict future climate conditions instead of maximum air temperature, since the hybrid delta approach cuts off the annual extremes that better represent maximum air temperature.

Fig. 6. Time evolution of mean annual ice cover thickness, mean annual air temperature and mean summer air temperature at WLB. Beginning in 2012, the three hybrid climate deltas (10th, 50th and 90th percentile) of the CESM-LE output are shown.

The ice cover of WLB exhibited different trends during the future simulations. First, the ice thickened between the years 2012 and 2025 at a rate of 0.22 m per decade (p<0.001), reaching a maximum thickness of 3.66 m (under the 10th percentile climate delta). This thickness increase occurred because of the decrease in mean annual air temperatures between the historical and the first years of the future values. From 2025 to 2054 the ice cover thinned at a rate of 0.07 m per decade (p<0.001). Since 2054, each simulated climate delta surpassed the mean summer air temperature of the ‘flood year’. After 2054, the ice cover thinned more quickly until 2098 (0.19 m per decade; p<0.001), dropping to a minimum of 2.34 m in the 90th percentile.

Monthly average surface temperatures of the ice cover, and monthly ice thickness change rates in WLB, for each projected climate delta, are shown in Figure 7. For each climate delta, the ice surface temperature peaks occurred in December (coinciding with the peak of the ablation rates), and the temperature valleys occurred in August. Ablation rates of the median projected climate delta decreased by ~0.86 cm per month in average between the historical and mid-term time periods, and they increased by about 1.52 cm per month in average between the mid- and long-term time periods. Conversely, melting rates at the ice–water interface under the median projected climate delta rose by~0.22 cm per month in average between the historical and mid-range time periods, and they stayed relatively constant between the mid- and long-term time periods. Among the present and mid-term time periods, the ice growth rates under the median projected climate delta decreased by 0.46 cm per month in average, and these rates increased to ~1.29 cm per month in average from the mid-term to the long-term time periods.

Fig. 7. Monthly average simulated surface temperatures (black dashed line) and monthly cumulative ice thickness change rates in WLB. (a) Historical time period: 1996–2011. (b) Mid-term time period: 2040–2055. (c) Long-term time period: 2084-2099. The projected future periods shows the three hybrid climate deltas considering the 40 ensemble members from the CESM-LE.

Historical and future simulations of mean total ice thickness change of the WLB ice cover; considering ablation, accretion and bottom melt for each month, are summarized in Table 3. Results show that the ice cover lost ice every austral summer (January, February, March) and part of spring (November, December), peaking in December when ice surface temperature and ablation rates were higher (Fig. 7). The exception is March, which exhibits an ice lost during the historical and mid-term periods, but in the long-term period its thickness increases. This shift in ice cover loss was associated with the increase of ice growth rates between the mid- and long-term time periods (Figs 7b, c). On the other hand, in autumn (April, May, June) and winter (July, August, September), the ice cover of WLB gained the majority of its ice, especially in midwinter (July, August).

Table 3. Monthly average total ice thickness change rate (cm month⁻¹) in WLB, defined as the ice growth minus the ablation and the bottom melt. To show interannual variability the monthly standard deviation is also shown

4.2.2. Crooked Lake

The evolution of the mean annual ice thickness between 2003 and 2099 at Crooked Lake is shown in Figure 8; the evolution of the mean annual air temperature, and the mean summer air temperature between 1996 and 2099 at Crooked Lake are also presented in Figure 8. Under future CESM-LE simulations, the mean summer air temperatures at Crooked Lake rose at a rate of 0.47^°C per decade (p<0.001). Maximum mean air summer temperatures occurred in different years, and were different for each climate delta. The 2098–2099 austral summer, under the 50th percentile climate delta, was the warmest with a mean air temperature of 2.97^°C. This maximum temperature would occur after the ice cover totally melted on 2069 when the mean summer air temperature is estimated to reach 0.61^°C.

Fig. 8. Time evolution of mean annual ice thickness cover, mean annual air temperature and mean summer air temperature at Crooked Lake. Beginning in 2004, the three hybrid climate deltas (10th, 50th and 90th percentile) of the CESM-LE output are shown. The ice cover completely melts between 2064 and 2069. The historical mean annual ice thickness cover evolution is not shown as it only corresponds to a 1-year period.

At Crooked Lake, the ice thinned at a rate of 0.07 m per decade (p<0.001) since 2004 until 2069 (under the 50th percentile climate delta), when the ice cover disappeared. The minimum ice thickness before 2069 was 0.89 m and occurred in 2068. The ice thickness evolution in Crooked Lake did not show an increase within the simulation period as in WLB. This was due to the low variability of the mean annual air temperatures for the historical period in Crooked Lake (see Fig. 8), compared to those in WLB (standard deviations of 0.72 and 1.08^°C, respectively)

Monthly average surface temperatures of the ice cover, and monthly ice thickness change rates in Crooked Lake, for each projected climate delta, are shown in Figure 9. For each climate delta, in the mid-term simulations, the ice surface temperature peak occurred in January (coinciding with the peaks of the ablation and bottom melt rates), and the temperature valleys occurred in July. In the historical case, the maximum monthly average surface temperature occurred in December (the information from January was not available). Ablation rates of the median projected climate delta decreased by ~0.63 cm per month in average during February–December, between the historical and mid-term time periods. Similarly, melting rates at the ice–water interface under the median projected climate delta decreased by ~1.33 cm per month in average during February–December, between the historical and mid-range time periods. Oppositely, within the present and mid-term time periods, the ice growth rates under the median projected climate delta increased by 1.93 cm per month in average in February–December. Comparisons with the long-term simulations are not presented because of the totally absence of ice.

Fig. 9. Monthly average simulated surface temperatures (black dashed line) and monthly cumulative ice thickness change rates in Crooked Lake. (a) Historical time period: 2003. There are no data for January owing to the lack of information (the measurements of ice thickness began on 18 January). (b) Mid-term time period: 2040-2055, showing the three hybrid climate deltas. As the ice cover of Crooked Lake melted completely on 2065, the long-term time period is not included.

The historical and mid-term simulations of the mean total ice thickness change for each month in Crooked Lake are summarized in Table 4. Results showed that this ice cover lost ice during January and February (austral summer) and every spring (October, November, December); including March of the historical period, peaking in January when ice surface temperature and ablation rates were higher (Fig. 9). Steep changes in ice cover loss occurred in February and March, between the historical and mid-term time periods. These steep changes were related to the increase of ice growth rates among these periods (Figs 9a, b). Furthermore, in autumn (April, May, June) and winter (July, August, September), the ice cover in Crooked Lake gained the majority of its ice, with a maximum growth in May.

Table 4. Monthly average total ice thickness change rate (cm month⁻¹) in Crooked Lake, defined as the ice growth minus the ablation and the bottom melt. To show interannual variability, the monthly standard deviation is also shown. As the historical data from Crooked Lake contemplates only 1-year period, its standard deviation is not shown. January has no data. Long-term simulations are not shown due to the absence of ice in that period

5. Discussion

Mean annual ice cover thickness of WLB responded to climate change temperatures under a RCP8.5 scenario. The ice cover thickness reaches a minimum of 2.34 m in the 90th percentile climate delta on 2098, at a thinning rate of 0.19 m per decade. If the thinning continues at the same rate, the ice cover of WLB will extinguish around the year 2220 becoming an ephemeral ice-covered lake for the first time in at least 320 years (Scott, 1905). In this way, WLB will be exposed to wind-driven turbulence, which will enhance mixing of the shallow water and may perturb the permanent vertical stratification that exists nowadays (Hoare and others, 1964; Obryk and others, 2016). This effect, along with the change in absorbed radiation, most likely will modify primary production rates and phytoplankton dynamics (Fritsen and Priscu, 1999); forcing a total adaptation to the new WLB ecosystem. Nonetheless, these results must be considered carefully. For WLB, we adjusted F _w to represent in the best possible way the few available data of measured ice thickness. We obtained two fluxes: F _w = 3.8 and 5.9 W m⁻² before and after the ‘flood year’, respectively. Future simulations held the last F _w value until 2099. However, as mentioned before, in 2054 the mean summer air temperatures reached the historical unusual values of the ‘flood year’. Hence, the sensible heat flux from the lake should increase to be consistent with our model formulation. Even when the F _w should increase in the future simulations, we kept this heat flux constant to avoid increasing the uncertainty of model results. Nevertheless, it is likely that an increase in F _w will lead to an important reduction of ice growth rates in autumn and winter. In consequence, the ice cover of WLB will thin faster than our predictions.

The most likely physical sources of F _w in WLB are three: the heating of the upper layers of the water column due to solar penetration through the ice, the meltwater fluxes from Taylor Glacier (Fig. 1a)and the heating of the the upper layers of the water column by the stored heat of deep waters, which is conducted by a thermal gradient (Obryk and others, 2016). Particularly, in the WLB simulations, both β and κ are constants, thereby the sensible heat flux from the lake is not driven by solar radiation, which is the main difference with the model of Obryk and others (2016, 2019). Furthermore, given that the thermal structure of the water was not considered in the simulations, F _w is not associated to the stored heat of deep waters. Therefore, our hypothesis is that changes in meltwater fluxes from the glacier (connected to abrupt changes in air temperature and solar radiation, as stated by Doran and others, 2008 and Gooseff and others, 2017) produce a response on the sensible heat flux of the lake. This is why in our simulations, the F _w was forced to change in the period 2001–2002.

The incoming heat flux from the lake, combined with the conductive heat flux at the ice bottom, regulates the latent heat of freezing released to the atmosphere during periods of ice growth through ice conduction (Lepparanta, 2015), as shown in Eqn (6). For each modeled case, the long-term period exhibited a significant increase in ice growth rates of April and March (Fig. 7c). The reason for this change is that the thinner the ice, the smaller the distance available for heat conduction. Thus, the larger the ice growth rate of the ice cover. The yearly average distribution of simulated F _cb for the three evaluation time periods is shown in Figure 10a. It can be seen that long-term F _cb reached the value of F _w = 5.9 W m⁻² in the first days of March, contrary to the historical and mid-term simulations (mid March). Then, from April to mid-September, conductive heat flux at the ice bottom grew 17% in average (respect to mid-term simulations), under the median climate delta. Also, this effect is reflected in the March total ice thickness change rate results (Table 3), which in the long-term time period under each climate delta, the ice cover of WLB gained ice instead of losing it as the trend of other periods.

Fig. 10. Yearly average of daily F _cb in: (a) WLB and (b) Crooked Lake; for historical, mid-term and long-term simulations, considering each hybrid climate delta. In Crooked Lake, the results are shown since 18 January for comparison, due to the absence of data before that day in the historical time period. The difference between F _cb values of the ice covers from the two lakes is an indicator of how distinct are their water thermal structures.

The increase of F _cb is clearly a modeling issue, closely related to the decision of maintaining F _w constant, which is the main limitation of this work. Studies in an Arctic lake showed that bottom melt is only relevant when the ice thins to <1 m (Heron and Woo, 1994; Dugan and others, 2013), implying that in the ice cover of WLB, the bottom melt rates should be negligible (Dugan and others, 2013). Therefore, in our model formulation for WLB, bottom melt may be augmented by numerical conditions (F _w used as a parameter). This is one of the reasons why our model validation in WLB presents a good fit with ice thickness measurements (Fig. 3b), but did not represent well the ablation measurements (Fig. 4) with an RMSE of 0.28 m, compared to the 0.11 m calculated by Obryk and others (2016). In this basis, ice thickness change rates found in this study should be considered as a total, and not individually (ice growth, ablation, bottom melt).

It is also interesting to compare our results at WLB with those obtained by Obryk and others (2019). They established an envelope of climate possibilities based on climate extremes from a 16-year meteorological record, as opposed to our methodology that projects future climate using the CESM-LE output (Kay and others, 2015) under a business-as-usual scenario of climate change (RPC8.5), combined with the hybrid delta approach (Hausner and others, 2014). At one extreme, Obryk and others (2019) found that the ice cover of WLB vanished between ~2025 and ~2055 when maximum solar radiation values were simulated. Greater values of solar radiation result in an increased temperature of the upper layers of the water column, increasing the sensible heat toward the ice (Obryk and others, 2016), which quickly reduces the ice thickness. At the opposite extreme, Obryk and others (2019) found persistent ice cover thicknesses of ~4–6 m at 2060 in their scenarios with minimum solar radiation, air temperature, wind speed and relative humidity. The results presented in this work are more conservative than those of Obryk and others (2019), predicting that the ice cover of WLB will maintain its perennial behavior until at least 2100, for two reasons. In our hybrid delta method, the 10th, 50th and 90th percentiles of GCM projections are used as future climate drivers, reducing the influence of the most extreme events. Second, as explained above, we kept F _w constant to avoid increasing the uncertainty of our model results. In our model, F _w also accounts for processes such as additional penetration of solar radiation into the upper layers of the water column and the melted ice in the cover. As a modeling exercise, we performed additional simulations in which F _w was increased between 3 and 10% each year (data not shown). These results were more consistent with Obryk and others (2019), forecasting that the WLB ice cover could shift from perennial to seasonal between ~2085 and ~2040, respectively. Therefore, it is likely that the timing of shift between perennial and seasonal ice cover at WLB is between the predictions of Obryk and others (2019) and those reported in this work.

In Crooked Lake, February and March ice growth rates showed a great increase; whereas April and August ice growth rates decreased (Figs 9a, b) between the historical and the mid-term simulations. The bottom melt ceased in March. As in WLB, these particular changes were, most likely, due to the parametrization of F _w used in Crooked Lake. As indicated before, the linear regression established by Reid and Crout (2008) was considered to determine T _w and to calculate F _w. Nevertheless, the linear regression had a r ² = 0.69, being the lowest (July, August) and highest air temperatures (January, February) the values with the greatest dispersion respect to the linear tendency (see Fig. 11 of Reid and Crout, 2008). This shift is also appreciated in Table 4, in which the increment of total ice thickness change is noticeable, especially in February and March. It can be seen in Figure 10b, that between February and April, the conductive heat flux at bottom increased, enhancing the ice growth.

Crooked Lake loses its ice cover under each hybrid climate delta simulated: (i) in the year 2064 for the 90th percentile climate delta, when mean summer and annual air temperatures are of −0.3 and − 7.6^°C, respectively; (ii) on the year 2067 under the 10th percentile climate delta, when mean summer and annual air temperatures are of −0.1 and −8.3, respectively; and (iii) on the year 2069 under the 50th percentile climate delta, with mean summer and annual air temperatures of 0.1 and −7.5^°C, respectively. Will this ice cover come back? Our model did not consider ice formation in the complete absence of it. This is why the ice cover of Crooked Lake did not recover after it vanished (Fig. 8). However, it is very likely that this lake will recover its ice cover in autumn or winter, exhibiting an ephemeral behavior. This behavior can be supported by analyzing the projected DDAF, under each climate delta, as shown in Figure 11.

Fig. 11. Historical and projected degree days above freezing (DDAF) for WLB (a) and Crooked Lake (b), under each hybrid climate delta. The dashed line indicates the limit between a perennial and a seasonal behavior of the ice cover (based on Crooked Lake projected results).

In the median projected simulation, before the ice melted completely in 2069, Crooked Lake experimented the austral summer (2068/69) with more projected DDAF (= 206) (Fig. 11b). Hence, more days of surface melting occurred, and maximum values of sensible and latent heat flux until that summer were observed: 248 and 62 W m⁻², respectively. For these reasons, the ice cover of Crooked Lake could not stop melting, until it disappeared in the first days of February. When the ice cover thickness reached values <0.2 m (mid January), it stopped being stable. Hence, other aspects of ice, such as mechanics and dynamics, should be taken into account (Lepparanta, 2015). Also, in some future years, under every percentile climate delta, the projected DDAF are less than the summer when the ice vanished. Therefore, the ice cover of Crooked Lake may become perennial again. Nonetheless, under the median and 90th percentile climate delta, the ice cover would not exhibit a perennial behavior after 2080, while under the 10th percentile climate delta, the ice cover will be ephemeral after 2084. Note that the model implemented here corresponds to a spatial average, i.e. there were no considerations of the differences in the surface ice cover such as melt ponds, lateral melting and ice ridges. Thus, it is likely that the ice cover of Crooked Lake will melt in a greater extension than that predicted in our analysis. On the other hand, WLB (Fig. 11a) will reach the projected DDAF of the perennial–ephemeral transition of Crooked Lake in ~2091 (under the 90th percentile climate delta). This could be the onset of the ice cover disintegration in WLB.

The methodology employed in this study showed that using a one-dimensional ice-lake model combined with the hybrid delta approach can describe how Antarctic lakes will respond to climate change. Nonetheless, our estimates still have a relevant component of uncertainty. This uncertainty is associated with the lack of information regarding the thermal structure of the lakes, the inherent uncertainty associated to climate projectionsand the changes in interannual ice properties that are not taken into account, such as the ice albedo that depends on the complex surface evolution of the ice, which still is a challenge in climate projections (Hohenegger and others, 2012). Therefore, the future projections presented in this work should be interpreted with caution.

6. Conclusions

An ice-lake model was developed to determine the ice thickness evolution of two Antarctic lakes: WLB and Crooked Lake; under three different future climate deltas determined by the hybrid delta approach. The model is forced by radiative fluxes, atmospheric turbulent fluxes and a turbulent sensible heat flux from the water body beneath the ice cover. The model solves the 1-D heat equation within the ice considering a set of moving boundary conditions. In this way, the model predicted the evolution of the ice cover thickness of each lake. We found that under the median climate delta, WLB will exhibit a dramatic thinning since 2091, and Crooked Lake most likely will be ephemeral from 2069.

Even though the model validation in both lakes was satisfactory, future results should be interpreted with caution. As this work pointed out, ice-covered lakes respond to change in air temperatures, thinning their ice covers toward the future. However, this ice reduction was biased by the assumptions made to determine a turbulent heat flux from the lakes (constant flux in WLB, and a linear regression in Crooked lake), which caused an increase in ice growth rates. Therefore, the thermal structure of the water column in ice-covered lakes is a relevant physical feature that must be investigated, especially in a climate change analysis. Nonetheless, this work provides a first-order analysis that can be used to comprehend the ice cover dynamics of these lakes, and to recognize how Antarctic ice-covered lakes will respond to a particular climate change scenario.