Modeling Terrestrial Dissolved Organic Carbon Loading to Western Arctic Rivers.

The mobilization and land-to-ocean transfer of dissolved organic carbon (DOC) in Arctic watersheds is intricately linked with the region's climate and water cycle, and furthermore at risk of changes from climate warming and associated impacts. This study quantifies model-simulated estimates of runoff, surface and active layer leachate DOC concentrations and loadings to western Arctic rivers, specifically for basins that drain into coastal waters between and including the Yukon and Mackenzie Rivers. Model validation leverages data from other field measurements, synthesis studies, and modeling efforts. The simulations effectively quantify DOC leaching in surface and subsurface runoff and broadly capture the seasonal cycle in DOC concentration and mass loadings reported from other studies that use river-based measurements. A marked east-west gradient in simulated spring and summer DOC concentrations of 24 drainage basins on the North Slope of Alaska is captured by the modeling, consistent with independent data derived from river sampling. Simulated loadings for the Mackenzie and Yukon show reasonable agreement with estimates of DOC export for annual totals and four of the six seasonal comparisons. Nearly equivalent loading occurs to rivers which drain north to the Beaufort Sea and west to the Bering and Chukchi Seas. The modeling framework provides a basis for understanding carbon export to coastal waters and for assessing impacts of hydrological cycle intensification and permafrost thaw with ongoing warming in the Arctic.

Introduction

Arctic rivers export disproportionately large amounts of freshwater to coastal waters, with marked seasonal variations serving as a defining feature of each river's discharge regime. Rivers draining into Arctic coastal waters are rich in organic matter (McClelland et al., 2012), with Arctic watersheds contributing 32 Tg of dissolved organic carbon (DOC) to fluvial networks annually (Kicklighter et al., 2013). Across Alaska, riverine carbon export to the coast composes a large fraction of total (lateral and vertical) aquatic fluxes (Stackpoole et al., 2017). Rivers serve as a conduit for transferring terrestrial materials to the ocean, and in doing so they integrate hydrologic and biogeochemical processes across the contributing watershed that control how and in what form material and nutrients arrive to coastal waters (Abbott et al., 2016; Tank et al., 2016). Terrestrially derived DOC and other nutrients are a dominant component of river export to the Arctic Ocean, which strongly influence coastal food web dynamics (Connelly et al., 2015; Kicklighter et al., 2013; McMeans et al., 2013). Thawing permafrost that contains ancient soils is particularly concerning because mobilized organic matter may be readily decomposed, thereby causing a positive feedback to climate change (Gao et al., 2018). Ancient DOC delivered to fluvial networks under projected thawing conditions (Mann et al., 2015; Spencer et al., 2015) may alter aquatic food web structure and function in the decades to come (O'Donnell et al., 2021).

Manifestation of climate warming, including hydrological cycle intensification (Rawlins et al., 2010) and permafrost thaw (Biskaborn et al., 2019; Schuur et al., 2009; Turetsky et al., 2020) have the potential to alter the storage and movement of water‐borne constituents in northern high latitudes watersheds. A significant increase in warm season rainfall over recent years across Arctic watersheds (Arp et al., 2020) is altering hydrological connectivity and influencing terrestrial fluxes of organic and inorganic materials (Beel et al., 2021). Model projections point to future increases in convective storms and rainfall (Bintanja & Andry, 2017; Bintanja, 2018; Poujol et al., 2020, 2021). Upward trends in river discharge (Arp et al., 2020; Peterson et al., 2002; Rood et al., 2017; Shiklomanov & Lammers, 2009) in particular will directly impact organic matter fluxes to Arctic coastal environments, as water yield is a major control on fluxes of DOC in Arctic rivers (Holmes et al., 2012; Raymond et al., 2007). Studies have pointed to sustained long‐term increases in DOC delivery to the ocean (McGuire et al., 2010; Tank et al., 2016). Modeling results suggest that changes in DOC loading over the twentieth century are primarily a result of climate‐induced increases in runoff (McGuire et al., 2010).

Permafrost strongly influences soil moisture by controlling vegetation and water drainage, especially in near‐surface soils close to the land‐atmosphere interface (Streletskiy et al., 2017; Walvoord et al., 2012; Woo et al., 2008). While permafrost presents a barrier that generally limits deep surface‐groundwater transfer, shallow groundwater in permafrost regions flows through seasonally thawed (active layer) soils. This supra‐permafrost groundwater contributes DOC and nutrients to fluvial networks and Arctic coastal waters. Contributions from supra‐permafrost groundwater are particularly important during late summer (Connolly et al., 2020; Neilson et al., 2018). Recent studies have connected a warming climate and permafrost thaw to changes in the biogeochemistry and export characteristics of Arctic rivers (Coch et al., 2020; Frey et al., 2007; Frey & McClelland, 2009; Raymond et al., 2007; Striegl et al., 2005). Thus, observable changes in river‐borne constituent concentrations and fluxes are useful markers of environmental change (McClelland et al., 2008; Striegl et al., 2005; Walvoord & Striegl, 2007). Future change to watershed DOC concentrations and loading to Arctic rivers may be substantial. Simulations from regional and global‐scale biogeochemical models for the period 2010–2299 indicate losses of permafrost of between 3 and 5 million for climate of RCP4.5 and between 6 and 16 million for RCP8.5 (McGuire et al., 2018).

