Towards modeling data-poor lakes at the regional scale using parameters from data-rich lakes and relationships to lake characteristics

ABSTRACT

Lakes pivotal for recreation and economically relevant activities are often remote and not well studied, which hinders the application of predictive lake models for their management. Here, we provide an approach to simulate—by means of the process-oriented model MyLake—water temperature, ice cover duration, dissolved oxygen, and light attenuation in 198 data-poor lakes based on parameters obtained for a subgroup of 12 data-rich lakes and morphometric data. Specifically, the model is first calibrated using a genetic algorithm on well-studied lakes. Simple relationships between the fitted parameters and lake-catchment morphometric properties are then derived, and the results of simulations using fitted and derived parameters are compared. The loss in goodness-of-fit, expressed as root mean square error (RMSE) incurred by using estimated rather than calibrated parameters, is 0.17 °C for water temperature and 0.82 mg L⁻¹ for dissolved oxygen. These general relationships are then used to provide the model parameters for 198 data-poor lakes distributed throughout Sweden and to model these lakes. Overall, this proof of concept allows simulating lakes selected based on their relevance for lake management rather than based on the availability of extensive field datasets.

KEYWORDS:

• climate change impact
• data-poor lakes; lake modeling
• lakes
• model calibration
• oxythermal habitats

Introduction

Ecosystem services provided by lakes (Dodds et al. 2013) are defined by a wide array of morphological variables such as area, bathymetry, and fetch (Zentner et al. 2019) and physicochemical variables such as water temperature, dissolved oxygen (DO) concentration, and water clarity (Adrian et al. 2009, Stasko et al. 2012, Holubová et al. 2021). These variables respond to external drivers such as air temperature, ice cover, and water colour, which are closely linked to both weather and climate. Rapid changes in these variables have been observed, particularly in the boreal lakes (Adrian et al. 2009, Stasko et al. 2012). As such, climate change leads to warmer, oxygen-depleted (Schindler 2009, O’Reilly et al. 2015), and darker lakes (Van Dorst et al. 2020).

Warmer air temperature increases heating of the surface water and causes stronger thermal stratification, with a sharper thermocline reducing heat transfer to the hypolimnion and resulting in a warmer epilimnion and possibly a cooler hypolimnion (Magee et al. 2018). An increasing disparity in temperature between the epilimnion and hypolimnion has thus been reported (Kraemer et al. 2015, Richardson et al. 2017, Bartosiewicz et al. 2019, Pilla et al. 2020), along with changes in the timing and duration of seasonal stratification, namely earlier spring and later autumn turnovers (Magee et al. 2018). Finally, warmer temperatures also impact ice cover duration (Bartosiewicz et al. 2020), with a cascading effect on light and DO (discussed later).

The ice break-up and freeze-up dates influence the heat regulation of the water column and the timing of spring and autumn turnover (Beier et al. 2012). The duration of ice cover and timing of turnover also influence light penetration in spring as well as oxygen availability from changes in diffusion and photosynthesis (Bouffard et al. 2019). Ice cover duration is thus linked to the ventilation of the water column in early spring and the availability of oxygen after autumn turnover (Pilla and Williamson 2021). Longer open-water periods may increase oxygen input from the atmosphere (e.g., Couture et al. 2015) and promote autotrophy. Because ice phenology is strongly related to air temperature, cloud cover, wind, and precipitation, it is therefore also sensitive to climate change. A noticeable reduction in ice cover duration has been reported for lakes in the Northern Hemisphere over the last century (Woolway et al. 2020).

The increase of coloured dissolved organic matter (CDOM) in lakes reduces light availability (Thrane et al. 2014), which is antagonistic to the effect of shorter ice cover duration. As such, increasing CDOM also affects heat content and distribution in the lake (Solomon et al. 2015, Knoll et al. 2018, Bartosiewicz et al. 2019, Pilla et al. 2020) because more heat is absorbed in a darker epilimnion, favouring a sharper thermocline. The observed increase in CDOM in boreal lakes, referred to as “browning,” has been shown to thin the photic zone and shift predatory fish optical habitat (Weidel et al. 2017, Van Dorst et al. 2020). Finally, increasing dissolved organic carbon (DOC) loads along with CDOM influences DO concentration (Couture et al. 2015). DOC, comprising labile organic matter, fuels microbial metabolism and increases DO consumption in both the water column and the sediment. Longer stratification increases anoxia as the presence of the thermocline hinders bottom water ventilation, highlighting a synergy between air temperature and DOC in enabling oxygen depletion.

