Quantifying Dissolved Organic Carbon Dynamics Using a Three‐Dimensional Terrestrial Ecosystem Model at High Spatial‐Temporal Resolutions

Arctic terrestrial ecosystems are very sensitive to the global climate change due to the large storage of soil organic carbon and the presence of snow, glacier, and permafrost, which respond directly to near surface air temperature that has warmed in the Arctic by almost twice as much as the global average. These ecosystems play a significant role in affecting regional and global carbon cycling, which have been traditionally quantified using biogeochemical models that have not explicitly considered the loss of carbon due to lateral flow of water from land to aquatic ecosystems. Building upon an extant spatially distributed hydrological model and a process‐based biogeochemical model, we have developed a three‐dimensional terrestrial ecosystem model to elucidate how lateral water flow has impacted the regional dissolved organic carbon (DOC) dynamics in the Tanana Flats Basin in central Alaska. The model explicitly simulates the production, consumption, and transport of DOC. Both in situ observational data and remote sensing‐based products were used to calibrate and validate the model. Our simulations show that (1) plant litter DOC leaching exerts significant controls on soil DOC concentration during precipitation and snowmelt events, (2) lateral transport plays an important role in affecting regional DOC dynamics, and (3) DOC export to the Tanana River is approximately 9.6 × 106 kg C year−1. This study provides a modeling framework to adequately quantify the Arctic land ecosystem carbon budget by considering the lateral transport of carbon affected by permafrost degradation. The quantification of the lateral carbon fluxes will also improve future carbon cycle modeling for Arctic aquatic ecosystems.


Introduction
Arctic terrestrial ecosystems are sensitive to global climate change due to the presence of snow cover, mountain glacier, and permafrost as well as vast soil organic carbon store that responds directly to temperature changes.These ecosystems also affect the global climate system through both biogeophysical and biogeochemical feedbacks.The observed warming of near surface air temperature in the Arctic has been almost twice as large as the global average (Hartmann et al., 2013).Many studies have suggested that Arctic terrestrial ecosystems have already undergone rapid changes in recent decades (Liao & Zhuang, 2017a;Sturm et al., 2001;Tanski et al., 2016).As the Arctic continues to warm faster than the global mean, improving our understanding in the water and carbon dynamics is critical for predicting how Arctic terrestrial ecosystems will change in the future (Bintanja & Selten, 2014;Pachauri et al., 2014;Vaughan et al., 2013).
Arctic land ecosystems store about 1,600 Pg organic carbon in soils, which is twice as large as the global total carbon in the atmosphere (about 800 Pg C) (Tarnocai et al., 2009;Zimov et al., 2006).Most organic carbon is stored within the permafrost, the frozen ground that remains at or below 0 °C for at least two consecutive years (Woo, 2012).There is evidence that permafrost has started to thaw due to warming (Marchenko et al., 2007).As temperature continues to increase in the Arctic, the rate of permafrost degradation will increase, and more carbon is expected to be released into the atmosphere (Frey & Smith, 2005) through vertical fluxes across the soil-air interface (Davidson et al., 2006) and into the ocean through lateral flow across the soil-water and landocean interfaces (Guo et al., 2007;Ping et al., 2011).Nevertheless, specific mobilization pathways and mechanisms of soil organic carbon from permafrost to waterways remain elusive (Gao et al., 2018).
The Arctic terrestrial ecosystems have complex landscape.Driven by glacier dynamics, seasonal snow cover, and soil thawing and freezing dynamics, the Arctic landscape features various aquatic and land ecosystems (Hinzman et al., 2005).In some regions, lake can occupy more than 30% of the land area (Zimov et al., 1997).The Arctic landscape can shift from one state to another within a very short period of time.For example, during snowmelt onset, land snow cover decreases, and stream discharge can reach its peak value within one week (Liao & Zhuang, 2017a).
To date, few studies have quantified the water and carbon dynamics associated with the complex landscape of the Arctic where lateral flow may exert a large influence.Traditionally, studies of land and aquatic ecosystems are conducted by different research communities.Hydrological models often focus on quantifying water cycle and lateral water flow (e.g., overland runoff) to form streamflow (Gao et al., 2010;Neitsch et al., 2005).In contrast, land surface models focus on both water and carbon cycles on large spatial domains without explicitly considering the role of lateral flow/fluxes.These land surface models assume that the land surface is "flat" within each relatively large grid cell, so horizontal interactions between adjacent land units are often ignored (e.g., Zhuang et al., 2002;Govind et al., 2009;Bocaniov & Scavia, 2018).Aquatic ecosystem models, on the other hand, usually focus on biogeochemical processes within aquatic ecosystems (e.g., lake and stream) (e.g., Tan & Zhuang, 2015).
Most existing land surface models also generally ignore the impacts of loss of dissolved organic carbon (DOC) in their carbon budget quantifications (e.g., Kalbitz & Kaiser, 2008;Neff & Asner, 2001;Wickland et al., 2007).It was not until the recent decades that the importance of DOC was beginning to be recognized.DOC participates not only in various microbial activities in terrestrial ecosystems, it also serves as an important nutrient source for aquatic ecosystems (Gao et al., 2018;Rier & Stevenson, 2002;Yurova et al., 2008;Znachor & Nedoma, 2009).To date, only a few land surface model studies have explicitly simulated DOC dynamics (e.g., Fan et al., 2010;Kicklighter et al., 2013;Lauerwald et al., 2017;Nakhavali et al., 2018;Tang et al., 2018;Wu et al., 2014;Yurova et al., 2008).However, these studies have not explicitly considered the impacts of lateral flow of water.Meanwhile, many models were developed to estimate DOC dynamics at site levels or regional scales (Camino-Serrano et al., 2018;Fan et al., 2010;Grieve, 1991;Lambert et al., 2014).Similar to land surface models, some of these DOC models do not consider lateral flow.Others used statistical and spatially distributed methods to estimate DOC dynamics (e.g., Futter et al., 2007;Futter et al., 2011;Lessels et al., 2015).However, large uncertainties remain in their DOC quantifications.
In the present study, a process-based three-dimensional ECOsystem model (ECO3D) was developed to quantify DOC dynamics in the Tanana Flats Basin in central Alaska.Special attention has been paid to the production, consumption and transport of DOC, by explicitly modeling the lateral flow and its impacts on DOC dynamics.We first developed a framework to simulate the water and carbon cycle dynamics using a three-dimensional modeling approach.We then quantified the role of lateral transport in affecting regional DOC dynamics.Finally, we quantified the DOC dynamics in both land and aquatic ecosystems in the region.