The largest fraction of annual DOC export from Arctic rivers occurs during the spring freshet, an approximate 2–4 week period typically starting in May when snowmelt‐driven runoff leaches decaying tundra vegetation and organic‐rich shallow soils above frozen ground, thereby mobilizing massive quantities of DOC that is shunted through Arctic fluvial networks (Holmes et al., 2012; McClelland et al., 2012; Spencer et al., 2008). Measurements from the outlet of Arctic rivers also show that concentrations of DOC are highest during the spring freshet period (Connolly et al., 2018; Holmes et al., 2012; McClelland et al., 2012). Studies across a wide range of spatial scales, including the North Slope of Alaska (McClelland et al., 2007, 2014) and the broader pan‐Arctic drainage region (Holmes et al., 2012; Kicklighter et al., 2013; McClelland et al., 2012, 2016; Raymond et al., 2007) have demonstrated the prominent role that the spring freshet plays in the export of water and water‐borne constituents from rivers to the Arctic Ocean. Less attention has been given to the role of seasonal changes in sub‐ and supra‐permafrost groundwater flow and soil organic carbon interaction (Connolly et al., 2020; Walvoord & Kurylyk, 2016).

Warming and associated climate‐driven changes to Arctic watershed loadings and river export motivate the need for modeling studies that advance understanding of processes influencing the flows of water and carbon across the critically important regions. Vertical fluxes (i.e., land‐atmosphere) of carbon are much larger compared to lateral fluxes, and have received much attention given the considerable influence on atmospheric . Assessments of ecosystem carbon stocks and fluxes require accurate estimates of the lateral transport of DOC to the coast. Few land surface models account for lateral carbon fluxes to aquatic ecosystems (McGuire et al., 2016), although approaches to model carbon mobilization, leaching, and loading to Arctic rivers have been developed and incorporated into ecosystem and land surface models (Kicklighter et al., 2013; McGuire et al., 2010; Nakhavali et al., 2018; Wen et al., 2020). Simulations of river DOC dynamics can be derived by coupling of a DOC production module with a hydrological model (Lessels et al., 2015). Model schemes for DOC production and loading typically employ a soil organic carbon pool and process representations that account for the influence of soil temperature and moisture on organic matter decomposition and leaching.

This study describes the implementation and validation of a model that estimates DOC loading to watersheds spanning the western Arctic. Descriptions of the model algorithms and validation, which includes relevant processes and temporal variation, is followed by discussion of seasonal and annual total loadings. Finally, we describe uncertainties in the model estimates and highlight the study's conclusions.

Methodology

Study Domain

The spatial domain includes all land areas that drain river basins with outlets between (and including) the Yukon and Mackenzie Rivers (Figure 1, hereafter the “Y2M region”). We make use of validation data, namely observations of river discharge and associated DOC export, for the large Mackenzie and Yukon Rivers, and two on Alaska's North Slope—the Colville and Kuparuk Rivers. Besides northern AK, the Mackenzie and Yukon basins, the remainder of the Y2M region is occupied by what we term “western Alaska” for partitioning the total DOC loading, as described in Section 3.3. Relatively small rivers span this western Alaska region, and no observations of DOC export exist with which to validate the model simulation. To better understand spatial variability in DOC concentration of river waters draining the North Slope, its domain is resolved into 24 basins. The full Y2M region spans roughly 3,018,675 and is resolved into 4804 25 × 25 km grid cells on the Northern Hemisphere EASE‐Grid (Brodzik & Knowles, 2002).

Figure 1

Y2M study domain encompassing drainage basins Mackenzie (black), Yukon (red), and regions western Alaska (green) and northern Alaska (blue). Two river basins on the North Slope, the Colville and the Kuparuk, shown in light blue. The Colville, Kuparuk, Mackenzie, and Yukon comprise the four focus basins.

Numerical Model

Numerical modeling and analysis of measured data for river discharge and DOC concentration, loading, and export were used to quantify daily DOC loading from each grid. Simulations of DOC production and loading to rivers across the Y2M region were made with an updated version of the Permafrost Water Balance Model (PWBM), achieved by adding new algorithms associated with evapotranspiration (ET) and snow sublimation, and a new DOC module, to the PWBM v3 as described in Rawlins et al. (2013, 2019). An overview of relevant model components and new updates to the model, hereinafter PWBM v4 (Figure 2), are described in the sections that follow.

Figure 2

(a) Diagram of the vertical profiles for the Permafrost Water Balance Model (PWBM) showing characteristic water storage and flux elements for winter, late spring when active layer is developing, and a thawed profile consistent with summer in an area absent of permafrost. (b) Flowchart showing model parameterizations, forcings, outputs, and submodels used to estimate DOC loading and export to coastal zones. River routing of DOC mass loading is not employed in the present study.

Water Balance Modeling and Updates

The PWBM v4 simulates all major elements of the water cycle, including snow storage, sublimation, transpiration, and surface evaporation (Rawlins et al., 2003, 2013). The model includes a multi‐layer soil module with representation of unfrozen water dynamics and phase change, as well as specification of the thermal and hydrological properties of organic soils (Nicolsky et al., 2007). The snow model simulates the effects of seasonal changes in snow density and, in turn, snow thermal conductivity (Liston et al., 2007; Sturm et al., 1995). Input and output fields are resolved at a daily time step. Simulations are commonly performed at an implicit daily time step and forced with meteorological data. Daily air temperature, precipitation, and wind speed data from the Modern‐Era Retrospective Analysis for Research and Applications (MERRA; Rienecker et al., 2011) were used as forcing. The MERRA data were recently used as forcing for simulations that helped quantify baseline conditions and investigate the changing characteristics of hydrological elements across Alaska's North Slope (Rawlins et al., 2019). Daily runoff across a model grid cell is the sum of surface and subsurface runoff from one or more soil layers. A routing procedure can be applied to propagate runoff through a digital representation of river networks to estimate discharge at the coast (Rawlins et al., 2019). The distribution of tundra and forests across the study domain are parameterized with vegetation classes constructed from a number of sources for simulations with the Terrestrial Ecosystem Model (TEM) at 0. resolution (Melillo et al., 1993). Each grid cell in the simulation described here is assigned a single vegetation class.

