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

Earth’s lakes vary greatly in size and morphometry, from small circular lakes of hundreds of m2 to large and deep fractal systems of several thousands of km2. Previous research has demonstrated a link between the size of lakes and their carbon dynamics. However, the influence of lake morphometry on lake carbon biogeochemistry remains largely unexplored. Here, we analyze the morphometry and carbon concentrations of more than 250 lakes across boreal Quebec, encompassing a wide range in lake size from 0.002 to 4300 km2. We show that, in addition to lake size, the biogeochemistry of lake carbon is influenced by the circularity, shoreline complexity and vertical profile of the lake. Yet the type and degree of influence vary among the different carbon species. A comparative exercise shows that taking into account the morphometry of lakes moderately increases the predictive power of empirical models of carbon concentration across lakes. Therefore, future studies might benefit from adding lake morphometry metrics to the empirical rules used for prediction and upscaling.


Introduction
Lakes cover only a small portion of the Earth's surface, but they are essential to many biogeochemical and ecological processes in landscapes as well as a source of freshwater and food for a multitude of species, including humans. Lake ecology and responses to anthropogenic pressures have thus been at the forefront of ecological research for decades (e.g. Schindler 1980). More recently, the recognition of inland waters as an active component of the global carbon (C) cycle (Cole et al 2007, Tranvik et al 2009 has placed the dynamics of lake C under the spotlight. This is because lakes receive large amounts of C derived from terrestrial ecosystems that, once in the lake, undergoes a multitude of physical and biogeochemical processes. Some of these processes include the bio-and photodegradation of dissolved organic C (Stubbins et al 2010, Guillemette andDel Giorgio 2011), the emission of CO 2 and CH 4 to the atmosphere (Bastviken et al 2011, Raymond et al 2013), or the burial of particulate C in lake sediments (Ferland et al 2014, Mendonça et al 2017, where it can be preserved from degradation for millennia. There is today no question that lakes play a key role in the processing of C on its journey from land to ocean, and this role cannot be understood without considering the distribution of lakes in the landscape.
Lakes can range in size from very small systems of a few hundreds of m 2 to vast water bodies of more than 100 000 km 2 , as for example the Great Laurentian Lakes or the Caspian Sea. Most of the lakes on Earth are located in the Northern latitudes (Feng et al 2016, Messager et al 2016, particularly across the boreal biome where also the larger stocks of soil organic C are (Köchy et al 2015). Estimates of the total number of lakes vary widely, from 21 to 304 millions (Downing et al 2006, Verpoorter and Tranvik 2014, Messager et al 2016 and strongly depend on the methods used, especially on the resolution and treatment of satellite images. In contrast, there is consensus on the fact that the size-distribution of Earth's lakes is strongly skewed towards small systems: there are just a few large lakes but millions of small lakes (Downing et al 2006, Verpoorter and Tranvik 2014, Cael and Seekell 2016, Messager et al 2016, Steele and Heffernan 2017.
At the catchment and regional scale, the C biogeochemistry of lakes is very sensitive to the regional climate and geology as well as to land cover and anthropogenic pressures (Pace and Cole 2002, Larsen et al 2011, Lapierre et al 2015, Hastie et al 2017, Seekell et al 2018. Yet once these effects are accounted for, the size of lakes can strongly influence the dynamics of lake C. For example, a negative relationship between the size of lakes and the concentrations of CO 2 and CH 4 in surface water has been widely reported in the literature (Kelly et al 2001, Sobek et al 2003, Humborg et al 2010Lapierre and del Giorgio 2012, Rasilo et al 2015, Rocher-Ros et al 2017, a pattern that could be linked to a stronger connectivity with terrestrial ecosystems, a higher influence of sediment activity, and lower gas exchange coefficients in smaller lakes. A lower air/ water gas exchange in smaller lakes has also been linked to a positive relationship between lake area and the 13 C signature in dissolved inorganic C ( 13 C-DIC, Bade et al 2004). Likewise, a negative correlation between dissolved organic carbon (DOC) concentrations and lake size has been observed across boreal lakes of Canada and Finland (Kelly et al 2001, Kortelainen 2009, Einola et al 2011, again highlighting the role of lake size in determining connectivity with land C sources. Beyond the wide variability in size, lakes also vary in their morphometry. Lakes can have very different morphometries in terms of contour (circular vs. elongated), shoreline complexity (smooth vs. fractal) or depth profiles (bowl vs. dish-shaped). The morphometry of lakes clearly influences how lakes interact with the surrounding land and the atmosphere, and also influence some key aspects of their internal functioning. For example, the depth profile of a lake has implications for its stratification dynamics (Imboden and Wüest 1995), the internal movement of water (Håkanson 2005), the proportion of the water column that is exposed to sunlight (Wetzel 2001), and the importance of benthic versus pelagic processes as well as the relative importance of littoral areas. The lake contour and shoreline complexity, on the other hand, have implications on the intensity of contact with the surrounding land and the overall interaction with surface winds. Since all of the lake properties and processes that are influenced by the shape of lakes are directly or indirectly related to lake C sources and sinks, it is reasonable to expect that the shape of lakes may as well influence lake C dynamics beyond the known effects exerted by lake size. Yet, to our knowledge, the link between lake morphometry and C biogeochemistry remains largely unexplored. What complicates matters further is that each of the major C species in lakes (CO 2 , CH 4 , DOC, DIC) likely has different responses to these morphometric variables, if at all, because their sources and transformation pathways differ. Thus, identifying the influence of lake shape on a given C species does not necessarily allow inference on the behavior of the other species. Since the potential influence of lake shape on C dynamics is likely subtle relative to the more notorious drivers of C biogeochemistry discussed above, it is necessary to assess this question over a large number of lakes that span not only a wide range in ambient concentrations of the four C species, but also a wide range in both lake size and morphometry.
In this study, we have assessed the influence of lake morphometry on the concentrations of CO 2 , CH 4 , DOC and DIC in the surface water of 265 boreal lakes of varying size and shape. These lakes span several orders of magnitude in surface area (0.002-4300 km 2 ) and are located in six distinct regions that represent wide environmental gradients across the boreal biome of Québec, Canada (figure 1(a)). We hypothesize that, once lake size is controlled for, the morphometry of the lake influences the concentration of the four C species differently through their effect on the degree of connectivity with terrestrial ecosystems and the degree of influence of sediment activity. For example, we hypothesize that lakes with more fractal shorelines should have higher concentrations of C species that derive from surrounding soils, such as CO 2 and DOC (Aitkenhead-Peterson et al 2003, Maberly et al 2012. Among all C species, we expect CH 4 to show the strongest link to lake morphometry, as it is locally produced and is overwhelmingly driven by in-lake processes such as methanogenic activity in littoral anoxic sediments (Bastviken et al 2004, DelSontro et al 2016. In particular, we expect CH 4 concentrations to respond to the lake depth profile, since dish-like lakes should be more influenced by sediment activity than its bowl-like counterparts. It is important to note that ours is a large-scale comparative study that focuses on the main carbon pools in the surface waters of lakes during the late spring and summer period. It is therefore not meant to recreate an annual cycle of these C pools, but rather to detect a signal of morphometry during this period of the annual cycle.

Study region and field sampling
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 2 and CH 4 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.

Lake carbon species
Sampling and analysis of DIC, DOC, CO 2 and CH 4 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 2 (pCO 2 ) 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 4 (pCH 4 ) 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 size and morphometry
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 Dynamic ratio (Håkanson 1982) is an area (in km 2 ) to mean depth (Zmean; m) ratio of the lake. Low values indicate lakes that are bowlshaped, 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 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.
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 4 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.

Statistical analyses
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 2 , pCH 4 , 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 collinearityinduced 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.
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.

Allometric scaling of lake morphometry
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).

Relationship between C concentrations and lake size
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 2 and pCH 4 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 2 and pCH 4 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 2 and pCH 4 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 2 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.

The influence of lake morphometry on lake C biogeochemistry
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 2 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 2 concentrations in lakes. Lakes with higher shoreline complexity tend to have higher pCO 2 , 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 2 via direct groundwater and runoff inputs. In addition to shoreline complexity, the circularity of the lake seems to also slightly influence lake CO 2 dynamics, such that more circular lakes show higher surface water pCO 2 . 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 2 (Vachon et al 2013). We therefore suggest that the influence of circularity is related to the input of CO 2 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 2 evasion.
Among the four C species, CH 4 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 2 , a potentially larger fraction of the CH 4 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 4 . The elastic net model for Symbols are colored from orange to blue proportionally to the first axis' score of the environmental PCA (see section 2), which is highly correlated to the mean annual temperature and annual net primary productivity in the lake catchment ( figure 1(b)).
pCH 4 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 4 . 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 4 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 4 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 4 , which may reflect a higher degree of connection with terrestrial CH 4 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 4 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 4 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 2 and pCH 4 . 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). The size of the symbol is proportional to the absolute value of the coefficient, while the color indicates whether it is positive (blue) or negative (red). The magnitude of the standardized coefficients reflects the relative predictive importance of each explanatory variable. Lake volume and the two first axes of the environmental PCA ( figure 1(b)) were also included as explanatory variables to account for potential dilution effects as well as climatic and catchment-scale drivers of lake C concentrations (see section 2). Detailed coefficients of all explanatory variables can be found in supplementary table 1. The R 2 of the models are 0.33, 0.50, 0.45 and 0.12 for CO2, CH4, DOC and DIC, respectively.
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 2 , CH 4 , and DOC. Lake DIC, composed of CO 2 , HCO 3 − and CO 3 2− , 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 Figure 5. Explained variability (R 2 ) of three different models for CO2, CH4, DOC and DIC that have (a) climate and land cover (represented by the two first axes of the environmental PCA, figure 1(b)), (b) climate, land cover and lake area, and (c) climate, land cover, lake area and lake morphometry as explanatory variables. 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).

Implications for upscaling and regional patterns
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 2 and pCH 4 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 2 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).
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.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// figshare.com/s/7a00c11757ea28cac31a.