The role of lake morphometry in modulating surface water carbon concentrations in boreal lakes

For this exploration we are using data that were generated in a large-scale comparative study of boreal lakes that our group carried out between 2006 and 2014. The study focused on the main carbon pools in the surface waters of 265 lakes distributed across six regions in northern Québec (Canada; figure 1(a)), encompassing a wide range of climatic, topographic, and environment gradients. Various aspects of this project have been previously published, including the main limnological features of these lakes (Roehm et al 2009, Lapierre and del Giorgio 2012), the ambient DOC and DIC concentrations (Lapierre et al 2015), and CO₂ and CH₄ partial pressures (Roehm et al 2009, Rasilo et al 2015). Here we provide only a summary description of the main approaches used, further details can be found in those publications. All lakes were sampled during the late spring and summer season (from late May to September; supplementary figure 1 (available online at stacks.iop.org/ERL/16/074037/mmedia)), and the timespan of the sampling campaigns were 2006–2008 for the Eastmain region, 2010 for the Abitibi and James Bay regions, 2011 for Chibougamau and Chicoutimi, 2012 for Schefferville, and 2013–2014 for lakes in the Cote Nord. Most lakes were accessed by truck and sampled by boat when possible, and some were accessed by hydroplane. Water samples for analysis were taken at 0.5 m depth over the deepest part of the lake. Some sites were sampled more than once during the same year, generally in late spring, summer and early fall, and in these cases concentration values were averaged before data analysis. Since most of the lakes were sampled only once, note that temporal variations in lake C concentrations throughout the late spring and summer periods might not be fully represented in our analyses.

Sampling and analysis of DIC, DOC, CO₂ and CH₄ are described in detail in Roehm et al (2009), Lapierre and del Giorgio (2012) and Rasilo et al (2015). Briefly, the in situ partial pressure of CO₂ (pCO₂) was measured in the early years using a PP-Systems EGM-4 infrared gas analyzer and a membrane that equilibrates the air with the water sample (Roehm et al 2009). The partial pressure of CH₄ (pCH₄) was measured by triplicate using the headspace method with 60 ml polypropylene syringes. Thirty milliliters of lake water were equilibrated with 30 ml of ambient air by vigorously shaking the syringes for 1 min. The equilibrated headspace was then injected to a 30 ml glass vial, stored in cold and dark conditions, and then analyzed in a gas chromatograph in the lab. In the later years (starting in 2010), both methods have been replaced by a new headspace method combined with a portable greenhouse gas analyzer (Los Gatos Research (LGR), California). In this case, the headspace technique was performed in a 1 l glass bottle equipped with to valves, with a 1:1 volume of water and ambient air. Bottles were vigorously shaken for 3 min, and the headspace was then extracted by pumping more of the same lake water into the bottle, pushing the headspace gases into a dark, vented gas container. The equilibrated headspace was then analyzed using the LGR portable greenhouse gas analyzer.DOC and DIC were measured in lake water samples filtered through 0.45 μm cartridges and later analyzed on an OI 1010 TOC analyzer after sodium persulfate digestion.

