Sensitivity of lake ice regimes to climate change in the Nordic region

A one-dimensional process-based multi-year lake ice model, MyLake, was used to simulate lake ice phenology and annual maximum lake ice thickness for the Nordic region comprising Fennoscandia and the Baltic countries. The model was first tested and validated using observational meteorological forcing on a candidate lake (Lake Atnsjøen) and using downscaled ERA-40 reanalysis data set. To simulate ice conditions for the contemporary period of 1961–2000, the model was driven by gridded meteorological forcings from ERA-40 global reanalysis data downscaled to a 25 km resolution using the Rossby Centre Regional Climate Model (RCA). The model was then forced with two future climate scenarios from the RCA driven by two different general circulation models (GCMs) based on the Special Report on Emissions Scenarios (SRES) A1B. The two climate scenarios correspond to two future time periods namely the 2050s (2041–2070) and the 2080s (2071–2100). To take into account the influence of lake morphometry, simulations were carried out for four different hypothetical lake depths (5 m, 10 m, 20 m, 40 m) placed at each of the 3708 grid cells. Based on a comparison of the mean predictions in the future 30year periods with the control (1961–1990) period, ice cover durations in the region will be shortened by 1 to 11 weeks in 2041–2070, and 3 to 14 weeks in 2071–2100. Annual maximum lake ice thickness, on the other hand, will be reduced by a margin of up to 60 cm by 2041–2070 and up to 70 cm by 2071–2100. The simulated changes in lake ice characteristics revealed that the changes are less dependent on lake depths though there are slight differences. The results of this study provide a regional perspective of anticipated changes in lake ice regimes due to climate warming across the study area by the middle and end of this century.


