Rising methane emissions from boreal lakes due to increasing ice-free days

The Arctic Lake Biogeochemistry Model (ALBM) is a one-dimensional process-based lake biogeochemical model designed for predicting both thermal and carbon dynamics of aquatic ecosystems. It mainly consists of several modules including those for the radiative transfer, the water/sediment thermal circulation, the water/sediment biogeochemistry, and the gas diffusive and ebullition transportation. Although this model was originally developed for Arctic lakes (Tan et al 2015a), it can achieve comparable or even better performance than those widely used one-dimensional lake models in representing the physical and biogeochemical processes of other northern lakes (Guseva et al 2020). Detailed information about ALBM can be found in Tan et al (2015a, 2017, 2018). We introduce the key governing equations of methane processes in ALBM below.

Methane production rate P (μmol m⁻³ s⁻¹) in sediments is calculated from labile carbon content and temperature:

$$\begin{equation}P = {R_c}{C_{labile}}PQ_{10}^{\left( {T - {T_{pr}}} \right)/10}\end{equation} \tag{ 1 }$$

where ${R_c}$ is the fraction of carbon decomposed per second (s⁻¹), ${C_{labile}}$ is the labile carbon content (μmol m⁻³), PQ₁₀ is the factor by which the production rate increases with a 10 °C rise in temperature, and T_pr is a reference temperature (°C). Methane can be oxidized after released into the water and the oxidation rate V_oxid (μmol m⁻³ s⁻¹) is described by

$$\begin{equation}{V_{oxid}} = {Q_{C{H_4}}}OQ_{10}^{\left( {T - {T_{or}}} \right)/10}\frac{{{C_{{O_2}}}}}{{{k_{{O_2}}} + {C_{{O_2}}}}}\frac{{{C_{C{H_4}}}}}{{{k_{C{H_4}}} + {C_{C{H_4}}}}}\end{equation} \tag{ 2 }$$

where ${Q_{C{H_4}}}$ is the maximum oxidation potential (μmol m⁻³ s⁻¹), OQ₁₀ is the factor by which the oxidation potential increases with a 10 °C rise in temperature, T_or is a reference temperature (°C), ${C_{{O_2}}}$ and ${C_{C{H_4}}}$ and gas concentrations (μmol m⁻³), and ${k_{{O_2}}}$ and ${k_{C{H_4}}}$ are the Michaelis-Menten constants (μmol m⁻³). Together, the modeled methane concentration in water columns is calculated by

$$\begin{equation}\frac{{\partial {C_{C{H_4}}}}}{{\partial t}} = \frac{\partial }{{\partial z}}\left( {{D_{C{H_4}}}\frac{{\partial {C_{C{H_4}}}}}{{\partial z}}} \right) - {V_{oxid}} \pm {L_{C{H_4}}}\end{equation} \tag{ 3 }$$

where ${D_{C{H_4}}}$ is the diffusivity of methane (m² s⁻¹), and ${L_{C{H_4}}}$ is the gas exchange between bubbles and the ambient water (μmol m⁻³ s⁻¹). Finally, methane within water is transported to the atmosphere. The diffusive transfer velocity k (m s⁻¹) is defined as

$$\begin{equation}k = \sqrt {{{\left( {0.00015U} \right)}^2} + {{\left( {0.07{{\left( { - \beta {Z_{AML}}} \right)}^{\frac{1}{3}}}} \right)}^2}} \times S{c_m}^{ - 0.5}\end{equation} \tag{ 4 }$$

where U is the wind speed at 2 m (m s⁻¹), β is buoyant flux (β < 0 if losing heat and vice versa), z_AML is the depth of the actively mixing layer (m), and Sc_m is the Schmidt number of methane. Since we lacked ebullition flux observations and therefore, we were unable to validate the modeled ebullition emissions, we only quantified diffusive emissions in this study.

The numerical experiment consists of three steps: (1) the model calibration using observations of diffusive emissions during 1998–1999 from 39 individual lakes; (2) regional simulations of 1990–1999 by applying ALBM to the Ranta10 data product; (3) regional simulations of 2010–2019 and 2090–2099 under the representative concentration pathway (RCP) scenarios of RCP4.5 and RCP8.5. For all the simulations, a spin-up period of two years was run.