The response of water temperature, DO, ice cover duration, and water colour to ongoing environmental and climate change is a current preoccupation of managers who need to focus limited resources on lakes most likely to respond positively to conservation efforts. The preservation of cold-water fish habitats, for instance, provides an impetus for management efforts (Stefan et al. 2001, Kraemer et al. 2021). Such fish are found in boreal lakes and seek cold, well-oxygenated (Jacobson et al. 2010), and well-lit waters (Stasko et al. 2012). During summer, the habitat of those fish may be markedly reduced as the temperature preference limits the habitat to the hypolimnetic layers where the DO concentration is likely to diminish with increasing warming. Ultimately, the volume of suitable habitat may shrink to lethal levels (Nürnberg 1995). Reduced fish habitat impacts not only their biomass, but also the ecosystem components that depend on top-down control, such as community composition, prey abundance, and the biomass of primary producers (Karlsson et al. 2009, Vasconcelos et al. 2019).

Assessing the impact of climate changes on the availability of oxythermal (Stefan et al. 1995, Magnuson et al. 1997, Jacobson et al. 2010) and optical (Stasko et al. 2012, Weidel et al. 2017, Keyler et al. 2019, Van Dorst et al. 2020) habitats helps inform science-based management decisions relevant for both recreational and economical purposes. This assessment thus calls for accessing multiple lakes simultaneously, especially when considering the ongoing impact of climate change.

Process-oriented lake models are powerful tools to forecast the evolution of lake dynamics under climate change. Modeling of lakes at the regional scale has been performed using 1D lake models (Moore et al. 2021). The hydrodynamic modules of such models can be transferred from site to site without extensive site-specific calibration (Winslow et al. 2017), which likely is not the case for the prediction of other variables such as DO because their calibration generally requires lake-specific field data (Fang and Stefan 2009, Fang et al. 2012). These data are typically not available at the regional scale (Robson 2014), where dataset availability ranges from minimal, such as only morphology and sparse chemistry surveys (Jacobson et al. 2010, Van Zuiden et al. 2016), to comprehensive, where moored instruments provide high-density time series (e.g., Fang and Stefan 2009, Fang et al. 2012, Couture et al. 2015).

We propose an approach to model these variables at the regional scale using the sparse and/or meagre available data to provide a proof-of-concept to enable lake managers and governmental agencies to utilize datasets comprising lakes with sparse data and lakes with extensive time series. We describe the approach enabling the modeling of a large number of lakes, not only for water temperature and ice cover duration but also for state variables requiring lake-specific model calibration, namely DO and light, using a group of data-rich and data-poor lakes from throughout Sweden. First, the data from 12 data-rich lakes were used for model calibration, and the model parameters for 198 data-poor lakes were then derived from a simple statistical method, enabling the simulation of their physicochemical dynamics under climate change scenarios.

Methods

Dataset

We chose data from 210 dimictic, oligotrophic, boreal lakes distributed across Sweden. These lakes were selected based on the coexistence between Arctic char (Salvelinus alpinus) and pike (Esox lucius) as well as their relevance for recreational fisheries. They span a large range of areas, from 0.4 to 1900 km². The morphometric and catchment data, such as maximum and average lake depths, lake and catchment area, and residence time, were retrieved from the Swedish Meteorological and Hydrological Institute (SMHI 2017) database. For modelling purposes, the area at each 1 m depth increment (A_z, m²) was calculated after by multiplying the lake surface area (A_L, m²) with the proportional area at the depth (z, m) using a coefficient estimating the basin shape (S) based on the maximum (z_max, m) and mean depth (z_mean, m) (Lester et al. 2004): (1) ${\mathrm{A}}_{\mathrm{z}}={\mathrm{A}}_{\mathrm{L}}\times (1-(\mathrm{i}/{\mathrm{z}}_{\max }{)}^{\mathrm{S}}{)}^2,$(1) where (2) $\mathrm{S}={\displaystyle \frac{3\mathrm{r}+{\left({\mathrm{r}}^2+8\mathrm{r}\right)}^{0.5}}{4\left(1-\mathrm{r}\right)}\mathrm{and}\mathrm{r}=\left({\mathrm{z}}_{\mathrm{mean}}/{\mathrm{z}}_{\max }\right).}$(2) Finally, sediment area (A_S, m²) was estimated as the sum of the difference in the horizontal area between 2 consecutive depth levels (Saloranta and Andersen 2007). Within the 210 lakes, only 12 had temperature, light, and oxygen profiles or time series. Model calibration was thus constrained by the availability of field data. The 12 lakes with data are hereafter referred to as the “focus lakes,” which are distributed throughout the study area (Fig. 1) and present a range of lake shape, depth, area, and residence time (Table 1).

Field and hydrographic data