In Arctic rivers, terrestrial hydrology, in particular, runoff, is a dominant control on DOC loading to fluvial networks (Holmes et al., 2012; Kicklighter et al., 2013; Raymond et al., 2007). In this study, process representations in the PWBM involving ET and snow sublimation were modified to improve the simulated runoff estimates. Potential ET is estimated with the Hamon method (Hamon, 1963), an empirical function based on air temperature. While the Hamon method is a relatively simple empirical relation, a study comparing 11 potential ET methods suggested that this approach exhibited an appropriate empirical response to the interaction of vegetation type and climate, and provided satisfactory results for studies of contemporary climate (Vörösmarty et al., 1998). In the PWBM, the soil column is discretized into 23 layers, with ET calculated as the sum total of water extracted from each layer within the rooting zone (Figure 2). Daily simulated ET (mm ) for each grid cell in the model domain depends on atmospheric demand, and surface and soil conditions

where is unadjusted ET (mm ), is potential ET (mm ), is the proportion of roots in the soil layer, and is the proportion of water above field capacity relative to the available storage in the soil layer. Here a calibration coefficient () was used to account for spatial variations in canopy conductance and to scale the daily grid ET, with the assumption that canopy conductance is equal to surface conductance:

Daily snow sublimation (mm ) is computed as

where is the model simulated value for snow water equivalent (mm) and is potential sublimation (mm ) that takes the form

where is potential snow‐atmosphere flux from the Hamon method (mm ), is monthly average leaf area index (), is the calibration coefficient (dimensionless) for tundra, and (dimensionless) is the calibration coefficient for forest environments. Canopy interception and sublimation of snow is an important process in forest environments (Hedstrom & Pomeroy, 1998; Pomeroy & Schmidt, 1993; Storck et al., 2002), and transpiration is closely related to Leaf Area Index (H. Yan et al., 2012; Zhang et al., 2008). These differences motivated our inclusion of different calibration coefficients for forest and tundra grids (Table 1). The variation imposed by changing LAI in forest areas when LAI 1.0 accounts for the influence of canopy interception on sublimation, with higher LAI resulting in less snow mass at the surface. Monthly long‐term (1981–2015) mean LAI were drawn from data derived from the Advanced Very High Resolution Radiometer (AVHRR) Global Inventory Modeling and Mapping Studies (GIMMS) LAI3g version 2 (Mao & Yan, 2019).

Table 1

Major Parameters Related to Hydrological Flows in the Model Simulations

Hydrology parameters
Parameter	Value	Description
k_ET tundra	0.38,0.45	Evapotranspiration coefficient, tundra grid
k_ET forest	1.48,2.00	Evapotranspiration coefficient, forest grid
k_subl tundra	0.06,0.10	Sublimation coefficient, tundra grid
k_subl forest	0.90,2.36	Sublimation coefficient, forest grid

Note. Values for and are listed for the Y2M region excluding the Mackenzie River basin, with values for the Mackenzie following the comma. The parameters are dimensionless scalars.

Process Modeling of DOC Loading

Watershed characteristics exert controls on the mobilization and export of riverine DOC (Connolly et al., 2018; Larouche et al., 2015; Neilson et al., 2018). Simulations with the PWBM provide the framework for modeling DOC production and loading to rivers. Production of DOC is assumed to occur through incomplete decomposition of soil organic matter (SOM), with the rate of decomposition, and thus DOC production, influenced by soil temperature, soil moisture, and the amount of SOM (Kicklighter et al., 2013; McGuire et al., 2010). In our simulations, DOC production occurs as follows