Model forcing data include air temperature, surface pressure, 10 m wind, relative humidity, precipitation, snowfall, downward short-wave radiation and downward long-wave radiation. The historical simulation was driven by the climate data retrieved from the European Center for Medium-Range Weather Forecasts (ECMWF) Interim re-analysis (ERA-Interim, Dee et al 2011) with a resolution of 0.75° × 0.75°, and organized into daily datasets. For future climate scenarios, we used a down-scaled bias-corrected dataset of the Intersectoral Impact Model Intercomparison Project (ISIMIP) output from HadGEM2-ES (Frieler et al 2017) that is set on a 0.5° × 0.5° global grid and is divided into daily time steps. This climate dataset is bias-corrected based on ERA-Interim, which guarantees the consistency of historical and future simulation results.

Lakes with an area smaller than 200 m² were omitted in our simulations due to the large uncertainties in the mapping of these lakes, leaving 176 876 lakes covering 36 690 km². In general, the region north of 67 °N has the highest lake density but relatively sparse observations of thermal or carbon dynamics (figures 1(a) and 1(b)). Over 90% of the lakes are smaller than 0.1 km² (figure 1(c)), which are not included in the GLWD-3 database. Depth information was lacking for over 90% of the lakes even by combining the Ranta10 and the GLWD-3 database. As such, we applied a statistical approach to construct the full lake depth dataset. We first grouped the lakes into 10 bins bounded by areas of 0, 0.01, 0.05, 0.1, 0.5, 1, 2, 3, 4, 5, 10, 100, 1500 km², respectively. We then generated a histogram of depths in each group. We randomly assigned the depths following the fitted probability distribution in each group. By following this approach, we aimed to construct lakes profiles that match the diversity of the real lake system in Finland. In terms of lake bathymetry, we assume a linear decrease of the cross-section area with increasing depth.

Zoom In Zoom Out Reset image size

Water temperature, nutrient and methane concentrations were measured at four levels of depth four times (before and after ice melt, late summer and autumn circulation) either in the 1998 or 1999. The diffusive methane fluxes were calculated following equation (4) (see also Juutinen et al 2009). For the simulation purpose, several assumptions and approximations were made on for carbon and phosphorous input that is required by the model, including (1) dissolved organic carbon (DOC) concentration is equal to total organic carbon (Mattsson et al 2005, Kortelainen et al 2006, Mopper et al 2006); (2) dissolved inorganic carbon (DIC) concentration is calculated from pH and alkalinity; (3) soluble reactive phosphorus (SRP) concentration is equal to PO₄ and (4) particular organic carbon (POC) concentration is approximately 1/5.1 of DOC concentration (Rachold et al 2014). We produced an input map of DOC, DIC, PO₄ and SRP at 1° × 1° with observations averaged at each grid and filled with nearest-neighbor interpolation (figure S1 (stacks.iop.org/ERL/15/064008/mmedia)).

Since the lake shape has large impacts on the thermal dynamics (Mazumder and Taylor 1994, Woolway et al 2016, Woolway and Merchant 2017) and carbon dynamics (Schilder et al 2013, Pighini et al 2018), we decided to conduct calibration and thus simulations by lake groups. The simulated lakes were firstly divided by the surface area into six groups bounded by 0, 0.05, 0.1, 0.2, 1, 10, 1500 km², respectively and then by the shape factor, defined by $\sqrt {{\text{Area}}} /{\text{Depth}}$, into two groups bounded by 0, 0.1 and 10, respectively. Thirty-nine lakes that represent various depths and shape factors were used for calibration.