Introduction
According to the European Environmental Agency (EEA), there are more than 500 000 natural lakes larger than 0.01 km 2 (1 ha) in Europe, and three quarters of the lakes are located in Norway, Sweden, Finland and the Karelia-Kola part of Russia (EEA database, 2012).In northern Europe, most lakes are characterized by extended periods of winter ice cover each year (Blenckner et al., 2010).The ice season can be up to 7 months long and the thickness of the ice can reach 100 cm (Leppäranta, 2010).Lake ice phenology (freeze-up date, break-up date and ice cover duration) and ice thickness are regarded as a good proxy for the past and present climate (Magnuson et al., 2000); in some cases it can be considered a more robust measure of climate variability and change than air temperature (Livingstone, 1997).The various forms and processes of freshwater ice such as lake ice phenology and ice thickness are directly controlled by atmospheric conditions that influence the radiative, sensible and latent heat fluxes (Dibike et al., 2011a) and hence their spatial and temporal trends can be used as indicators of climate variability and change (Prowse et al., 2011).Changes in lake ice cover characteristics due to anthropogenic causes will have significant ecological, hydrological and socioeconomic impacts (Prowse et al., 2011).This concern has increased as the analysis of historical lake ice observations indicates significant changes in the timing and duration of the ice season that is commensurate with global warming.Magnuson et al. (2000) investigated the freeze-up and break-up trends for 26 rivers and lakes in the Northern Hemisphere.Accordingly, from 1846 to 1995, average freeze-up dates were delayed by 0.58 days per decade, while break-up occurred 0.65 days earlier per decade.Benson et al. (2011) updated the study by Magnuson et al. with more lakes and more recent data  and found that freeze-up occurred 0.3-1.6 days later per decade whereas break-up occurred 0.5-1.9days earlier per decade.
Climatic variables that are used to compute the heat balance on a lake surface (temperature, precipitation, relative humidity, wind speed, cloudiness and air pressure) are strongly affected by climate change.Air temperature is the most important factor of all influencing the thermal structure and ice cover regime of lakes.The magnitude of regional climate change signals varies depending on the chosen general circulation models (GCMs).The general trend however is that high-latitude regions will experience greater temperature increases than elsewhere and that these changes are likely to occur ever faster because of the positive feedback effects of melting snow and ice (Quesada et al., 2006;MacKay et al., 2009).Projections from GCMs have coarse spatial resolution (∼100 km at best), and these projections are particularly inappropriate for lakes which are highly dynamic and respond in complex ways to short-term changes in weather (Samuelsson, 2010).Hence, GCM projections shall be downscaled to higher spatial resolutions (25 km or less) to capture the spatial variability in meteorological forcing.The ENSEMBLES project is one such initiative which was supported by the European Commission's 6th Framework Programme as a 5-year integrated project from 2004-2009.One of the objectives of the ENSEMBLES project was to develop an ensemble prediction of climate change using high-resolution regional climate models (RCMs) (van der Linden P. and Mitchell, 2009).
Several one-dimensional thermodynamic lake models have been developed over the years, such as the Dynamic Reservoir Simulation Model -DYRSMICE (Gosink, 1987), Minnesota Lake Model -MINLAKE (Stefan and Fang, 1997), Lake Ice Model Numerical Operational Simulator -LIMNOS (Walsh et al., 1998), Canadian Lake Ice Model -CLIMo (Duguay et al., 2003), Multi-year Lake Simulation Model -MyLake (Saloranta and Andersen, 2007), Freshwater Lake model -FLake (Mironov, 2008) and high resolution thermodynamic snow and ice model -HIGHTSI (Yang et al., 2012).The models have been successfully applied to evaluate the influence of climate and lake depth on the lake thermodynamic balance and ice-phenology.There have been, for example, quite a number of studies of model applications to characterize at-site lake ice conditions (e.g.MacKay et al., 2009;Fang et al., 1996;Pour et al., 2012).A number of studies have also used one-dimensional lake models for assessing the impact of anticipated climate change on lake thermal regimes (Fang and Stefan, 1999;Stefan et al., 1998;Fang and Stefan, 2009;Saloranta et al., 2009) and ice cover characteristics (Fang and Stefan, 1998;Stefan and Fang, 1997;Fang and Stefan, 2009).
The coupling of one-dimensional lake models and gridded climate data over large areas is quite a recent phenomenon.Walsh et al. (1998) used a modified version of the lake ice model LIMNOS and produced ice phenology over the Northern Hemisphere on a 0.5 • by 0.5 • latitude-longitude grid for the period 1931-1960.The study was conducted on hypothetical lakes with mean depths of 5 m and 20 m using average monthly climate data as input.Some studies have reported modelling results of future lake ice conditions using deterministic, process-based models with input forcing from GCMs or RCMs.Dibike et al. (2011a) examined possible changes to the lake ice regime in North America from 40 • N to 75 • N under future climate.They found that break-up will advance by 10-20 days, while freeze-up will be delayed by 5-15 days, resulting in a reduction of the ice cover duration by 15-35 days in 2041-2070 compared to 1961-1990 under the Special Report on Emissions Scenarios (SRES) A2.Dibike et al. (2011b) carried out a similar exercise, but for the whole of the Northern Hemisphere between 40 • N and 75 • N (excluding Greenland), and came up with a similar conclusion: lake ice freeze-up timing will be delayed by 5-20 days and break-up will be advanced by approximately 10-30 days, thereby resulting in an overall decrease in lake ice duration by about 15-50 days, for the period 2040-2079compared with 1960-1999. Brown et al. (2011) ) carried out a grid-based simulation of lake ice conditions for Arctic North America above 58 • N and found that projected changes to the ice cover using the 30-year mean data between 1961-1990 and 2041-2070 suggest a shift in break-up and freeze-up dates for most areas ranging from 10-25 days earlier (break-up) and 0-15 days later (freeze-up).Kourzeneva et al. (2012) simulated a set of model lakes with the Flake model using a 1 • gridded global climate forcing data set for the period 1986-2006 to investigate lake parameterization for numerical weather prediction models.They found reasonable agreement with longterm temperature, but for boreal lakes, a shift in the annual cycle and cold spring caused problems with the timing of the ice period.
The above examples show offline applications of lake models, but since lakes are important in the atmospheric boundary layer, there is ongoing work to improve lake parameterization in RCMs (e.g.Eerola et al., 2010;Rontu et al., 2012) and thereby provide the opportunity to extract climatic impacts on ice directly from the RCM.For example, a coupled RCM and lake model were used to estimate the lake impacts on the climate over Europe, showing lake induced warming over most of the area (Samuelsson et al., 2010).
Hitherto, no detailed physically based modelling of lake ice conditions has been carried out at a spatially high resolution over the Nordic region and the Baltic countries which is the focus of this study.The objective of this research is to evaluate the sensitivity of lake ice conditions to climate change in the future in the study region, particularly focusing on ice formation, break-up, duration of the ice cover and maximum annual ice thickness.Through grid-based modelling of four different hypothetical lake depths (5 m, 10 m, 20 m and 40 m) using an offline lake model and applied to a 25 km spatial resolution, this study is meant to give an indication of the anticipated changes in the future of lake ice phenology, ice thickness and thermal stratification.In addition,