Lake geometry was calculated in ArcGIS v.10.4 based on lake polygons obtained from the National Hydro Network, which is provided by Natural Resources Canada (https://open.canada.ca/data/en/dataset/a4b190fe-e090-4e6d-881e-b87956c07977). The polygons are derived from Landsat 7 (30 m) or higher resolution source data along with provincial NRCan 1:10 000–1:50 000 datasets, and have been rigorously evaluated using automatic or interactive validation tools at the production stage (GeoBase^® 2010). The basic morphometry of lakes was characterized through four metrics: circularity, shoreline complexity, dynamic ratio and littoral fraction, and they all involve combinations of lake area, mean depth, and perimeter. The description of each metric and their calculation are summarized in table 1. For the majority of these lakes there were no actual morphometric measurements, so volume and depth had to be modeled. The volume of each lake was modeled based on the elevation change within a system-specific buffer around the lake, following the models developed by Heathcote et al (2015) for boreal lakes. The mean depth was then calculated as the volume divided by the area of the lake. Most CH₄ ebullition in boreal lakes occurs in littoral areas shallower than 3 m (DelSontro et al 2016). Therefore, for the purpose of this study, lake littoral fraction was defined as the proportion of lake area that is shallower than 3 m (see table 1 for calculation). Most of the metrics used in our work imply a reduction from three-dimensional lake basins to a two-dimensional space. This is a powerful way to synthesize complex 3D information, but it is important to note that we also lose information that could be relevant for further understanding the role of lake morphometry in modulating lake carbon concentrations.

Table 1. Description and calculation of the lake morphometry metrics.

Metric	Calculation	Description
Circularity	$\frac{{{\text{Area}}}}{{{\text{MBC}}}}$	Circularity compares the area of the lake (m2) with the area of the minimum boundary circle (MBC; m2) around it. It differentiates between circular (value approaches 1) and elongated lakes (value approaches 0). Unitless.
Shoreline complexity	$\frac{{{\text{Perimeter}}}}{{\sqrt {{\text{Area}}} }}$	Shoreline complexity is a perimeter (m) to area (m2) ratio. Increasing values express increasing complexity. Area is square-rooted to account for size dependency. Unitless. This metric can have different names such as shoreline development or shoreline convolution.
Dynamic ratio	$\frac{{\sqrt {{\text{Area}}} }}{{{Z_{{\text{mean}}}}}}$	Dynamic ratio (Håkanson 1982) is an area (in km2) to mean depth (Zmean; m) ratio of the lake. Low values indicate lakes that are bowl-shaped, whereas high values are associated to dish-like lakes. Area is square-rooted to account for size dependency. Units are km m−1.
Littoral fraction	$1 - {\left( {1 - \frac{3}{{{Z_{{\text{max}}}}}}} \right)^{\left( {\frac{{{Z_{{\text{max}}}}}}{{{Z_{{\text{mean}}}}}} - 1} \right)}}$	Littoral fraction is the relative proportion of lake surface area that is shallower than 3 m. Zmax and Zmean are the maximum and mean depth of the lake in m.

One of the main challenges in this exercise is to isolate the influence of lake morphometry from other environmental drivers of lake C biogeochemistry, as for example the terrestrial productivity and land cover in the lake catchment. To simplify the integration of multiple environmental drivers in our exercise, we ran a principal component analysis (PCA) based on several climatic and catchment-scale variables (figure 1(b)). The first axis explained 50% of the total variance across our dataset and was linked to climatic variables such as mean annual temperature and evapotranspiration as well as annual terrestrial net primary productivity, whereas the second axis (24%) was linked to catchment properties such as the organic carbon stock and the areal proportion of wetlands in the lake catchment. In our subsequent analyses, the first two axes of the PCA are used as an integrative proxy for the different climatic and catchment drivers of lake carbon concentrations across Québec. The PCA analysis was performed in R using the package vegan (Oksanen et al 2013).

Pair-wise relationships between variables were assessed through Spearman correlation coefficients, and a correlation map was used to assess the relationship between the different size and morphometry metrics. Since the littoral fraction metric was extremely correlated with the mean depth of the lake (Spearman rho = −0.99), we excluded littoral fraction from subsequent analyses.

Elastic net regression analysis (Zou and Hastie 2005) was used to identify and rank morphometry predictors of pCO₂, pCH₄, DOC and DIC variability across lakes. We chose elastic net over stepwise ordinary least square regression because the penalization implemented in the elastic net protects the estimates of the coefficients against collinearity-induced enhancement. This is important because lake morphometry metrics are generally all correlated, albeit with different strengths (supplementary figure 2). Lake area and volume were also included in the elastic net models. Lake area was included because there is a size-scaling of lake morphometry, such that larger lakes are more fractal, less circular and more dish-like (figure 2). Thus, by including lake area in the models we can separately assess the influence of lake size and morphometry on the concentration of the different C species. Lake volume was included because larger lakes have a higher dilution capacity and generally present longer water residence times. In addition, we included the first two axes of the environmental PCA (figure 1(b)) to account for the influence of climatic and catchment-scale drivers on lake carbon concentrations. To fit the elastic net models, the response and explanatory variables were standardized to zero mean and unit variance. Thus, the absolute size of the estimated coefficients we report reflects the relative importance of each explanatory variable. Since accurate standard errors are not available for regularized regression models, it is not advisable to estimate p-values in Elastic Net models. Instead, it is recommended to evaluate the effectiveness of each variable based on whether its coefficient was shrunk to zero or not by the model. Predicted versus measured values for the four models are shown in supplementary figure 3. The elastic net model was fitted using the glmnet package (Friedman et al 2010) in R (R Core Team 2020), with a ten-fold cross-validation to avoid overfitting.

Zoom In Zoom Out Reset image size

As mentioned above, in the vast majority of cases the volume and mean depth of lakes were calculated using the published empirical models in Heathcote et al (2015), which are inevitably associated with a prediction error that could potentially impact our results. We therefore evaluated the sensitivity of our elastic net model results to variations in volume and mean depth using a Monte Carlo simulation with 1000 iterations. We used the models' prediction error (5% for volume, 49% for mean depth; Heathcote et al 2015) to generate a Gaussian distribution of possible lake volume and mean depth values from which to randomly assign a value for each iteration. The distribution of mean depth values was then propagated to the dynamic ratio using the formula shown in table 1. The coefficient of variation of the different model coefficients after the 1000 iterations was generally very low (median = 1.5%, interquartile range = 0.4%–4.7%; supplementary table 1), indicating that our results are robust against potential errors derived from the calculation of lake volume and mean depth.

The comparison of the different morphometry metrics with lake size shows that there is an allometry of lake morphometry across the boreal region of Québec. For example, we observed a high degree of covariation between lake size and shoreline complexity, such that the shoreline of lakes become more fractal or convoluted with increasing size (figure 2(a)). This is in line with previous research by Cael and Seekell (2016), who reported a similar scaling of lake shoreline complexity with lake size globally. Likewise, our dataset shows that smaller lakes tend to be more circular and present a low dynamic ratio, whereas larger lakes are characterized by more elongated shapes and higher dynamic ratios (figures 2(a) and (b)). Note however that there is still a large degree of uncoupling between lake size and morphometry, in the sense that there is a wide range of potential lake morphometries within any given category of lake size (figure 2).

Consistent with previous studies, our data shows that the C biogeochemistry of boreal lakes is influenced by the size of the lake, yet the degree of influence strongly varies among the different C species (figure 3). Among the four C species, DOC is the only one that did not significantly covary with lake size (Spearman rho= −0.05, p = 0.41). In contrast, the pCO₂ and pCH₄ were significantly negatively correlated with lake area (Spearman rho= −0.40 and −0.61, respectively; p < 0.001 in both cases), a pattern that has been previously described in other regions of the boreal biome (Sobek et al 2003, Humborg et al 2010, Rasilo et al 2015, Rocher-Ros et al 2017), and that becomes even more apparent when climatic effects are accounted for (figure 3). This negative size scaling of pCO₂ and pCH₄ has been related to a decrease in the relative importance of terrestrial C inputs relative to lake area and especially to lake volume (DelSontro et al 2018). In addition, the influence of wind on the gas exchange between the lake surface water and the atmosphere tends to increase with lake size (Vachon et al 2013), hence lower pCO₂ and pCH₄ in larger lakes may be also explained by higher outgassing rates. DIC covaried weakly but positively with lake area (Spearman rho= 0.21, p = 0.003), which is somewhat counterintuitive considering that CO₂ is one of the constituents of DIC. However, we suggest that such a pattern may in reality reflect a covariance between geology and average lake size across regions rather an effect of lake area per se. For example, the highest DIC concentrations were recorded in the flat region of Abitibi, which is characterized by very large, shallow lakes and is also characterized by a volcanic and sedimentary geology that contributes high alkalinity to aquatic networks.

Zoom In Zoom Out Reset image size

Our results show that lake morphometry does influence the biogeochemistry of C in boreal lakes, albeit with clear differences among C species both in terms of which morphometry metrics are linked and of the strength of influence (figure 4). The model for pCO₂ indicates that, once environmental and size-volume effects are accounted for (see section 2 for details), the complexity of the lake shoreline plays an important role in shaping surface water CO₂ concentrations in lakes. Lakes with higher shoreline complexity tend to have higher pCO₂, which is most likely linked to a higher contact surface with terrestrial environments per unit lake area and therefore higher areal inputs of terrestrial CO₂ via direct groundwater and runoff inputs. In addition to shoreline complexity, the circularity of the lake seems to also slightly influence lake CO₂ dynamics, such that more circular lakes show higher surface water pCO₂. Previous research have shown that circular lakes have a stronger biophysical coupling between the littoral sediment and water column compartments (Dolson et al 2009). Also, more elongated lakes might have higher fetch, which facilitates the transfer of energy from wind to the lake surface and therefore the evasion of CO₂ (Vachon et al 2013). We therefore suggest that the influence of circularity is related to the input of CO₂ derived from sediment respiratory processes (Algesten et al 2005, Kortelainen et al 2006) as well as to the effect of fetch on lake surface water turbulence and CO₂ evasion.

Zoom In Zoom Out Reset image size

Among the four C species, CH₄ appeared to be the one most strongly related to lake morphometry, since it linked to virtually all the morphometry metrics we explored (figure 4). Compared to CO₂, a potentially larger fraction of the CH₄ in lakes is produced in situ, mostly through the anaerobic degradation of organic matter in bottom waters and sediments (Bastviken et al 2004; but see Grossart et al 2011, Bogard et al 2014). It is therefore not surprising that aspects of lake morphometry, particularly those related with the lake vertical profile, would influence surface water pCH₄. The elastic net model for pCH₄ identified dynamic ratio as the most influential morphometry metric, indicating that lakes with a dish-like profile (i.e. high dynamic ratio) tend to have higher surface water pCH₄. Higher dynamic ratio values imply a higher proportion of the water column that is in contact with the bottom of the lake. Thus, the positive effect of dynamic ratio on lake pCH₄ can probably be explained by a stronger influence of sediment and lake littoral processes such as methane ebullition (DelSontro et al 2016) and the activity of macrophytes (Bergström et al 2007). This would be in agreement with the work from DelSontro et al (2017), who recently demonstrated that a large fraction of the CH₄ in lakes is produced in littoral areas and transported to the lake center via lateral diffusion. Both shoreline complexity and circularity had a positive effect on lake pCH₄, which may reflect a higher degree of connection with terrestrial CH₄ sources in more fractal lakes, and higher influence of littoral sediments and less water turbulence in more circular lakes. Our analysis further shows that pCH₄ was slightly positively related to the mean depth of the lake. We do not have a clear explanation for this result; however, we suspect that it could reflect an increased contribution of algal-linked oxic CH₄ production in deeper lakes (Bogard et al 2014), perhaps through the more likely development and higher thickness of deep chlorophyll maxima with increasing lake depth (Leach et al 2018).

It is worth noting that, even when morphometry metrics are included in the models, lake area is still the most influential variable on lake pCO₂ and pCH₄. This indicates that the common negative size scaling of the two C gases (figures 2(a) and (b)) is probably not related to a decreasing relative influence of lake sediments and terrestrial C inputs, which is already accounted for by the shoreline complexity and dynamic ratio metrics, but rather to the higher turbulence and increased gas exchange coefficient in larger lakes (Vachon et al 2013).

DOC was positively linked to the dynamic ratio of the lake, such that dish-like lakes tend to have higher DOC concentrations (figure 4). This relationship with dynamic ratio may reflect on the balance between sedimentation and release of DOC from sediments to the water column (e.g. Dadi et al 2016, Peter et al 2016). Additionally, it is possible that dish-like lakes have increased macrophyte coverage, which has been shown to release large amounts of DOC to the surrounding waters (Demarty et al 2009). The elastic net model also revealed a weak negative influence of lake area and volume on surface water DOC (supplementary table 1), which only shows up if environmental drivers are accounted for. A tendency to lower DOC concentrations in larger lakes could be explained by dilution but also by the fact that larger volumes generally correspond to longer water residence times and therefore increased exposure to degradation processes (Catalán et al 2016, Vachon et al 2017).

The model for DIC explained only 10% of the variability across lakes, suggesting that the patterns in lake DIC concentrations are driven by an entirely different set of drivers than those of CO₂, CH₄, and DOC. Lake DIC, composed of CO₂, HCO₃ ⁻ and CO₃ ²⁻, partly derives from mineral weathering in the catchment that is routed through drainage networks into lakes (Marcé et al 2015). Therefore, lake DIC is more likely to reflect the geology of the catchment. Supporting this, we found a positive correlation between lake DIC and the percentage of volcanic and sedimentary rocks in the lake catchment (Spearman rho = 0.38 and 0.25, p < 0.001 and p = 0.006, respectively). In this regard, the positive effect of lake area and negative effect of mean depth on lake DIC (figure 4) is probably a coincidence derived from the fact that Québec regions underlain by metamorphic rocks (slow weathering; low alkalinity) tend to be steeper, and are dominated by small and canyon-like, deep lakes. In contrast, regions underlined by intrusive and volcanic rocks that produce higher alkalinity levels are much flatter and dominated by large and shallow lakes.

Here we linked lake morphometry metrics, which are rather static variables, to lake C concentrations that may vary over shorter time scales from days to seasons. This implies that the role that lake morphometry plays in modulating lake C dynamics may vary throughout an annual cycle. Therefore, it should be kept in mind that our results are dependent on the period of sampling (late spring and summer), and that they do not capture potential additional roles that lake morphometry might play in other seasons, for example through an influence on the timing and characteristics of the ice cover during the winter period (Denfeld et al 2018).

Large-scale estimates of lake C emissions have traditionally relied on averaged-based approaches, that is, multiplying an average emission rate by the regional or global surface areas covered by lakes (Cole et al 2007, Tranvik et al 2009, Bastviken et al 2011). More recently, the negative scaling of C gas concentrations with lake size has proven useful to improve the prediction of lake pCO₂ and pCH₄ in lakes. Thus, lake area is now commonly taken into account in regional and global upscaling exercises of lake C emissions (e.g. Hastie et al 2017, DelSontro et al 2018), leading to more realistic estimates. Our results raise the possibility that including lake morphometry within empirical models will further improve the predictability of lake C concentrations. To evaluate this, we compared the R² of three different models for each of the C species: (a) a model that only has the two first axes of the environmental PCA (figure 1(b)) as explanatory variables, thus accounting only for climate and catchment drivers; (b) a model that in addition includes lake area; and (c) the elastic net model shown in figure 4 that also includes the different morphometry metrics. The comparison of the three models shows that lake morphometry has only a moderate potential to improve the prediction of lake C concentrations at the regional scale (figure 5). This is probably because of the observed covariation between lake size and morphometry (figure 2), in the sense that lake area alone implicitly accounts for a significant proportion of the influence of morphometry on the biogeochemistry of lake C. Regardless, lake morphometry metrics such as shoreline complexity or circularity can be easily obtained through GIS analysis. Therefore, even though they might only moderately improve empirical models of lake C biogeochemistry, we suggest their incorporation into future upscaling exercises is nevertheless warranted (see Casas-Ruiz et al (2021) for an example).

Zoom In Zoom Out Reset image size

The different lake morphotypes that result from the various morphometry and size combinations are not randomly distributed across the boreal landscape. Rather, certain lake morphotypes are more recurrent within distinct geomorphological regions (Jakobsson 2018; Håkanson and Karlsson 1984), defining a regional 'morphogeography' of lakes across the boreal biome (Jakobsson 2018). Therefore, beyond their application for upscaling exercises, our results would imply that patterns in lake morphometry may impose a certain degree of spatial structure in boreal lake C dynamics that is independent of other catchment and climatic drivers.