Volumetric inflow rates and temperature of inflowing waters for all lakes were obtained from the SMHI databases of simulated runoff into Swedish lakes. Vertical profiles of water temperature and DO concentrations, along with Secchi depth (SD), were measured as part of the monitoring of the 12 focus lakes. The data cover at least 10 years from 1989 to 2017, with datasets ranging from 48 days of data (lake #1) to 7017 days (lake #12). The other lakes did not have sufficient coverage of temperature, oxygen, and SD data for this period. Ice-on and ice-off dates were available for 52 of the 210 lakes, including the 12 focus lakes, for different durations between 1866 and 2011 (SMHI 2019). For most of the lakes, observations of temperature and DO were available only for the upper part of the water column, which is above the estimated seasonal thermocline.

Figure 1. Location of the 210 study lakes throughout Sweden. Focus lakes and data-poor lakes are identified with black and white diamonds, respectively.

Table 1. Morphometric characteristics of the 12 focus lakes. The coordinates indicate the main inflow of the lakes.

Lake #	Name	Location		zmax m	zmean m	AL km^2	VL km^3	AS km^2	ASC km^2	AC km^2	C:L	SC:L	μice d	RT yr
1	Ralången	57°54′ N,	4°49′ E	5	2.86	4.94	0.017	4.97	20.07	591.30	119.81	4.07		0.102
2	Bösjön	61°19′ N,	14°16′ E	17	4.28	1.18	0.0057	1.21	9.39	9.39	7.97	7.97		0.959
3	Västra Skälsjön	59°56′ N,	15°33′ E	18.9	3.33	0.392	0.0015	0.41	1.71	1.71	4.37	4.37		2.88
4	Solgen	57°32′ N,	15°06′ E	20	7.16	21.62	0.17	21.80	86.85	621.50	28.75	4.02	110	0.933
5	Viaredssjön	57°41′ N,	12°44′ E	28.63	12.4	3.86	0.050	3.97	31.21	68.66	17.80	8.09		1.08
6	Storvindeln	65°38′ N,	17°27′ E	35.9	15.3	52.16	0.82	52.64	218.91	2901.63	55.63	4.20	197	0.399
7	Betarsjön	63°44′ N,	16°46′ E	45.08	10.6	34.09	0.38	34.48	123.91	478.16	14.03	3.64	182	2.13
8	Sommen	58°08′ N,	15°02′ E	53	13.8	130.30	1.86	131.25	92.68	808.93	6.21	0.71	100	4.42
9	Åsunden	58°08′ N,	15°41′ E	53.53	14.3	52.58	0.78	53.19	259.80	1925.15	36.61	4.94		2.26
10	Storsjön	63°17′ N,	14°27′ E	74	19.9	456.33	9.31	458.84	1000.53	12061.40	26.43	2.19	140	1.22
11	Stor-Ullen	59°58′ N,	13°44′ E	74	20.5	7.42	0.16	7.75	28.80	38.78	5.22	3.88		9.54
12	Vättern	58°31′ N,	15°02′ E	128	40.3	1885.98	76.95	1895.58	1757.99	5733.68	3.04	0.93	64	56.00

z_max = maximum depth; z_mean = mean depth; A_L = lake surface area; V_L = lake volume; A_S = sediment area; A_sc = sub catchment area; A_c = catchment area; C:L = catchment area/lake area; SC:L = sub catchment area/lake area; μ_ice = mean ice cover duration; RT = residence time

Meteorological forcing data

Meteorological forcing data were retrieved from the World Climate Research Program Coordinated Regional Downscaling Experiment (EURO-CORDEX) project (Jacob et al. 2014) at 0.11° resolution, centred on the European domain of ∼27–72°N, ∼2°W–45°E and available at the daily time scale directly usable by the lake models. We used time series from the CORDEX data for each lake and its catchment. If the lake or catchment area was entirely contained within one of the CORDEX pixels, then the model used the data from that pixel; if the lake or catchment geometry covered multiple CORDEX pixels, we used the fractional area of the lake in each pixel, ignoring islands, and calculated a weighted mean.

We used 6 global climate models (GCMs) and corresponding regional circulation models (RCMs) for Europe that comprise a subset of CORDEX data (Supplemental Table S1), omitting 1 or 2 simulation(s) for a given period and greenhouse gas (GHG) scenario. For instance, daily near-surface relative humidity was unavailable for the Max Planck Institute (MPI) model and was set to 50% throughout, while cloud cover and near-surface humidity were not available for the Institut Pierre-Simon Laplace (IPSL) model and were set to 0.65% and 50%, respectively (Karlsson 2003; description of the GCM–RCM combinations in Jacob et al. 2014). We used 5 decades of data covering historical periods (1971–1980 and 2001–2010) and free-running simulated decades (2031–2040, 2061–2070, and 2091–2010). The global climate models were driven by 2 separate trajectories of GHG concentrations, Representative Concentration Pathways RCP4.5 and RCP8.5 representing the moderate and extreme GHG concentration trajectories adopted by the United Nations Inter-governmental Panel on Climate Change (IPCC) for its 5th Assessment Report (Moss et al. 2008, Meinshausen et al. 2011, Stocker 2014). These RCPs correspond respectively to the stabilization of radiative forcing after the 21st century at 4.5 W m⁻² or a rising radiative forcing crossing 8.5 W m⁻² by the end of the 21st century.

Lake modeling

The 1D process-based lake model MyLake (Saloranta and Andersen 2007) was used to simulate time series of daily water temperature, light penetration, and DO at 1 m vertical resolution for each study lake as well as daily ice cover thickness. The model focuses on in-lake physical processes (diurnal vertical distribution of thermal energy, surface heat balances, shortwave radiation penetration, wind mixing, ice formation, and break-up) and a simplified oxygen module comprising water column and sediment oxygen demand equations (Stefan and Fang 1994, Kiuru et al. 2018), with the latter replacing the sediment diagenesis module (Couture et al. 2015). This version of the model uses meteorological and hydrological input variables (Table 2) and simulates time series of the following state variables: water temperature, ice and snow thickness, DO and DOC concentrations, and photosynthetically active radiation (PAR) irradiance (400–700 nm).

Table 2. Daily input variables used by the MyLake model.

Variables	Units	Variable names
Cloud fraction (0–1)	—	Ctl
Surface radiation: shortwave downwelling	MJ m^−2	rsds
Near-surface relative humidity	%	hurs
Precipitation	mm d^−1	Pr
Near-surface air temperature	°C	Tas
Near-surface atmospheric pressure	hPa	Ps
Near-surface wind speed	m s^−1	sfcWind
Inflow volume	m^3 d^−1	Inflow_V
Inflow temperature	°C	Inflow_T
Inflow DOC concentration	g m^−3	Inflow_DOC
Inflow DO concentration	mg m^−3	Inflow_DO

Calibration of the focus lakes

The calibration procedure allowed the following model parameters to vary from their initial values: PAR light attenuation coefficient (SWA_b1, 1 m⁻¹), non-PAR light attenuation coefficient (SWA_b0, 2.5 m⁻¹), wind sheltering coefficient (c_shelter, 0.7), and oxygen demand coupled with the decomposition rate of organic carbon in the water column (k_BOD, 0.1 d⁻¹) and in the sediment (k_SOD, 500 mg m⁻² d⁻¹) (equations pertaining to sensitive parameters in Supplemental Material). Further, the inflow scalers were also allowed to vary to approximate inflow volume (I_scV, 1), temperature (I_scT­, 0), DO (I_scO, 1), and DOC (I_scDOC, 1). Note that values for the diffusion parameter a_k during open water and ice cover periods, which is used with the stratification stability to estimate the vertical diffusion coefficient (K; m² d⁻¹), were not obtained from the calibration but calculated from the lake surface area using the optional default parameterisations included with MyLake.

The model calibration, or determination of the model parameter values, for the 12 focus lakes was done using the genetic algorithm (GA) tool provided by the Matlab global optimization toolbox (MATLAB 2019) with a population size of 48 parallel model runs and a maximum of 60 generations, leading to 2880 model runs per lake. The parameter values were initialized with values from the previously calibrated MyLake application at the boreal lake Langtjern (De Wit et al. 2018) and then allowed to vary with the automated calibration procedure. Calibration was done using the historical period 2001–2010, with the normalized RMSE used as a performance metric (Los and Blaas 2010). Performance against ice observations was considered for validation only using the additional historical period of 1971–1980. This calibration method provided the values for the parameters mentioned earlier (Table 3).

Table 3. Relations between parameters and predictors: the equations used for each parameter with R² and p-value.

Parameters	=	Predictors	R^2	p-value
Log (SWAb1)	=	Log (C:L) + log (longitude)	0.566	0.233
Log (SWAb0)	=	Sqrt (SC:L) + log (area)	0.821	<0.001
c_shelter	=	Log (C:L) + log (area) + log (latitude)	0.683	0.0216
Log (IscV)	=	Log (area) + log (longitude) + log (latitude) + log (mean depth)	0.857	<0.01
IscT	=	SC:L + longitude + latitude	0.697	0.0195
Ln (IscO)	=	Sqrt (SC:L) + ln (longitude) + log (area)	0.684	0.0178
Log (kSOD)	=	Log (C:L) + sqrt (SC:L)	0.706	0.0491
kBOD	=	C:L + longitude + mean depth	0.830	0.0364
IscDOC	=	C:L + IscT + SWAb0	0.858	<0.01

SWAb1 = PAR light attenuation coefficent; SWAb0 = non-PAR light attenuation coefficient; c_shelter = wind sheltering coefficient; IscV = inflow volume; IscT = inflow temperature; IscO = inflow DO; IscDOC = inflow DOC; kSOD = sediment oxygen demand; kBOD = water column oxygen demand; C:L = catchment area/lake area; SC:L = subcatchment area/lake area.

Derivation of model parameters for the data-poor lakes

The parameter values for the lakes with no or sparse data were determined using the following procedure. First, we determined relationships between lake and catchment characteristics (Table 1), which included the predictors and the values of the optimized parameters of the 12 focus lakes (i.e., the response variable). To select the potential predictor candidates, the variance inflation factor (VIF) was used to evaluate the possible collinearity between predictors. Six of the 10 parameters listed earlier were considered potential predictors with a threshold of VIF < 7. A stepwise regression approach (SR) was done by using the step function in base R (R Core Team 2020) to generate generalized linear models. Predictors were either log-transformed, taken as a square root (for areas), or left as-is when needed for normalization and to reduce heteroscedasticity (Table 3). We ranked relationships using the Akaike information criterion (AIC) and the Bayesian information criterion (BIC), metrics developed to provide an objective estimate of the model error (lower values are better model fit). In addition to R² and RMSE, the AIC and BIC criteria include a penalty that increases the error when including additional terms in the regression. Thus, they help choose the more parsimonious model, which is a sensible approach (Hastie et al. 2018). According to the handbook of statistics (Narisetty 2020), AIC is a function of residual sums of square of the corresponding model and of a term that penalizes the model having a higher complexity. BIC builds on AIC but imposes a larger penalty for higher complexity models (formulation details in Narisetty (2020). For each parameter, we selected the models that ranked first according to both criteria, as calculated in R. The final equation for each parameter was finally selected considering the RMSE, p-value, and R² for each comparison. This exercise yielded equations relating lake and catchment characteristics with model parameters obtained for the 12 focus lakes. Except for the predictor of c_shelter, none of the models had strong indications of violating the standard assumptions of linearity (Supplemental Fig. S1). We then used those equations, with catchment characteristics as input parameters, to generate lake-specific values for the chosen model parameters (Table 3).

Simulations

Using the parameter values obtained from the SR equations, state variables for the 210 lakes were simulated at a daily time step for 5 decades (2 of which were historical periods and 3 were future periods) and the 2 climate scenarios RCP 4.5 and RCP 8.5 using projections of the 6 CORDEX climate models (Supplemental Table S1). The calculations led to 10 080 decadal simulations comprising the following outputs: vector of ice cover presence or absence, matrices (day, depth) of water temperature (°C), DO concentrations (mg L⁻¹), and total photon flux (µmol m⁻² s⁻¹). Prescribed initial conditions were set to isothermal water temperature and DO at saturation (Golub et al. 2022).

Results

Water temperature, DO concentration, and SD in the 12 focus lakes

In all 12 focus lakes, water temperature followed a dimictic seasonal pattern, with the warmer temperature at the bottom of the lake during winter, an isothermal profile during seasonal mixing, and warmer temperatures at the surface during summer. Smaller lakes had larger seasonal variations in bottom temperature than deeper lakes. Field data from lake #2 (Fig. 2) show that over the whole dataset, lake temperatures varied from 3.7 °C at the bottom to 22.9 °C at the surface in summer and from 4.2 °C at the bottom to 0.2 °C at the surface in winter. Overall, DO concentration was higher in winter than summer (Fig. 3), varying between 4.37 and 12.62 mg L⁻¹ in summer and 5.35 and 15.13 mg L⁻¹ in winter. The 12 focus lakes encompass both clear and coloured lakes, with the deepest SD at 17.5 m in lake #12 and the shallowest at 0.5 m in lake #1.

Figure 2. Observed daily temperature above (red squares) and below (blue squares) the thermocline at lake #2, along with simulated daily temperature using parameters determined by the genetic algorithm (GA; solid line) or estimated by stepwise regression (SR; dashed lines). The grey zones represent simulated periods of ice cover.

Figure 3. Observed daily DO concentrations above (red circles) and below (blue circles) the thermocline at lake #2, along with simulated daily DO using parameters determined by the genetic algorithm (GA; solid line) or estimated by stepwise regression (SR; dashed lines). The grey zones represent simulated periods of ice cover.

Calibration performance

Calibration performance for temperature, DO, and light (SD) (Table 4) shows that the model performs better with more data, such as lake #2 (Fig. 2 and 3). For all lakes and calibrations, the best fit for temperature was lake #12 (RMSE_Temp = 1.72 °C), for DO was lake #9 (RMSE_DO  = 0.76 mg L⁻¹), and for SD was lake #1 (RMSE_SD = 0.45 m). Conversely, the worst performance for temperature was lake #1, for DO was lake #4, and for SD was lake #8 (Table 4). The SR approach found significant relationships (p < 0.05) between the parameters and the lake characteristics (i.e., the predictors), which varied in terms of predictive power, with R² values ranging from 0.57 for SWA_b1 to 0.89 for the inflow scalers I_scV and I_scDOC (slopes of estimated relationships in Table 3 and summary in Fig. 4; performance of GA and SR summary in Fig. 5 and 6)

Figure 4. Relationships between parameter values estimated with genetic algorithm (GA) and optimized using linear regression models obtain with stepwise regression (SR) approach (Table 3). Parameters influencing water temperature are identified with squares and those influencing DO with circles. SWA_b1 = PAR light attenuation coefficient, SWA_b0 = non-PAR light attenuation coefficient, c_shelter = wind sheltering coefficient, Isc_V = inflow volume, Isc_T = inflow temperature, Isc_O = inflow DO, IscDOC = inflow DOC, k_SOD = sediment oxygen demand, k_BOD = water column oxygen demand.

Figure 5. Target diagrams summarizing model performance for temperature (squares, upper panel) and oxygen (circles, lower panel) using parameters constrained by the genetic algorithm (solid symbols) or estimated by stepwise regression (open symbols) for the 12 focus lakes. The blue squares and red circles show the performance for Lake #2. The diagrams used the normalized bias and normalized centered RMSE between the value observed and simulated from 2001 to 2010 for all depths measured.

Table 4. Model performance summary for temperature, dissolved oxygen (DO), and Secchi depth (SD) for the focus lakes with RMSE, normalized RMSE, and R² of the accepted parameter set for genetic algorithm (GA) and stepwise regression (SR).

Lake	Max depth (m)	Number of sampling day	Calibration method	RMSE			Normalized RMSE			R^2
				Temp. (°C)	DO (mg L^−1)	SD (m)	Temperature	DO	SD	Temperature	DO	SD
1	5	48	GA	4.97	2.39	0.45	0.21	0.19	0.26	0.87	0.34	0.05
			SR	4.99	3.08	0.50	0.21	0.25	0.29	0.87	0.14	0.01
2	17	700	GA	3.08	2.17	1.18	0.14	0.17	0.33	0.72	0.29	0.0002
			SR	3.38	2.48	0.93	0.16	0.19	0.26	0.69	0.18	0.00
3	18.9	345	GA	3.99	1.67	1.28	0.18	0.15	0.23	0.82	0.24	0.14
			SR	3.97	1.69	5.91	0.17	0.15	1.05	0.82	0.34	0.25
4	20	245	GA	3.46	3.89	1.60	0.16	0.30	1.60	0.49	0.15	0.50
			SR	3.46	4.94	0.40	0.16	0.37	0.40	0.49	0.00	0.33
5	28.63	174	GA	3.75	2.12	0.80	0.18	0.18	0.36	0.80	0.41	0.0005
			SR	3.71	2.86	0.93	0.17	0.24	0.42	0.80	0.31	0.01
6	35.9	57	GA	2.95		1.71	0.17		0.23	0.95		0.36
			SR	3.25		4.11	0.18		0.55	0.92		0.00
7	45.08	131	GA	4.07		0.63	0.19		0.23	0.88		0.01
			SR	4.50		1.96	0.21		0.72	0.89		0.08
8	53	134	GA	3.42	2.41	2.20	0.14	0.18	0.29	0.85	0.26	0.01
			SR	3.66	4.66	3.02	0.14	0.34	0.39	0.77	0.06	0.00
9	53.53	123	GA	4.40	0.76	0.84	0.19	0.16	0.23	0.96	0.82	0.23
			SR	4.49	5.68	1.76	0.19	1.24	0.49	0.96	0.43	0.24
10	74	275	GA	2.47	2.62		0.13	0.42		0.80	0.07
			SR	2.50	2.67		0.13	0.43		0.81	0.04
11	74	256	GA	2.57	1.24	1.68	0.12	0.17	0.46	0.85	0.31	0.53
			SR	2.57	1.86	1.18	0.12	0.25	0.32	0.85	0.26	0.00
12	128	7017	GA	1.72	2.17	1.65	0.09	0.23	0.24	0.89	0.25	0.03
			SR	2.01	2.79	2.41	0.10	0.30	0.36	0.84	0.11	0.00

Simulated versus observed ice cover durations were plotted to validate the ice module, with ice cover data available for 47 of the 210 lakes. The comparison considers both parameterization methods (i.e., GA for the 12 focus lakes and SR for the other lake; Fig. 7). The results show that the 2 methods present a similar duration of the ice cover, with an overestimation for the lakes with the shortest ice cover duration.

Figure 6. Simulated (y-axis) versus observed (x-axis) Secchi depth range, water temperature (square), and oxygen (circle) in the 12 focus lakes. Simulations were performed using parameters determined with genetic algorithm (left panels, solid symbols) or estimated using stepwise regression (right panels, open symbols). Sampling depths for observations are represented by the colour scale from shallow (red) to deep (blue), with relative depth ranging from 0 (Surface) to 1 (z_max). The dashed line shows the 1:1 line and the solid line a linear regression, and the light grey zones encompass 95% of the simulations.

Figure 7. Observed versus simulated annual average ice cover duration for 1971–1980 for 47 lakes with available ice observations. Simulations were performed using parameters determined with genetic algorithm (black diamonds) or estimated using stepwise regression (white diamonds). The dashed line shows the 1:1 relationship.

Projected changes over a 120-year period

The 210 lakes present a large range of future surface temperatures, ice cover duration, and DO concentration (Fig. 8). Within the limitations of the current model performance during calibration on the 12 focus lakes, and neglecting the possible evolution of phosphorus loads and thus of primary production over the next decades, such a framework at the regional scale can be used to explore the response of those 210 lakes to ongoing climate change. Under the Danish Meteorological Institute (DMI) model climate projections, when comparing past (1971–1980) and future (2091–2100) decadal periods, the model predicts on average (1) an increase of 0.75 °C (RCP 4.5) or 2.1 °C (RCP 8.5) for the surface temperatures, (2) a decrease of ice cover duration of 12.6 d (RCP 4.5) or 48.4 d (RCP 8.5), and (3) and a decrease of 0.31 mg L⁻¹ in bottom water DO under RCP 4.5 and, conversely, an increase of 0.12 mg L⁻¹ under RCP 8.5. The density plots also reveal a unimodal distribution for the surface temperatures and multimodal distribution for the ice-cover period, with larger lakes changing less than the smaller lakes. DO concentration at the bottom of the lakes presents a unimodal distribution, with slight changes in most lakes. Results generated from the RCP 8.5 scenario show a much more pronounced change than from the RCP 4.5 scenario, especially for temperature and ice cover duration.

Figure 8. Density curves showing the distribution of change in simulated surface temperature, bottom DO concentration, and ice cover duration between past (1971–1980) and future (2091–2100) decadal periods. Change is calculated as the difference between daily averages simulated for the historical 1971–1980 and future 2091–2100 period, using either the RCP 4.5 (grey line) or the RCP 8.5 (coloured line) scenario. The change in ice cover is the difference between the average ice cover duration of each period divided by the number of decades between past and future periods. Surface temperatures were calculated using simulations of the water layer at 0.5 m while bottom temperatures used the simulations of the last water layer before maximum depth.

Discussion

Performance of the model and comparison with previous studies

The model performance reported here is lower than in previous studies. RMSE for temperature and DO are at the higher end of values reported previously. In particular, we found that DO is less well predicted at depth (Markelov et al. 2019, Pilla and Couture 2021). MyLake has historically performed well when applied to single lakes. For example, low RMSE ranging from 0.02 to 1.20 °C for temperature and 1.67 to 1.70 mg L⁻¹ for DO concentration were reported for the modeling of Lake Langtjern, Norway (Couture et al. 2015, 2018). Applications of the lake version of the General Ocean Turbulence Model (GOTM-Lake) at Lake Langtjern (Clayer et al. 2021) reported RMSE values of 1.33 °C and 1.9 mg L⁻¹ for water temperature and DO concentration, respectively, in line with the values obtained using MyLake at the same site. Low RMSE for temperatures ranging from 0.63 to 1.00 °C were also reported using MyLake at Lake Kuivajarvi, Finland (Kiuru et al. 2018). For 1D lake models (Golub et al. 2022), performance has been satisfactory in multi-model studies using LakeEnsembleR, which included the MyLake model among others, with RMSE ranging from 0.5 to 1.3 °C depending on the model evaluated (Moore et al. 2021).

Past modeling experiments suggest, however, that the calibration performance tends to decrease when an increasing number of lakes are calibrated together. For instance, in a multi-lake study of 26 Minnesota lakes, the MINLAKE2012 model yielded RMSE values of 1.57 °C and 1.72 mg L⁻¹ for temperature and DO, respectively (Jiang and Fang 2016). An RMSE of 2.79 °C was reported for the analysis of 10 774 lakes using uncalibrated iterations of the GLM model (Winslow et al. (2017). Both studies reported RMSE higher than reported for single-lake applications but lower than those reported here. The lack of lateral homogeneity among the lakes selected here may play a role in the reported model performance. Nevertheless, despite high RMSE values of the GA on the 12 focus lakes, simulations using parameters from SR for those 12 lakes yield RMSE values close to those derived from GA. The fit indicates that the parameter inference technique works as intended.

Relationships between parameters and lake morphometry

The SR analysis yielded coherent results, in line with relationships previously observed in boreal lakes at a large scale (Hastie et al. 2018). All parameters were derived in relation to 2 or more predictors (Table 3). The target graph (Fig. 5) shows that calibrating the model based on SR does not decrease performance when compared to a GA calibration except for light, where GA yields better results (Fig. 6). While this relationship is not a strict validation of the model for the data-poor lakes, it tends to support the idea that the determination of the parameters from morphometric characteristics is reasonable. The model shows a cold, low-DO bias for most lakes regardless of the calibration method, with normalized RMSE (i.e., with no units) varying between 0.1 and 0.21 for the temperature and 0.15 and 1.24 for DO concentration.

Projections for data-poor lakes

On average, surface water temperature increases faster in shallow lakes than in deep lakes, consistent with shallow lakes being more responsive to air temperature and therefore more sensitive to climate (De Stasio et al. 1996, Adrian et al. 2009, Fang and Stefan 2009, Dibike et al. 2011, Missaghi et al. 2017). Ice cover diminution is also in line with trends observed in Sweden from historical data (Hallerbäck et al. 2021) and forecasted by previous studies (Guzzo and Blanchfield 2017, Woolway et al. 2020). Our results project a reduction of ice cover duration of ∼4.03 (median 5.1) days per decade between 1971 and 2100 under the RCP 8.5 scenario (Fig. 8), consistent with previous estimates (Bartosiewicz et al. 2016) and within the range of 4–12 days per decade projected for subarctic lakes for 2040–2079 (Dibike et al. 2011). Our projected decrease of DO under RCP 4.5 (Fig. 8) is in agreement with that observed globally (De Stasio et al. 1996, Jane et al. 2021). Conversely, RCP 8.5 leads to a slight increase in DO concentrations, reflecting a combined response to both disappearing ice cover and enhanced summer metabolism.

Sources of uncertainty in our projections arise from the lack of future projections for DOC loads in the lakes. Changing DOC would affect all 3 key state variables (Toming et al. 2020). Uncertainties also arise from the dataset we used for calibration (Fang and Stefan 2009, Missaghi et al. 2017). Further work is needed to quantify the data frequency that maximizes model performance (Jackson-Blake and Starrfelt 2015) and to quantify the amount of focus lakes necessary for large-scale studies. Specifically, the equations used to estimate our parameters are only based on the variability within a small sample of lakes that cannot exhibit all the variability in characteristics that continental-scale studies have achieved. Building on the presented method, robust relationships could be fitted between characteristics and parameters drawing on a larger lake dataset. Finally, because DO modeling is more complex than modeling water temperature only, it bears more parameter uncertainty.

Conclusions

The amount and nature of data available for each lake are likely to vary widely, ranging from minimal, such as morphology and sparse chemistry data (Jacobson et al. 2010, Van Zuiden et al. 2016), to comprehensive, such as locations with intensive time-series monitoring stations (Fang and Stefan 2009, Fang et al. 2012, Couture et al. 2015). Real-world datasets likely to be encountered by lake managers and governmental agencies are expected to comprise both lakes with sparse data and lakes with extensive time series. We used a technique in which the parameters for a group of data-poor lakes could be inferred from the parameters of a smaller group of data-rich lakes and morphometric properties (predictors) of these lakes. This technique provided the values of the parameters suitable to correctly simulate water temperature, DO concentration, ice cover duration, and light conditions in 198 data-poor lakes based on the parameters obtained from the calibrations of 12 data-rich lakes.

Acknowledgements

We thank Koji Tominaga (Nanyang Technological University, Singapore) and Benjamin Laken (Cervest Inc., London, United-Kingdom) for the retrieval and preparation of the climate data.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

RMC acknowledges funding from the Sentinel North program of Université Laval, made possible in part thanks to funding from the Canada First Research Excellence program. Support from the Natural Sciences and Engineering Research Council of Canada, through the Discovery Grant program, from the Advancing climate science in Canada project “Changing carbon sinks in subarctic Canada” and from the Institut nordique du Québec (INQ) is also acknowledged. GE, DOH, TA and AGF acknowledge support from the Research Council of Norway projects #224779 and #221410.