Study Area
The Tanana Flats Basin (TFB), located in interior Alaska near Fairbanks, drains an area of approximately 16,000 km 2 (Figure 1).This basin is a typical glacier and snowmelt-fed catchment with extreme relief from Mount Hayes in the eastern Alaska Range.Surface elevation decreases from about 4,000 m at Mount Hayes to around 120 m at the Tanana Flats in less than 50 km.In the southern part, the eastern Alaska Range forms the southern barrier of the glacier outwash.In the central part, glacial outwash and alluvial fan from the Alaska Range occupy the lowlands.The northern part of the alluvial fans is commonly recognized as the Tanana Flats near Fairbanks.The average elevation gradient in the Tanana Flats is only 1 m km −1 .Only a few well-developed stream channels including the Clear Creek and Wood River run through the flats prior to joining the Tanana River.In the northern part, the Tanana River, which is the largest tributary of the Yukon River, runs from the east near Big Delta to the west near Nenana.
Based on the latest geologic map of Alaska, approximately 68% of the study area is covered by unconsolidated sediments (e.g., soil and surficial deposits) of Quaternary age.Most sediments are deposited along the alluvial fans of the Alaska Range.The major soil type is "Gelisols" (United States, Soil Survey Staff, 1975).At the foot of the Alaska Range, sedimentary rocks, mainly Nenana gravel and coal-bearing rocks with ages from late Miocene to Pliocene, cover up to 11% of the study area.At the Alaska Range, land surface and bedrock are mainly composed of igneous rocks (Wilson et al., 2015).Mountain glaciers near the Alaska Range cover up to 1% of the basin based on the World Glacier Inventory data (Figure 2) (WGMS and NSIDC, 1999).Recent studies suggested that the volumes of these glaciers are decreasing due to the warming climate (Arendt et al., 2002;Hill et al., 2015).Seasonal snowmelt patterns are also changing in recent decades (Liao & Zhuang, 2017a).
The study area lies entirely within discontinuous permafrost zones under current climate.Approximately 60% of the study area is underlain by permafrost (GIPL, 2011).Permafrost-free zones are the major pathway for groundwater discharge to the land surface.Both groundwater and surface water (GW-SW) near  1).Purple circles are soil carbon sites.Cyan circles are soil dissolved organic carbon sites.Yellow circles are stream dissolved organic carbon sites (Table 2).Indices within circles (e.g., 1, 2, and 3) are IDs in corresponding figures and tables.Fairbanks have been studied intensively for decades (Anderson, 1970;Viereck et al., 1993;Williams, 1965).
Ecosystem of the Tanana Flats has also drawn lots of attentions, with an emphasis on permafrost degradation and vegetation dynamics (Jorgenson et al., 2001;Viereck et al., 1993).Currently, evergreen forest and shrub are the dominant vegetation types occupying the lower and higher elevation areas, respectively (Fry et al., 2011).Many studies have suggested that warming-induced permafrost degradation has gradually led the ecosystem shifting from boreal forests to fens and bogs in this region (Jorgenson et al., 2001;Viereck et al., 1993).
Multiple data sets were used to conduct our model simulation and evaluation.Details of these data sets are listed in Table 3.

Model Description 2.2.1. Overview
The ECO3D model was developed based on a spatially distributed hydrological model, Precipitation Runoff Modeling System (PRMS) (Leavesley et al., 1983), and a process-based biogeochemistry model, Terrestrial Ecosystem Model (TEM) (Zhuang et al., 2003).The PRMS model is a spatially distributed and process-based modeling system to quantify the impacts of climate and land use on streamflow and watershed hydrology (Markstrom et al., 2015).TEM is a process-based ecosystem model driven by air temperature, precipitation, and radiation to estimate terrestrial ecosystem carbon and nitrogen dynamics (Chen et al., 2011).To simulate DOC dynamics, a process-based DOC module was developed based on earlier work (e.g., Fan et al., 2010;Yurova et al., 2008).Because both PRMS and TEM were described in earlier studies (Liao & Zhuang, 2015, 2017a), only the core algorithms are provided below.Additional algorithms are presented in the supporting information.
In the ECO3D, water and carbon cycles are seamlessly coupled.First, for an individual unit, water and carbon cycles are always coupled similar to most land surface models (Running et al., 2010;Zhuang et al., 2003).For example, soil heterotrophic respiration is directly influenced by soil moisture.Second, water and carbon cycles are coupled in the spatial domain through an interaction module, in which each computational unit is connected with its neighbors.During simulation, different interactive processes are simulated at each time step.For example, when two land units are neighbors, both surface and subsurface runoff along with corresponding DOC fluxes are transported from upslope to downslope.Another example is that when a land unit  2).Red dots represent the glacier from World Glacier Inventory.
and a stream unit are connected, surface and subsurface runoff as well as DOC fluxes are transported from land to stream.
Below we first introduce the water cycle algorithms that quantify water pools and fluxes.Second, we describe the carbon cycle algorithms to model the changes of various carbon pools and fluxes.Finally, we illustrate the numerical method to simulate carbon and water dynamics in a 3-D manner.

Water Cycle
Most hydrological modules in the PRMS model were implemented or revised within the ECO3D model.New modules were added to consider additional processes.In PRMS, hydrological modules and algorithms including canopy interception, snow accumulation and ablation are used to estimate the water stored or exchanged among various hydrologic response units (Markstrom et al., 2015).The newly added and revised modules from PRMS are briefly described below.

Canopy Interception
In the PRMS model, canopy density is defined using summer/winter canopy density, which may introduce great uncertainty due to coarse temporal resolution.In ECO3D, we revised the interception algorithms using the Moderate Resolution Imaging Spectroradiometer (MODIS) leaf area index product to consider the impacts of temporal variations of canopy density (supporting information Text S1) (Myneni et al., 2015).
Because canopy interception was revised, the evapotranspiration (ET) algorithm has been updated.The new ET module redistributes the potential evapotranspiration (PET) to plant canopy, snow cover, and soil zone.