where (t) is the rate of DOC production (g C ) on day , represents the production rate coefficient (), SOM is the density of soil organic matter (g C , and and are the rate dependence on soil temperature () and moisture (), respectively. Table 2 lists the adjustable parameters for DOC production and decomposition. The parameters were set separately for the Mackenzie, Yukon, and remainder of study domain to help to account for the differing DOC yields expressed across the large region. The influence of soil temperature on DOC production in each soil layer is modeled through the commonly used relationship: , where estimates the rate increases with soil temperature. The influence of soil moisture is modeled with , where takes the value of 1.0, which is within the typical range of 0.75–3.0 for most soils (Z. Yan et al., 2018). For each model soil layer, SOM densities were assigned from the 0–100, 100–200, or 200–300 cm layer value in the Northern Circumpolar Soil Carbon Database (NCSCD v2.2) (Hugelius et al., 2013).

Table 2

Parameters for DOC Production and Loading From Surface and Within Soils

DOC parameters
Parameter	Value (day^−1)	Description
k_prod	0.11 × 10^−6, 0.17 × 10^−6, 0.43 × 10^−7	DOC production rate coefficient, subsurface
k_prod	0.44 × 10^−11, 0.15 × 10^−10, 0.59 × 10^−11	DOC production rate coefficient, surface
k_decomp	0.60 × 10^−2, 0.83 × 10^−2, 0.50 × 10^−1	DOC decomposition rate coefficient, subsurface
k_decomp	0.81 × 10^−1, 0.97 × 10^−1, 0.92 × 10^−1	DOC decomposition rate coefficient, surface

Note. The three values shown are applied for each grid in Mackenzie, Yukon, and remainder of study domain. The surface are scaled to account for enhanced production and leaching during and after snowmelt as described in Section 2 and Equation 8.

The decomposition of DOC is assumed to be lost to carbon dioxide and/or sorbed to the mineral soil as expressed by

where is the decomposition rate coefficient () and () is soil water DOC storage. Estimates of daily DOC production (Equation 5) and loss to decomposition (Equation 6) are used to update the soil DOC pool

where is soil DOC storage on day , before leaching (Equation 9) and is soil DOC storage before production and decomposition occurs.

In north‐draining Arctic rivers, peak DOC concentration at the coastal outlet occurs during the spring freshet (Connolly et al., 2018; Holmes et al., 2012; McClelland et al., 2014), when high flow produces a large fraction of the total annual DOC export. During this time, melt water accumulating at the base of warming snowpacks in late spring (late‐April to early‐May) interacts with decaying plant litter and shallow organic‐rich soil layers, including biomass produced during the prior growing season (Guo & Macdonald, 2006; Neff et al., 2006; Spencer et al., 2008; O'Donnell et al., 2021). To account for the high DOC concentration during peak flows, the SOM value for surface leaching was scaled from the NCSCD 0–100 cm data value when there is water in the snowpack or meltwater at the surface. The was also modified on days with snowmelt water. The scaling factor () was set under the following conditions when snowpack liquid water is present, and for the days when snowmelt water is present:

The modifications to surface SOM and production reflect the relatively high DOC mass flux that occurs during the spring freshet. No vertical movement of soil water DOC takes place in the model simulation described here. Transfer of DOC from a soil layer to the stream and river network takes place whenever surface or subsurface runoff occurs where DOC(t) is mass load (g C ) and Q is runoff (m ). Leaching of DOC occurs when runoff is generated from the surface water pool or any soil layer, and the term “leachate” refers to the DOC‐laden water loaded to stream and river networks. The soil DOC pool is updated to reflect the loss to the stream and river network

Because loading is the product of runoff volume and DOC concentration, no DOC leaves the soil pool on days with no runoff. Since the model explicitly simulates soil temperature, the direct effects of seasonal active layer development on hydrological flow paths and DOC production and mobilization are quantified. Soil leaching from the surface carbon pool is dependent on the surface runoff while subsurface leaching is a function of subsurface runoff. As with the hydrological variables, simulated concentrations or loadings for specific basins were computed as averages and/or totals over the contributing drainage area. The (Equation 5) and (Equation 6) were manually adjusted to minimize MBE relative to data from published studies (Section 2.3).

Given uncertainties with the level of in‐stream processing of DOC, model estimates of runoff and DOC loading and concentration in this study represent an average or total integrated across the respective river basin. A 50‐year spin‐up over the year 1980 was performed to stabilize soil temperature, moisture, and soil DOC pools, followed by the transient run. Additional background information on model configurations are described in Rawlins et al. (2003, 2013, 2019); Yi et al. (2015) and appendices within.

Model Calibration

Given differences in surface properties, parameters and were determined separately for grids with tundra and forest vegetation classes. Initial analysis revealed large differences in runoff‐precipitation ratio across the Mackenzie River basin relative to several other study basins. Therefore and were also selected separately for the Mackenzie basin and for the remainder of the domain. During calibration, the and values (Table 1) were manually adjusted across a range of values to minimize the mean bias error (MBE, Section 2.4) in long‐term average (1981–2010) annual runoff for four river basins (the Colville, Kuparuk, Mackenzie, and Yukon) with the most complete discharge records in the Y2M region. While the model effectively captures the seasonal cycles of sublimation and ET, calibration and validation targeted annual runoff given its prominent role in DOC loading (Kicklighter et al., 2013; McGuire et al., 2010; Tank et al., 2016). The and parameters were stratified for the Mackenzie, Yukon, and remainder of the domain (Table 2) to better account for differences in DOC yield and concentration based on data in published studies.

Historical river discharge data was retrieved from U.S. Geological Survey (USGS) data repositories for the Colville River at Umiat

(https://waterdata.usgs.gov/ak/nwis/uv/?site_no = 15,875,000), Kuparuk River near Deadhorse

(http://waterdata.usgs.gov/nwis/uv?15,896,000) and Yukon River at Pilot Station

(https://waterdata.usgs.gov/ak/nwis/uv?site_no = 15,565,447). Discharge data for the Mackenzie River at Tsigachik were obtained from the Water Survey of Canada

(https://wateroffice.ec.gc.ca/).

Model Performance Assessment

Simulated runoff and DOC loading and concentration are evaluated against other measured and modeled data (Connolly et al., 2018; Holmes et al., 2012; Kicklighter et al., 2013; Raymond et al., 2007; Spencer et al., 2009) using common measures of average model performance, including mean‐bias error (MBE), mean‐percentage error (MPE), and explained variance derived from Pearson's correlation coefficient (Pearson R). Data for the Yukon and Mackenzie Rivers are archived in the Arctic Great Rivers Observatory (https://arcticgreatrivers.org/data/). Model estimates of DOC concentration of runoff leachate for 24 North Slope of Alaska basins (which can contain multiple smaller rivers) draining into coastal waters are validated through comparisons to independent estimates derived from river water samples and statistical relationships with watershed characteristics (Connolly et al., 2018).

The MBE and MPE are computed as where is the number of model estimates and observations, is the model estimate and are the thought‐to‐be reliable and pairwise matched observations. These metrics are more natural measures of average error, and unlike RMSE, are unambiguous (Willmott & Matsuura, 2005; Willmott et al., 2015). Model validation centers on four focus river basins, the Kuparuk, Colville, Yukon, and Mackenzie, which have drainage areas ranging from 8000 to 1,800,000 (Table 3). We employ these common metrics of model performance while noting that the observations () are not identical quantities to the model estimates. In particular, it is important to keep in mind that we are comparing simulated DOC loadings of river networks to in‐stream estimates of concentrations and fluxes that have been modified to some degree by in‐stream processes. While our river network loading estimates do not account for in‐stream processing, this subject is addressed in the discussion section. Nonetheless the differences from the various independent data sets are reported here with the common error measures shown in Equations 11 and 12.

Table 3

River Basin Area, Annual Precipitation (MERRA), Observed Runoff, Simulated Runoff, and Mean Percentage Error (MPE) for the Colville, Kuparuk, Mackenzie, and Yukon Rivers

Basin	Area (km^2)	Precipitation (mm yr^−1)	Basin average runoff
			Observed runoff (mm yr^−1)	Simulated runoff (mm yr^−1)	MPE (%)
Colville	64,095	399	240	232	−3.3
Kuparuk	8169	319	144	126	−11.3
Mackenzie	1,783,972	473	170	162	−4.7
Yukon	854,700	509	241	250	3.7

Note. Observed runoff in unit depth (mm ) was derived by distributing the river discharge flux across the respective basin area.

Model Validation and Results

Runoff

Simulated runoff is compared with runoff calculated from historical discharge data for four river basins in the region with the longest records; the Colville and Kuparuk in northern Alaska and the large Mackenzie and Yukon rivers. Monthly climatologies were examined given the overarching goal of estimating and understanding patterns in annual and seasonal DOC loading. The magnitude of annual total runoff and seasonal changes from high freshet runoff through recession and winter low flows are generally well captured (Figure 3). Mean percent error (MPE) is 3.3, 11.3, 4.7, and 3.7% for the Colville, Kuparuk, Mackenzie, and Yukon, respectively (Table 3), suggesting relatively well constrained annual total freshwater export. Peak monthly runoff (June) aligns with peak observed runoff for the Colville, Mackenzie, and Yukon. The model performs well in capturing high freshet runoff for the Colville basin (defined by the gauge site at Umiat, Alaska), attributable to higher snowfall rates across the Brooks Range. For the Kuparuk the model peak runoff leads peak discharge by 1–2 weeks, resulting in an overestimation for May and underestimation for June, with freshet (May–June) simulated runoff amounting to 71% of the observed total. The simulations capture the higher ET rates and relatively lower runoff‐precipitation ratios across the Mackenzie River basin, attributable to more extensive forests (with deeper rooting depths) and warmer temperature in the southern parts of the basin. The basin also has less permafrost coverage, which leads to less runoff per unit of precipitation. Model performance demonstrated here is similar to recent results for the North Slope of Alaska from a simulation forced with the MERRA reanalysis data (Rawlins et al., 2019).

Figure 3

Monthly precipitation (mm ), runoff (mm ), and air temperature (c) as climatological averages over the period 1981–2010 from the baseline simulation for the (a) Colville, (b) Kuparuk, (c) Mackenzie, and (d) Yukon River basins. Runoff is expressed in unit depth equivalent. Observed runoff in unit depth was obtained by distributing the measured river discharge flux volumes across the contributing basin area.

Analysis of Relevant Processes and Temporal Variations in DOC Concentrations and Loadings

Model simulated estimates of leachate DOC concentrations and loading were quantified to understand model behavior and agreement with the various independent validation data. Simulated DOC concentrations for the Kuparuk River basin vary over two orders of magnitude for surface runoff and nearly four orders of magnitude for subsurface runoff (Figure 4). These variations are, in part, related to spatial variations in SOC content: Across the Kuparuk basin, 0–100 cm SOC density in the NCSCD database varies from 3 kg C in the foothills of the Brooks Range to 45 kg C near the coast. DOC concentrations also vary as a function of water discharge rate and season. Simulated DOC leachate concentrations in the Kuparuk are highest during spring in association with surface snowmelt runoff. Within seasons a dilution effect is evident in both surface and subsurface runoff (DOC concentrations generally decreasing with increasing runoff rate). However, peak concentrations occur at flow rates between 0.3 and 1 mm for the surface waters and at less than 0.05 mm for subsurface waters. Daily basin averages (Figure 5b) are broadly similar to observed DOC concentrations of approximately 5–20 mg along the North Slope of Alaska (Connolly et al., 2018, 2020; Neilson et al., 2018).

Figure 4

Simulated leachate DOC concentration (mg ) of snowmelt runoff, surface runoff, and subsurface runoff loaded to streams and rivers of the Kuparuk River basin. Points plotted for days with runoff, at the surface or from a soil layer, for each basin grid cell. Symbol size is scaled by the magnitude of DOC loading, which is the product of simulated runoff and DOC concentration.

Figure 5

Daily simulated runoff (mm ), DOC concentration of leachate (mg ), and DOC loading (Gg C) for the four focus basins during 2006–2010.

Time series of basin average runoff, leachate DOC concentration, soil water DOC storage, and loading further illustrate model behavior for the four focus basins (Figure 5) The highest DOC leachate concentrations occur near the time of maximum runoff for each basin (see Supporting Information S1). Maximum DOC concentration is seen shortly after peak runoff for the two North Slope rivers and shortly before peak runoff for the Yukon and Mackenzie. For the Kuparuk, basin average DOC concentrations reach a maximum (20–25 mg ) at or soon after peak discharge, when runoff emanates from saturated near‐surface soils with relatively high SOM. The maximum runoff and DOC concentration timing differences are attributable largely to the extent, location, and concentration of DOC leachate in surface runoff. For example, for the Colville, much of the basin experiences surface runoff when basin average runoff peaks in late May, with concentrations of leachate loaded to river networks of 2–20 mg . Less area experiences surface runoff when basin average DOC concentration reaches a maximum several days later, and less water interacts with available surface carbon, resulting in higher concentration leachate as the freshet wanes. Subsurface runoff emanates from a similar proportion of the basin during this period from maximum runoff to maximum DOC concentration, and subsurface leachate concentrations are roughly similar (20–40 mg ). For the Mackenzie, DOC concentration of surface runoff is higher at the time of maximum basin average concentration in early April. More total runoff from snowmelt is carried to river networks when the runoff maximum occurs several weeks later, and it originates further north, primarily along the eastern slopes of the Rocky Mountains. Also, surface leachate concentrations loaded from those areas at this time are much lower. DOC concentrations loaded from the subsurface runoff are only slightly higher (56 vs. 43 mg ) during this time, and do not offset the influence of more surface water and much lower associated DOC concentrations. Rising DOC concentrations in the Colville in late fall occur from subsurface runoff, in the center of the basin, where soil carbon content is relatively high (40–50 kg C ), at extremely low flow rates. For the two North Slope basins, the increase from summer to fall is greater in the Colville (from 1–10 mg ) than the Kuparuk (from 4–9 mg ), owing to more variability in SOM. Maximum daily loading is nearly coincident with maximum runoff in most years for the Colville, Kuparuk, and Yukon.

For the Yukon, simulated DOC concentrations during winter low flow of around 5 mg agree well with data collected at Pilot Station in 2009 and 2010 (Wickland et al., 2012). However, our model estimates in late spring before and during the freshet (13–30 mg ) are notably higher than the observations (13–17 mg ). As expected, simulated DOC concentrations are consistently lower in autumn and winter, with relatively low temporal variability. Overall, the simulations point to high DOC concentrations during the spring freshet, followed by relatively lower concentrations in summer that originate primarily through subsurface runoff from the thawed active layer or in non‐permafrost soils. Production in soil water adds DOC to the soil water pool during autumn, winter, and spring, while mobilization during runoff events reduces the storage pool. Loading occurs primarily through surface runoff during the freshet.

The percent of DOC loading from subsurface flow (1981–2010) is 56%, 48%, 37%, 34% for the Mackenzie, Yukon, Colville, and Kuparuk respectively (Figure 6a). These fractions are generally consistent with the amount of thawed soil, both vertically (i.e., thawed depth) and horizontally (i.e., permafrost coverage) in each basin. Differences in the timing of maximum loading, relative to peak runoff among the four basins, is likely the result of a convergence of several factors including surface and soil carbon content, snowmelt runoff magnitude, and the form of the surface production, decomposition, and leaching algorithms.

Figure 6

(a) Fraction of DOC loading from subsurface runoff as a percent of total loading across the Y2M study domain. (b) Average (1981–2010) annual total DOC yield (g C ).

Leachate DOC concentration and loading was also derived for 24 drainage areas that span the North Slope of Alaska (Figure 7). Simulated leachate DOC concentrations in spring (May–June) and summer (July–September) for these basins show a distinct east‐west gradient that is broadly consistent with independent estimates derived from empirical relationships between measured fluvial DOC concentration and remotely‐sensed watershed slope (Connolly et al., 2018). Slope‐based river DOC concentrations at the outlet of each drainage basin were quantified by applying the empirical relationships defined in Connolly et al. (2018) to slope estimates for delineated watershed boundaries (using a 100 m digital elevation model) that align with the simulated basin boundaries of this study. Simulated spring and summer DOC concentrations for the Kuparuk average 11.4 and 5.3 mg , respectively. These estimates are consistent with the 8.4 and 6.0 mg averages derived from the slope‐based estimates. Likewise, the simulations captures lower concentrations in the adjacent Sagavanirktok River, with spring and summer averages of 3.4 and 1.5 mg . The east‐west gradient is generally attributable to shifts in topographic relief, subsurface hydrology, and soil organic carbon content that correspond to variations in the relative amounts of organic matter‐rich coastal plain terrain (higher in the west) versus organic matter‐poor mountainous terrain (higher in the east) in their watersheds (Connolly et al., 2018). In a few instances (particularly in the spring), the model overestimates DOC concentrations for the more western drainage areas. One possible explanation is that simulated DOC is precluded from in situ photo‐ and bio‐degradation, which is expected to lower DOC concentrations in river water samples that make up the slope‐based estimates (Cory et al., 2014). Another explanation is that the model uncertainty (e.g., in soil organic matter) is greater for several of these very small lake drainage basins that are resolved by one or two grid cells. Nonetheless, the overall agreement between simulated and slope‐based estimates of river DOC concentrations reveals that the processes driving fluvial DOC concentrations and loadings to north‐flowing Arctic rivers reflect seasonal changes in climate and surface/subsurface hydrology as well as regional variations in watershed landscape characteristics.

Figure 7

Simulated (sim) and empirically‐derived (obs) (a) spring and (b) summer DOC concentration (mg ) for 24 North Slope rivers arranged west (left) to east (right). The easternmost basin (Kongakut) is not adjacent to the Mackenzie delta, as there is a small area of the Northern AK region (Figures 1 and 9) just west of the delta, but east of the Kongakut. The empirically‐derived estimates were calculated from watershed slope‐concentration relationships (Connolly et al., 2018) and reported as averages over May–June and July–September, respectively. Simulated concentrations loaded to rivers represent the median from each seasonal distribution.

Figure 9

(a–e) boxplots of monthly total DOC loading (Gg C) for the Y2M region, Mackenzie, Yukon, western, and northern Alaska regions shown in Figure 1. Northern AK region includes a small area just west of the Mackenzie delta that is not part of the 24 North Slope basins examined in Section 3.2 (Figure 7). Each boxplot depicts the distribution of monthly total loading across each regions for the 30 months over years 1981–2010. Boxplot rectangles bracket the 25th and 75th percentiles. Whiskers extend to the 5th and 95th percentiles. Thick and thin horizontal lines mark the distribution mean and median respectively.

Seasonal and Annual Total DOC Loadings

Simulated DOC loadings (see Supporting Information S2) for the Mackenzie and Yukon were also compared to annual export data from river sampling. The annual total loadings agree well with the river data reported in other studies, and are well within the range of those export totals (Figure 8). Annually the majority of loading in the Mackenzie basin comes from subsurface runoff, while in the Yukon basin nearly half arises from surface runoff, including surface runoff attributed to snowmelt. In spring, loading originates primarily from surface runoff in both basins, though nearly half of the DOC loading emanates from subsurface flow in the Mackenzie. All loading in summer comes from subsurface runoff. In the Mackenzie, winter sees a significant amount of loading from snowmelt runoff, which mostly arises during April across southern parts of the basin. In the Yukon loading in winter comes primarily from subsurface runoff.

Figure 8

Annual, spring (May–Jun), summer (Jul–Oct), and winter (Nov–Apr) total DOC export and loadings (Gg C ) derived from observations and model simulations for the Mackenzie (a) and Yukon (b) River basins. Observed total export (Obs) from Holmes et al. (2012) for years 2001–2008. Simulated loadings from this study (Sim) are indicated as totals from leachate associated with subsurface, surface (non‐snowmelt), and surface snowmelt runoff for the period 1981–2010. Symbols mark the annual totals from other studies. For the Mackenzie: Annual export , Raymond et al. (2007); , Kicklighter et al. (2013). For the Yukon other estimates include Striegl et al. (2007); Spencer et al. (2009); and a range from Guo et al. (2012) shown by the vertical line. These totals differ slightly from those shown in Table 4 comparisons to annual and seasonal totals from Holmes et al. (2012), over different time intervals.

Table 4

Annual and Seasonal Total Observed DOC Export, Simulated DOC Loading and Percent Difference (%) for the Mackenzie and Yukon Rivers

River basin, season	DOC flux
	Export (Gg C)	Loading (Gg C)	Difference (%)
Mackenzie, annual	1377	1589	11.4
Mackenzie, spring	493	583	18.2
Mackenzie, summer	610	659	8.0
Mackenzie, winter	275	347	26.1
Yukon, annual	1472	1820	23.6
Yukon, spring	783	924	18.0
Yukon, summer	508	669	31.7
Yukon, winter	182	226	24.2

Note. Observed values (Gigagrams C) from Holmes et al. (2012). Consistent with that study, simulated loadings shown here are calculated as averages for 1999–2008 and 2001–2008 for the Mackenzie and Yukon, respectively, and thus differ slightly from those shown in Figure 8 for the full simulation period 1981–2010.

For the Mackenzie and Yukon, simulated annual loading shows reasonable agreement (25%) with observed DOC export totals reported by Holmes et al. (2012), with percent differences of 11.4% and 23.6%, respectively (Table 4). For seasonal totals, the MPE suggests reasonable agreement for four of the six seasonal comparisons. Consistent with the export totals, the model simulates the highest loadings in spring (May–June) and summer (July–October) for the Yukon and Mackenzie, respectively. Overestimates are not unexpected, as a fraction of the DOC loaded to river networks is assumed to be lost to in‐stream processing (Cory et al., 2014; Lauerwald et al., 2012). Much lower simulated loadings in winter are also consistent with the observations of riverine DOC export.

Unit area yields (1981–2010) range from 1 to 3 g C (Figure 6b), and are broadly consistent with estimates from ecosystem modeling (Kicklighter et al., 2013). Yield for the Mackenzie and Yukon basins average 0.95 g C and 2.1 g C , respectively. These simulated yields are in very good agreement (20%) with estimates derived from river sampling, which average 0.82 and 1.8 g C for the Mackenzie (1999–2008) and Yukon (2001–2008), respectively (Holmes et al., 2012).

Simulated DOC loading for the Y2M region averages 3907 Gg C annually for the simulation period 1981–2010. Loading to rivers draining to the Beaufort Sea (the combined Mackenzie and North Slope rivers) is 1943 Gg C. This amount is roughly equivalent to the loading to the Bering and Chukchi Seas from the Yukon and western Alaska rivers, which totals to 1964 Gg C. Loading reaches a maximum in May each year across the Mackenzie and Yukon basins, rivers across both western and northern Alaska, and for the full Y2M region (Figure 9). In light of water travel times to the coast, the timing in peak simulated loading agrees well with the June peak in DOC export reported in field studies using river sampling (Holmes et al., 2012; McClelland et al., 2014). Higher variability in loading during spring, as estimated by the coefficient of variation, is present for northern and western Alaska relative to the larger Mackenzie and Yukon rivers. Virtually no loading takes place in western and northern Alaska watersheds during January–April each year.

Uncertainties in Simulated Hydrology and DOC Loading

Confident assessment of model simulated DOC loadings is challenged by several factors, most notably the paucity of observations across the region at frequent intervals through the year for medium and smaller watersheds. Validation of riverine DOC loading and mass flux in the Arctic is hindered by the limited number of observations in smaller rivers. Flux estimates also require coincident discharge and DOC concentration measurements, limiting the utility in the latter for estimating mass exports. This lack of information limits our ability to understand the magnitude and fate of terrigenous DOC exported to coastal aquatic environments. Uncertainties in mass exports arise due to logistical constraints that limit sampling coverage across the entire flow season. In a meta‐analysis of 204 peer‐reviewed studies Shogren et al. (2020) estimated that 30% of the annual DOC export from the Kuparuk River occurs during the “shoulder seasons” (e.g., May and late August and September), noting that June, July, and early August are the most intensively studied times. The study also suggested that such sampling biases likely result in underestimates of material transported by Arctic rivers.

The flux of DOC in Arctic rivers is strongly dependent on water yield (runoff) from the contributing watershed (Holmes et al., 2012; Raymond et al., 2007). While simulated annual runoff aligns well with discharge measurements of the Colville, Kuparuk, Mackenzie, and Yukon rivers, incongruencies may be present for the other, smaller drainage basins. This may be an issue, in particular, where the reanalysis precipitation is not well constrained, such as river basins where orographic effects enhanced seasonal snow accumulation. Soil moisture variability is strongly influenced by subgrid‐scale topographic variations, which are not accounted for in our simulations at the moderate 25 25 km grid size. While the model includes thermal and hydrological effects of organic soils (Nicolsky et al., 2007), water phase change influence on soil thermal dynamics, and soil water drainage redistribution through numerical solution of the Richard's equation (Rawlins et al., 2013), the model lacks local effects on carbon mobilization and transport, such as thermokarst lake drainage (Jones & Arp, 2015), land subsidence (Coch et al., 2020; Streletskiy et al., 2017), stream‐lake biogeochemical interactions, rapid surface‐subsurface exchanges (Neilson et al., 2018), and hillslope water tracks (Evans et al., 2020; Rushlow & Godsey, 2017), and vertical movement of soil water DOC. While including these influences could potentially improve simulated DOC estimates at the local scale, process representations are not feasible to implement at the regional scale investigated here.

The DOC loading estimates in this study, though correlated with mass fluxes, do not represent coastal DOC export magnitudes. Losses of DOC in transit to carbon dioxide and sequestration in lake sediments, as well as other types of storage, are not explicitly modeled in this study. For a number of reasons, it is important to note that alignment of peak model simulated runoff and peak discharge at the river mouth are not expected, as surface impoundments of snowmelt due to vegetation and other local influences can slow water movement into river networks, and travel times in the large rivers can be a week or more. It is plausible that lags in water transience throughout a given watershed, particularly the larger ones in our study domain, can explain a portion of the error in timing of peak DOC concentration.

While estimates of losses and rates of in‐stream processing in aquatic systems are limited around the Arctic, studies of DOC processing for the Kuparuk River (Cory et al., 2014; Holmes et al., 2008; Rocher‐Ros et al., 2021) and regional‐scale estimates of carbon fluxes for inland aquatic ecosystems in Alaska (Stackpoole et al., 2017) suggest that lateral exports are less than the loadings for each river. However, lateral DOC export and DOC loadings may be more similar during the spring freshet, when cold temperatures, turbidity, and decreased residence times limit in‐stream decomposition. Moreover, in‐stream DOC production (Zou et al., 2006) offsets in‐stream losses of terrestrially loaded DOC. Studies continually put into question the degree of bio‐lability of riverine DOC and watershed inputs (Abbott et al., 2014; Behnke et al., 2021; Holmes et al., 2012; Mutschlecner et al., 2018; Spencer et al., 2008; Textor et al., 2019). River‐borne DOC may undergo complete or partial oxidation and subsequent export to the Arctic Ocean where it has important consequences for coastal ecosystems (Vonk & Gustafsson, 2013). Enhanced decomposition could potentially lead to increased DOC mobilization and loading, coincident with rising land‐atmosphere fluxes of carbon dioxide and methane, with some of those fluxes arising from decomposition of soil and riverine DOC pools (Xu et al., 2009). Significant amounts of DOC may reach Arctic waterways should permafrost soils be thawed and flushed into aquatic environments. In turn, some amount of riverine carbon may undergo photochemical decomposition (Cory et al., 2014). Considerable uncertainty exists in how decomposition of soil carbon, and in particular soil water DOC, will impact aquatic environments. Measurement programs that help elucidate the influence of environmental conditions such as temperature, moisture, substrate quality and nutrient availability are needed to better understand the mobilization and fate of soil DOC on arctic aquatic environments in a warming arctic. Future studies should also investigate whether ongoing warming will likely result in shifts to greater riverine DOC exports from groundwater sources (Walvoord & Kurylyk, 2016). Key calibration and validation data to assist model development include river chemistry and discharge measurements during the shoulder seasons (Shogren et al., 2020) and observations of the quantity and quality of DOC in subsurface flow (Neilson et al., 2018; Connolly et al., 2020).

Conclusions

This study investigates the loading of DOC to rivers of the western Arctic derived from process modeling that involves dynamics of terrestrial hydrology and DOC production and mobilization. Results suggest that simulated runoff reflects the characteristic pattern of high flow following snowmelt, subsequent recession, and the relatively low flows during the cold season. Simulated annual totals for the four large basins examined agree well with observed runoff, and the seasonal cycles, based on monthly climatological runoff, are well simulated.

Our simulated leachate DOC concentrations for Alaska's North Slope drainage areas show a distinct east‐west gradient owing to spatial variations in soil organic carbon content and watershed hydrology (Connolly et al., 2018). While basin averages for the Kuparuk are similar in magnitude to available observations, the simulated DOC concentrations for Kuparuk basin grids exceed the range in observations, suggesting that loading per unit of soil organic matter may be less sensitive than the model predicts. DOC yields for the Mackenzie and Yukon show reasonable agreement with annual and seasonal total loading for spring and winter, although overestimates do occur for summer loading for the Yukon. The modeling points to relatively equivalent loadings to watersheds draining north to the Beaufort Sea, and west into the Bering and Chukchi Seas. Overestimates relative to flux data derived from river measurements near the coast are expected, as processing within streams and rivers through biological degradation and photo‐oxidation reduces, in effect, the total mass transported downstream in the DOC pool (Lauerwald et al., 2012; Cory et al., 2014). In a scenario where losses for the North Slope rivers average 10%–20% and the large rivers 15%–25%, estimates in this study become closely aligned with river export mass fluxes described in recent studies.

Our modeling approach captures the spatial variability in DOC export, the characteristic seasonal variations in leachate DOC concentrations, and differences in concentrations of surface and subsurface runoff observed in field studies. Process modeling of the lateral transfer of carbon and other nutrients in Arctic rivers provide a contemporary baseline for future model refinements and related studies of the biogeochemistry of Arctic coastal waters.

Supporting Information S1

Click here for additional data file.

Movie S1

Click here for additional data file.

Movie S2

Click here for additional data file.

Movie S3

Click here for additional data file.