We calibrated the parameters related to water temperature (Ks, Cps, Pi, Rous, Roun, Feta, Wstr and Ktscale) and methane diffusive emission (OQ10, QCH4, KCH4, PQ10 n, Rcn and Rca as in equations (1) ∼ (4)). The detailed descriptions and corresponding ranges of the parameters are listed in table S1. Since the sensitive parameters of the water temperature and the methane diffusive emission simulations were different, the calibration was conducted separately. We first applied a Monte Carlo calibration method for water temperature calibration using 6000 sets of parameters for each lake. For the methane emissions, a 'history matching' method (Mcneall et al 2013, Williamson et al 2013, 2015, 2016) was adopted for higher accuracy and efficiency. This method requires less computation time than a full Monte Carlo simulation (Williamson et al 2013) and has been shown to largely reduce simulation biases (Salter et al 2019). It was carried out in the following steps: (1) use the Sobol sequence sampling method (Sobol' 1967) to generate a perturbed parameter ensemble (PPE) by sampling from the parameter space; (2) run simulations for all perturbed parameter (PP) sets; (3) rule out regions of the parameter space based on the outputs where a predefined metric exceeds the threshold; (4) repeat step 1–3 until a certain number of iterations or the desired outcome is achieved. In our study, 1200 PP sets were generated for each round and the metric used was the root-mean-square error (RMSE). Parameter spaces resulting in RMSEs over the 50% of the observations at each site were ruled out at each round, until 3 rounds were finished.

For parametric sensitivity analysis, running the model for the whole region over a 10-year time period takes about five days, and thus it would be rather time-consuming to run full PPE simulations. Instead, we run short PPE simulations for a single year from 1998 to 1999 with 20 PP sets sampled from the remaining parameter space after history matching. It has been proved that the results from short simulations match well to the longer-term simulations (Yun et al 2016), especially when methane emission response to air temperature is a relatively well-defined process (Sanches et al 2019) and thus can be captured within a one-year range.

The model overall showed good performances on reproducing both water temperature and methane emissions (table 1) with an average RMSE of 1.59 °C for full-profile water temperature simulations and a correlation coefficient (r) of 0.69 for methane emissions. There was a lake with an estimated methane flux over 200 mmol m⁻² yr⁻¹ while the corresponding modeled value was only one-third of the observed value. This underestimation was possibly due to the missing of DOC concentration measurement at this lake. Since it is a very humic shallow lake (color >90 Pt mg l⁻¹, mean depth <3 m, see Juutinen et al 2009), taking the average value of other lakes in the same grid likely results in a much lower DOC concentration. The metrics without considering this lake are also listed in table 1.

Note that the uncertainty range can be wider using the history matching approach than using the Bayesian method, which is expected because the former focuses on confining the output whereas the latter on confining the input, i.e. the parameters themselves. Therefore, the PPE by the latter would be more representative of the parameter distributions and thus results in smaller uncertainties.

Simulations indicated that the methane emissions from Finnish lakes were 0.12 ± 0.03 Tg CH₄ yr⁻¹ in the 1990 s. There was only a 4% and 6% increase during the period 1990 ∼ 2019 under RCP4.5 and RCP8.5, respectively (figure 2) when less than 2 °C of warming has occurred. Walter et al (2007) estimated that the current total diffusion from all lakes north of 45 °N is 1.12 ± 0.22 Tg CH₄ yr⁻¹. If assuming the same total lake area over the Arctic and the same lake size distribution as the Finnish lakes (Text S1), we estimated the emission to be 3.65 ± 1.06 Tg CH₄ yr⁻¹. The difference can be due to: (1) Walter et al (2007) simply used a constant flux calculated from several glacial lakes and thermokarst lakes for all other lakes, which may underrepresent the variation; (2) Based on several Siberian thaw lakes, an ice-free period of 120 d was assumed in the calculation for all lakes, leading to underestimation in warmer regions, for example, the mean ice-off periods in Finland is about 170 d.

Zoom In Zoom Out Reset image size

We estimate that the Finnish lake diffusive methane emission will increase by 25.8% from 0.12 ± 0.04 to 0.16 ± 0.05 Tg CH₄ yr⁻¹ under the RCP4.5 scenario while it will increase by about 58.9% from 0.13 ± 0.04 to 0.20 ± 0.06 Tg CH₄ yr⁻¹ under the RCP8.5. The magnitude of the growth is relatively small compared to Tan et al (2015a) and Tan and Zhuang (2015b)) who predicted an 80% increase in Northern Europe even under the RCP2.6 scenario. It is likely because we did not model ebullition emission. It has been found that future warming may have much larger effects on ebullition even altering diffusive-emission dominant lakes to ebullition-dominant ones (Aben et al 2017). Therefore, the amount of extra increase is likely due to the enhanced ebullition processes.

Based on the analysis of 297 lakes worldwide, Sanches et al (2019) found that considering only diffusive emissions would cause an average underestimation of 277%. By taking into account the ebullition emissions, the current annual methane emission for all northern lakes would be 13.76 ± 3.99 Tg CH₄ yr⁻¹, in the lower range of the previous estimation, 24.2 ± 10.5 Tg CH₄ yr⁻¹ by Walter et al (2007) and about the same as the 16.5 ± 9.2 Tg CH₄ yr⁻¹ by Wik et al (2016). Furthermore, Aben et al (2017) predicted that ebullition would increase faster than diffusive emissions, by 51% with 4 °C warming. Similarly, Thornton et al (2015) predicted an increase of around 56% from the present to the 2040–2079 period using observations from three subarctic lakes in Sweden. If, based on the warming magnitude, we estimate a 50% and 100% growth of ebullition emissions under the RCP4.5 and RCP8.5, respectively, this will result in a total methane emission of 20.02 ± 5.8 Tg CH₄ yr⁻¹ and 27.97 ± 8.11 Tg CH₄ yr⁻¹, respectively from the whole Arctic lakes. This projection is slightly lower than the previously modelled 32.7 ± 5.2 Tg CH₄ yr⁻¹ from lakes north of 60 °N using the GLWD map by Tan et al (2015b). Their potential overestimation can result from: (1) They calibrated the model using only five lakes, four of which are thermokarst lakes that are found to have higher ebullition to diffusive emission ratio than other types of lakes (Wik et al 2016). According to their observation, this ratio can be over 30 which is much higher than we assumed; (2) They also assigned a shallow depth of 3 m to all lakes missing depth information, leading to even more overestimation of the ebullition emission (Joyce and Jewell 2003, Bastviken et al 2004, Del Sontro et al 2016).