Glacier Dynamics
In the PRMS model, glacier was treated as normal snowpack.In ECO3D, a glacier module was developed to account for glacier dynamics.The mass balance of the glacier thickness is expressed as where GWT t is glacier water equivalent (m) at time step t (from now on t represents time step unless otherwise stated); P is precipitation in forms of rain or snow (m day −1 ); E is glacier sublimation (m day −1 ), which is modified from snow sublimation algorithm; SM is surface glacier melting, which is estimated using classical energy-based snowmelt algorithm (m day −1 ); and GM is glacier melting at the glacier-rock interface (m day −1 ).
Due to the regelation effect at the glacier-rock interface, glacier melting occurs even when air temperature is below 0 °C (Hock, 2005;Van der Veen, 2013).The term GM is used to account for this effect according to where K gm is the linear melting coefficient (fraction), T air is near surface air temperature (K), and T base is the reference temperature, which is an elevation-dependent parameter (K).

Surface and Subsurface Cascade
To route both surface and subsurface water flow, we developed a cascade module based on the USGS Cascading Routing Tool, which is used in PRMS (Henson et al., 2013;Markstrom et al., 2008).This routing module is also the core component of the interaction module in our ECO3D model.In this new cascade module (Figure 3), the fraction of flow from upslope to downslope is calculated according to where Fraction i,j is the fraction of flow from upslope grid i to downslope grid j, Elev i and Elev j are the elevations of upslope and downslope grid cells, respectively (m), and N is the total number of downslope directions of each grid cell.

Stream Flow Routing
Because the PRMS model uses a segment-based approach to simulate the stream flow routing, it is impossible to simulate the solute concentration at different locations within the same stream segment (Markstrom et al., 2015).Therefore, we developed a new reach-based stream flow routing module (Figure 4).For each stream reach, the incoming flow is calculated as 10.1029/2019MS001792 Journal of Advances in Modeling Earth Systems where Inflow reach is the total inflow (m 3 min −1 ), Inflow upstream is inflow from upstream (m 3 min −1 ), Inflow land is inflow from land surface (m 3 min −1 ), Inflow soil is inflow from subsurface (m 3 min −1 ), and Inflow external is inflow from external source (m 3 min −1 ).The same Muskingum routing method from PRMS is used to calculate the outflow from each reach.

Carbon Cycle
Carbon cycle was developed based on TEM.In addition, we developed the DOC model to account for DOC pools and fluxes.The ECO3D model mainly simulates three carbon pools including vegetation carbon, plant litter carbon, and soil carbon.
Many studies suggested that DOC dynamics are associated with various biogeochemical and biophysical processes.For example, litter-derived DOC plays an important role in soil DOC dynamics (Froberg et al., 2005;Hafner et al., 2005;Scheibe & Gleixner, 2014;Uselman et al., 2007).Therefore, we explicitly simulate the production, consumption, and transport of DOC using three modules: litter DOC, soil DOC, and stream DOC modules.Below, we first introduce the core carbon pool and flux algorithms in vegetation, litter, and soil with littler DOC and soil DOC included in litter carbon and soil carbon, respectively.Then, we introduce the stream DOC dynamics modeling.Because there are few lakes in the study area, lake DOC is not considered here.

