Scaling waterbody carbon dioxide and methane fluxes in the arctic using an integrated terrestrial-aquatic approach

In the Arctic waterbodies are abundant and rapid thaw of permafrost is destabilizing the carbon cycle and changing hydrology. It is particularly important to quantify and accurately scale aquatic carbon emissions in arctic ecosystems. Recently available high-resolution remote sensing datasets capture the physical characteristics of arctic landscapes at unprecedented spatial resolution. We demonstrate how machine learning models can capitalize on these spatial datasets to greatly improve accuracy when scaling waterbody CO2 and CH4 fluxes across the YK Delta of south-west AK. We found that waterbody size and contour were strong predictors for aquatic CO2 emissions, attributing greater than two-thirds of the influence to the scaling model. Small ponds (<0.001 km2) were hotspots of emissions, contributing fluxes several times their relative area, but were less than 5% of the total carbon budget. Small to medium lakes (0.001–0.1 km2) contributed the majority of carbon emissions from waterbodies. Waterbody CH4 emissions were predicted by a combination of wetland landcover and related drivers, as well as watershed hydrology, and waterbody surface reflectance related to chromophoric dissolved organic matter. When compared to our machine learning approach, traditional scaling methods that did not account for relevant landscape characteristics overestimated waterbody CO2 and CH4 emissions by 26%–79% and 8%–53% respectively. This study demonstrates the importance of an integrated terrestrial-aquatic approach to improving estimates and uncertainty when scaling C emissions in the arctic.