Methane emissions are generally highest in regions dense with small lakes (figures 1(a) and 3), a pattern which was also found in previous studies (Bastviken et al 2004, Saarnio et al 2009, Del Sontro et al 2016, Wik et al 2016). This finding can be explained with two mechanisms. First, methane can be oxidized along diffusive transportation, and thus deeper lakes usually mean more loss by oxidation. Second, it was found that smaller lakes are more likely to have abundant organic substrates in their sediments, and thus they are potentially more productive for methane.

Zoom In Zoom Out Reset image size

Interestingly, the southern and the northern part of Finland show different trends. By the same 4 °C of air temperature warming (figure S2), the south experiences a much more severe increase in emissions than the north (figure S3). Here we define the south and north be distinguished by latitude 67.5 °N, which north and south mean annual air temperature is above or below 0 °C in the 2010s, respectively. Therefore, the absolute increase in air temperature is not a complete proxy for emission. Instead, we found that the length of ice-free days plays a key role. By the end of the century, the mean ice-free days in the north are extended by 10.8% and 23.5% under RCP4.5 and RCP8.5, respectively while they are extended by 13.2% and 34.1%, respectively in the south (table S3). Methane fluxes can be blocked by ice cover and then oxidized in the water column during ice-on seasons. Therefore, a large amount of methane trapped by ice cover currently could be emitted into the atmosphere in the future. Generally, ice covers of lakes in the north are thicker and thus would take more energy input to melt before methane can be released in winter.

The influence of ice-on days is further verified by the seasonality of emission (figure 4). From May to November, the warming climate affects emission in a similar pattern in both regions. However, the difference is shown in the south from December to April when we hardly see any emissions in the 2010s, but we expect much higher emissions in the future for the period due to warming. Such a shift may not be initiated in the north within this century and thus, the increase of emission in the north is much less severe. Wik et al (2014, 2016) also indicated that the ice-free days could potentially be the primary proxy for emission based on the observation of three lakes. Here, we confirm the finding with all the lakes in Finland and further verify it from a seasonality perspective. They predicted a 30% increase in total emission from northern lakes within the century by assuming a 20-day increase of ice-free days for all lakes. However, we emphasize that differentiating assumptions on ice-free day increase among regions is necessary for better precision in the estimation of regional budgets and future projections. Additionally, our simulated increase of ice-free days under RCP8.5 is two to three times higher than their assumption, indicating an underestimation in the previous study.

Zoom In Zoom Out Reset image size

Our model simulations did not consider the impacts of DOC dynamics on lake methane emissions. Boreal lakes, in addition to the climate change, have experienced moderate to severe browning over the last decades, probably due to the increased import of DOC from soils (Roulet et al 2006, Monteith et al 2007, Haaland et al 2010, Larsen et al 2011, Seekell et al 2015, Isidorova et al 2016). If DOC concentrations increase at the current rate, they can be doubled by the end of this century (de Wit et al 2016). It is found that methane diffusive methane emission is positively related to DOC concentrations (Sanches et al 2019). However, the relationship can involve several processes. More nutrients will be available for microbes and primary producers. However, the resulting increases in turbidity will weaken photosynthesis and thus oxygen production which limits methane oxidation. The model will need to be improved to predict the effect of lake browning. So far, no study has included this effect when estimating future emissions.

We assumed a constant landscape in our simulations that no lake expansion or drainage was considered until 2100. This is because boreal lakes in Finland are not formed over permafrost and thus not affected by the active response of thawing and groundwater penetration processes in responding to climate as thermokarst lakes (van Huissteden et al 2011, Katey et al 2018). Wik et al (2016) predicted that with a 20-day increase of ice periods, even the total lake area decreases by 30%, the total methane emission can still grow by 20% which is 10% less than assuming constant lake area. Therefore, we will still expect the boreal lake methane emissions will be affected by considering the lake area changes. Incorporating lake areal dynamics into future quantification is necessary.