Vegetation Carbon Dynamics
Carbon fluxes including photosynthesis and respiration are simulated.Most algorithms of these processes are from TEM and MODIS Gross Primary Production (GPP) algorithms (Chen et al., 2011;Heinsch et al., 2003).A few new modules were added to account for additional processes.Wood mortality is currently not considered in the model.More details of these algorithms are discussed in the (supporting information Text S2.  10.1029/2019MS001792

Journal of Advances in Modeling Earth Systems
For each vegetation type, the changes in total vegetation carbon (C V ) is expressed as (Running et al., 2010) where C t V is the total vegetation carbon (kg C m −2 ), GPP is gross primary production (kg C m −2 day −1 ), R A is autotrophic respiration (kg C m −2 day −1 ), C lf is canopy litterfall (kg C m −2 day −1 ), and C dead ¯root are dead root decomposition (kg C m −2 day −1 ).GPP was estimated based on the MODIS GPP algorithm (Heinsch et al., 2003).Shortwave radiation is estimated using the PRMS radiation module (Markstrom et al., 2015).At each time step, GPP is allocated to canopy, stem and root.The R A was estimated as the sum of maintenance respiration (R M ) and growth respiration (R G ) (Zhuang et al., 2003).Similar to GPP, R A is also allocated to vegetation canopy, stem, and root.Net primary production (NPP) is estimated as the difference between GPP and R A , and net ecosystem production (NEP) of each land unit is then modeled as where R H is total heterotrophic respiration estimated using the classical Q 10 function (kg C m −2 day −1 ).DOC net is the total net DOC flux from DOC routing from one pixel to others (kg C m −2 day −1 ).It can be either positive or negative (see supporting information).Starting from here, all DOC net from routing can be either positive or negative unless otherwise specified.

Litter Carbon Dynamics
Because litter DOC can only be released from litter soluble organic carbon (LSOC), multiple pools are required to account for the transitions between LSOC and insoluble organic carbon (LIOC).Besides as many studies suggested, litter decomposition rate varies with substrate quality at different stages (Berg & McClaugherty, 2003).Taken together, we modeled the litter substrate using a three-pool model and litter decomposition is modeled with three phases.In the early phase, litter mass loss rate decreased with time due to rapid consumption of available carbon by microorganisms, with high rates of carbon mineralization (CO 2 production by heterotopic respiration).Soluble organic matter, which ultimately turns into DOC, also has the highest production rate.In the intermediate phase, as decomposition became carbon-limited, lignin degradation is hampered.Previous incompletely decomposed organic carbon continues to decompose.This is also the phase when a fraction of LIOC is converted to LSOC.Both CO 2 and DOC production rates decrease.In the later phase, CO 2 production is limited by DOC and they are positively correlated.Highly decomposed litters enter into soil in the form of humus.This model also considers continuous fresh litterfall inputs and water flow (surface runoff and leaching).
For the litter as a whole, the change of total litter carbon was calculated as where C t L is the total litter carbon (kg C m −2 ); C lf is carbon flux from canopy fresh litterfall (kg C m −2 day −1 ), which is calculated in vegetation carbon dynamics module (process 1 in Figure 5); DOC net is overland net DOC flux from DOC routing (kg C m −2 day −1 ) (process #4); DOC leaching is vertical DOC flux due to surface leaching (kg C m −2 day −1 ) (process #3); C humus is vertical carbon flux that leaves the litter in the form of humus (kg C m −2 day −1 ) (process #2); and R H ¯litter is litter heterotrophic respiration (kg C m −2 day −1 ) (process #9 and #10).Because litter respiration occurs throughout all decomposition phases, term R H ¯litter is the total respiration from these phases.Other studies also suggested that litters, regardless of whether they are decomposed or not, can be transported during intensive surface runoff events (Meyer et al., 1998).
The kinetic method was used to simulate different decomposition phases: (1) litter organic carbon to LIOC (particulate organic carbon (LPOC)) (process #6) and LSOC (process #7), (2) LIOC to LSOC (Berg & McClaugherty, 2003) (process #8), and (3) LIOC and LSOC to humus (process #2).During all phases, some organic carbon is consumed through heterotrophic respiration (R H ¯litter ).When surface runoff or precipitation passes through plant litters, LSOC is carried away through water flow in the forms of surface runoff and leaching (process #4) (Froberg et al., 2005;Hafner et al., 2005;Meyer et al., 1998).This is also the only 10.1029/2019MS001792 Journal of Advances in Modeling Earth Systems process that produces DOC.While LPOC can also be transported through water flow, it is not currently considered in the ECO3D model.
The mass balance changes in insoluble organic carbon and soluble organic carbon are calculated as where C t ioc is the total insoluble organic carbon (kg C m −2 ), C ioc ¯new is the newly decomposed LIOC from undecomposed litter (kg C m −2 day −1 ), C ioc ¯2 ¯soc is carbon flux from LIOC to LSOC (kg C m −2 day −1 ), C ioc ¯2 ¯humus is carbon flux that leaves LIOC to soil O horizon (kg C m −2 day −1 ), R h ¯ioc and R h ¯soc are respiration from LIOC and LSOC, respectively (kg C m −2 day −1 ), C t soc is the total soluble organic carbon (kg C m −2 ), DOC soc ¯2 ¯runoff is the DOC flux leaving SOC (kg C m −2 day −1 ), DOC soc ¯2 ¯leaching is vertical DOC flux that leaves LSOC due to surface leaching (kg C m −2 day −1 ), and C soc ¯2 ¯humus is carbon flux that leaves LSOC in the form of humus (kg C m −2 day −1 ).
The total amount of LSOC that can be carried away through surface runoff and leaching is calculated based on the total water flow, maximum DOC concentration, available LSOC and DOC from upslope.The outflow DOC flux is partitioned between downslope surface runoff and surface leaching according to

Soil Carbon Dynamics
Soil carbon module was modified based on the TEM soil model (Zhuang et al., 2002).In addition, we explicitly considered the impacts of DOC fluxes: where C t soil is total soil carbon (kg C m −2 ), C humus is vertical carbon flux that enters into soil O horizon in the form of humus from litter carbon model (kg C m −2 day −1 ) (process 2 in Figure 6), C root is carbon flux from dead root decomposition from vegetation carbon model (process #3) (kg C m −2 day −1 ), DOC net is subsurface net DOC flux from DOC routing (process #5 and #6) (kg C m −2 day −1 ), DOC leaching is vertical DOC flux from litter surface leaching (process #4) (kg C m −2 day −1 ), and R H ¯soil is soil heterotrophic respiration (process #1) (kg C m −2 day −1 ).More details were discussed in the (supporting information S4. The soil DOC module was developed based on previous studies (Fan et al., 2010;Kalbitz et al., 2000;Yurova et al., 2008).In this module, soluble organic carbon is divided into two states: dissolved (DOC) and sorbed, which is potentially mobile, but currently part of the particulate phase (SPSOC).The balance between the two states changes with time due to changes in first-order kinetic in adsorption and desorption.A convection-dispersion equation was originally introduced in the one-dimensional model without lateral and vertical water flow.ECO3D considers both vertical and horizontal water fluxes (Figure 6).
P and μ 1 are calculated according to (Yurova et al., 2008) T−T basal 10 10 ( 16) where P basal is reference microbial DOC production rate (mg g −1 hr −1 ); μ basal is reference microbial DOC mineralization rate (hr −1 ); T and T basal are temperatures (K).For anaerobic conditions, both production and mineralization rates are adjusted.
The DOC fluxes term ∇ • (q × c) is calculated using lateral DOC flow and vertical DOC leaching according to where DOC upslope and DOC downslope are lateral DOC fluxes from upslope and to downslope (kg C m −2 day −1 ), and DOC leaching is vertical DOC flux due to litter surface leaching (process 4) (kg C m −2 day −1 ).
DOC adsorption and desorption are calculated according to (Yurova et al., 2008): where τ des is kinetic rate constant of desorption (hr −1 ), K d is distribution coefficient, which characterizes an equilibrium between the sorbed and dissolved phases (L g −1 ), and s is the mass of SPSOC (fraction).More details of the soil carbon model are given in the supporting information.
Because the groundwater module in PRMS uses an empirical method to simulate the interactions between soil water and shallow groundwater, it cannot simulate deep groundwater systems.As a result, the current ECO3D model cannot simulate the lateral DOC fluxes from groundwater systems.Considering the fact that groundwater flow in the study area is very slow due to the presence of permafrost, we assume this impact is negligible.
To simulate the DOC transport in hydrological networks, we developed a mass balance solute transport module based on the USGS PRMS stream flow routing module and Groundwater Solute Transport Simulator (MT3D-USGS) (Zheng & Wang, 1999).The mass balance equation of DOC of each reach is calculated according to where A is the cross-section area of reach (m 2 ), x is length of reach (m), C is reach DOC concentration (mg C L −1 ), t is (inner) time step (minute) (Table 4), Q is stream flow rate (m 3 min −1 ), D is dispersion coefficient (fraction), DOC lateral ¯land is lateral DOC fluxes from surface runoff (kg C min −1 ), DOC lateral ¯soil is lateral DOC fluxes from soil zone subsurface flow (kg C min −1 ), DOC net is net DOC flux from routing (kg C min −1 ), DOC source is external DOC sources flux (kg C min −1 ), and DOC sink is external DOC sinks flux (kg C min −1 ).The current ECO3D does not simulate stream biogeochemical processes (e.g., POC degradation and DOC decomposition).However, we believe DOC photodegradation can be negligible because boreal rivers are normally turbid and light penetration is limited.

Numerical Method 2.3.1. Spatial and Temporal Discretization
In ECO3D, we introduced a term called "Column Unit" (CU), which is an extension of hydrologic response unit and grid cell.Horizontally, similar to most hydrological models, ECO3D uses the grid-based spatial 10.1029/2019MS001792 Journal of Advances in Modeling Earth Systems discretization (Markstrom et al., 2015).Each CU occupies a grid cell, which has the same resolution as the DEM data.As a result, the spatial domain is discretized into a 2-D matrix of CUs (Figure S1).The spatial extent and resolution are defined using the typical watershed delineation process (Figure 2).The hydrological networks were retrieved from watershed delineation using the 1:100,000-scale (medium resolution) NHD flowline and the downscaled DEM data (Table 3) (Gesch et al., 2002;United States Geological Survey, 2013).The System for Automated Geoscientific Analyses and Arc Hydro tool were used to delineate the watershed following several steps: (1) "burn" the NHD flowline through DEM reconditioning, (2) fill the depressions in DEM, (3) calculate the flow direction and accumulation, (4) define the drainage line and subbasin, and (5) correct potential errors (e.g., spikes on the edges) (Esri Water Resources Team, 2011; Olaya, 2004).
Vertically, each CU is discretized into a number of layers, which depends on the type of CU.For example, a land CU with vegetation coverage has canopy layer, litter layer, and soil layer (Figure 7), whereas a stream CU has stream layer and groundwater layer (Figure 8).Atmosphere layer is shared among all CUs.Currently, a total of five CU types, including glacier, lake, land, stream, and swale, are defined.
ECO3D has a default daily (outer) time step for all water and carbon cycle processes.For individual process, a subdaily (e.g., half daily or hourly) (inner) time step is used when necessary (Table 4).For example, as stream flow is simulated at reach level, the travel time within each reach decreases significantly.Based on Manning's equation, the inner time step is set to minute based in the new stream flow routing module (Te Chow, 1964).Consequently, the solute transport module also has a minute-based temporal resolution.

Time-Variant Data and Parameters
We explicitly considered the ecosystem dynamics through time-variant data and parameters.First, time series climatic data including temperature, precipitation, and vapor pressure were used at daily time step (Liao & Zhuang, 2017a).Second, time series land use and land cover data were used to account for vegetation dynamics.In addition, we used high spatial-temporal resolutions leaf area index data to consider the vegetation canopy dynamics.Third, we considered the spatial and temporal variations in model parameters.For example, subsurface soil water flow parameters were adjusted when soil is partially frozen in winter.Vegetation parameters were updated whenever there is a land cover change.

Computational Sequence and Implementation
A computational sequence was designed to control the model simulation.First, a number of system state variables are initialized or read in as the initial state.Second, time-invariant data sets including DEM and  ECO3D was written in C ++ 11 and the Object-Oriented Programming method was used exclusively.For example, each CU is treated as an object during simulation period.OpenMP was used to improve computational efficiency.As a three-dimensional land surface model, ECO3D is highly nonlinear, which makes the model parameterization process highly computationally expensive.Model parameterization was calibrated with two steps.First, because the ECO3D model is developed based on the PRMS and TEM models, calibrated parameters from previous studies were imported directly (Heinsch et al., 2003;Liao & Zhuang, 2017a;Running et al., 2010;Zhuang et al., 2003).Second, we used a model-independent Parameter ESTimation and uncertainty analysis package (PEST) to automate the model calibration (Doherty et al., 1994).

Model Calibration
In PEST, the cost function Φ is defined as where c is the observation vector, c 0 is the "linearized" system simulated vector, J is the Jacobian matrix of the linearization function M which maps n-dimensional parameter space into m-dimensional observation space, b is the system parameter vector, and b 0 is the corresponding system parameter vector using M. To minimize the cost function, the "linearized" system parameter vector b 0 is updated through iterations until Φ is reduced to certain criteria.We used the Message Passing Interface (MPI) version of PEST, BeoPEST, to conduct parameter estimation (Hunt et al., 2010).To further reduce the number of parameters, the spatial domain was  Journal of Advances in Modeling Earth Systems divided into a number of zones, within which some parameters are assumed constant.These zones are defined using vegetation and soil types, etc. Special attention has been paid to topography effects because some parameters are elevation-dependent (Bell & Moore, 1999;Molotch et al., 2005;Winstral et al., 2002).
To calibrate the model, we mainly used hydrological data of stream discharge and snow water equivalent to calibrate the hydrological processes using PEST (Table 5).Because in situ carbon measurements were rare in this remote region, they were only used for evaluation purposes.
We ran the ECO3D model until it reached steady state, in which the changes in major carbon pools (vegetation carbon, litter carbon, and soil carbon) are within predefined thresholds.After that, we continued the simulation from preindustrial conditions (the start of year 1850) to the end of year 2015.
The model captured the time series data of stream discharge and snow water equivalent dynamics at the basin outlet and SNOTEL site, which is consistent with earlier studies (Figures S2-S5) (Liao & Zhuang, 2017a).

Model Evaluation
After ECO3D was calibrated, we evaluated the model performance of carbon cycle using both remote sensing-based products and in situ measurements.
First, we evaluated the simulated GPP with MODIS products (MOD17A2).Although the GPP algorithm used in ECO3D is the same as the MODIS GPP algorithm, differences between our estimates and MODIS GPP are still visible (Figure 10).On average, simulated GPP is 5% lower than MODIS GPP, possibly due to the overestimation in MODIS GPP (Turner et al., 2006).Spatial patterns of the simulated GPP and MODIS GPP are similar (Figures 10a and 10b).We also evaluated the simulated NPP with MODIS NPP products (MOD17A3).The simulated annual NPP is also lower than MODIS NPP.Similar to GPP, the simulated and observed NPP spatial patterns are consistent (Figures S6a and S6b).We further evaluated the simulated soil carbon, soil DOC concentration, and stream DOC concentration using in situ measurements from the Bonanza Creek Long Term Ecological Research Program (BNZ-LTER) (Figure S7).Details of these sites are described in Table 2.Because in situ soil carbon data only include those from the top 15 cm, an empirical relation was developed to estimate total soil carbon from a nearby site (Hicks-Pries & Schuur, 2009;Hollingsworth & Schuur 2007).The comparisons show that the differences between simulated soil carbon and in situ measurements are less than 5% at two sites (TKN0034 and TKN0133) and are about 20% at the other site (TKN0102).The distribution boxplots also show that soil carbon varies as much as 50% between different seasons.The range (15-26 kg C m −2 ) of soil carbon is also consistent with other studies with similar environments.
In situ soil DOC concentrations show significant variations within samples at the same location (Figure 11).
For example, at the PL site, the measured minimal and maximal soil DOC concentrations are 0.7 and 42 mg C L −1 , respectively (Kane et al., 2006, Petersen et al., 2012).The distribution of corresponding simulated soil DOC concentrations near the date of collection also has a large spread (0.1-40 mg C L −1 ).This variability is likely caused by local hydrological conditions within soil profile.
Last, we evaluated the simulated stream DOC concentration using observed DOC concentration (Moran, 2007;Rinella et al., 2008) (Table 2).Within the TFB, the comparisons show that the model performs reasonably well with a relatively large bias near glacier-fed headwater (Figure 12).The differences may also be caused by distances from in situ sites and actual stream channels.
Because the Tanana River receives inflow from upstream and continuous inflow DOC observations are thus far unavailable, no direct comparisons were made for the Tanana River.Instead, we compared its tributaries' DOC concentration with nearby observations.Comparisons show that the model captured both the magnitude and temporal variations.The highest simulated DOC concentration exceeds 30.0 mg C L −1 in several stream channels, whereas the highest observed concentration is 19.1 mg C L −1 in the Salcha River and Chena River (Cai et al., 2008).The trends of the simulated and observed DOC concentrations are also similar (Figure 13).The simulated stream DOC concentrations show more temporal variations when compared with in situ measurements.This could be due to the overestimates in litter DOC fluxes, affecting both lateral litter DOC and soil DOC fluxes.More details of model evaluation are given in the (supporting information Text S5.

Water Dynamics
In general, simulated water fluxes and reservoirs are close to previous PRMS simulation results (Liao & Zhuang, 2017a).For example, the Nash-Sutcliffe efficiency coefficients of simulated snow water equivalent and stream discharge are 0.7 and 0.8, respectively.
The simulated stream discharge showed significant spatial-temporal variations.First, the highest stream discharge occurs in spring due to snowmelt for most stream channels.Because this basin also receives upstream inflow, the highest stream discharge in the Tanana River (more than 1.5 × 10 8 m 3 day −1 ) is in summer.On an annual basis, the TFB contributes approximately 10% to the total stream discharge.Second, simulated stream discharge gradually increases from upstream to downstream.For example, during snowmelt/freshet period, the stream discharge of most stream headwater is less than 8.0 × 10 4 m 3 day −1 , whereas it is more than 5.0 × 10 5 m 3 day −1 near the subbasin outlet (Figure S8).

Carbon Dynamics
The simulated GPP shows significant spatial-temporal variations.First, the spatial patterns of GPP and vegetation type are similar.In regions near the Alaska Range, GPP is relatively low (less than 1.5 g C m −2 day −1 ) in growing season, but GPP is much higher (more than 6.5 g C m −2 day −1 ) in the Tanana Flats (Figure 10a).In nongrowing season, the simulated GPP is generally negligible.
Similar to GPP, the spatial pattern of the simulated NPP is comparable to that of vegetation type.NPP is close to zero in regions near the Alaska Range and it is much higher (about 3.5 g C m −2 day −1 ) in the Tanana Flats during growing seasons.No significant trend was detected in the annual NPP during the simulation period, and the average annual NPP is 2.8 × 10 6 kg C m −2 year −1 .
The simulated NEP, which indicates whether the ecosystem is a carbon source or a sink, also shows significant spatial-temporal variations (Figure S9).On an annual basis, approximately 30% of the study area has a NEP less than 15 kg C m −2 year −1 .Only less than 10% of the study area is a carbon source with NEP less than −50 kg C m −2 year −1 .The total NEP is positive during the growing season whereas its negative in  Journal of Advances in Modeling Earth Systems nongrowing season (Figure 17).Overall, the study area is close to neutral.However, the time series Mann-Kendall trend test shows that the annual NEP has increased by 2.6 g C m −2 year −1 (two-sided P value less than 0.001).
The simulated soil carbon storage shows significant spatial variations that are correlated with surface elevation.In regions near the Alaska Range, the average soil carbon is less than 10 kg C m −2 , whereas it is approximately 30 kg C m −2 in the Tanana Flats (Figure S10).Our estimates are very close to both in situ measurements and model simulations in similar environments (Johnson et al., 2011;Jorgenson et al., 2013;Michaelson et al., 1996).The simulated litter carbon storage is relatively small (less than 0.1 kg C m −2 ) when compared with vegetation (≈9.6 kg C m −2 ) and soil carbon storage.However, its importance is manifested by the litter DOC transport.Compared with other carbon pools and fluxes, DOC is relatively small in magnitude (Figure 17).

DOC Dynamics
The simulated litter-derived DOC shows significant spatial-temporal variations.First, litter-derived DOC is insignificant in winter but in other seasons, litter-derived DOC often increases significantly after precipitation events or snowmelt onset, which is also observed in other studies (Froberg et al., 2007;Hafner et al., 2005).Second, the spatial pattern of simulated litter-derived DOC is highly correlated with the cascade pathway, which is defined as the overland and subsurface flow path with high flow accumulation (Esri Water Resources Team, 2011).For example, on day 130 of year 2010, the spatial distribution of the simulated litter-derived DOC shows a visible channel shape pattern throughout the basin.The magnitude of litterderived DOC is also much higher (more than 100 mg C m −2 day −1 ) within the cascade pathway (Figure 14a).The simulated litter-derived DOC accounts for the abrupt changes in soil DOC concentration, which is correlated with the maximum litter-derived DOC concentration.Our estimates of litter-derived DOC concentration (more than 50 mg C L −1 ) are also in agreement with other studies (Deng et al., 2017;Froberg et al., 2007).
Lateral litter DOC flux contributes approximately 54% to the annual total DOC export, mainly during snowmelt onset.When lateral carbon flow is disabled, the simulated litter-derived DOC changes dramatically.Not only is its spatial pattern no longer correlated with the cascade pathway, but its magnitude is low near stream channels (Figure 14b).

Journal of Advances in Modeling Earth Systems
The simulated soil DOC concentration also shows significant spatial-temporal variations in the study region.Mediated by SPSOC adsorption and desorption, soil DOC concentration is relatively stable.However, the simulated soil DOC concentration exhibits a strong seasonality.In early spring, soil DOC concentration usually changes significantly due to abrupt increases in litter-derived DOC (Guo & Macdonald, 2006).In growing season, the average soil DOC concentration is approximately 30.0 mg C L −1 in the Tanana Flats, whereas it is less than 5.0 mg C L −1 in regions near the Alaska Range.The highest soil DOC concentration is usually near the headwaters and the cascade pathway (Figure 15a).Our estimates are consistent with previous results (Fan et al., 2010;Yurova et al., 2008).
Lateral soil DOC contributes approximately 46% to the total DOC export.When lateral carbon flow is disabled, the spatial distribution of soil DOC concentration is very different.For example, the simulated soil DOC concentration changes from the highest (with lateral DOC flux) to the lowest (without lateral DOC flux) near the cascade pathway (Figure 15b).The domain-average soil DOC concentration also becomes much lower (about 10.0 mg C L −1 ).
There are large spatial-temporal variations in the simulated reach DOC concentrations.In general, reach DOC deceases with increasing stream discharge.For example, the simulated reach DOC concentration at the Clear Creek outlet often reaches its lowest (2.0 mg C L −1 ) in early summer when its discharge rate is the highest.Reach DOC concentration usually reaches its maximum value in early spring due to snowmelt.Within the study region, the simulated reach DOC concentration decreases gradually from upstream to downstream (Figure 16).
Time series of total lateral DOC export to the streams indicate that DOC export reaches its maximum (2.5 × 10 6 kg C day −1 ) in spring.In some years, it may have more than one peak value.On an annual basis, approximately 9.6 × 10 6 kg C year −1 of DOC is exported to the hydrological networks (Figure 17).If this rate holds for the whole Arctic basin and considering the total Arctic land area is about 1, 000 times of the TFB, the total DOC export to Arctic river systems will be 10 Tg C year −1 , which is about one third of the total Arctic river DOC export reported by other studies (Holmes et al., 2012;Kicklighter et al., 2013;Raymond et al., 2007).Without lateral carbon flow, the ECO3D model cannot produce reach DOC dynamics because lateral DOC flux is the only DOC source other than DOC fluxes from upstream.Our study also confirmed that the total lateral DOC export is relatively small when compared to the total NPP/NEP in the region.On average, the annual total NPP is approximately 300 times the total DOC export (Figure 17).

The Importance of Lateral Flow
Our fully spatially distributed three-dimensional model has significantly extended the capability of existing hydrologic and ecosystem models.Unlike traditional 1-D approach, which omits or "teleports" water and carbon to aquatic ecosystems, ECO3D explicitly models the horizontal water and carbon flows between adjacent grid cells.
Because of that, we were able to trace the DOC all the way from plant litter to soil, and eventually to the stream networks.We were also able to improve the estimate of spatial distribution of carbon and water cycle as a whole.The ECO3D model can also be used to model the hillslope and riparian zone dynamics of carbon and water.

Challenges in DOC Modeling
At site level, DOC dynamics are directly affected by microbial activities, which are controlled by several environmental factors including soil moisture and nutrients availability.Spatially, DOC dynamics are strongly influenced by hydrological processes.Taken together, DOC dynamics are highly nonlinear in the three-dimensional domain, especially with snow melting and permafrost thawing.Therefore, it is imperative to use a 3D approach to model DOC dynamics under a changing climate.
Our simulations have shown that our modeled system is very sensitive to environmental changes because DOC pool and fluxes are orders of magnitude smaller than soil carbon and other fluxes.The mass balance of DOC requires modeling all DOC components including litter DOC, soil DOC, and others in the whole system.Failure to consider any component will likely to change the whole spatial pattern of DOC dynamics.

Implications for the Carbon Cycle
Although DOC flux is only a small fraction in the total carbon budget, it is a highly reactive carbon pool and could determine whether an area is a carbon source or sink (Figure 17).If we consider the changing climate in this area, the role of DOC is even more important.First, the increased litterfall resulting from increased productivity or landscape changes will increase DOC availability for leaching and transport.Second, ancient carbon and microbes stored within the permafrost will participate in carbon cycling once permafrost degradation occurs, resulting in a spatial expansion in DOC dynamics domain.Third, hydrological processes, which control DOC transport, will change significantly in both spatial and temporal patterns.Taken together, it is necessary to consider DOC, as well as dissolved inorganic carbon, POC, in future carbon budget study for a region that include both land and aquatic ecosystems across the landscape.

Model Limitations
There are a number of limitations in the model.First, groundwater and surface water interactions were simulated using an empirical method, which does not consider deep groundwater below the permafrost table.Consequently, possible DOC transport from the groundwater systems (e.g., Liao & Zhuang, 2017b;Walters et al., 1998) was not considered.Second, we did not explicitly simulate permafrost dynamics.Therefore, the estimated soil carbon only represents the carbon in the active layer.Third, we did not consider the impacts of natural disturbances (e.g., wildfires) in the region.Under a warming climate, groundwater flow, permafrost distribution, and wildfires are expected to affect both carbon and water dynamics significantly (e.g., Yoshikawa et al., 2002).Therefore, these processes and impacts should be factored into future analysis.

Conclusions
We have developed a ECO3D based on a spatially distributed hydrological model and a process-based ecosystem model.The model was calibrated and evaluated using both in situ observations and remote sensing-based products and applied to the Tanana Flats Basin in central Alaska.Our simulations show that the model satisfactorily captured the spatial-temporal variations of the dynamics of major carbon pools and fluxes.Although our simulations show much less carbon stored in plant litters, plant litter carbon provides an important DOC source for soil DOC and stream DOC during snowmelt onset and precipitation events.Lateral litter DOC flux and soil DOC flux contribute 54% and 46% to the total DOC export (9.6 × 10 6 kg C year −1 ), respectively.The developed model can be applied to the whole Arctic region to quantify the land ecosystem DOC dynamics.

Figure 1 .
Figure 1.The spatial location and major hydrological networks of the study area from Google Earth.The blue icon near the city of Fairbanks represents the snow telemetry (SNOTEL) site.Red circles with indices are United States Geological Survey (USGS) gage stations (Table1).Purple circles are soil carbon sites.Cyan circles are soil dissolved organic carbon sites.Yellow circles are stream dissolved organic carbon sites (Table2).Indices within circles (e.g., 1, 2, and 3) are IDs in corresponding figures and tables.

Figure 2 .
Figure 2. Hydrologic networks of the study area.HUC819040507 is the hydrologic unit code (HUC) from the USGS Watershed Boundary Dataset (WBD).Colored lines are stream segments.The Tanana River is represented by a number of connected stream segments, and the basin outlet is at the last stream reach.Green stars represent the five DOC concentration sites (Table2).Red dots represent the glacier from World Glacier Inventory.

Figure 3 .
Figure 3.The spatial distribution of cascade parameters of cascade module.Grey dots represent the stream channels.Arrows with different colors and thicknesses represent the directions and percentages of water flow from each CU to its downslope CU.

Figure 4 .
Figure 4. Schematic diagram of the reach-based stream routing algorithm.Grids of the same colors are reaches within the same stream segment.Blue and yellow arrows are lateral surface runoff and subsurface runoff, respectively.Red arrows are in-reach flow.

Figure 5 .
Figure 5. Schematic diagram of the modeled pools (square boxes) and fluxes (green arrows with indices) that control litter DOC production and transport.Litter consists of both undecomposed and decomposed solid litter.The latter consists of insoluble and soluble organic carbon (OC).The insoluble OC can be further decomposed to soluble OC.During all decomposition process, OC is consumed through respiration and leaves in gaseous form.Fully decomposed litters enter into soil as humus.Water fluxes (rain, runoff, etc.) carry in/out DOC.
downslope ¼ DOC avail −DOC soc¯2¯leaching (13) where DOC cap is the maximum DOC flux leaving the litter (kg C m −2 day −1 ), DOC max is the maximum DOC concentration (mg C L −1 ), DOC avail is the total DOC flux (kg C m −2 day −1 ), Infil is surface infiltration (m day −1 ), Runoff is surface runoff to downslope (m day −1 ), and DOC downslope is DOC flux to downslope (kg C m −2 day −1 ).More details of the litter carbon model are discussed in the (supporting information Text S3.

Figure 6 .
Figure 6.Schematic diagram of the modeled pools (square boxes) and fluxes (yellow arrows with indices) that control soil carbon, DOC production, consumption, and transport.

Figure 7 .
Figure 7. Water and carbon cycle processes simulated by the ECO3D model for a land column unit.Green lines and curves represent the vegetation.Blue arrows with indices represent key carbon and hydrological processes listed in the rightside table.
10.1029/2019MS001792 Journal of Advances in Modeling Earth Systems LIAO ET AL. hydrological networks are assigned.Third, ECO3D starts to run in either steady state or transient simulations.During each time step, water and carbon cycle processes are simulated within all CUs.After that, the interaction module simulates the horizontal interactions between adjacent CUs.At the end of each time step, system state variables of all CUs are updated.The general work flow of the model simulation is illustrated in Figure 9.

Figure 8 .
Figure 8. Water and carbon cycle processes simulated by the ECO3D model for a stream column unit.Blue arrows with indices represent key hydrological processes listed in the rightside table.

Figure 9 .
Figure 9. Computational sequence scheme of the ECO3D model.SS/TR represents steady state or transient simulation."Read" represents reading time-variant input data."Interaction" represents lateral flow interactions among Column Units."Update" represents updating both system states and model parameters.

Figure 10 .
Figure 10.Spatial distributions of the simulated GPP and MODIS GPP on day 185 of year 2010 (g C m −2 day −1 ).

Figure 11 .
Figure 11.Comparisons between observed and simulated soil dissolved organic carbon (DOC) concentration at two BNZ-LTER sites.The red dots represent the mean simulated soil DOC concentration near the collection date.Horizontal and vertical boxplots represent the distribution of observed and simulated soil DOC concentration, respectively (mg C L −1).

Figure 12 .
Figure 12.Comparisons between observed and simulated stream dissolved organic carbon (DOC) concentration at four Bonanza Creek Long Term Ecological Research Program (BNZ-LTER) sites.The red dots represent the mean simulated time series stream DOC concentration near the sample collection date.Boxplots represent the distribution of simulated stream DOC concentration (mg C L −1 ).

Figure 13 .
Figure 13.Time series of the simulated and observed stream dissolved organic carbon (DOC) concentration at different locations from 2002 to 2005.Line features represent simulated stream DOC concentrations at major river outlets.Symbols represent observed stream DOC concentrations along the Tanana River (Table 2) (mg C L −1 ) (USGS National Research Program and Cai et al., 2008).

Figure 14 .
Figure 14.Spatial distributions of the simulated litter-derived DOC with and without lateral carbon flow on day 130 of year 2010 (g C m −2 day −1 ).

Figure 15 .
Figure 15.Spatial distributions of the simulated soil DOC concentration with and without lateral carbon flow on day 130 of year 2010 (mg C L −1 ).

Figure 16 .
Figure 16.Spatial distribution of the simulated stream reach DOC concentration on day 130 of year 2010 (mg C L −1 ).

Figure 17 .
Figure 17.Time series of the simulated daily total lateral DOC export to the Tanana River and NPP/NEP (kg C day −1 ).The black, blue, and red lines and axis represent the daily DOC export, NPP and NEP, respectively.

Table 1
List of United States Geological Survey Gage Stations Used Journal of Advances in Modeling Earth SystemsLIAO ET AL. 4496 19422466, 2019, 12, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2019MS001792,Wiley Online Library on [11/11/2023].See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions)on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

Table 2
List of In Situ Sites Used for Evaluation a WGS84.bThese two sites are at unnamed tributaries.10.1029/2019MS001792

Table 3
List of Data Sets Used

Table 5
Estimated Values of the ECO3D Model Parameters for Major New Modules LIAO ET AL.