Model description
The lake model used is MyLake (Multi-year simulation model for Lake thermo-and phytoplankton dynamics) which is a one-dimensional process-based model.It simulates, at a daily time, step the vertical distribution of lake water temperature, evolution of seasonal lake ice and snow cover thickness and phosphorus-phytoplankton dynamics (Saloranta and Andersen, 2007).The most important physical, chemical and biological processes are included in a relatively simple and transparent model structure which is coded in Matlab.The model has been described previously in the literature (Saloranta and Andersen, 2007;Dibike, et al, 2011a, b) and hence we only include a brief description hereunder.A more detailed description is included for the interested reader as Supplement.
The model includes computational methods for simulating heat transfer at the air-water interface, advective heat transfer of inflow and outflow and the internal vertical diffusion of thermal energy.The vertical profile of the lake is conceptualized as horizontal layers of equal height ( Z), each of them having a defined area and hence volume.The onedimensional heat conservation equation for the temperature distribution in a horizontally homogeneous, vertically stratified lake is solved in each layer.The equation can be written as follows (Saloranta and Andersen, 2007): where T is laterally averaged water temperature [ For the top layer of the lake the net surface energy flux is calculated as follows: where H SW is the incoming solar radiation absorbed by the layer, H LW is the net long-wave radiation, H Sen and H Lat are the sensible and latent heat fluxes, H P R is heat advected due to precipitation and H Sed is the heat flux at sediment-water interface (all in J m −2 d −1 ).For the subsurface layers, only H SW and H Sed contribute to the local heating rate H * .During the ice cover period, only the sediment water flux and the shortwave radiations penetrating through snow and ice contribute to the local heating.The sediment heat flux module was not used in the present study.In order to compute the heat fluxes the following daily meteorological variables should be provided as input series: air temperature at 2 m height ( • C)C), cloud cover (0-1), relative humidity (%), solar radiation (J m −2 d −1 ), wind speed at 10 m height (m s −1 ), air pressure at station level (hPa) and precipitation (mm d −1 ).
If solar insolation observations are not available, MyLake computes H SW using the MATLAB Air-Sea Toolbox (Beardsley et al., 1998) as a function of Julian day, latitude and longitude.The Air-Sea Toolbox is also used for the calculation of the radiative and turbulent heat fluxes, and surface wind stress.Wind has a strong influence on the lake thermal stratification in the open water period as the surface water is forced to move due to wind stress.Depending on the wind strength, this movement is transmitted to the underlying water layers up to a certain depth.Wind-induced effects in MyLake are handled by checking if the total kinetic energy accumulated over a daily time step is sufficient to induce mixing in the water column.
The model triggers ice formation when water layer temperature drops below the freezing point.The sensible heat deficit in the super-cooled layer is turned into a latent heat of freezing and an initial ice-layer is created.The initial solid ice cover, associated with the freeze-up date predicted by the model, only appears when the accumulation of ice reaches a threshold thickness of 3 cm.Once solid ice cover has been formed, additional ice thickness due to congelation ice growth is calculated whenever the air temperature T a ( • C) is below the freezing point.When the weight of snow cover exceeds the buoyancy capacity of the ice layer, the ice surface submerges and water floods the top of the ice.This water mixes with the lower layer of the snow cover, forms slush and becomes "snow-ice" when it freezes.Snow-ice properties are assumed to be the same as congelation ice, and thus the newly formed snow-ice layer is subtracted from the snow cover and added to the ice layer.
Lake ice decay or melt is computed by considering the net heat flux at the air-snow or air-ice interface.The net heat flux is used to melt the snow cover before ice melt starts.Snow and ice albedo are respectively assigned to default values of 0.77 and 0.30 in the melting algorithm.The model also takes into account lake inflow and outflow.The inflow either flows over the reservoir surface if the inflow density is less than the density of the surface layer (i.e. the top 1 m layer of the water column) or plunges beneath the surface to be placed over a layer having density higher than the inflowing water.The sediment heat flux has been neglected in this study in light of computational time.
An important aspect for the model selection is the ability in MyLake to simulate snow accumulation and snow impact on the ice formation.This is important for insulation that hinders ice thickening on the one hand and the formation of snow-ice on the other, and it has been shown to influence the accuracy of ice simulations (Pour et al., 2012;Semmler et al., 2012;Brown and Duguay, 2011).3 Data and methods