Introduction
Regions of permafrost soils, perennially frozen ground, store approximately twice as much carbon as is currently in the entire atmosphere [1,2]. With accelerated warming in high latitudes, permafrost is thawing across the Arctic, leading to increasing emissions of CO 2 and CH 4 [3][4][5][6] as it is decomposed directly to CO 2 and CH 4 or transported through landscapes to inland waterbodies [4,[7][8][9]. Globally, inland waterbodies receive 2-3 Pg-C yr −1 from terrestrial landscapes, most of which is emitted as CO 2 fluxes to the atmosphere (0.5-2.1 Pg-C yr −1 ) [10][11][12][13][14], an amount comparable to the global net terrestrial carbon sink [15][16][17]. Inland waterbodies are also a globally significant source of CH 4 (70-150 Tg-C yr −1 ), accounting for up to half of all CH 4 emissions from natural sources and close to a quarter of global CH 4 emissions [13,14,[18][19][20][21][22][23]. Lateral carbon transport accounts for ∼20% of terrestrial net ecosystem productivity in high latitude ecosystems (compared to 1% for temperate and tropical ecosystems) due to the abundance of lakes and ponds, highlighting the importance of CO 2 and CH 4 fluxes from inland waters [11,[24][25][26]. Despite their importance in global carbon cycling, inland waterbody CO 2 and CH 4 budgets remain uncertain due to the wide range in fluxes reported and the uncertainty in waterbody areal estimates, particularly for small lakes and ponds [12,16,27].
Observations of inland aquatic CO 2 and CH 4 fluxes are highly variable, caused by differing contributions from processes across a hierarchy of scales from watershed transport, carbon cycling within waterbodies, to microcosms of microbial productivity. The drivers of inland aquatic carbon dynamics are often complex and non-linear. Observations of these drivers lack the spatial representation necessary for scaling, e.g. water temperature, dissolved organic carbon concentration and lability [28][29][30][31][32]. Traditional bottom-up scaling estimates of inland waterbody CO 2 and CH 4 fluxes applies an average or median flux to an estimated total water surface area on local, regional, or global scales [10,11,14,18,21,26]. Several scaling studies have shown improved estimates of inland waterbody CH 4 emissions by including lake size and landscape history as categorical drivers [14,19,21,23,33]. Lake productivity (calibrated from a remotely sensed analog) has also been used as a linear predictor of lake dissolved CO 2 concentrations [14,34]. However, a large amount of variation in CO 2 and CH 4 waterbody fluxes remains unexplained in scaling studies, resulting in high uncertainty in regional and global carbon budget estimates [10,11,14,19,21]. Applying an average flux to waterbodies within a region or lake size-class could also create a biased carbon estimate, for example if smaller high-flux lakes are more abundant in the observation dataset than they are in the landscape [22,35,36].
Top-down carbon estimates from inversion studies rarely consider inland waterbodies. Most often, inversion studies mask inland water, functionally attributing a zero flux, or categorize waterbodies as wetlands (but see Tan et al [37]). In lake-rich regions this mis-attribution will cause either over or under estimation of fluxes in wetland or terrestrial environments to compensate, a source of uncertainty in the attributed carbon fluxes that is rarely quantified or discussed. A recent top-down and bottom-up comparison of CO 2 fluxes from the North Slope of Alaska found that inland waters were likely a significant source of CO 2 during the early cold season, and attributing a flux to waterbodies in bottom-up estimates was necessary to match airborne observations [38]. Correctly attributing waterbodies in top-down inversion analyses requires a gridded waterbody carbon flux map for use as a flux prior, which are not often available or produced in bottom-up studies (but see Tan and Zhuang [39]).
Advances in remote sensing and computational abilities have led to steady improvements in inland waterbody areal estimates in recent years [19,40,41], but accurately mapping small lakes and ponds still remains a challenge [13,42]. Delineating open waterbodies from vegetated wetland is particularly important because both are critical ecosystems for carbon emissions but with differing governing processes [22]. Vegetated wetlands are often narrow features bordering waterbodies along shorelines or channel networks that cannot be detected without high resolution imagery (<30 m) [43][44][45][46]. Furthermore, wetlands and lakes are often mapped separately (but see Olefeldt et al [47]), which can lead to double counting. This uncertainty around the waterbody and wetland area could explain why bottom-up CH 4 budgets in the Arctic are twice as high as top-down atmospheric inversion estimates [20,48].
Here we create accurate scaling models of waterbody CO 2 and CH 4 emissions that reflect the underlying processes driving carbon cycling in these ecosystems and reduce uncertainty and bias in carbon budget estimates. We use an integrated terrestrialaquatic approach by combining high-resolution remote-sensing imagery of watershed-level and waterbody drivers from an object-based imagery analysis. We include waterbody size and contour characteristics and the surrounding landcover, hydrology, terrain, and landscape heterogeneity as possible variables affecting carbon fluxes. We train boosted regression tree models, a type of machine learning, to predict waterbody CO 2 and CH 4 diffusive fluxes. While ebullition is a large component of aquatic methane fluxes, we do not include it here as ebullitive fluxes are unlikely to be driven by watershed level processes and are stochastic in nature, and therefore remain a unique challenge to scaling. We demonstrate this scaling technique using the Yukon-Kuskokwim (YK) Delta of Alaska as our study region. The YK Delta is subarctic tundra abundant in lakes and wetlands and underlain by discontinuous permafrost [49]. Atmospheric inversion models of high latitudes have shown the YK Delta to be a regional hotspot of CO 2 and CH 4 emissions [50][51][52][53]. Despite this, the YK Delta has been historically understudied with few ground-based observations of carbon fluxes [54,55]. We leverage recent datasets of high-density observations of CO 2 and CH 4 measurements from waterbodies in the YK Delta to train scaling models and map waterbody CO 2 and CH 4 emissions [56, 57].

CH 4 and CO 2 observations
This study uses a dataset of surface water samples (n = 364) analyzed for dissolved CO 2 (n = 235) and CH 4 (n = 294), collected from various waterbodies in the central-interior of the YK Delta from the first half of July 2016-2019 [56,57]. The majority (>85%) of the samples in these datasets were from Figure 1. Schematic diagram of the geospatial analysis and remote sensing imagery used in this study. The thumbnail imagery examples were created from remote sensing layers used in the study but are not presented at accurate scales. The layers in the grey box were used to create the landcover map via k-means. The raw bands and derived layers (including landcovers) were averaged over each set of basins to create landscape-level drivers. An object-based image analysis of the waterbody product from the landcover map was used to create waterbody size, perimeter, and contour drivers, as well as waterbody reflectance, elevation, and maximum flow accumulation. lakes and waterbodies within wetlands. All waterbodies were sampled at the water surface, either from shore for small waterbodies or from a boat. While all waterbodies were sampled in triplicate for dissolved gases with the average reported, the largest few waterbodies were sampled multiple years and at multiple locations. The average of all observations was used for those waterbodies with multiple samples. Variation in concentrations between waterbodies was far greater than interannual variability or spatial variation within waterbodies. The waterbodies in this region are uniformly shallow (less than 2 m deep) and flat, and consequently well mixed [32]. As a result, we do not consider lake depth as driving variable of CO 2 and CH 4 concentrations. A further description of the site, dataset, sample processing, and size distribution of sampled waterbodies can be found in the supplement (figure S1) and Ludwig et al [32].

Geospatial waterbody and sub-basin analyses
We use remote sensing imagery to: (1) identify waterbodies and quantify waterbody contour characteristics, (2) quantify watershed characteristics that might be related to landscape-level drivers, hydrology, or indirectly affect waterbody CO 2 or CH 4 biogeochemistry, and (3) scale results to map diffusive fluxes of CO 2 and CH 4 (figure 1).

Remote sensing imagery
We use Google Earth Engine to select imagery in the YK Delta from 2016 through 2019 to coincide with the timing of the samples in the dissolved CO 2 and CH 4 observation dataset. We used a composite of cloud-free level-2A Sentinel-2 multispectral images (red, green, blue, near-infrared (NIR), and short-wave infrared (SWIR) bands) from July 2019 (within 1 week of water sample collection timing). Level-2A are surface reflectance Sentinel-2 products, with cloud removal, orthorectification, and sen2cor atmospheric corrections applied using the Sentinel-2 Toolbox. We use Sentinel-1 C-band short-aperture radar (SAR) images, pre-processed using the Sentinel toolbox [58,59]. We created four SAR images; mean composites of July and December, VV and VH backscatter (10 m resolution). We use elevation from the mosaiced ArcticDEM (2 m resolution), after filling and detrending the digital elevation model (DEM) using System for Automated Geoscientific Analyses (SAGA) [60]. From these remote sensing layers, we derive further layers such as slope and flow accumulation (from the DEM before detrending), the normalized difference vegetation index (NDVI), and normalized difference water index [61,62].

Landcover map
We use a 5 by 10 m resolution landcover map created for the region (https://doi.org/10.3334/ORNLDAAC/ 2178) [63]. For detailed methods, see Ludwig et al [32]. We used Google Earth Engine's 'entropy' function to create a spatial texture layer from the landcover map. Higher entropy values occur at borders and transitions between landcover types, and we interpret the average entropy in an area as a metric of landscape heterogeneity. We used Google Earth Engine's 'fastDistanceTransform' function to create a gridded layer of distance to nearest water for the study region, using the landcover category identified as 'surface water' and excluding waterbodies >10 km across. While the surface water landcover was highly accurate (balanced accuracy >0.95; Ludwig et al [32]), the validation procedure was not stratified by waterbody size. We compared the size distribution of waterbodies used here to a higher resolution, independent waterbody product validated by size against ground truth points in an overlapping region of the YK Delta (Mullen et al [64]). The mapped waterbodies were remarkably similar despite differences in seasonality and interannual variation in the underlying remote sensing imagery used in the two maps ( figure  S2). Most notably, the waterbody product used here did not under-sample small ponds despite approaching the limit of detection. Therefore, we chose to include waterbodies <0.001 km 2 though they become pixelated, as omitting them would lead to an underestimate of carbon emissions.

Waterbody object-based image classification and shape analysis
We used an object-based imagery analysis of the surface water classification from the landcover map to identify and then quantify aspects of waterbodies in the region. We used the 'reduceConnectedComponents' algorithm in Google Earth Engine, thresholding to exclude waterbodies larger than 10 km in width (due to lack of in situ data), to calculate the area, perimeter, and ratio of area:perimeter of each uniquely identified waterbody. We specified a maximum neighborhood of 1024 pixels in any dimension using eight-way connectedness with a scale of 10 m. We created gridded layers of the average red, green, blue, NIR, and SWIR reflectance from each waterbody by reducing the Sentinel-2 composite imagery in section 2.2.1 over the waterbody objects. We similarly created an average waterbody elevation layer, using the detrended DEM, which reflects each waterbody's relative landscape position on (higher elevation) or between (lower elevation) peat plateaus. In QGIS (Quantum Geographic Information System, we used the 'catchment area' algorithm to create a flow accumulation layer from the DEM, and then reduced this over the waterbody objects to create a maximum flow accumulation per waterbody gridded product. A higher flow accumulation within a waterbody indicates more pixels 'pouring' into it, which correlates with a larger watershed. All layers were reprojected to 10 m resolution while reduced to averages over waterbody objects.

Sub-basin analyses
We split the study region into non-nested contiguous sub-basins that are distinct hydrologic units. We used the 'channel network' SAGA algorithm in QGIS inputting the filled-DEM and flow accumulation map as the channel initialization to construct these subbasins. We created three sets of sub-basins using different channel thresholds: (1) >1 × 10 5 m 2 created 16 036 small sub-basins with an average area of 0.15 km 2 , (2) > 1 × 10 6 m 2 created 1988 medium sub-basins with an average area of 1.2 km 2 , and (3) > 1 × 10 7 m 2 created 278 large sub-basins with an average area of 8.5 km 2 . These sizes were chosen to span the range of watershed sizes from the waterbodies in the observation dataset [32]. For each set of sub-basins, we masked out waterbodies and then calculated the sub-basin average values for the remote sensing imagery and derived indices described in section 2.2.1, the sub-basin average entropy (which we term 'landscape heterogeneity'), and average distance to nearest waterbody. We also calculated the percent cover of each landcover type within each of the three sets of sub-basins. Finally, each set of these sub-basin metrics were assigned to the waterbodies located within those sub-basins, to create images of pseudo-watershed level drivers. We used this subbasin approach to enable all the potential drivers to be organized as a stack of images, which is efficient to operate over for predicting and scaling. The alternative approach, using actual watersheds delineated for each waterbody such as in Ludwig et al [32], would be computationally prohibitive. We identified 17 071 distinct waterbodies used in this study, and even if watersheds for all 17 071 could be derived efficiently, the nested nature of most of the watersheds would prevent operating over watershed-level variables as images. Not all waterbodies will exhibit similar connectivity to the sub-basins used to aggregate remotesensing imagery. While flow accumulation is related to connectivity, better metrics, particularly of lateral or sub-surface flow, could improve predictability in future scaling models. We further expect connectivity between waterbodies and sub-basins to vary seasonally: for example, different landcover variables might be retained when using late season flux observations and imagery where deeper thaw depths change the hydrologic regime. Future scaling studies using seasonally representative datasets will likely need to take this seasonality into account.

Statistical modeling and scaling
We use boosted regression tree models to predict and scale waterbody CO 2 and CH 4 fluxes. Recent studies have used machine learning to accurately model and predict aquatic carbon cycling [31,32,65]. Machine learning methods are particularly useful when using surface reflectance over inland waters, where atmospheric corrections can sometimes cause artifacts that lead to nonsensical results in approaches that predict using band ratios, band linear combinations, or other empirical or analytical models [66][67][68]. For example, gradient boosting models, such as used here, were effective in predicting Chl-a using Sentinel 2 and 3 imagery of inland waters [69,70]. We use a similar approach to Ludwig et al [32], using the gradient boosting machine ('gbm') package in R v.3.6.1 [71]. Both CO 2 and CH 4 values were log-transformed to achieve normality. Our potential drivers include all waterbody-specific variables (section 2.2.3), sub-basin averages and percent landcover areas (section 2.2.4). We sampled this stack of image layers at every waterbody observation point to create a tabular dataset of drivers for model training. We used 'gbm.step' as described in Elith et al [71] to tune the number of trees and drop variables to avoid overfitting, using ten-fold cross validation, with a bag-fraction of 0.65. We tested the learning rate (lr = 0.005) and tree complexity (tc = 2) manually to optimize. We used percent deviance improvement over null model from cross-fold validation to determine model predictability, and regression between observations and fitted values of dissolved gases to determine model fit ('lm' function in R). We repeated this analysis with ten random seeds for bag-fraction to quantify data uncertainty in model training. We calculated the relative influence of each predictor variable, which were scaled to sum to 100 [71,72]. We used partial dependence plots (the mean and standard deviation of the ten model runs with different random seeds) to investigate the average predicted response across all observations for a given predictor variable [72]. All partial dependence plots were centered on zero µM predicted CO 2 or CH 4 .
We used the 'predict' function in the 'raster' package in R to apply our boosted regression tree scaling models to the study region using the stack of imagery described in sections 2.2.3 and 2.2.4. These predictions were then back-transformed to convert into µM, with a prediction-based backtransformation bias correction applied [73]. Dissolved gas concentrations were converted to diffusive fluxes (mg C m −2 d −1 ) using the relationship in equation (1) [74,75], where C aq is the surface dissolved gas concentration, C atm is the atmospheric concentration, F is the diffusive flux, and k is the gas transfer velocity, Gas transfer velocity values were obtained from paired observations of chamber-based diffusive fluxes and dissolved surface water concentrations from waterbodies in the YK Delta (n = 55 for CO 2 , n = 65 for CH 4 [56,57]). Gas transfer velocities were calculated separately for CO 2 and CH 4 and normalized to k 600 using Schmidt numbers [75]. High outliers of k 600 were deemed to be contaminated by ebullitive fluxes and removed [33]. Since there was no detectable relationship to lake size or wind speed we used the mean, 25th, and 75th percentiles of gas transfer velocity from the observations in the region.

Model performance
Our models were able to accurately fit and predict dissolved CO 2 and CH 4 observations using only remotely sensed drivers. Model fit for dissolved CO 2 and CH 4 was not significantly different from a slope of one based on 95% confidence intervals, with high coefficients of determination (CO 2 ; r 2 = 0.76, RMSE = 175 µM, CH 4 ; r 2 = 0.75, RMSE = 4.8 µM, figure 2). The average of residuals from both the model predictions (figures 2(a) and (b)) and the back-transformed results (figures 2(c) and (d)) were zero (p-value ≫ 0.05 in one-sample t-tests). Both scaling models performed equally well across waterbody sizes, with no relationship between model residuals and waterbody area (figure S3). Models that also consider biogeochemical and mechanistic drivers (e.g. dissolved oxygen, dissolved organic carbon), which cannot be detected by remote sensing, outperform our scaling models (e.g. Ludwig et al [32]: CO 2 ; r 2 = 0.94, CH 4 ; r 2 = 0.88).
We measured the predictive strength of our models using ten-fold cross-validation using percent reduction in deviance as the metric of success. Our CO 2 model had the best predictability with 63% reduction in predictive deviance, while our CH 4 model had 39% reduction in predictive deviance. This is only slightly reduced predictive ability compared to the models used to describe dissolved CO 2 and CH 4 in unburned watersheds in the YK Delta (79% and 52% respectively) that incorporated numerous biogeochemical but non-scalable drivers, and is better predictability than the equivalent models for burned watersheds (61% and 36% respectively) [32].

Figure 2. Modeled and observed dissolved gas concentrations for CH4 (a) and (c) and CO2 (b) and (d). Model goodness-of-fit is shown with log-transformed data (a) and (b) and un-transformed and back-transformed data with bias corrections (c) and (d).
Solid lines and shading depict the fitted slope and 95% CI from linear regressions, with the one-to-one lines as dashed lines.

Spatial drivers of dissolved CO 2 and CH 4
Dissolved CO 2 was primarily driven by variables related to waterbody contour, with waterbody perimeter, area, and the ratio of area:perimeter accounting for more than two-thirds of the total explanatory power of the model ( figure 3(a)). In contrast, the influence of waterbody contour drivers on dissolved CH 4 was similar to watershed landcover and watershed hydrology drivers ( figure 3(b)). For both CO 2 and CH 4 , smaller waterbodies and those with more complex contours (smaller area:perimeter) exhibited higher dissolved gas concentrations (figures 5(a), (b), S4(a) and S6(d)), consistent with patterns observed globally and in the Arctic [13,19,21]. Higher carbon emissions from waterbodies with complex contours could be caused by more relative abundance of lake-edge landcovers: the transitions between wetland and aquatic ecosystems have long been recognized as biogeochemical hotspots [76].
Wetland landcover was the most important driver of dissolved CH 4 ( figure 3(b)), with a threshold effect where more CH 4 was predicted as wetland percent area rose above 15% (figures S4(f) and (j)). Basin averaged NDVI, SWIR reflectance, and SAR backscatter were also significant drivers of dissolved CH 4 in waterbodies and likely also depict the importance of wetlands on downstream CH 4 . Wetland areas were distinctly visible in NDVI, SWIR, and SAR imagery, identifiable as the greenest areas (NDVI effect), while SAR backscatter and SWIR reflectance are commonly used in wetland mapping due to relationships to soil water content and canopy structure [59,[77][78][79][80]. Basin averaged NDVI, SAR, and SWIR correlated with basin wetland percent area in the study region (Pearson's r = 0.81, 0.56, −0.58 respectively). Wetlands are often a significant, if not the dominant, source of natural CH 4 emissions in ecosystems. Our results show wetlands also impact CH 4 concentrations in nearby and downstream waterbodies ( figure 3(b)). This downstream effect could be caused by groundwater flow transporting dissolved CH 4 from where it was produced in a wetland to an open waterbody [81]. Wetland-related drivers were also significant but less important for predicting CO 2 ( figure 3(a)). Wetlands may encourage conditions towards CH 4 and CO 2 production downstream in the watershed through longer water residence times, depletion of oxygen, and increased dissolved organic carbon inputs [82]. Regardless of the mechanism, cohesive high-resolution mapping of waterbodies and wetlands is important for scaling CH 4 emissions.
Basin hydrology plays an important role for both dissolved CO 2 and CH 4 . Predicted CH 4 concentrations peak when surface waters are present in ∼10% of the surrounding basin, but decline if there is either less or more surface water present ( figure S4(c)). Error bars indicate standard deviation of relative influence from ten boosted regression tree model runs with different random seeds. Relative influences are scaled to 100%, with similar variables combined (e.g. watershed SAR backscatter from VV, VH, July, and December composites). The full listing of variable relative influence can be found in tables S1 and S2 in the supporting information.
Predicted CO 2 strictly declined with increasing surface water area in the surrounding basin ( figure  S6(e)). Predicted CH 4 and CO 2 concentrations both increased as a function of increasing basin-average distance to nearest water (figures S4(g) and S6(l)). Flow accumulation had an overall negative effect on CH 4 concentrations ( figure S4(h)). Surface water area, distribution, and flow accumulation relate to water residence times in the surrounding landscape, and could indicate basins with longer water residence times can promote higher waterbody CH 4 , possibly through increased interaction with soil pore water, more oxygen depletion, or more leaching of carbon substrates [83,84].
Waterbody surface reflectance in blue, red, and NIR bands contributed significantly to predicting dissolved CO 2 and CH 4 (figure 3). Combinations of these bands from Landsat surface reflectance have been used to remotely sense chromophoric dissolved organic matter (CDOM) in rivers and inland waters in other environments [85][86][87]. Previous results from the YK Delta have indicated that CDOM is an important driver of both dissolved CO 2 and CH 4 [32]. Small waterbodies (<0.01 km 2 ) may have land-adjacency effects in surface reflectance that are complicating signals from waterbody color. It is possible for these smaller waterbodies that the role of surface reflectance instead indicates increasing edge-effects through land-adjacency, which is a demonstrated driver of both CO 2 and CH 4 . To test how small-waterbody reflectance values, land-adjacency, and size metrics affect our scaling models, we re-ran them using the same hyper-parameters, drivers, and random-seeds while replacing waterbody reflectance, size, contour, and perimeter variables with NA's for those waterbodies under 0.01 km 2 in area. There was a slight decrease in predictive performance (percent deviance explained decreased by 3% and 2% for CH 4 and CO 2 ). We regressed the predicted results from our original models against those with small lake data withheld, with little discernable difference (slopes of 1, intercepts of 0, and r 2 of 0.97 and 0.99). The remote sensing-based models in this study could be capturing both CDOM and edge-effects (for the smallest waterbodies) indirectly through waterbody surface reflectance.

Lake size effects
Using the boosted regression tree models developed from observations of dissolved CO 2 and CH 4 and remotely sensed drivers (figure 1), we created maps of waterbody CO 2 and CH 4 fluxes for our study region in the YK Delta ( figure 4). The effect of waterbody size on CH 4 and especially CO 2 fluxes can be seen clearly in the mapped predictions (figures 5(a) and (b)). The smallest waterbodies account for the highest fluxes of both CO 2 and CH 4 , and fluxes decline precipitously with increasing waterbody size for CO 2 . The smallest waterbodies (<0.001 km 2 ) have disproportionately high fluxes and comprise 38% of the 17 071 uniquely identified waterbodies in the region, but account for only 1% of the overall surface area of water (table 1). Because of this, the smallest waterbodies do not contribute significantly to the total CO 2 and CH 4 fluxes (figures 5(c) and (d)), accounting for only 3% and 2.8% of the total fluxes of CO 2 and CH 4 in the region respectively (table 1). In contrast, Holgerson and Raymond [13] estimated very small ponds (<0.001 km 2 ) account for ∼15.1% and 40.6% of global CO 2 and CH 4 emissions from lentic inland waters. Holgerson and Raymond [13] estimated the distribution of very small ponds by extrapolating the Pareto distribution, a log-abundance logsize regression relationship [88], yielding a much greater area of very small ponds relative to total water area. While our waterbody mapping could be missing sub-pixel very small ponds (<0.00005 km 2 ), our lake distribution are similar to higher resolution products from the same region [42,64]. Muster et al [42] found that waterbody distributions in permafrost lowland regions were not representable by a power-law relationship. Consequentially, the Pareto distribution would overestimate very small pond abundance. Small lakes (0.001-0.1 km 2 ) are the most abundant in the study area, have high CO 2 and CH 4 emissions relative to their area, and contribute the majority of total CO 2 and CH 4 emissions (table 1). Medium and large lakes are the least abundant, and while they are still significant contributors to CH 4 emissions, they have a small or negative contribution to CO 2 emissions (table 1). Our results suggest that very small ponds in the Arctic need to be explicitly mapped using high-resolution techniques to avoid underestimating aquatic C emissions from their exclusion, or overestimating aquatic C emissions from inaccurate areal estimates.

Scaling CO 2 and CH 4
Carbon emissions from inland waterbodies remains one of the least certain portions of the global C cycle. We compared the total fluxes when scaled using our boosted regression models to two simpler, more traditional approaches. The simplest scaling method multiplies the average areal flux rate by the total waterbody surface area [10,11,18].
The other approach applies the average observed flux in a waterbody size-class to the total area of water in that size-class [13,14,21,26,33]. With the simplest approach, CO 2 total fluxes are overestimated by 79% and CH 4 total fluxes are overestimated by 53% ( figure 6) compared to our approach. When scaling using the average flux by size-class, CO 2 and CH 4 fluxes were overestimated in the three smallest size-classes of waterbodies ( figure 6). The largest lakes (>1 km 2 ) were underestimated for CO 2 and CH 4 in the size-class average scaling, as well as the mid-sized lakes (0.1-1 km 2 ) for CO 2 ( figure 6). Overall the scaling approach using the average flux by lake size-class overestimated the diffusive waterbody fluxes of CO 2 and CH 4 from the region by 26% and 8% respectively. This bias in simpler scaling results is likely due to a lack of spatial representativeness in the observation dataset, a common problem in arctic flux datasets. Synthesis studies in particular will lack spatially representative flux observations. For example, Kuhn et al [21] attribute the relatively poor fit of their aquatic CH 4 model compared to their terrestrial models' performance to spatial under-sampling.
Recent efforts to improve arctic CH 4 estimates include the Boreal-Arctic Wetland and Lake methane Dataset (BAWLD), a circumpolar database of wetland landcover jointly mapped with lake size and abundance [22,47]. The BAWLD database solves issues of double-counting and the representativeness of empirical flux observations and is well disposed for large-scale modeling of circumpolar CH 4 emissions. However, the random forest models used in BAWLD has the lowest predictive power when mapping small lakes (<0.1 km 2 ), and did not differentiate waterbody sizes below 0.1 km 2 . Our results are complimentary to the BAWLD approach and demonstrate how highresolution remote sensing, machine learning models, b Figure 5. Box plots of predicted fluxes of CO2 (a) and CH4 (b) as a function of waterbody size. The lower and upper hinges correspond to the first and third quartiles, the whiskers extend to 1.5 times the interquartile range, with outliers indicated as points. Cumulative CO2 (c) and CH4 (d) fluxes as a function of waterbody size, with cumulative waterbody area as a secondary y-axis. Note the log10 scale x-axis. Shading indicates the predicted flux using the 75th and 25th percentile of k600 from observations in the region. Table 1. Waterbody distribution and area by size-class. CO2 and CH4 emissions are the predicted dissolved concentrations from our scaling models, converted to diffusive fluxes using the mean gas transfer velocity (k in equation (1)) from observations, multiplied by each waterbody's surface area, and summed over all the waterbodies within each size-class. The range in parentheses for CO2 and CH4 emissions is calculated using the 25th and 75th percentile of gas transfer velocities to calculate fluxes. and extensive field observations can be leveraged to explicitly and accurately map open-water diffusive C fluxes at the regional scale. Our models do not account for temporal variation in fluxes, but rather provide a snapshot map of peak growing season fluxes. Typically, waterbody scaling studies extrapolate an average flux to a seasonal or annual budget estimate. Recent syntheses have improved carbon budget estimates by accounting for ice-free days and large fluxes from ice-off events [19]. Studies of specific lakes have documented diurnal patterns in dissolved CO 2 and CH 4 in surface waters, though sampling regimes rarely capture this temporal variation [89]. Similarly, variable gas evasion from wind and convection affecting turbulence is an important source of temporal variability in lake fluxes [89][90][91]. While our models improve the spatial accuracy of waterbody fluxes, we recommend considering seasonal and diurnal variation in dissolved gas concentrations and gas transfer velocities Figure 6. Total carbon emissions from each waterbody size-class using different scaling methods: this study, using an integrated terrestrial-aquatic approach (green), using the average flux by size-class (purple), and the overall average flux (orange). Scaled CO2 total fluxes are in the top panel (a) and scaled CH4 total fluxes are in the bottom panel (b). Carbon fluxes from the smaller water bodies are overestimated and largest lakes are underestimated by the simpler scaling methods. Error bars depict the range flux totals derived using the 25th-75th percentile of k600 from observations in the region. and seasonal variability in water body size when using our flux maps.

Conclusion and implications
Inland aquatic carbon emissions remain one of the most uncertain components of the global carbon budget. We can reduce this uncertainty using scaling models driven by watershed and waterbody processes. Aquatic CO 2 emissions can be predicted well using lake size and contour as continuous variables, which could be applied to larger regions with suitably high-resolution waterbody maps. Very small ponds in the Arctic need to be explicitly mapped using high-resolution techniques to avoid biasing aquatic carbon emissions from their exclusion or inaccurately extrapolated areal estimates. Waterbody size metrics are insufficient for scaling CH 4 emissions, which were primarily driven by wetland landcover and wetland-related variables within watersheds. Our terrestrial-aquatic integrated approach using watershed landcovers and hydrology-related remote sensing drivers improved our ability to scale both CO 2 and CH 4 waterbody emissions. Concurrent wetland landcover and waterbody mapping is necessary to avoid double-counting open water areas and for integrating terrestrial effects on aquatic emissions. As the Arctic warms with climate change, new waterbodies will form from thawing permafrost while others will drain and be replaced by wetlands. Our results imply the increased abundance of small ponds and replacement of large lakes with wetlands would lead to higher emissions of CO 2 and CH 4 for the YK Delta.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.18739/A23775V7T and https://doi.org/10. 18739/A22804Z8M. Maps of CO 2 and CH 4 produced in this study can be found at https://doi.org/10.3334/ ORNLDAAC/2178.