Model testing
The lake ice model, MyLake (Saloranta and Andersen, 2007) was first tested and validated using in situ observations on Lake Atnsjøen in south-central Norway (Fig. 1) which has long-term measured temperature profiles and ice phenology data.Atnsjøen has a maximum depth of 80.2 m, a mean depth of 33 m and a surface area of 5.1 km 2 .The model requires daily average input data for six meteorological variables, including 2 m air temperature ( • C), precipitation (mm), 2 m relative humidity (percent), 10 m wind speed (m s −1 ), surface air pressure (hPa) and cloud cover (-).To discriminate between rain and snow, a fixed threshold temperature, set in this study to the default value of 0 • C, is applied.The shortwave radiation (W m −2 ) reaching the surface is computed within the model as a function of solar altitude and cloud cover.The meteorological data at nearby stations (within a distance of 20-40 km) were obtained from the Norwegian Meteorological Institute (met.no).Daily inflow to the lake and inflow temperatures are also required to produce an ac-curate thermodynamic simulation of the lake profile.Outflow is computed within the model based on daily water balance and a fixed maximum water surface elevation.Finally, we need to input initial temperature profile, ice and snow thicknesses.Typically, simulations start in the ice-and snow-free season and these values are set to zero.Inflow discharge and temperature, lake temperature profiles and ice phenology data were obtained from the Norwegian Water Resources and Energy Directorate (NVE).The model simulates a number of lake characteristics, notably water temperature profile, freeze-up and break-up dates as well as ice and snow thicknesses.Model parameters are either fixed (e.g.dz = 1 m, snow albedo = 0.77, ice albedo = 0.30) or computed from lake surface areas (e.g.vertical heat diffusivity coefficient (K) and wind sheltering coefficient (W str )).Sensitivity analyses for model boundary conditions (input meteorology), model parameters and initial conditions are all important tasks in model applications.In this study, we considered only model sensitivity to input meteorology and hydrology as model parameters were left to default values and no calibration was carried out.The sensitivity analysis involves changing one variable at a time and applying global perturbations of ±10 %, ±20 % and ±30 % on each input variable and evaluating the changes in the mean values of freeze-and break-up dates.The Atnsjøen Lake model was also driven with the Regional Climate Model (RCA) downscaled reanalysis data to evaluate how well the observed lake temperature and ice characteristics are reproduced.
Forcing data for the regional simulation was obtained from the online data portal of the EU ENSEMBLES project (http://ensemblesrt3.dmi.dk/),maintained at the Danish Meteorological Institute.The control period ) lake ice conditions were simulated using the ERA-40 global reanalysis data (Uppala et al., 2005) dynamically downscaled to a 25 km resolution by using the Rossby Centre Regional Climate Model (RCA) developed by the Swedish Meteorological and Hydrological Institute.ERA-40 is a reanalysis of meteorological observations produced by the European Centre for Medium-range Weather Forecasts (ECMWF) and integrates measurements of varying accuracy from a wide variety of in situ and remote sensing instruments to produce a comprehensive global data set of climate data for the period 1957-2002 at 2.5 • spatial resolution (Uppala et al., 2005).RCA is a coupled regional climate model which has been developed at the Swedish Meteorological and Hydrological Institute to perform optimally in a spatial resolution range of 10 to 50 km.The model was originally developed from the high-resolution weather prediction model HIRLAM; however, the majority of the physical parameterization schemes have been replaced or substantially redeveloped to allow for accurate operation of the model in climate mode (Jones et al., 2004).A comparison of RCA output driven by ERA-40 boundary conditions and a high-quality observational data set from the Climate Research Unit (CRU) showed that the model provides an accurate simulation of the seasonal and interannual evolution of the near-surface temperature in Europe in general and the Nordic region in particular (Jones et al., 2004).Hence, the MyLake model results based on the reanalysis data are believed to represent the seasonal evolution of lake temperature profile and ice cover growth and ablation.In addition, as our main objective is to look at changes in ice regimes, it is expected that the effects of the biases in meteorological forcings will cancel out when considering differences in simulated lake ice characteristics.
Hypothetical lakes of 5 m, 10 m, 20 m and 40 m depths were placed at every grid cell (25 × 25 km) which gave 3708 grid cells over the study area; and MyLake simulations were carried out for every grid cell and the respective lake depths.The four hypothetical depths considered are thought to represent a large percentage of the lakes in the region and show the variation of anticipated changes with lake depths.The gridded control period simulation was also compared with historical lake ice observations (with lakes of similar depth) from the Global Lake and River Ice Phenology Database maintained at http://nsidc.org/data/lake_river_ice/(Benson and Magnuson, 2000) which are shown in Fig. 1.In addition, linear trends in simulated (using the downscaled ERA-40 reanalysis data set) lake ice phenology (freeze-up and breakup dates, and ice cover duration) as well as annual maximum lake ice thickness for the contemporary period of 1961-2000 are determined at each grid cell via linear regression.
To evaluate the agreement between measured and simulated values, two basic statistical measures of accuracy have been used, as appropriate: the mean absolute error (MAE) and the mean bias error (MBE): where P i and O i are respectively predicted and observed values.The MAE measures the average magnitude of the errors in a set of forecasts without considering the signs and hence is unambiguous and the most natural measure of average error magnitude (Fekete et al., 2004;Willmott and Matsuura, 2005).Hence, it is the average over the "total error" that is obtained by summing the absolute values of the errors where all the individual differences are weighted equally in the average.The MBE, on the other hand, is usually intended to indicate average model "bias", that is, average over-or underprediction, and can convey useful information in relation to model performance (Fekete et al., 2004).A negative value means the model systematically over-predicts and a positive value means the reverse.

Future scenario modelling
The future climate corresponding to the Inter-Governmental Panel on Climate Change (IPCC) SRES A1B (Nakićenović and Swart, 2000) was derived from two different GCMs (ECHAM5, Max Planck Institute, Germany; HadCM3Q3, Hadley Centre, UK) downscaled using the RCA RCM at 25 km resolution.The two GCMs were chosen because of their wide application in various climate change impact studies in the Nordic region (Bergström et al., 2007;Beldring, et al., 2006).The SRES A1 group represents a future world of very rapid economic growth, global population that peaks in mid-century and declines thereafter and the rapid introduction of new and more efficient technologies (IPCC, 2001).
The SRES A1B describes a technological emphasis leading to a balance across all sources of energy.This emissions scenario was the one chosen in the EU-funded ENSEMBLES project for inter-comparison of RCMs which is the source of forcing data for this study.The mean annual increase in temperature for 2041-2070 is between 1.5 The delta-change approach was used to perturb the control period meteorological data using the mean monthly climate change signals.The delta-change method is simple to implement and has been widely applied in climate impact research worldwide (Lawrence and Hisdal, 2011;Teutschbein and Seibert, 2010;Hay et al., 2000).Climate models commonly exhibit biases in simulating the present-day climate.To minimize the influence of biases on estimates of future impacts of climate change, it is recommended to use the delta-change approach (Jylhä et al., 2004).The method assumes that future model biases for both mean and variability will be the same as those in present-day simulations (Bader, 2008).As our main objective is to look at mean changes in ice and thermal regimes rather than extremes, the delta-change approach is well suited for the task.Monthly delta changes, m , (in • C for temperature and in percent for the five other elements) are derived as the difference between the mean monthly values for modelled 30-year future climate and the ones for the current climate .The daily values for the future climate for an element X are then computed as follows: where i is the day number and m is the month.Equation ( 5) is used for air temperature and Eq. ( 6) is used for the other five elements.A smoothing routine according to Sheng and Zwiers (1998) is used to eliminate sharp and abrupt changes between days at the beginning and end of each month while maintaining the mean change values of each month.In addition, care was taken so that the perturbed new series lie within the physical limits of the variables (such as maximum cloud  cover = 1, or maximum relative humidity = 100 %, etc.).
The potential changes to lake ice regimes on a regional basis are evaluated for two time periods in the future: mid-century (2041-2070) and end-of-century (2071-2100) by comparing model simulations of the future periods with the control period.

Morphometry of the hypothetical lakes
We identified lakes from the Norwegian Lakes Database which have a maximum depth between 18 and 24 m (33 lakes) and between 38 and 44 m (27 lakes).These depth ranges were chosen to represent hypothetical lakes with maximum depths of 20 m and 40 m respectively.The areaelevation maps were downloaded from the Norwegian Water Resources and Energy Directorate online database (http: //atlas.nve.no) and digitized to be able to extract points at 2.5 m depth intervals.The relative area, A(z)/A s , versus relative depth, z/z m , curves were computed, where A(z) is lake area at a given depth z, z m is the lake maximum depth and A s is the lake surface area.Finally, the typical lake geometries were created by taking the average value of A/A s at a given depth.MyLake computes two model parameters from lake surface area, namely, the vertical heat diffusion coefficient, K, and the wind sheltering coefficient, W str .The median lake surface areas were computed for the two lake depth categories representing 20 m and 40 m deep hypothetical lakes as 0.52 km 2 and 1.48 km 2 respectively and used as lake surface areas.For 5 m and 10 m hypothetical lakes, the z/z m versus A(z)/A s relationship and the lake surface area of the 20 m deep hypothetical lakes have been used.A summary of lake morphometry data for the four lake depths considered is given in Table 1.

Model testing on Lake Atnsjøen
Detailed model tests on Lake Atnsjøen using observational data revealed that MyLake can simulate temperature profiles with good accuracy.The simulation results show a very small systematic bias with MBE = −0.04• C, and a MAE of 0.76 • C. Typical simulated and observed temperature profiles are shown in Fig. 2. Similarly, Fig. 3 shows the comparison of observed and simulated lake ice phenology as well as lake ice thickness.Further, the simulations were also carried out using RCA downscaled ERA-40 reanalysis data focussing on ice phenology and ice thickness.The modelling accuracy for freeze-up date using the two data sources (Table 2) is an MBE of −10 days (earlier) for freeze-up for both data sets.For break-up, the local climate data provides a very accurate prediction (MBE = −1 day) whereas the RCA downscaled ERA-40 data predicts with an MBE of +8 days.The delay in break-up using the ERA-40 data can be explained in part by the cold bias of the ERA-40 data compared to observed temperatures.The biases in seasonal temperatures amount to −1.6 • C for winter (DJF) and summer (JJA), −1.9 • C for spring (MAM) and −0.6 • C for autumn (SON).The overall errors of prediction as measured by the MAE are, however, of the same order of magnitude for freeze-up, break-up and ice cover duration.It is also important to note that the simulations were carried out keeping all MyLake model parameters at default values (without calibration) to ascertain the applicability of the model for region-wide simulation.

Sensitivity to meteorological and hydrological forcing
To compare the relative influence of the variation in meteorological (air temperature, precipitation, relative humidity, wind speed, cloudiness and air pressure) and hydrological The  (inflow and inflow temperatures) forcings, a simple sensitivity analysis that involves changing one variable at a time has been carried out.We conducted the sensitivity analysis by applying global perturbations of ±10 %,±20 % and ±30 % on each input variable and evaluating the changes in the mean values of freeze-and break-up dates.Some of the model inputs are interdependent, and significant correlations (correlation coefficient, r 2 > 0.30) exist between cloud cover and relative humidity (r 2 = 0.34), air temperature and inflow temperature (r 2 = 0.71) and air temperature and inflow volume (r 2 = 0.32).However, these interdependencies were not taken into account during the sensitivity analysis.From the results of the sensitivity analysis (Fig. 4a, and b), it can be seen that air temperature has, by far, the strongest influence on computed freeze-up and break-up dates.The relative humidity and cloudiness, which influence the heat fluxes exchanges, are the two other main influential climate factors.The inflow temperature and volume influence the lake ice timing only to a certain degree.In addition, it should be noted that the inflow characteristics show more influence on the computed break-up dates than on the computed freeze-up dates.

Validation of the grid-based simulation
For practical reasons, the grid-based simulation of hypothetical lakes was carried out assuming the lakes were at steady state with no inflow and outflow.The region-wide gridded simulations (using the downscaled ERA-40 reanalysis data) were compared to observations at 16 sites obtained from the Global Lake and River Ice Phenology Database (GLRIPD) (Fig. 1 and Table 2).These results also proved that the model reproduces freeze-up and break-up dates reasonably well with the gridded reanalysis data, especially when considering the uncertainty in the reanalysis data as well as ice observations.Seven lakes (1-7 in  depths ranging from 19 m to 26 m were compared with the hypothetical 20 m deep lake simulations.For each ice phenology indicator, the observation was compared to the average value of the four closest RCA cells.There was, in general, close correspondence between simulated and observed mean freeze-up dates with MAE of 3.9 days, the predictions being earlier than the observations in all but one lake.The MAE of the predicted mean break-up date was 13.1 days, the predictions being always later than observations.Overall, RCA-based simulations over-estimated the ice cover duration by 17 days.For the 9 lakes (8-16 in Table 2) used to compare the 40m deep lake simulations, the MAE of freezeup simulations was 5.3 days while for break-up it was 11.3 days.Overall, the mean ice cover duration is over-estimated with a MAE of 12.5 days.It must be noted that the predicted mean break-up date is always later than the observed one by an average of 13.1 days for 20 m deep lakes and 11.3 days for 40 m deep ones.Cold biases in the downscaled ERA-40 data set may contribute a significant portion of the prediction error of break-up.In addition, lake ice observations are usually carried out near the shores of lakes and may not be representative of the lake as a whole (IPCC, 2007).Lake ice near lake shores where observations are typically carried out may break up earlier whereas model predictions represent ice break-up for the entire lake; this may contribute a part in the almost 2-week later prediction errors.Moreover, ice events (freeze-up and break-up) are not abrupt events and require observer judgment, thus a part of the uncertainty perhaps can be assigned to the in situ observations (Latifovic and Pouliot, 2007).Another factor we investigated was the uncertainty related to lake surface area (A s ) representation.Lake surface areas were parameterized as 0.52 km 2 for the 5 m, 10 m and 20 m deep lakes and 1.48 km 2 for 40m deep lakes in the gridded simulation.These areas were determined using the median values of 60 actual lakes used to determine the shape of the lakes.We carried out the simulation for the same period  by using the actual surface areas (Table 2) while maintaining the shape of the lake as in the hypothetical characterization.It has been found that the observed freezeup dates show a much stronger dependence on lake surface areas than do the break-up dates.Similar observations have been made in relation to lake depths by (Vavrus et al., 1996).Using the actual surface areas, which are much larger than the hypothetical ones, delays freeze-up by an average of 7.4 days and break-up by an average of 2.3 days hence reducing the mean ice cover duration by 5.1 days (Fig. 5).However, the hypothetical representation of lake surface areas in the gridded simulation may not significantly affect our results especially when considering the fact that our main interest is to find an indication of the changes in lake ice conditions in the future climate.At least a significant portion of these errors are expected to cancel out when comparing the future simulations with that of the contemporary period simulations, as the errors occur (though possibly with different magnitudes) in both the present and future simulations.

Contemporary trends of simulated lake ice characteristics
We performed linear trend analysis of the contemporary lake ice simulations  of ice phenology (freezeup date, break-up date and ice cover duration) as well as annual maximum lake ice thickness at each of the 3708 (25 km by 25 km) grid cells.On a regional basis, freezeup shows nearly no trend (0.0 ± 0.6 days decade −1 ) whereas break-up has shown a regional mean advance of 0.9 day earlier per decade (−0.9 ± 2.4 days decade −1 ).The regional mean advance of break-up is slightly higher than the value reported by Magnusson et al. (2000) which was 0.63 days earlier per decade.Ice cover duration showed a weak regional trend of 0.4 days (−0.4 ± 2.9 days decade −1 ) shorter per decade whereas ice thickness reduced by a regional average of 1.1 cm per decade (−1.1 cm ± 2.1 cm decade −1 ). Figure 6 shows the lake ice phenology (Julian day) and ice thickness (cm) simulations for the control period  for the four hypothetical lake depths.

Changes in meteorological forcing
There is considerable confidence that climate models provide credible quantitative estimates of future climate change, with the confidence being higher for some climate variables (e.g.temperature) than for others (e.g.precipitation) (Randall et al., 2007).As the dominant climatic variable influencing ice cover characteristics and lake water thermal structure is the near-surface air temperature, such confidence adds to the reliability of using climate change signals for assessing future conditions.The future scenarios for the A1B emissions scenario by two GCMs (ECHAM5 and HADCM3Q3) dynamically downscaled to a 25 km resolution using the RCA RCM (from the Rossby Centre in Sweden) are used in this study.
The two GCMs are among a suite of GCMs used in the EU FP6 ENSEMBLES project (van der Linden P. and Mitchell, 2009), and have been widely used for impact studies in the Nordic region (Bergström et al., 2007;Beldring et al., 2006).
The two models are in agreement that winter warming is the highest both in the 2050s and the 2080s, and that the warming increases with latitude.The overall mean changes in the six meteorological forcings are summarized as cumulative density function (CDF) plots in Fig. 7.A more comprehensive discussion on anthropogenic climate change impacts in the study area can be found in the BACC report, Assessment of Climate Change for the Baltic Sea Basin (Phil Graham et al., 2008).

Lake ice phenology and ice thickness
Figures 8 and 9 respectively show expected future changes in lake ice phenology and ice thickness for the mid-century (2041-2070) and end of the century (2071-2100) compared to the control period .For the period 2041-2070, freeze-up will be delayed by 2 to 22 days with the largest changes occurring at lower latitudes.Break-up timing, on the other hand, will advance on a larger margin of 8 to 74 days.The advance in break-up will also be larger at lower latitudes.The delay in freeze-up and the advance in breakup will reduce the duration of ice cover by 12 to 68 days.The mean annual maximum lake ice thickness will show a marked reduction ranging from 3 cm to 59 cm and the largest diminution occurring along the Norwegian mountain ranges.
For the end-of-century period (2071-2100), freeze-up dates are delayed by between 5 and 26 days while break-up dates advance by between 11 and 73 days leading to a reduction in ice cover duration of between 21 to 85 days (compared to 12 to 68 days for 2041-2070).The sensitivity of the expected changes to various lake depths is depicted in Fig. 10    shows that the expected changes do not show clear and significant dependency with lake depth.We further compared our results with a large-scale study of the Northern Hemisphere above 40 • N by Dibike at al. (2011b), where they evaluated changes for the period 2040-2079 versus 1960-1999 considering 20 m deep hypothetical lakes placed at 2.5 • spatial grids.The same MyLake model was employed in their study.The comparison is made with our simulation of the changes between the future period 2041-2070 and the control period 1961-1990.The change in freeze-up dates is closely comparable; where the hemispherical study gives changes of the order of 0 to 20 days later for our study area, whereas our study at 25 km grid shows a range from 2 to 22 days later.Break-up, on the other hand, is expected to advance by 0 to 30 days in the case of the hemispherical study, whereas in our simulation it will advance by 7 to more than 70 days.Ice thickness reduction ranged from 0 to 40 cm in the hemispherical study whereas the present study showed values in the range of 2 to 60 cm.The comparison shows that the changes are more pronounced in the highresolution study we carried out compared to the large-scale study.In addition to the spatial resolution, the studies used different climate forcing both for the contemporary period (ERA-40 data versus downscaled ERA-40 data) and for the future climate where different GCMs and RCMs have been used for input forcing for the future period.Hence, the uncertainty associated with input forcing both due to spatial resolution and differing GCM-RCM combination imparts uncertainties on the expected changes in the future.

Lake water thermal structure
Changes in water temperature profiles are also expected as a result of the changes in meteorological forcing data.The mean annual temperature profiles were simulated for the current and the mid-century (2041-2070) future time period (average of simulations using the two GCMs chosen for the study) along a latitudinal gradient varying from 56.25 to 66.25 • N at 2 • intervals and along two longitudinal axes, 15 • E and 25 • E. The computations were carried out for 20 m deep hypothetical lakes.The expected changes are then computed as a difference between the future and current period simulations.According to the results (as shown in Fig. 11), the summer stratification is expected to start earlier and last longer under the future climate conditions.The largest changes will arise during the open water period since when an ice cover occurs, water is isolated from the air.Consequently, southern lakes will experience changes in water temperature most of the year while changes in water temperature for northern lakes will be confined to summer and autumn.It is also observed that upper regions of the water column will have larger temperature increases than the lower regions.This will lead to generally steeper vertical temperature gradients and enhanced thermal stability similar to what has been reported by a number of other modelling studies (Hondzo and Stefan, 1993;Stefan et al., 1998;Fang and Stefan, 2009).The largest change in water temperature arises during the ice break-up period.In fact, since under future climate conditions lakes are expected to be free of ice earlier, air will be able to warm the lake surface sooner.To a smaller extent water is also cooled down slower in the autumn period.Slight decreases in water temperature during winter at the lake bottom are probably due to the reduction of the ice cover and its insulating effect.Since the summer stratification lasts longer, it delays the autumn turnover of the lake, which is why slight decreases in water temperature are also experienced in autumn under future climate conditions.

Potential implications of findings
The IPCC Working Group II, Third and Fourth Assessment Reports (Anisimov et al., 2007) have stated that there is very high confidence that components of the terrestrial cryosphere and hydrology are increasingly being affected by climate change.A number of studies have demonstrated that the freshwater cryosphere and especially lake ice duration and thickness will be considerably reduced in the future (Dibike et al., 2011a, b;Brown and Duguay, 2011;Walsh et al., 1998).Changes in thermodynamic characteristics of lakes due to diminished ice cover duration (such as the one predicted in this study) will lead to a number of physical, biological and chemical changes (Corell and Cleveland, 2010;ACIA, 2005).Summer stratification is expected to start earlier due to the advance of break-up and last longer due to the delay in freeze-up as shown in our results.Decreases in lake ice duration combined with higher temperatures during the increasingly long open-water period will lead to increased evaporation and lowering of lake levels (AMAP, 2011).The largest changes in water temperature are expected in the icefree period resulting in increased evaporation from freshwater lakes.Periods that normally have ice cover will be ice-free due to climate change impacts and this will increase biological activity in the lakes.A great economic impact is also likely to arise from a reduced ice thickness and bearing capacity which could restrict the size and load limit of traffic (ACIA, 2005;Prowse et al., 2011).A shorter ice season and thinner ice cover could have some desirable consequences for hydroelectric power operation.It can lead to reduced static ice loads on dams (Prowse et al., 2011), and more water could be available during winter due to a reduction of immobilized frozen water (Seidoua et al., 2007).

Conclusions
This study has examined the effects of climate change on lake ice phenology and lake ice thickness in the Nordic and Baltic regions of northern Europe.A one-dimensional, process-based model of lake water temperature, ice cover growth and ablation (MyLake) was used to simulate lake ice phenology and ice thickness on a spatial grid of 25 km.The model underwent a thorough testing and validation process using a candidate lake that provided reasonable evidence that the model accurately captured the essential elements of the physical system.The results showed that freeze-up will be delayed by an average of 1 to 3 weeks in 2041-2070 and 1 to 4 weeks in 2071-2100.Break-up, on the other hand, will advance by between 1 to 10 weeks in 2041-2070 and by 2 to 10 weeks in 2071-2100.The combined effects of a delay in freeze-up and an advance in break-up will shorten the ice duration by 1 to 11 weeks in 2041-2070 and 3 to 14 weeks in 2071-2100.Ice thickness shows a reduction between 1 and 60 cm in 2041-2070 and 1 to 70 cm in 2071-2100.Climate model predictions show a consistent enhanced warming at higher latitudes, the winter season showing the larger percentage of the changes.A physics-based modelling using the change signals have shown that this enhanced warming, together with changes in other atmospheric variables, will result in a significant decline in lake ice duration and thickness at the end of the 21st century in the Nordic-Baltic region.In addition, there will also be significant changes in the thermal structure of lakes, with lakes in lower-latitudes warming much more than lakes at higher latitudes.These changes may have consequences on the socioeconomic and ecological functions of the lakes.
The Supplement related to this article is available online at doi:10.5194/tc-8-1589-2014-supplement.

Figure 1 .
Figure 1. Figure showing the study area domain, Lake Atnsjøen (white circle) used for detail model testing and the location of an additional 16 lakes (black circles) used for validation of the gridded simulation.

Figure 2 .
Figure 2. A partial view of simulated (solid lines) and observed (dots) temperature profiles for Lake Atnsjøen using observational forcing data.The bottom graph shows the overall performance of the model in predicting water temperature profiles.

Figure 3 .
Figure 3. Results of model tests on Lake Atnsjøen based on observed meteorological forcing data; also shown are the mean bias error (MBE) values: (a) freeze-up date, (b) break-up date, (c) ice cover duration, (d) total lake ice thickness.

Figure 4a .Figure 4b .
Figure 4a.Sensitivity of ice phenology (freeze-and break-up dates) simulation results to changes in input meteorological and hydrological forcing data for the case of Lake Atnsjøen.

Figure 5 .
Figure 5. Evaluation of lake ice simulation using downscaled ERA-40 gridded data set and observational data for 20 m and 40 m deep lakes; the legend shows the simulated values using actual surface areas.Circles represent simulations based on the hypothetical lake characterization whereas the cross symbols are the respective observations.

Figure 6 .
Figure6.Lake ice phenology and ice thickness simulation for the control period for the four hypothetical lakes as simulated using dynamically downscaled ERA-40 reanalysis data with a spatial resolution of 25 km.

Figure 7 .
Figure 7. Cumulative density function (CDF) plot of mean delta changes for the two future periods in the six meteorological forcings (TM is 2 m air temperature, PR is precipitation, WS is 10 m wind speed, RH is 2 m relative humidity, CL is total cloudiness, and AP is air pressure).The same legend is used for all plots; H refers to HadCM3Q3 and E refers to ECHAM5.

Figure 10 .
Figure 10.Box plot showing the variability of the future changes (a) for 2041-2070 and (b) for 2071-2100 for the various lake depths (box plots show the minimum, maximum, 25th and 75th percentile and median values).

Figure 11 .
Figure 11.Modelled mean daily lake water temperature profiles during the 1961-1990 (a and b) and 2041-2070 (c and d) time periods as well as the corresponding changes between these periods (e and f) at 2 • latitude intervals along 15 • E (a, c and e) and 25 • E (b, d and f) longitudes.Note that the temperature changes have been constrained to a maximum of 4 • C for better representation.
• C to 3.0 • C for ECHAM5 GCM whereas HadCM3Q3 gives a warming between 2.0 • C and 4.4 • C. For the 2080s, ECHAM5 gives a warming range of between 2.2 • C and 4.7 • C whereas the corresponding increase for HadCM3Q3 lies between 2.8 • C and 5.7 • C. Based on the seasonal distribution of the mean change signals it is observed that the warming is much more intensive during the winter (DJF) season.

Table 1 .
Lake morphometric characteristics for gridded simulation.

Table 2 .
Lake ice phenology results for Lake Atnsjøen using various sources of input meteorological forcing.
FU is freeze-up date, BU is break-up date, ICD is ice cover duration.

Table 3 .
Characteristics and locations of lakes used to validate grid-based lake ice simulation using RCA downscaled ERA-40 data set (taken from GLRIPD)