Linking global terrestrial and ocean biogeochemistry with process-based, coupled freshwater algae–nutrient–solid dynamics in LM3-FANSY v1.0

Estimating global river solids, nitrogen (N), and phosphorus (P), in both quantity and composition, is necessary for understanding the development and persistence of many harmful algal blooms, hypoxic events, and other water quality issues in inland and coastal waters. This requires a comprehensive freshwater model that can resolve intertwined algae, solid, and nutrient dynamics, yet previous global watershed models have limited mechanistic resolution of instream biogeochemical processes. Here we develop the global, spatially explicit, and process-based Freshwater Algae, Nutrient, and Solid cycling and Yields (FANSY) model and incorporate it within the Land Model (LM3). The resulting model, LM3-FANSY v1.0, is intended as a baseline for eventual linking of global terrestrial and ocean biogeochemistry in next-generation Earth system models to project global changes that may challenge empirical approaches. LM3-FANSY explicitly resolves interactions between algae, N, P, and solid dynamics in rivers and lakes at 1° spatial and 30 min temporal resolution. Simulated suspended solids (SS), N, and P in multiple forms (particulate or dissolved, organic or inorganic) agree well with measurement-based yield (kg km⁻² yr⁻¹), load (kt yr⁻¹), and concentration (mg L⁻¹) estimates across a globally distributed set of large rivers, with an accuracy comparable to other global nutrient and SS models. Furthermore, simulated global river loads of SS, N, and P in different forms to the coastal ocean are consistent with published ranges, though regional biases are apparent. River N loads are estimated to contain approximately equal contributions by dissolved inorganic N (41 %) and dissolved organic N (39 %), with a lesser contribution by particulate organic N (20 %). For river P load estimates, particulate P, which includes both organic and sorbed inorganic forms, is the most abundant form (64 %), followed by dissolved inorganic and organic P (25 % and 11 %). Time series analysis of river solid and nutrient loads in large US rivers for the period ∼ 1963–2000 demonstrates that simulated SS and N loads in different N forms covary with variations of measurement-based loads. LM3-FANSY, however, has less capability to capture interannual variability of P loads, likely due to the lack of terrestrial P dynamics in LM3. Analyses of the model results and sensitivity to components, parameters, and inputs suggest that fluxes from terrestrial litter and soils, wastewater, and weathering are the most critical inputs to the fidelity of simulated river nutrient loads for observation-based estimates. Sensitivity analyses further demonstrate a critical role of algal dynamics in controlling the ratios of inorganic and organic nutrient forms in freshwaters. While the simulations are able to capture significant cross-watershed contrasts at a global scale, disagreement for individual rivers can be substantial. This limitation is shared by other global river models and could be ameliorated through further refinements in nutrient sources, freshwater model dynamics, and observations. Current targets for future LM3-FANSY development include the additions of terrestrial P dynamics, freshwater carbon, alkalinity, enhanced sediment dynamics, and anthropogenic hydraulic controls.

1 Introduction

Dramatic increases in fossil fuel combustion, deforestation, agriculture, fertilizer use, and sewage outflows have increased loadings of terrestrial sediments and nutrients (e.g., nitrogen – N, phosphorus – P) to rivers and coastal waters and changed N:P ratios (Cordell et al., 2009; Fowler et al., 2013; M. Lee et al., 2019; Sytvitski et al., 2005). These changes in sedimentary and nutrient loadings have altered turbidity and biogeochemistry in many freshwater and coastal ecosystems. The consequences include (1) changes in ecosystem productivity and carbon (C) exports (Liu et al., 2021), (2) increases in frequency, duration, and severity of harmful algal blooms (HABs) (Anderson et al., 2002, 2021; Heisler et al., 2008; Paerl et al., 2018) and hypoxic dead zones (Diaz and Rosenberg, 2008), and (3) perturbations of aquatic plant, seagrass, and coral reef ecosystems, incurring substantial socioeconomic costs (Lacoul and Freedman, 2006; McLaughlin et al., 2003; Restrepo et al., 2006).

Resolving prominent drivers of the aforementioned aquatic ecosystem consequences requires a freshwater biogeochemistry model that captures intertwined algae, nutrient, and solid dynamics. In general, strong positive relationships have been observed between P and phytoplankton production in freshwaters, while N increases have been linked with the development of large algal blooms and hypoxic events in estuarine and coastal waters (Howarth and Marino, 2006; Smith, 2003). In particular, excessive inorganic nutrients, which are generally more bioavailable than organic forms (Sipler and Bronk, 2014), have been recognized as critical drivers of algal blooms (including non-HABs) and hypoxic events (Kemp et al., 2005). Shifts in community composition towards more toxic or harmful algal species have furthermore been attributed to changes in nutrient supply ratios, including N:P (Anderson et al., 2002; Heisler et al., 2008) and relative abundance of different N and P forms (e.g., nitrate – NO${}_{\mathrm{3}}^{-}$, Parsons and Dortch, 2002; ammonium – NH${}_{\mathrm{4}}^{+}$, Trainer et al., 2007, Leong et al., 2004; urea, Glibert et al., 2001, Glibert and Terlizzi, 1999; dissolved inorganic N and P – DIN and DIP, Glibert et al., 2008). Such shifts can be explained by differences in algal species-specific nutrient acquisition pathways that are controlled by nutritional status and preferences, uptake capability, and physiological status (Anderson et al., 2002). Furthermore, nutrient and algae dynamics are strongly linked with solid dynamics, for example, through phosphate (PO${}_{\mathrm{4}}^{\mathrm{3}-}$) sorption–desorption interactions with solid particles (McGechan and Lewis, 2002) and algae growth reduction due to light shading by suspended solids (SS) (Di Toro, 1978). Estimating river solids, N, and P in both quantity and composition resulting from intertwined algae, nutrient, and solid dynamics is thus necessary for understanding the development and persistence of many HABs and hypoxic events.

Building confidence in projected global freshwater biogeochemistry changes rests in part on the development of process-based models that are robust under unprecedented conditions expected in the next century. Process-based freshwater biogeochemistry models of the N cycle alone (LM3-TAN, Lee et al., 2014; DLEM, Tian et al., 2020; IMAGE-DGNM, Vilmin et al., 2020; INCA, Wade et al., 2002) or the P cycle alone (DLEM, Bian et al., 2022; IMAGE-DGNM, Vilmin et al., 2022) have been widely applied across scales. However, prior applications of models that capture coupled algae, multi-nutrient, and/or solid cycles, such as RIVE (Billen et al., 1994) and QUAL2K (Pelletier et al., 2006), have generally been limited to relatively small watersheds. Modeling river nutrient yields and loads on a global scale, in both magnitude and form, has been challenged by the difficulty of balancing desired details of instream biogeochemical processes with limitations imposed by available knowledge, input, and validation datasets.

Global NEWS (Mayorga et al., 2010) and IMAGE-GNM (Beusen et al., 2015, 2016) are global watershed models that have been widely used for pioneering simulations of river nutrient loads on global scales. Global NEWS estimates have been shown to be consistent with measurement-based estimates across a globally distributed set of major rivers and provided important nutrient inputs for global and regional ocean biogeochemistry model simulations. Global NEWS and IMAGE-GNM, however, do not resolve coupled algae, nutrient, and solid dynamics in freshwaters despite their intertwined dynamics. Global NEWS, representing a hybrid of empirical, statistical, and mechanistic components, formulates and implements different elements and their chemical forms independently based on basin-averaged properties. IMAGE-GNM applied at a global scale does not differentiate dissolved, particulate, inorganic, and organic nutrient forms, though such differentiation has been considered at regional scales. Global applications of both models do not mechanistically resolve instream biogeochemical processes.

Prior global watershed models are also limited in their capacity to represent process-based nutrient storage in terrestrial plants and soils. Global NEWS assumes that nutrients are in steady state and do not accumulate on land. IMAGE-GNM does not explicitly simulate terrestrial nutrient dynamics, such as vegetation growth, leaf fall, natural and fire-induced mortality, and soil microbial processes, but takes a mass balance approach to calculate dynamic soil nutrient budgets. Organic nutrient delivery to rivers, however, depends on changes in vegetation and soil organic nutrient storage in response to the aforementioned terrestrial dynamics under long-term (multiple decades to centuries) historical climate and land use changes. Terrestrial storage changes, for example, have been shown to significantly alter multi-decadal river nutrient trends (Van Meter et al., 2018; M. Lee et al., 2019) and seasonal to multi-year river nutrient extremes (Kaushal et al., 2008; Lee et al., 2016, 2021).

Here we develop the global, spatially explicit, and process-based Freshwater Algae, Nutrient, and Solid cycling and Yields (FANSY) model and incorporate it within the National Oceanic and Atmospheric Administration (NOAA)/Geophysical Fluid Dynamics Laboratory (GFDL) Land Model (LM3), which is capable of resolving coupled water, C, and N dynamics and storage changes in a vegetation–soil system (M. Lee et al., 2014, 2019). The resulting coupled terrestrial–freshwater ecosystem model LM3-FANSY constitutes a significant step toward a more process-based representation of coupled freshwater algae, nutrient, and solid dynamics. LM3-FANSY v1.0 is aimed at linking global terrestrial and ocean biogeochemistry towards next-generation Earth system models. Here we provide a detailed model description, performance assessment against measurement-based estimates of solids and nutrients across major world rivers, and sensitivity evaluation to a range of components, parameters, and inputs.

2 Model description

2.1 LM3-FANSY framework

LM3-FANSY is an expansion of NOAA/GFDL LM3-Terrestrial and Aquatic Nitrogen (TAN) (M. Lee et al., 2014, 2019) to include a terrestrial soil erosion process and comprehensive freshwater sediment and biogeochemical dynamics (Sect. 2.2). The terrestrial component LM3, which has been described in detail elsewhere (Gerber et al., 2010; Milly et al., 2014; Shevliakova et al., 2009), captures coupled water, C, and N dynamics within a vegetation–soil system. LM3 simulates transfers and transformations of three N species (i.e., organic N, NH${}_{\mathrm{4}}^{+}$, and NO${}_{\mathrm{3}}^{-}$) for vegetation and soil systems, considering the effects of anthropogenic N inputs, land use, atmospheric CO₂, and climate over timescales of hours to centuries. LM3 simulates the distribution of five vegetation functional types (C₃ and C₄ grasses, temperate deciduous, tropical, and cold evergreen trees) based on prevailing climate conditions and C–N storage in vegetation including leaves, fine roots, sapwood, heartwood, and labile storage. There are four soil organic pools (fast or slow litter and slow or passive soil) and two soil inorganic pools (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$). Scenarios of land use states and transitions are used to simulate four land use types (primary lands – lands effectively undisturbed by human activities, secondary lands – abandoned agricultural land or regrowing forest after logging, croplands, and pastures). LM3 captures key terrestrial dynamics that affect the state of vegetation and soil C–N storage, such as vegetation growth, leaf fall, natural and fire induced mortality, deforestation for agriculture, wood harvesting, reforestation after harvesting, and various soil microbial processes. LM3 extended to include a global river-routing and lake model (Milly et al., 2014) is thus well suited to simulate the delivery of terrestrial N to rivers and coastal waters.

The terrestrial component LM3, including the newly introduced terrestrial soil erosion process (Sect. 2.2.1), receives N inputs of fertilizer applications and atmospheric deposition, simulates biological N fixation, and estimates N outputs including net harvest (N in harvested wood, crops, and grasses after subtracting out internally recycled inputs, e.g., manure applied to croplands and sewage), emissions to the atmosphere, and eroded sediment and N fluxes from terrestrial to river systems. In addition to terrestrial runoff of three N species (dissolved organic N – DON, NH${}_{\mathrm{4}}^{+}$, and NO${}_{\mathrm{3}}^{-}$) introduced in our previous study (Lee et al., 2014), here we have added particulate organic N (PON) fluxes from terrestrial litter and soils (Sect. 2.2.1). M. Lee et al. (2014, 2019) provide further details on the terrestrial model.

The freshwater component FANSY receives N, P, and solids in multiple forms either from LM3 or from prescribed inputs (Sect. 3.1). It then simulates biogeochemical transformations and transport of each form of the nutrients and solids within streams, rivers, and lakes (Sect. 2.2).

2.2 Freshwater component FANSY

FANSY constituents of algae, nutrients, and solids in rivers and lakes are listed in Table 1 and described in Fig. 1. FANSY has 13 prognostic state variables and 5 diagnostic state variables. Inorganic suspended solid (ISS) is delivered from the terrestrial soil erosion process (Sect. 2.2.1). ISS dynamically interacts with benthic sediment inorganic solid (Sed) through deposition and resuspension processes. Primary interactions between SS and other model components are through the shading effect of turbidity on algae growth (Sect. 2.2.2) and the sorption of PO${}_{\mathrm{4}}^{\mathrm{3}-}$ to ISS as particulate inorganic P (PIP) (Sect. 2.2.4). Algae take up N and P, which is subsequently partitioned between organic and inorganic N and P pools via algae mortality (Sect. 2.2.2). Algae C (C_al) and algae P (P_al) are diagnosed from algal N (N_al) assuming the Redfield $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio (Chapra, 1997; Redfield et al., 1963). Chlorophyll a (Chl) is derived using the photoacclimation model of Geider et al. (1997) to predict a Chl-to-C ratio (r_ChlC), and Chl is calculated from r_ChlC and C_al. The five prognostic N variables contain reduced and oxidized dissolved inorganic forms (NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$), as well as dissolved and two particulate (suspended and sedimentary) organic forms (DON, PON, and benthic sediment N – SedN; Sect. 2.2.3). The five prognostic P variables include the same organic forms as for N variables (dissolved organic P – DOP, particulate organic P – POP, and benthic sediment P – SedP; Sec. 2.2.4). The other two prognostic P variables are dissolved and particulate inorganic forms (PO${}_{\mathrm{4}}^{\mathrm{3}-}$ and PIP). FANSY does not distinguish between PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DIP, and soluble reactive phosphorus (SRP). The subsections that follow (Sect. 2.2.1–2.2.4) provide a detailed description of each variable and associated processes.

Table 1 FANSY prognostic and diagnostic state variables.

Figure 1 FANSY structure with arrows depicting fluxes of constituents of algae, nutrients, and solids in rivers and lakes. The constituents are listed in Table 1.

Added solids and nutrients to streams and rivers are subject to retention within rivers and lakes or transformed during transport to the coastal ocean. Each model grid cell contains one river reach and/or one lake. Water containing solids and nutrients in each river reach or lake flows to another river reach in the downstream grid cell following a network that ultimately discharges to the ocean (Milly et al., 2014). Freshwater physics, hydrology, and hydrography are described in detail elsewhere (Milly et al., 2014). In each river reach or lake, for each prognostic variable i, settling–resuspension dynamics and/or biogeochemical reactions (SB_i) are calculated according to the process-based formulations described in Sect. 2.2.1–2.2.4. For example, if i is PON, SB_i is Eq. (28). If i is SedN, SB_i is Eq. (29). A general mass balance for a variable i in a river reach or lake at each computation time step (30 min in this study) is written as

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{X}_i}{\mathrm{d}t}}& =F_i^{\mathrm{in}}-F_i^{\mathrm{out}}+I_i+{\mathrm{SB}}_i,\text{if}i\text{is a river/lake water}\\ & \text{column variable, i.e.,}i=\text{all variables except}\\ & \text{Sed, SedN, or SedP},\end{array}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{X}_i}{\mathrm{d}t}}& ={\mathrm{SB}}_i,\text{if}i\text{is a river/lake benthic sediment}\\ & \text{variable, i.e.,}i=\text{Sed, SedN, or SedP},\end{array}\end{array}$$

where i is a prognostic variable listed in Table 1, X_i is the amount of variable i (kg), $F_i^{\mathrm{in}}$ and $F_i^{\mathrm{out}}$ are inflow and outflow of variable i (kg s⁻¹) through the river or lake network, I_i is input of variable i from terrestrial systems and/or the atmosphere (kg s⁻¹), and SB_i is settling–resuspension dynamics and/or biogeochemical reactions of variable i (kg s⁻¹).

2.2.1 Solid dynamics

In LM3-FANSY, terrestrial soil erosion is controlled by land surface slope, rainfall, and leaf area index (LAI) based on Pelletier (2012), as described in Eq. (3). N fluxes from terrestrial litter and soils, in the form of PON, in Eq. (4) are based on the simulated soil erosion fluxes and litter–soil N concentrations. This approach is consistent with that employed by several previous modeling studies (Tian et al., 2015; Zhang et al., 2022). The litter–soil concentrations for this purpose are estimated by using litter–soil contents and effective soil depths simulated by LM3 (Gerber et al., 2010).

Inorganic soil inputs to rivers are derived from the simulated soil erosion fluxes by subtracting the PON contribution, as described in Eq. (5). This requires an assumed ratio of POM:PON in eroded soils (i.e., 13.9, Table 2). Previous studies have shown a wide range of C content in tree biomass (∼ 42 %–61 %, Thomas and Martin, 2012) and of C:N ratios in litter and soils (∼ 5–500, Gerber et al., 2010, and references in Gerber et al., 2010's Table S1). This implies that POM:PON ratios in soil erosion fluxes can also vary significantly. We have found, however, that predicted river SS loads are insensitive to an order of magnitude variation in the ratio (i.e., 1.39 vs. 139, see Sect. 4.4) because organic contents in eroded soils are generally small. We thus used the POM:PON ratio of 13.9. The same ratio has been used to estimate the contribution of PON to SS in freshwaters, again noting that it is generally a small fraction of SS.

$$\begin{array}{}\text{(3)}& {\displaystyle }E=C_{\mathrm{1}}\cdot {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{b}}}{{\mathit{\rho }}_{\mathrm{w}}}}\cdot S^{\mathrm{5}/\mathrm{4}}\cdot r\cdot e^{-L},\text{(4)}& {\displaystyle }{E}^{\mathrm{PON}}=E\cdot \left({\displaystyle \frac{{\mathrm{N}}_{\mathrm{FL}}+{\mathrm{N}}_{\mathrm{SL}}+{\mathrm{N}}_{\mathrm{SS}}}{h_{\mathrm{s}}\cdot {\mathit{\rho }}_{\mathrm{b}}}}\right),\text{(5)}& {\displaystyle }{E}^{\mathrm{ISS}}=E-r_{\mathrm{DN},\mathrm{Ero}}\cdot E^{\mathrm{PON}},\end{array}$$

where E, E^PON, and E^ISS represent the terrestrial soil erosion flux (kg dry matter (D) m⁻² s⁻¹), terrestrial PON flux (kgN m⁻² s⁻¹), and terrestrial inorganic soil erosion flux (kgD m⁻² s⁻¹), respectively; C₁ is a free parameter of terrestrial soil erosion (unitless); ρ_b is soil bulk density (kgD m⁻³); ρ_w is water density (kg m⁻³); S is slope tan θ, with θ as the hillslope angle (unitless); r is rainfall (kg m⁻² s⁻¹); L is LAI (unitless); N_FL, N_SL, and N_SS represent the N content in the fast litter, slow litter, and slow soil pool, respectively (kgN m⁻²); h_s is effective soil depth (m); and r_DN,Ero is a POM-to-PON ratio in eroded fluxes (kgD kgN⁻¹).

Table 2 Model parameters along with their descriptions, values, units, and references and/or rationale.

Terrestrial soil erosion is known to be scale-dependent because it could be dominated by different spatial-scale processes (e.g., interrill, rill, and gully erosion, landsliding; Poesen et al., 1996; Renschler and Harbor, 2002). We thus adapt from Pelletier (2012) a degree of freedom via the coefficient of C₁ to account for the spatial resolution of input data (e.g., slope at 1° resolution). C₁ is a single global value scale parameter coarsely calibrated to reduce prediction errors of SS loads across major world rivers (see Sect. 2.2.5 for details in model calibration). Sensitivity of the model to C₁ is addressed in Sect. 4.4. It has been suggested to model soil erosion at event scales (daily or subdaily time steps) to account for episodic, substantial mass transport (Tan et al., 2017). We calculate soil erosion rates at the model time step (30 min).

Once introduced from the land model to the river and lake systems, particulate solids and nutrients (i.e., ISS, PON, PIP, and POP) are subject to either deposition or suspension based on a Rouse-number-dependent criterion, defined as settling velocity divided by the von Kármán constant and shear velocity (Pelletier, 2012).

$$\begin{array}{}\text{(6)}& R_{\mathit{\#}}={\displaystyle \frac{w_{\mathrm{S}}}{\mathit{\kappa }\cdot u_{\ast }}}={\displaystyle \frac{w_{\mathrm{S}}}{\mathit{\kappa }\sqrt{g\cdot z\cdot S}}},\end{array}$$

where R_# is the Rouse number (unitless), w_S is settling velocity (m s⁻¹), κ is the von Kármán constant (unitless), u_∗ is shear velocity (m s⁻¹), g is acceleration due to gravity (m s⁻²), and z is river or lake depth (m). If the Rouse number is less than 1.2 (a reported value in Pelletier, 2012), any newly introduced particulate matter, as well as that already in benthic sediments (i.e., SedN, SedP, and Sed), is suspended into the water column and subject to transport through the river network. Otherwise, the particulate matter is deposited to the benthic sediments. While organic material deposited to the sediments is assumed to be remineralized over time (i.e., no net long-term burial), PIP is subject to long-term burial in areas where resuspension does not occur.

Settling velocity (w_S) is estimated as a function of grain diameter, fluid viscosity, and fluid and solid density (Ferguson and Church, 2004).

$$\begin{array}{}\text{(7)}& w_{\mathrm{S}}={\displaystyle \frac{R\cdot g\cdot d^{\mathrm{2}}}{C_{\mathrm{2}}\cdot \mathit{\upsilon }+{\left(\mathrm{0.75}\cdot C_{\mathrm{3}}\cdot R\cdot g\cdot d^{\mathrm{3}}\right)}^{\mathrm{0.5}}}},\end{array}$$

where R is submerged specific gravity (unitless), d is particle diameter (m), υ is kinematic viscosity of the fluid (m² s⁻¹), and C₂ and C₃ are reported constants (unitless). For this initial FANSY implementation, a characteristic grain diameter (d) is assumed for all particulate material sinks.

For a batch river and lake system, settling and resuspension dynamics for ISS and Sed are written as

$$\begin{array}{}\text{(8)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{ISS}}=\left\{\begin{array}{ll}\frac{\mathrm{Sed}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ -\frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\frac{\mathrm{ISS}}{\mathrm{d}t}& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\text{(9)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{Sed}}=\left\{\begin{array}{ll}-\frac{\mathrm{Sed}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ \frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\frac{\mathrm{ISS}}{\mathrm{d}t}& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\end{array}$$

where ISS is inorganic suspended solid (kgD) and Sed is benthic sediment inorganic solid (kgD). ISS is gained by Sed resuspension and lost by deposition. The opposite holds for Sed. ISS deposition is modeled by implicitly solving for the ISS mass flux to Sed via w_S divided by a river or lake depth and multiplied by an ISS mass in the water column. This implicit scheme reduces the numerical burden and improves stability.

The conversion of PON to POM in freshwaters for the purpose of calculating SS to compare with observations is

$$\begin{array}{}\text{(10)}& {\displaystyle }\mathrm{POM}=r_{\mathrm{DN}}\cdot \mathrm{PON},\text{(11)}& {\displaystyle }\mathrm{SS}=\mathrm{ISS}+\mathrm{POM},\end{array}$$

where POM, PON, and SS are particulate organic matter (kgD), particulate organic N (kgN), and suspended solid (kgD), respectively, and r_DN is a POM-to-PON ratio in freshwaters (kgD kgN⁻¹).

2.2.2 Algae dynamics

Algae dynamics are governed by the balance of net growth (i.e., gross growth − respiration) and generalized mortality (i.e., non-predatory mortality + grazing + settling + excretion). The net growth is the difference between gross photosynthesis and respiration. The generalized mortality may include contributions from grazing, viruses, cell death, and excretion, though these diverse contributions are ultimately parameterized as a simple density-dependent loss term (Dunne et al., 2005). For a batch river and lake system, biogeochemical reactions for algae are written as

$$\begin{array}{}\text{(12)}& {\mathrm{SB}}_{{\mathrm{N}}_{\mathrm{al}}}=\mathit{\mu }\left(I_{\mathrm{av}},T,\mathrm{N},\mathrm{P}\right)\cdot {\mathrm{N}}_{\mathrm{al}}-m\left(T\right),\end{array}$$

where N_al is algae N (kgN), $\mathit{\mu }\left(I_{\mathrm{av}},T,\mathrm{N},\mathrm{P}\right)$ is algae net growth rate (s⁻¹) as a function of euphotic zone averaged irradiance I_av, temperature (T), N, and P, and m(T) is generalized algae mortality (kgN s⁻¹) as a function of T.

A dynamic regulatory model is adapted to predict a Chl-to-C ratio (r_ChlC) and net growth rate (μ) as a function of euphotic zone averaged irradiance, temperature, and nutrients (Geider et al., 1997). The μ is the difference between photosynthesis and respiration rates, as represented in Eq. (13). The r_ChlC is up- and down-regulated in accordance with light and nutrient conditions according to Eq. (14).

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathit{\mu }\left(I_{\mathrm{av}},T,\mathrm{N},\mathrm{P}\right)={\displaystyle \frac{P_m^{\mathrm{C}}}{\mathrm{1}+\mathit{\zeta }}}\cdot \left[\mathrm{1}-\exp \left({\displaystyle \frac{-{\mathit{\alpha }}^{\mathrm{Chl}}\cdot I_{\mathrm{av}}\cdot r_{\mathrm{ChlC}}}{P_m^{\mathrm{C}}}}\right)\right],\text{(14)}& {\displaystyle }{r}_{\mathrm{ChlC}}=max\left[r_{\mathrm{ChlC},\min },{\displaystyle \frac{r_{\mathrm{ChlC},\max }}{\mathrm{1}+\left(\frac{r_{\mathrm{ChlC},\max }\cdot {\mathit{\alpha }}^{\mathrm{Chl}}\cdot I_{\mathrm{av}}}{\mathrm{2}\cdot P_m^{\mathrm{C}}}\right)}}\right],\end{array}$$

where $P_m^{\mathrm{C}}$ is the C-specific light-saturated photosynthesis rate (s⁻¹), ζ is the cost of biosynthesis (unitless), α^Chl is the Chl-specific initial slope of the photosynthesis–light curve (gC m² gChl⁻¹ µmol photon⁻¹), I_av is euphotic zone averaged irradiance (µmol photon m⁻² s⁻¹), r_ChlC is the algae Chl-to-C ratio (kgChl kgC⁻¹), and r_ChlC,min and r_ChlC, max are minimum and maximum algae Chl-to-C ratios (kgChl kgC⁻¹).

The C-specific light-saturated photosynthesis rate ($P_m^{\mathrm{C}}$) is calculated as a function of temperature and nutrient limitation, also following the approach of Geider et al. (1997).

$$\begin{array}{}\text{(15)}& P_m^{\mathrm{C}}\left(\mathrm{T},\mathrm{N}\right)=P_{\max }^{\mathrm{C}}\left(T\right)\cdot min\left[\left({\lim }_{{\mathrm{NO}}_{\mathrm{3}}^{-}}+{\lim }_{{\mathrm{NH}}_{\mathrm{4}}^{+}}\right),{\lim }_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}\right],\end{array}$$

where $P_{\max }^{\mathrm{C}}\left(T\right)$ is temperature-dependent maximum photosynthesis rate (s⁻¹), and ${\lim }_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$, ${\lim }_{{\mathrm{NH}}_{\mathrm{4}}^{+}}$, and ${\lim }_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}$ are NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ limitations (unitless).

In LM3-FANSY, freshwater biogeochemical reaction rates approximately double for a temperature increase of 10 °C based on the Arrhenius equation, with a scaling factor θ (Chapra, 1997; Eppley, 1972). The simulated maximum and minimum water temperatures are limited to 30 and −3 °C, respectively.

$$\begin{array}{}\text{(16)}& P_{\max }^{\mathrm{C}}\left(T\right)=P_{\max }^{\mathrm{C}}\cdot {\mathit{\theta }}_{\mathrm{Al}}^{T-T_{\mathrm{ref}}},\end{array}$$

where T is temperature (°C), T_ref is reference temperature (°C), $P_{\max }^{\mathrm{C}}$ is the maximum photosynthesis rate at T_ref (s⁻¹), and θ_Al is the temperature correction factor for algae dynamics (unitless).

To combine the limiting effects of nutrients N and P, Liebig's law of the minimum is used. A NH${}_{\mathrm{4}}^{+}$ preference factor is used to account for inhibition of NO${}_{\mathrm{3}}^{-}$ uptake when NH${}_{\mathrm{4}}^{+}$ concentrations are high compared to an NH${}_{\mathrm{4}}^{+}$ half-saturation constant (Frost and Franzen, 1992). A saturating Monod relationship is used for handling the NH${}_{\mathrm{4}}^{+}$ and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ limiting effects. The maximum NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations are limited to 10 mol L⁻¹ to avoid numerical issues that can arise under extremely dry conditions.

$$\begin{array}{}\text{(17)}& {\displaystyle }{\lim }_{{\mathrm{NO}}_{\mathrm{3}}^{-}}={\displaystyle \frac{\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}{\left(k_{{\mathrm{NO}}_{\mathrm{3}}^{-}}+\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]\right)\cdot \left(\mathrm{1}+\frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{k_{{\mathrm{NH}}_{\mathrm{4}}^{+}}}\right)}},\text{(18)}& {\displaystyle }{\lim }_{{\mathrm{NH}}_{\mathrm{4}}^{+}}={\displaystyle \frac{\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{k_{{\mathrm{NH}}_{\mathrm{4}}^{+}}+\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}},\text{(19)}& {\displaystyle }{\lim }_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}={\displaystyle \frac{\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{k_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}+\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}},\end{array}$$

where $\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$, $\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]$, and $\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]$ are NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations (mgN L⁻¹ and mgP L⁻¹), and $k_{{\mathrm{NO}}_{\mathrm{3}}^{-}}$, $k_{{\mathrm{NH}}_{\mathrm{4}}^{+}}$, and $k_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}$ are NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ half-saturation constants for algae growth (mgN L⁻¹ and mgP L⁻¹).

Photosynthetically available, visible irradiance at the surface is used for algae growth dynamics. Light attenuation with depth is modeled by the Beer–Lambert law using an extinction coefficient (k_e, Chapra, 1997). The euphotic zone (depth where light intensity falls to 1 % of that at the surface) averaged light level (I_av) is used. The extinction coefficient is estimated dynamically to account for temporal and spatial variations in turbidity due to algae shading (Chapra, 1997; Riley, 1956), light extinction due to particle-free water and color, and variations in ISS and POM (Chapra, 1997; Di Toro, 1978).

$$\begin{array}{}\text{(20)}& {\displaystyle }{I}_z=I_{\mathrm{s}}\cdot e^{-k_{\mathrm{e}}\cdot z},\text{(21)}& {\displaystyle }{z}_{\mathrm{0.01}}=-{\displaystyle \frac{\ln \left(\mathrm{0.01}\right)}{k_{\mathrm{e}}}},\text{(22)}& {\displaystyle }{I}_{\mathrm{av}}=\left\{\begin{array}{ll}\frac{I_{\mathrm{s}}}{k_{\mathrm{e}}\cdot z}\left(\mathrm{1}-e^{-k_{\mathrm{e}}\cdot z}\right)& z\le z_{\mathrm{0.01}}\\ \frac{I_{\mathrm{s}}}{k_{\mathrm{e}}\cdot z_{\mathrm{0.01}}}\left(\mathrm{1}-e^{-k_{\mathrm{e}}\cdot z_{\mathrm{0.01}}}\right)& z>z_{\mathrm{0.01}}\end{array}\right.,\text{(23)}& {\displaystyle }\begin{array}{rl}{k}_{\mathrm{e}}& =k_{\mathrm{ew}}+s_{\mathrm{ISS}}\cdot \left[\mathrm{ISS}\right]+s_{\mathrm{POM}}\cdot \left[\mathrm{POM}\right]+s_{\mathrm{Chl}\mathrm{1}}\\ & \cdot \left[\mathrm{Chl}\right]+{\mathrm{s}}_{\mathrm{Chl}\mathrm{2}}\cdot [\mathrm{Chl}{]}^{\mathrm{2}/\mathrm{3}},\end{array}\end{array}$$

where I_z and I_s are irradiance at z and at the surface (µmol photon m⁻² s⁻¹), z_0.01 is river or lake depth where light intensity falls to 1 % of that at the surface (m), k_e is the light extinction coefficient (m⁻¹), k_ew is light extinction due to particle-free water and color (m⁻¹), s_Chl1 and s_Chl2 are algae self-shading factors (L µgChl⁻¹ m⁻¹ and L${}^{\mathrm{2}/\mathrm{3}}$ µgChl${}^{-\mathrm{2}/\mathrm{3}}$ m⁻¹), s_ISS and s_POM are constants accounting for the shading impacts of ISS and POM (L mgD⁻¹ m⁻¹), and [ISS], [POM], and [Chl] are ISS, POM, and Chl concentrations (mgD L⁻¹ and µgChl L⁻¹).

Biomass-specific algal mortality is assumed to increase nonlinearly with algae concentration, reflecting a presumed increase in predators with algal prey (e.g., Steele and Henderson, 1992).

$$\begin{array}{}\text{(24)}& m\left(T\right)=k_{\mathrm{m}}\left(T\right)\cdot {\mathrm{N}}_{\mathrm{al}}^{\mathrm{4}/\mathrm{3}},\end{array}$$

where k_m is temperature-dependent algae mortality rate (kgN${}^{-\mathrm{1}/\mathrm{3}}$ s⁻¹) reflecting the combined impacts of zooplankton grazing and other phytoplankton loss terms (e.g., viral-induced losses, cell death). Exponents between $\mathrm{4}/\mathrm{3}$ and 2 have been commonly applied in this relationship, with higher values corresponding to more tightly coupled top-down control (Dunne et al., 2005). We have adopted a value of $\mathrm{4}/\mathrm{3}$ to enable high biomass in nutrient-rich environments. The Arrhenius relationship with the same scaling as algae growth is applied to account for the effect of temperature on algae mortality (Eq. 16). The division of algal mortality between inorganic and organic and between dissolved and particulate nutrient pools is described in the following sections.

Algae P, C, and Chl are diagnosed from algae N using the Redfield $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio (Chapra, 1997; Redfield et al., 1963) and r_ChlC estimated above based on Geider et al. (1997).

$$\begin{array}{}\text{(25)}& {\displaystyle }{\mathrm{P}}_{\mathrm{al}}={\mathrm{N}}_{\mathrm{al}}\cdot r_{\mathrm{PN}},\text{(26)}& {\displaystyle }{\mathrm{C}}_{\mathrm{al}}={\mathrm{N}}_{\mathrm{al}}\cdot r_{\mathrm{CN}},\text{(27)}& {\displaystyle }\mathrm{Chl}={\mathrm{C}}_{\mathrm{al}}\cdot r_{\mathrm{ChlC}},\end{array}$$

where P_al, C_al, and Chl are algae P (kgP), C (kgC), and Chl (kgChl), and r_PN and r_CN are algae P-to-N (kgP kgN⁻¹) and C-to-N (kgC kgN⁻¹) ratios.

2.2.3 N dynamics

For a batch river and lake system, settling–resuspension dynamics and biogeochemical reactions are written for PON and SedN as

$$\begin{array}{}\text{(28)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{PON}}=\left\{\begin{array}{ll}{f}_{\mathrm{m},\mathrm{PON}}\cdot m\left(T\right)-k_{\mathrm{PON},\mathrm{d}}\left(T\right)& \\ \cdot \mathrm{PON}+\frac{\mathrm{SedN}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ f_{\mathrm{m},\mathrm{PON}}\cdot m\left(T\right)-k_{\mathrm{PON},\mathrm{d}}\left(T\right)& \\ \cdot \mathrm{PON}-\frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\frac{\mathrm{PON}}{\mathrm{d}t}& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\text{(29)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{SedN}}=\left\{\begin{array}{ll}-k_{\mathrm{SedN},\mathrm{d}}\left(T\right)\cdot \mathrm{SedN}& \\ -\frac{\mathrm{SedN}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ -k_{\mathrm{SedN},\mathrm{d}}\left(T\right)\cdot \mathrm{SedN}& \\ +\frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\frac{\mathrm{PON}}{\mathrm{d}t}& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\end{array}$$

where PON is particulate organic N (kgN), SedN is benthic sediment N (kgN), f_m,PON is the fraction of algae mortality which is deposited to the PON pool (unitless), and k_PON,d(T) and k_SedN,d(T) are temperature-dependent PON and SedN decomposition rates (s⁻¹).

In FANSY, PON is gained by algae mortality and benthic sediment N (i.e., SedN) resuspension and lost by deposition and decomposition. The same holds for SedN, except that it does not receive inputs from algae mortality. A fraction (f_m,PON) is adapted from Chapra et al. (2008) to represent the portion of algae mortality released as PON. First-order kinetics are used to describe various decay processes and transformations. PON and SedN are lost by decay processes that break down complex organic compounds into simpler organic N (i.e., DON) or into NH${}_{\mathrm{4}}^{+}$. Rate coefficients for these decay processes are thus much smaller than those for release of NH${}_{\mathrm{4}}^{+}$ due to DON decay (i.e., hydrolysis), oxidation of NH${}_{\mathrm{4}}^{+}$ to NO${}_{\mathrm{3}}^{-}$ (i.e., nitrification), and reduction of NO${}_{\mathrm{3}}^{-}$ to N₂ (i.e., denitrification) (Bowie et al., 1985; Chapra, 1997). The Arrhenius-based relationship (Eq. 16) is used to adjust the rate coefficients for temperature effects with different temperature correction factors of θ_Sed for sediment dynamics and θ for all dynamics except algae and sediment dynamics.

In FANSY, DON is gained by algae mortality and decomposition of PON and SedN and lost by hydrolysis. A fraction (f_m,DON) is adapted from Chapra et al. (2008) to represent the portion of algae mortality released as DON. Decomposition of PON and SedN releases both dissolved organic and inorganic N (i.e., DON and NH${}_{\mathrm{4}}^{+}$). Fractions (f_PON,DON and f_SedN,DON) are adapted from Bowie et al. (1985) to partition the fluxes to DON and NH${}_{\mathrm{4}}^{+}$.

$$\begin{array}{}\text{(30)}& \begin{array}{rl}{\mathrm{SB}}_{\mathrm{DON}}& =f_{\mathrm{m},\mathrm{DON}}\cdot m\left(T\right)+f_{\mathrm{PON},\mathrm{DON}}\cdot k_{\mathrm{PON},\mathrm{d}}\left(T\right)\\ & \cdot \mathrm{PON}+f_{\mathrm{SedN},\mathrm{DON}}\cdot k_{\mathrm{SedN},\mathrm{d}}\left(T\right)\cdot \mathrm{SedN}\\ & -k_{\mathrm{DON},\mathrm{d}}\left(T\right)\cdot \mathrm{DON},\end{array}\end{array}$$

where DON is dissolved organic N (kgN), f_m,DON is the fraction of algae mortality which is deposited to the DON pool (unitless), f_PON,DON and f_SedN,DON are fractions of PON and SedN decomposition fluxes which are deposited to the DON pool (unitless), and k_DON,d(T) is the temperature-dependent DON hydrolysis rate (s⁻¹).

In FANSY, NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ are removed by algae uptake during photosynthesis. NH${}_{\mathrm{4}}^{+}$ is returned to the water column through soluble excretions of algae (which is included in the generalized algae mortality term) and decomposition or hydrolysis of SedN, PON, or DON. Removal of NH${}_{\mathrm{4}}^{+}$ by nitrification generates NO${}_{\mathrm{3}}^{-}$, which is in turn lost by denitrification.

$$\begin{array}{}\text{(31)}& {\displaystyle }\begin{array}{rl}{\mathrm{SB}}_{{\mathrm{NH}}_{\mathrm{4}}^{+}}& =\left(\mathrm{1}-f_{\mathrm{m},\mathrm{PON}}-f_{\mathrm{m},\mathrm{DON}}\right)\cdot m\left(T\right)\\ & +\left(\mathrm{1}-f_{\mathrm{PON},\mathrm{DON}}\right)\cdot k_{\mathrm{PON},\mathrm{d}}\left(T\right)\cdot \mathrm{PON}\\ & +\left(\mathrm{1}-f_{\mathrm{SedN},\mathrm{DON}}\right)\cdot k_{\mathrm{SedN},\mathrm{d}}\left(T\right)\cdot \mathrm{SedN}\\ & +k_{\mathrm{DON},\mathrm{d}}\left(T\right)\cdot \mathrm{DON}-k_{\mathrm{nitr}}\left(T\right)\cdot {\mathrm{NH}}_{\mathrm{4}}^{+}\\ & -f_{{\mathrm{NH}}_{\mathrm{4}}^{+},\mathrm{up}}\mathit{\mu }\left(I_{\mathrm{av}},T,\mathrm{N},\mathrm{P}\right)\cdot {\mathrm{N}}_{\mathrm{al}},\end{array}\text{(32)}& {\displaystyle }\begin{array}{rl}{\mathrm{SB}}_{{\mathrm{NO}}_{\mathrm{3}}^{-}}& =k_{\mathrm{nitr}}\left(T\right)\cdot {\mathrm{NH}}_{\mathrm{4}}^{+}-k_{\mathrm{denitr}}\left(T\right)\cdot {\mathrm{NO}}_{\mathrm{3}}^{-}\\ & -(\mathrm{1}-f_{{\mathrm{NH}}_{\mathrm{4}}^{+},\mathrm{up}})\cdot \mathit{\mu }\left(I_{\mathrm{av}},T,\mathrm{N},\mathrm{P}\right)\cdot {\mathrm{N}}_{\mathrm{al}},\end{array}\text{(33)}& {\displaystyle }{f}_{{\mathrm{NH}}_{\mathrm{4}}^{+},\mathrm{up}}=\left({\displaystyle \frac{{\lim }_{{\mathrm{NH}}_{\mathrm{4}}^{+}}}{{\lim }_{{\mathrm{NO}}_{\mathrm{3}}^{-}}+{\lim }_{{\mathrm{NH}}_{\mathrm{4}}^{+}}}}\right),\end{array}$$

where NH${}_{\mathrm{4}}^{+}$ and NO${}_{\mathrm{3}}^{-}$ are ammonium and nitrate N (kgN), k_nit(T) and k_denit(T) are temperature-dependent nitrification and denitrification rates (s⁻¹), and $f_{{\mathrm{NH}}_{\mathrm{4}}^{+},\mathrm{up}}$ is the fraction of NH${}_{\mathrm{4}}^{+}$ uptake for algae growth (unitless).

2.2.4 P dynamics

Overall, P dynamics are similar to those of N, but with several differences. Because there are two suspended particulate forms of POP and PIP, SedP resuspension is divided into POP and PIP pools with a fraction of f_SedP,POP. Unlike N, P does not exist in a gaseous form, and FANSY includes no loss term for P to the atmosphere. Dissolved inorganic P sorbs strongly to solid particles. The exchange of PO${}_{\mathrm{4}}^{\mathrm{3}-}$ between the dissolved and particulate forms is modeled based on Freundlich kinetics (Garnier et al., 2005; Nemery, 2003), with the flux estimated as the disequilibrium between the two phases.

$$\begin{array}{}\text{(34)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{POP}}=\left\{\begin{array}{ll}{r}_{\mathrm{PN}}\cdot f_{\mathrm{m},\mathrm{PON}}\cdot m\left(T\right)-k_{\mathrm{POP},\mathrm{d}}\left(T\right)& \\ \cdot \mathrm{POP}+f_{\mathrm{SedP},\mathrm{POP}}\cdot \frac{\mathrm{SedP}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ r_{\mathrm{PN}}\cdot f_{\mathrm{m},\mathrm{PON}}\cdot m\left(T\right)-k_{\mathrm{POP},\mathrm{d}}\left(T\right)& \\ \cdot \mathrm{POP}-\frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\frac{\mathrm{POP}}{\mathrm{d}t}& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\text{(35)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{SedP}}=\left\{\begin{array}{ll}-k_{\mathrm{SedP},\mathrm{d}}\left(T\right)\cdot \mathrm{SedP}& \\ -\frac{\mathrm{SedP}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ -k_{\mathrm{SedP},\mathrm{d}}\left(T\right)\cdot \mathrm{SedP}& \\ +\frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\left(\frac{\mathrm{POP}}{\mathrm{d}t}+\frac{\mathrm{PIP}}{\mathrm{d}t}\right)& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\text{(36)}& {\displaystyle }{\mathrm{SB}}_{\mathrm{PIP}}=\left\{\begin{array}{ll}{F}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\text{\_to\_PIP}}+\left(\mathrm{1}-{\mathrm{f}}_{\mathrm{SedP},\mathrm{POP}}\right)\frac{\mathrm{SedP}}{\mathrm{d}t}& R_{\mathit{\#}}<\mathrm{1.2}\\ F_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\text{\_to\_PIP}}-\frac{w_{\mathrm{S}}}{z}{\left(\frac{\mathrm{1}}{\mathrm{d}t}+\frac{w_{\mathrm{S}}}{z}\right)}^{-\mathrm{1}}\frac{\mathrm{PIP}}{\mathrm{d}t}& R_{\mathit{\#}}\ge \mathrm{1.2}\end{array}\right.,\text{(37)}& {\displaystyle }\begin{array}{rl}{\mathrm{SB}}_{\mathrm{DOP}}& =r_{\mathrm{PN}}\cdot f_{\mathrm{m},\mathrm{DON}}\cdot m\left(T\right)+f_{\mathrm{POP},\mathrm{DOP}}\cdot k_{\mathrm{POP},\mathrm{d}}\left(T\right)\\ & \cdot \mathrm{POP}+f_{\mathrm{SedP},\mathrm{DOP}}\cdot k_{\mathrm{SedP},\mathrm{d}}\left(T\right)\cdot \mathrm{SedP}\\ & -k_{\mathrm{DOP},\mathrm{d}}\left(T\right)\cdot \mathrm{DOP},\end{array}\text{(38)}& {\displaystyle }\begin{array}{rl}{\mathrm{SB}}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}}& =r_{\mathrm{PN}}\cdot \left(\mathrm{1}-f_{\mathrm{m},\mathrm{PON}}-f_{\mathrm{m},\mathrm{DON}}\right)\cdot m\left(T\right)\\ & +\left(\mathrm{1}-f_{\mathrm{POP},\mathrm{DOP}}\right)\cdot k_{\mathrm{POP},\mathrm{d}}\left(T\right)\cdot \mathrm{POP}\\ & +\left(\mathrm{1}-f_{\mathrm{SedP},\mathrm{DOP}}\right)\cdot k_{\mathrm{SedP},\mathrm{d}}\left(T\right)\cdot \mathrm{SedP}\\ & +k_{\mathrm{DOP},\mathrm{d}}\left(T\right)\cdot \mathrm{DOP}-r_{\mathrm{PN}}\cdot \mathit{\mu }\left(I_{\mathrm{av}},T,\mathrm{N},\mathrm{P}\right)\\ & \cdot {\mathrm{N}}_{\mathrm{al}}-F_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\text{\_to\_PIP}},\end{array}\text{(39)}& {\displaystyle }\left[{\mathrm{PIP}}_{\mathrm{eq}}\right]=a\cdot [{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}{]}^{\mathrm{b}}\cdot [\mathrm{ISS}{]}_{\mathrm{2}},\text{(40)}& {\displaystyle }{F}_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\text{\_to\_PIP}}=\left(\right[{\mathrm{PIP}}_{\mathrm{eq}}]-[\mathrm{PIP}\left]\right)/\mathrm{d}t\cdot {\mathrm{H}}_{\mathrm{2}}\mathrm{O}\cdot {\mathrm{10}}^{-\mathrm{3}},\end{array}$$

where POP is particulate organic P (kgP), SedP is benthic sediment P (kgP), PIP is particulate inorganic P (kgP), DOP is dissolved organic P (kgP), PO${}_{\mathrm{4}}^{\mathrm{3}-}$ is phosphate P (kgP), H₂O is water volume (m³), k_POP,d(T) and k_SedP,d(T) are temperature-dependent POP and SedP decomposition rates (s⁻¹), k_DOP,d(T) is the temperature-dependent DOP hydrolysis rate (s⁻¹), f_SedP,POP is the fraction of SedP resuspension which is deposited to the POP pool (unitless), f_POP,DOP and f_SedP,DOP are fractions of POP and SedP decomposition fluxes which are deposited to the DOP pool (unitless), [PIP] and [PO${}_{\mathrm{4}}^{\mathrm{3}-}$] are PIP and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ concentrations (mgP L⁻¹), [ISS]₂ is ISS concentration (gD L⁻¹) (notice the concentration unit difference from [ISS] in Eq. 23), [PIP_eq] is PIP equilibrium concentration (mgP L⁻¹), $F_{{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\text{\_to\_PIP}}$ represents fluxes from PO${}_{\mathrm{4}}^{\mathrm{3}-}$ to PIP (kgP), and a (mgP gD⁻¹) and b (unitless) are reported empirical kinetic parameters.

2.2.5 Model calibration

Because many of the reported parameters required to simulate the coupled freshwater algae, nutrient, and solid dynamics within LM3-FANSY vary widely, it is difficult to assign a single global value for each parameter. Informed by parameter sensitivity analysis herein (Sect. 4.4), our approach was to coarsely calibrate a limited set of uncertain yet highly influential parameters within their broad observed ranges to reduce errors in simulated SS, N, and P loads in different forms across a globally distributed subset of large rivers.

First, as described above, the terrestrial soil erosion parameter C₁ was calibrated to reduce prediction errors of SS loads. We emphasize that this parameter is expected to be resolution-dependent.

Second, the generalized algal mortality constant k_m was tuned to produce reasonable chlorophyll a concentrations in globally distributed lakes, acknowledging the limitation of the present lake biogeochemistry in LM3-FANSY (see Sect. 4.5 for further discussion). We adapted an approach of Chapra et al. (2008) to partition nutrient fluxes from algae mortality to different pools of nutrients in different forms (i.e., particulate organic, dissolved organic, and inorganic) based on fixed fractions. Uncertainties in these fixed factions due to the lack of theoretical and empirical evidence are investigated in the sensitivity analysis (Sect. 4.4).

Third, we find that the rate coefficients for hydrolysis, nitrification, and denitrification are highly influential parameters for determining river nutrient loads in different forms relative to those for decay processes that break down complex organic compounds into simpler organic and inorganic compounds. These highly influential parameters have been calibrated to reduce prediction errors of nutrient loads in different forms. We adapted an approach of Bowie et al. (1985) to partition fluxes from complex organic nutrient decomposition to simpler organic vs. inorganic nutrient pools based on a uniform fraction. Uncertainties in the faction are investigated in the sensitivity analysis (Sect. 4.4).

3 Model forcing and simulations

3.1 Baseline simulations

LM3-FANSY was implemented globally at 1° spatial and 30 min temporal resolution with all inputs regridded to 1° resolution. Following ∼ 11 000 years of spin-up from M. Lee et al. (2019), the terrestrial component LM3 was run for the period 1700–1899 by recycling 30 years (1948–1977) of observation-based, historical climate forcing (Sheffield et al., 2006a) and Coupled Model Intercomparison Project (CMIP6) datasets for atmospheric CO₂ (Meinshausen et al., 2017a), atmospheric N deposition (Durack et al., 2024), and land use states and transitions (Hurtt et al., 2020a). Since the freshwater component requires a shorter time for equilibrium than vegetation and soil, the merged terrestrial and freshwater components for LM3-FANSY were run for only the 1900–2000 period using additional CMIP6 datasets for fertilizer N applications (Hurtt et al., 2020a) and reported N and P inputs to rivers (Beusen et al., 2015).

The observation-based, historical climate forcing data available for the period 1948–2010 (Sheffield et al., 2006a) include precipitation, specific humidity, air temperature, surface pressure, wind speed, and shortwave and longwave downward radiation at 1° and 3 h resolution. This forcing was cycled over a period of 30 years (1948–1977) to perform long-term simulations from 1700 to 1947, and the 1948–2000 forcing data were used for the simulations from 1948 to 2000. Annual atmospheric CO₂ estimates (Meinshausen et al., 2017a) available for the period 1–2500 are used for the corresponding period simulations from 1700 to 2000. The atmospheric N deposition data (Durack et al., 2024) include oxidized and reduced N (NO_y and NH_x) at 2.5° longitude by 1.9° latitude and 1-month resolution for the period 1850–2099. The NO_y and NH_x deposition for the year 1850 was applied to soil NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ pools, respectively, for the 1700–1849 simulation, and then the 1850–2000 deposition was applied for the 1850–2000 simulations.

The dataset of land use states and transitions as well as fertilizer N applications at 0.25° and 1-year resolution (Hurtt et al., 2020a) is available for the period 850–2100. The 1700–2000 land use state and transition data were used for the simulations from 1700 to 2000. Since the amount of fertilizer application in the dataset is zero until 1915, the 1916–2010 fertilizer N was applied for the simulations from 1916 to 2000. For land use and fertilizer applications, 12 land use types reported in the Land Use Harmonization (LUH2) (Hurtt et al., 2020a) were grouped into 4 types in LM3-FANSY: (1) primary land in LM3-FANSY is the sum of forested primary land and non-forested primary land in LUH2; (2) secondary land in LM3-FANSY is the sum of potentially forested secondary land, potentially non-forested secondary land, and urban land in LUH2; (3) cropland in LM3-FANSY is the sum of C₃ annual cropland, C₃ perennial cropland, C₄ annual cropland, C₄ perennial cropland, and C₃ N-fixing cropland in LUH2; and (4) pasture in LM3-FANSY is the sum of managed pasture and rangeland in LUH2. The sum of fertilizers allocated to the five croplands in LUH2 was applied to the cropland in LM3-FANSY.

The terrestrial soil erosion component requires slope, rainfall, and LAI as inputs. Rainfall was simulated by using 3-hourly precipitation and temperature from Sheffield et al. (2006a) and assuming that all of the precipitation falls as snow with the temperature less than 0 °C; otherwise, it is assumed to be rain. The slope input was derived from Danielson and Gesch (2011). The LAI input was from an observationally derived, monthly average global vegetation LAI dataset from the Global Inventory Modeling and Mapping Studies (GIMMS) Normalized Difference Vegetation Index (NDVI3g) for the period 1982–2010 (Zhu et al., 2013) to avoid potential errors that might be caused by using modeled LAI. The 1982 LAI was used to perform long-term simulations from 1900 to 1981 and the 1982–2000 LAI was used for the simulations from 1982 to 2000.

Solid and nutrient inputs from terrestrial systems and from the atmosphere to rivers are either simulated by LM3-FANSY or provided by Beusen et al. (2015) (Table 3). For solids, all inputs were simulated by LM3-FANSY. For N, all inputs were simulated by LM3-FANSY except aquaculture, wastewater, and atmospheric deposition, which were provided by Beusen et al. (2015). For P, which is not currently included in LM3, all inputs were provided by Beusen et al. (2015).

Table 3 Description of solid and nutrient inputs from terrestrial systems and the atmosphere to rivers, which were either simulated by LM3-FANSY or provided by Beusen et al. (2015). The numbers are the fractions of dividing TN and TP inputs provided by Beusen et al. (2015) into different N and P species. Organic N (ON) and organic P (OP) from Beusen et al. (2015) are considered to be mainly (70 %) particulate.

Beusen et al. (2015) provided 5-year interval data for the period 1900–2000 at 0.5° resolution. The data were regridded to our 1° resolution by summing up the values given in kilograms per year (kg yr⁻¹) and linearly interpolated across the 5-year intervals. The Beusen et al. (2015) wastewater N and P inputs were from the Morée et al. (2013) urban waste N and P discharge estimates to surface waters. Beusen et al. (2015) calculated aquaculture N and P inputs using the Bouwman et al. (2013b) finfish and the Bouwman et al. (2011) shellfish data. For atmospheric N deposition inputs, the Beusen et al. (2015) input for the year 2000 was from the Dentener et al. (2006) ensemble of reactive transport models and those for the years before 2000 were made by scaling the deposition with the Bouwman et al. (2013a) ammonia emissions. The Beusen et al. (2015) surface runoff P inputs included those leached from soil P budgets (i.e., the sum of fertilizer and animal manure minus crop and grass withdrawal) and those driven by soil erosion estimates based on Cerdan et al. (2010). The Beusen et al. (2015) P inputs of litter from floodplains were estimated as 50 % of total NPP with a C:P ratio of 1200. The Beusen et al. (2015) weathering P inputs were computed based on the Hartmann et al. (2014) chemical weathering P release estimates.

We divided yearly total N (TN) and total P (TP) inputs from Beusen et al. (2015) into different N and P forms (listed in Table 3) based on Vilmin et al. (2018). Vilmin et al. (2018) suggested fractions to divide TN and TP inputs into three N and P species according to the source. For sewage, fractions were given for three types (i.e., untreated, primary treated, and secondary/tertiary treated). The sewage inputs from Beusen et al. (2015), however, were aggregated. Thus, to divide the Beusen et al. (2015) sewage TN and TP inputs into three N and P species, we have taken a middle ground and used the fractions for primary treated sewage. Although we acknowledge that this is a simplification, we find that our results are relatively insensitive to alternatives, assuming that all sewage is untreated, all sewage has secondary treatment, and the two options are based on the fact that over 80 % of wastewater is discharged without “adequate treatment” (Environment and Natural Resources Department, 2022, Table 3). Specifically, the fractions are driven by assuming that (1) 80 %, 10 %, and 10 % of the Beusen et al. (2015) sewage is untreated, primary treated, and secondary/tertiary treated, respectively, and (2) 40 %, 40 %, and 20 % of the Beusen et al. (2015) sewage is untreated, primary treated, and secondary/tertiary treated, respectively. In all cases, the simulated river loads of each species change by ≤9 %, and the simulated total load does not change (see Sect. 4.4).

For TP fluxes from agricultural lands to rivers, distinct species fractions were given for two sources (i.e., surficial runoff and soil loss), while surface runoff TP inputs from Beusen et al. (2015) were aggregated. To divide the Beusen et al. (2015) agricultural surface runoff TP inputs into three P species, we assume nearly equal fractions. As is the case for sewage, perturbation experiments show that our results are relatively insensitive to a wide range of fractions (see Table 3 and the sensitivity analysis in Sect. 4.4). The six uncertainty simulations use (1) the fractions for surficial runoff for all agricultural fluxes, (2) the fractions for soil loss for all agricultural fluxes, and fractions driven by assuming that (3) 40 % and 60 % of the Beusen et al. (2015) agricultural surface runoff TP inputs are from surficial runoff and soil loss, respectively, (4) 60 % and 40 % of those are from surficial runoff and soil loss, respectively, (5) 20 % and 80 % of those are from surficial runoff and soil loss, respectively, and (6) 80 % and 20 % of those are from surficial runoff and soil loss, respectively. In all cases, the simulated river loads of each species change by ≤13 %, and the simulated total load changes by ≤1 %.

For aquaculture, two groups of fractions are given for particulate and dissolved sources. The fractions in Table 3 are driven by assuming that about 12 % and 31 % of aquaculture TN and TP inputs from Beusen et al. (2015) are particulate, approximating Fig. 3 of Bouwman et al (2013a). We did not conduct sensitivity simulations for aquaculture because aquaculture nutrient inputs are very small compared to the other sources (<1 % of total inputs for both N and P), and thus uncertainties associated with the fractions should be negligible for the global application herein.

3.2 Sensitivity simulations

Model sensitivities to the inputs, components, and parameters are examined by analyzing the responses of river solid and nutrient loads to their changes for the period 1948–2000. Model sensitivities to the nutrient inputs to rivers are examined by increasing each input source by 15 % or removing it. When the inputs were from LM3, the 15 % increase was applied to each input variable (i.e., NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, DON, and PON). When the inputs were externally prescribed, the 15 % increase was applied to each input source listed in Table 3. In both cases, the increases were applied uniformly over the grid and time.

One of the distinct features of LM3-FANSY is the capability to model interactions between algae, nutrient, and solid dynamics. Light shading by SS and algae themselves modulates the strength of algal productivity and, in turn, river solids and nutrients. In LM3-FANSY, the light extinction coefficient is dynamically simulated as a function of ISS, POM, and Chl (Eq. 23) instead of using a prescribed parameter. To evaluate how critical the dynamic light extinction component is for model predictive capacity, the component was replaced with a prescribed parameter value (k_e=0.15 or 0.45) and the river load responses to more active algal populations are examined.

Model sensitivities to the parameters are examined by decreasing each of the calibrated parameters listed in Table 2 by half or doubling them. For the parameters which have much smaller or much larger impacts compared to the other parameters, additional sensitivity tests have been performed to show responses of river loads to the parameter changes in their broad observed ranges.

3.3 Comparisons of measurement-based and modeled estimates

For cross-watershed evaluations, we compare LM3-FANSY results of river SS, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, DIN (the sum of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$), DON, total Kjeldahl N (TKN, the sum of NH${}_{\mathrm{4}}^{+}$, DON, and PON), PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DOP, and TP (the sum of PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DOP, PIP, and POP) yields (kg km⁻² yr⁻¹), loads (kt yr⁻¹), and concentrations (mg L⁻¹) with measurement-based estimates from 70 of the world's major rivers (Tables S1–S9, the GEMS-GLORI world river discharge database, Meybeck and Ragu, 2012). The 70 river basins cover 55 % of global land area (excluding the Antarctic) and are distributed globally across various climates and land use (Fig. S1). The LM3-FANSY performance is also compared with that of Global NEWS (Mayorga et al., 2010) by using the same measurement-based estimates, yet excluding a few unavailable rivers in Global NEWS (Tables S1, S4, S6–S9).

The 70 rivers were chosen by the following procedure. First, river basins with areas <100 000 km², about 10 grid cells in our 1° resolution, were excluded from the comparisons. Second, river locations were identified by matching latitudes, longitudes, and basin areas of the modeled and reported rivers, which are located either at the river mouths or near the river mouths. Third, rivers that were not properly represented within the LM3-FANSY river network were excluded. For example, the GEMS-GLORI database provides a Sanaga River SS concentration at 3.8° N and 10.1° E with its basin area of 119 300 km². LM3-FANSY, however, does not capture the Sanaga River at a comparable location (3.5° N and 10.5° E), where the river has a much smaller basin area (5607 km²). The Sanaga River was thus excluded. In total, nine rivers (i.e., Anabar, Brantas, Burdekin, Don, Hayes, Huai, Pyasina, Sanaga, and Sepik) were excluded. When a river in the LM3-FANSY river network captures two merged small rivers in GEMS-GLORI, the water-discharge-weighted mean concentration of the two rivers was used for analysis. When more than one data source was given for a river, the data selected for the first line were used, since they were considered to be “the most reliable and generally were obtained first-hand by local engineers or scientists (Meybeck and Ragu, 1997)”. When data in the first line were outdated (1970s–early 1980s) or reported as zero, the latest data were used if available in the next lines. When more than one data source was given for a river with different basin areas, the data monitored at the location with the largest basin area (i.e., nearest to the river mouth) were used.

Since hydraulic controls like damming, irrigation, and diversion affect many rivers, natural river water discharges are distinguished from actual, modified ones in the GEMS-GLORI database. LM3-FANSY does not resolve such hydraulic controls, and thus, if available, the natural discharges of GEMS-GLORI are used when calculating loads and yields from the GEMS-GLORI concentrations. Comparisons of the model results with loads and yields calculated by using the actual discharges are also presented in the Supplement. Since Global NEWS accounts for anthropogenic hydraulic controls, Global NEWS results are compared with the loads and yields calculated by using the actual discharges.

We report the Pearson correlation coefficient (r) and Nash–Sutcliffe model efficiency coefficient (NSE) between the log-transformed modeled and measurement-based estimates across the 70 rivers. We also report the prediction error computed as the difference between the modeled and measurement-based estimates of loads expressed as a percentage of the measurement-based load. For global evaluations, LM3-FANSY results are compared with reported global estimates from various references. For cross-watershed and global evaluations, annual results for the year 1990 are analyzed and presented. In parentheses, ranges of using annual results for the years 1990–2000 are also provided throughout the paper.

For time series evaluations, we used reported annual solid and nutrient loads across eight stations in large US rivers from the USGS National Water Quality Network (NWQN, Lee, 2022). We note that the robustness of evaluating simulated interannual variability against simple flow-weighted annual observations depends on the frequency and timing of chemical samplings (Lee et al., 2021). The reported annual loads from NWQN were estimated with the USGS's latest load estimation method, Weighted Regressions on Time, Discharge, and Season method with Kalman filtering (WRTDS-K), for all analyzed rivers, except for the Columbia River for which only the regression (REG) method estimates are available (Lee et al., 2017). The WRTDS-K method proved to be the most accurate annual load estimation method among methods studied by C. J. Lee et al. (2019).

The eight stations are located in the three largest river basins in NWQN (six stations in the Mississippi River Basin, one station in the St. Lawrence River Basin, and one station in the Columbia River Basin). The basin area of the Yukon River station is larger than those of the St. Lawrence and Columbia River stations, but the Yukon River station was excluded for this analysis because data are only available beyond 2000. The corresponding station locations in the LM3-FANSY river network were identified by matching latitudes, longitudes, and basin areas of the reported and modeled rivers. The yearly data provided by NWQN are based on “a water year” defined as the 12-month period from 1 October for any given year through 30 September of the following year. The corresponding water yearly results of LM3-FANSY for the period ∼ 1963–2000 were used for this analysis.

4 Results and discussion

4.1 Model performance analysis

Measurement-based and simulated annual SS estimates across 65 rivers are significantly correlated, with Pearson correlation coefficient (r) values equal to 0.65 (0.57–0.65) for yields, 0.76 (0.71–76) for loads, and 0.67 (0.67–0.69) for concentrations for the year 1990 (range for the years 1990–2000) (Fig. 2, Table 4). This result, which allows for a coarsely calibrated value of the one free terrestrial soil erosion parameter (C₁=0.012), demonstrates that LM3-FANSY reproduces the measurement-based SS estimates fairly well, especially given that the model contains only one calibrated parameter for SS. This model performance is competitive with other global model estimates using larger numbers of free parameters (Hatono and Yoshimura, 2020; Mayorga et al., 2010; Tan et al., 2017). For example, model performance of Global NEWS 2, when analyzed on the same dataset used for our model performance evaluation yet excluding a few unavailable rivers (Table S1), is slightly better than LM3-FANSY for yields and loads and slightly worse for concentrations based on correlations (Fig. S2). The total quantity of global river SS loads to the coastal ocean is 10 (10–11) Pg yr⁻¹ for the year 1990 (range for the years 1990–2010) for LM3-FANSY, at the lower bound of previous estimates (Table 5, Global NEWS estimate of 11–27 Pg yr⁻¹, Beusen et al., 2005; Discharge Relief Temperature sediment delivery model (QRT) estimate of 12.6 Pg yr⁻¹, Syvitski et al., 2005).

Figure 2 Pearson correlation coefficients (r) and p values (p) for the log-transformed measurement-based vs. simulated SS yields, loads, and concentrations across 65 rivers for the year 1990.

Table 4 Pearson correlation coefficients (r) and Nash–Sutcliffe model efficiency coefficients (NSE) for the log-transformed measurement-based vs. simulated loads across major world rivers for the year 1990 (range for the years 1990–2000 in parentheses). The prediction error is computed as the difference between the simulated and measurement-based estimates of loads expressed as a percentage of the measurement-based load.

Table 5 LM3-FANSY and published estimates of global river loads to the coastal ocean in Pg yr⁻¹ for SS and Tg yr⁻¹ for nutrients. LM3-FANSY results are for the year 1990 (range for the years 1990–2000).

Correlations between measurement-based and simulated NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, DON, and TKN yields, loads, and concentrations across 51, 37, 18, and 12 rivers, respectively (Fig. 3, Table 4), indicate that LM3-FANSY can also explain the observed spatial variations in river N in multiple forms and units to a reasonable extent. The fidelity of LM3-FANSY in terms of N is comparable to that of Global NEWS 2 (which does not estimate NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$ separately but only estimates their sum as DIN). Spatial DIN patterns evaluated by r values are better represented by Global NEWS 2, while LM3-FANSY estimates better spatial DON patterns (Tables S4, S9, Fig. S2). Measurement-based estimates of particulate nutrient compounds to evaluate our model implemented at 1° resolution herein are limited. The evaluation of modeled PON is limited to a few measurement-based TKN estimates that include PON, but these aggregate values match the model estimates reasonably well.

Figure 3 Pearson correlation coefficients (r) and p values (p) for the log-transformed measurement-based vs. simulated NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$, DON, and TKN yields, loads, and concentrations across 51, 37, 18, and 12 rivers for the year 1990.

Globally, TN inputs to rivers (N inputs hereafter) in LM3-FANSY are 85 (85–91) TgN yr⁻¹, about 59 (56–60) % of which is lost to the atmosphere or retained within freshwaters (N loss or retention hereafter) and the other 41 (40–43) % is exported to the coastal ocean (N exports hereafter) (Table 6). LM3-FANSY does not include net long-term N burial to bottom sediments, as all organic N delivered to the sediments is ultimately remineralized. While year-to-year sediment N inventories may vary, effectively all long-term N losses are to the atmosphere via freshwater denitrification. LM3-FANSY estimates that N loss or retention is 144 (130–148) % of the N exports, consistent with the 147 % and 143 % of Galloway et al. (2004) and Seitzinger et al. (2006) yet larger than the 73 % estimated by IMAGE-GNM (Beusen et al., 2016). LM3-FANSY estimates that 59 (56–60) % of the N inputs is lost or retained in freshwaters, consistent with the 60 % of Galloway et al. (2004) yet larger than the 42 % of IMAGE-GNM (Beusen et al., 2015). Recent estimates of global river TN loads to the coastal ocean vary widely, ranging from about 36.5 to 48 TgN yr⁻¹ (Table 5, Beusen et al., 2016; Boyer et al., 2006; Galloway et al., 2004; Green et al., 2004; Mayorga et al., 2010). Our global estimate of 35 (34–38) TgN yr⁻¹ for the year 1990 (range for the years 1990–2000) is thus at the lower bound of the published range. The simulated global river TN loads contain approximately equal contributions by DIN (the sum of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$, 14 (14–16) TgN yr⁻¹, 41 % of TN) and DON (13 (13–14) TgN yr⁻¹, 39 % of TN), with a lesser contribution by PON (7 (7–8) TgN yr⁻¹, 20 % of TN). The estimates of global river DIN loads are at the lower end of recent estimates, which range from 14.5 to 18.9 TgN yr⁻¹ (Mayorga et al., 2010; Green et al., 2004; Smith et al., 2003). This may be partly due to an overestimation of freshwater denitrification and/or algae-mediated transformations to organic forms. In contrast, our global river DON load estimate is slightly higher than a previous estimate of 10 TgN yr⁻¹ (Harrison et al., 2005). Our global river PON load estimate is considerably lower than a previous estimate of 13.5 TgN yr⁻¹ (Mayorga et al., 2010). See Sect. 4.5 for further discussion.

Table 6 LM3-FANSY and published estimates of global freshwater N and P budgets.

^a “Nr input to rivers” minus “Nr export to coastal areas” in Table 1 of Galloway et al. (2004). ^b The sum of denitrification in “rivers” and “lakes and reservoirs” in Table 1 of Seitzinger et al. (2006). ^c Inputs of P to dam reservoirs in Table 1 of Maavara et al. (2015). ^d Retention of P by dam reservoirs in Table 1 of Maavara et al. (2015).

Simulated river PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DOP, and TP yields, loads, and concentrations are in reasonable agreement with more limited measurement-based estimates across 47, 9, and 5 rivers, respectively (Fig. 4, Table 4). Global NEWS 2 has generally higher correlations for yields and loads and lower correlations for concentrations for the three species compared to LM3-FANSY (Tables S6–S8, Fig. S2).

Figure 4 Pearson correlation coefficients (r) and p values (p) for the log-transformed measurement-based vs. simulated PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DOP, and TP yields, loads, and concentrations across 47, 9, and 5 rivers for the year 1990.

Globally, TP inputs to rivers in LM3-FANSY are 8 (8–9) TgP yr⁻¹, about 9 (6–10) % of which is stored within freshwaters and the other 91 (90–94) % is exported to the coastal ocean (Table 6). IMAGE-GNM estimates that ∼ 56 % (5 of 9 TgP yr⁻¹) of the P inputs are stored within freshwaters (Beusen et al., 2016). This is a large difference from our estimate of 9 (6–10) % retention, but the difference is around a very uncertain number as the storage has not been directly measured. LM3-FANSY, which does not account for dams and reservoirs, likely underestimates global freshwater P retention by at least ∼ 12 % (Table 6, Maavara et al., 2015, see Sect. 4.5 for further discussion). The overall consistency of our SS, N, and P estimates with the observed cross-watershed constraints (Figs. 2–4, S6, Table 4), however, suggests that the bias introduced by the lack of dams and reservoirs may not be large. In contrast, underestimates of P exports to the coastal ocean from high-exporting basins such as the Amazon, Ganges, and Yangtze shown in Fig. 2 of Harrison et al. (2019) imply that IMAGE-GNM likely overestimates global freshwater P retention. Even though our freshwater P retention estimates are near the lower bound, our P retention estimates are far higher than those for N, mainly reflecting the sorption of PO${}_{\mathrm{4}}^{\mathrm{3}-}$ onto solids and its deposition to bottom sediments. Our estimate of global river TP loads to the coastal ocean (7, 7–8 TgP yr⁻¹) falls within the range of other estimates (9.04 TgP yr⁻¹ from Global NEWS 2, Mayorga et al., 2010; 4 TgP yr⁻¹ from IMAGE-GNM, Beusen et al., 2016, Table 5). LM3-FANSY estimates that, globally, rivers export 5 TgP yr⁻¹ as PP (64 % of TP), 2 TgP yr⁻¹ as PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (25 %), and 1 TgP yr⁻¹ as DOP (11 %). The global river PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DOP, and PP load estimates are consistent with previous estimates of 1.45–2.3 TgP yr⁻¹ for PO${}_{\mathrm{4}}^{\mathrm{3}-}$ (Harrison et al., 2010; Smith et al., 2003), 0.6 TgP yr⁻¹ for DOP (Harrison et al., 2005), and 6.6 TgP yr⁻¹ for PP (Mayorga et al., 2010).

Global watershed model performance of simulating N:P ratios has not been reported in prior publications. While our simulations are reasonably successful in capturing cross-ecosystem differences in both N and P species, variations in their ratios proved more challenging. The model captures the mean ratio of DIN and DIP but little of the variation (Fig. S3). One factor that likely contributes to this is the inconsistency between the N inputs to rivers, the majority of which (∼ 92 %) was simulated within LM3, and the P inputs, all of which were drawn from another model (IMAGE-GNM). Continued refinement is thus needed to reliably capture N:P ratio variations in rivers and their subsequent water quality implications (see Sect. 4.5).

Despite the relatively simple nature of lake biogeochemistry in LM3-FANSY (i.e., vertically unresolved mixed reactors, M. Lee et al., 2019), the model creates a reasonable range of chlorophyll a concentrations (Fig. S4) that generally fall within a range of in situ estimates from globally distributed lakes (Sayers et al., 2015). The in situ estimates in the compilation of Sayers et al. (2015) range from ∼ 0 to ∼ 100 mg m⁻³, mostly falling between 5 and 50 mg m⁻³ (Fig. 5 of Sayers et al., 2015, available at https://www.tandfonline.com/doi/full/10.1080/01431161.2015.1029099, last access: 20 June 2024).

Despite the significant correlation between the measurement-based and modeled estimates for each solid and nutrient form across various rivers, errors on a basin-by-basin scale are substantial, with high-load, large basins tending to have large absolute errors (Figs. 2–4). However, the ranges of prediction errors in our model simulation, as demonstrated in the interquartile range (IQR) and distribution of prediction errors (Table 4), are comparable to those of other models (Dumont et al., 2005; Harrison et al., 2005, 2010). This suggests that our model has a competitive correlation (r value), precision (IQR), and bias (median error) for each species compared to previous efforts even while including fewer observational constraints on the river sources and more explicitly parameterized freshwater biogeochemical processes.

4.2 Spatial pattern analysis

Spatial maps of river solid and nutrient yields and loads help identify global hotspots of water pollution and provide insight into which processes modulate the magnitudes and form of inputs. A global map of simulated terrestrial soil erosion rate from Eq. (3) is more strongly related to the basin slope map than to the rainfall or LAI maps, reflecting the prominent role of topographic steepness in controlling soil erosion (Fig. 5). This is consistent with previous studies (Pelletier, 2012; Syvitski et al., 2005). The eroded soil is transported as suspended load to rivers, some of which is stored within rivers and lakes, and the rest makes its way to large river outlets to the coastal ocean. Simulated river SS yields are high in mountains like the Andes, Rockies, and Himalayas and low in most gently sloping areas. The yields (i.e., loads per area) decrease with distance from mountains, as some of the soil is stored in lowland rivers and lakes and as basin areas (the denominator in yields) increase downstream. In contrast, simulated river SS loads tend to increase downstream because larger rivers carry more soils from many small streams and tributaries. The Ganges, Chang Jiang, Indus, and Huang He rivers in Asia, the Parana and Amazon rivers in South America, and the Mississippi and Columbia rivers in North America are among the largest river SS exporters (i.e., highest loads) in LM3-FANSY.

Figure 5 Global maps of the model inputs of LAI, rainfall, and basin slope and of the simulated soil erosion rate, river SS yields, and river SS loads for the year 1990.

The Mississippi, Chang Jiang, Ganges, Ob, Amazon, Parana, and Orinoco rivers are among the top exporters of all three N forms (DIN, DON, and PON) to the coastal ocean in LM3-FANSY (Fig. 6). These basins are characterized by tropical humid climates with high terrestrial productivity, high population and agricultural pressures, and/or high river water discharge. The highest river DIN yields and loads occurred in European and Asian rivers (e.g., Rhine, Elbe, Danube, and Zhujiang), despite their relatively low river water discharge and small basin areas, largely due to substantial anthropogenic N inputs (Dumont et al., 2005; Mayorga et al., 2010). In contrast, the lowest river DIN yields and loads are estimated for arid regions and most high-latitude basins with low population densities and less intensive agriculture. The Amazon, Parana, Orinoco, Zaire, Ganges, Zhujiang, Hong, Chang Jiang, Mississippi, Yukon, Ob, and Yenisey rivers are estimated to produce the largest river DON yields and loads. The largest river DON yields and loads are from humid tropical regions, despite lower human pressures, indicating a critical role of non-anthropogenic sources (i.e., terrestrial soil and litter fluxes from N-enriched natural forests) in exporting the dissolved organic form (Harrison et al., 2005). Low river DON yields and loads tend to occur in relatively dry regions with low anthropogenic pressures.

Figure 6 Global maps of the simulated river DIN, DON, and PON yields and loads for the year 1990.

The Mississippi, Chang Jiang, Ganges, Amazon, and Danube rivers are among the highest exporters of all three P forms (DIP, DOP, and PP – the sum of POP and PIP) to the coastal ocean (Fig. 7). Hotspots for river PO${}_{\mathrm{4}}^{\mathrm{3}-}$ yields and loads tend to occur in river basins including densely populated large urban centers, such as the Chang Jiang, Huang He, Mekong, Shatt el Arab, Ganges, Godavari, Narmada, and Danube. The critical role of urban areas with sewage effluents in producing high river PO${}_{\mathrm{4}}^{\mathrm{3}-}$ yields is consistent with previous studies (Harrison et al., 2010; Mayorga et al., 2010). High river PO${}_{\mathrm{4}}^{\mathrm{3}-}$ yields also occur in humid river basins characterized by high P weathering rates, such as the Amazon, Parana, Zaire, Niger, Ganges, Chang Jiang, and Mekong, or in river basins including intensively farmed areas like the Mississippi (Harrison et al., 2010). The highest DOP and PP yields and loads tend to follow a pattern similar to that of PO${}_{\mathrm{4}}^{\mathrm{3}-}$, but there are also differences in patterns of PO${}_{\mathrm{4}}^{\mathrm{3}-}$ yields from patterns of PP yields. The differences, in part, result from deforestation and agricultural expansion in river basins like the Columbia and Amur, demonstrating elevated PP yields in comparison to PO${}_{\mathrm{4}}^{\mathrm{3}}$ yields (Harrison et al., 2019).

Figure 7 Global maps of the simulated river PO${}_{\mathrm{4}}^{\mathrm{3}-}$, DOP, and PP yields and loads for the year 1990.

4.3 Time series analysis

Correlations between the measurement-based and simulated annual river solid and nutrient load time series across eight stations in large US rivers for the period ∼ 1963–2000 demonstrate that, while model results for individual rivers may be biased, the simulated solid and N loads in different N species covary with variations of the measurement-based loads (Table 7, Figs. 8, S5). LM3-FANSY, however, has less capability to capture interannual variability of the P loads. The prominent difference in the solid and N dynamics vs. the P dynamics in LM3-FANSY is that the large terrestrial litter and soil sources for N are simulated by LM3, while the corresponding P inputs are externally prescribed because LM3 does not yet include terrestrial P dynamics. The lack of terrestrial P dynamics in LM3 is one of the most plausible reasons for the lower capability to capture the P load interannual variability (see Sect. 4.5 for further discussion).

Table 7 Summary of the LM3-FANSY capacity to simulate observation-based estimates of time-varying solid and nutrient loads from eight stations in large US rivers with extended time series. Rows correspond to stations and columns to different nutrients or suspended solids. Each table entry contains multiple skill metrics. The first entry indicates when the time series starts, with all time series continuing to 2000. The second entry provides the Pearson correlation coefficient quantifying the correspondence between observed and simulated annual loads. A p value indicating the significance level of this correlation is also provided. The third entry gives the prediction error computed as the difference between the simulated and measurement-based load averages expressed as a percentage of the measurement-based load average over the time series. The fourth line indicates whether the observation-based and LM3-FANSY estimates have significantly positive (+) trends, negative (−) trends, or no trend (none) over the time series period. Thus, “+/+” indicates that both the observation-based and LM3-FANSY estimates have a significantly increasing trend, whereas “+/−” indicates that the observation-based estimates have a significantly increasing trend, while the LM3-FANSY estimates have a significantly decreasing trend. Trend significance was tested at the 5 % level using the p value of the t statistic for a linear trend.

Prediction errors of the average river solid and nutrient loads for the periods ∼ 1963–2000 across the eight stations (Table 7) are comparable to the errors shown in the cross-watershed analysis (Table 4), as well as to the errors of cross-watershed analyses in previous global watershed modeling studies (Harrison et al., 2010; He et al., 2011; Mayorga et al., 2010; Pelletier, 2012). The cross-watershed emphasis herein reflects the intended global application of the model, but significant biases for individual rivers are a natural consequence of the prioritization. For example, the Mississippi River N loads in the LM3-FANSY simulation herein are significantly underestimated despite the model's modest global load bias. Focused investigations of larger misfits will be pursued in future work to increase model robustness through targeted enhancements to better reflect regional variations. Alternatively, for regional applications, tuning a single parameter for N dynamics or solid dynamics allows LM3-FANSY to be calibrated for a specific river. For example, for the Mississippi River, reducing the freshwater denitrification rate coefficient from 0.15 to 0.05 d⁻¹ can reduce the errors from −60 % to −15 % for NO${}_{\mathrm{3}}^{-}$ and from −23 % to 8 % for TN, while increasing the correlation coefficients from 0.7 to 0.75 for NO${}_{\mathrm{3}}^{-}$ and from 0.65 to 0.67 for TN (Fig. 8). Reducing the free parameter of terrestrial soil erosion by half can reduce the error of SS from 150 % to 27 %. Lastly, the t statistic shows that 30 of the 37 measurement-based time series loads have no significant linear trends over the time periods at the 0.05 level (Table 7). LM3-FANSY captures 30 of the 37 measurement-based trends in terms of whether the trends are significant and, when significant, all of the slope signs demonstrating downward vs. upward trends.

4.4 Model sensitivity with changes in parameter settings and nutrient inputs

For solid dynamics, the free scale parameter of terrestrial soil erosion (C₁) plays a significant role in determining the overall quantity of river SS loads, with its decrease by half (doubling) leading to a 50 % decrease (99 % increase) in global river SS loads (Table 8). The decrease (increase) in C₁ also reduces (enhances) the erosion-associated fluxes from terrestrial litter and soils and, in turn, river PON loads. In addition, the decrease (increase) in C₁ reduces (enhances) sorption of PO${}_{\mathrm{4}}^{\mathrm{3}-}$ to solid particles, as reflected by DIP vs. PIP load changes. These results associated with terrestrial soil erosion are, however, found to be insensitive to 1 order of magnitude changes in POM-to-PON ratios in eroded terrestrial soils and/or freshwaters (r_DN,Ero and/or r_DN, Sect. 2.2.1).

Table 8 Model sensitivities to the changes in inputs, components, and parameters examined based on percentage (%) differences in global river loads between the sensitivity and baseline simulations for the year 1990. The changed parameter values other than a decrease by half or a doubling are given in parenthesis.

Figure 8 Pearson correlation coefficients (r) and p values (p) for the measurement-based vs. simulated annual loads from large US rivers for the period ∼ 1968–2000. Gray shows the load responses to the changes in the freshwater denitrification rate coefficient from 0.15 to 0.05 d⁻¹ for N dynamics and the free parameter of terrestrial soil erosion by half for solid dynamics (Table 8).

For nutrient dynamics, the responses of river loads to a removal of each nutrient input source or an increase of it by 15 % suggest that terrestrial runoff and erosion are the dominant drivers of river N loads, followed by wastewater (Table 8). For river P loads, terrestrial runoff and erosion are also the dominant drivers, but unlike for N, the second dominant driver is weathering, followed by wastewater. Terrestrial runoff and erosion also play a critical role in explaining model efficiency and spatial distribution of river nutrient loads (Table 9). A removal of them reduces NSE and r values substantially. Wastewater plays a relatively small role in explaining model efficiency and spatial distribution of river NH${}_{\mathrm{4}}^{+}$ and PO${}_{\mathrm{4}}^{\mathrm{3}-}$ loads. Weathering plays a modest role in explaining model efficiency and spatial distribution of P loads in all species. A removal of aquaculture or atmospheric deposition has almost no impacts on the quantity, model efficiency, and spatial variation of river loads, suggesting that inaccuracies in these inputs have minor impacts on regional and global model estimates relative to the inaccuracies associated with the other model inputs. However, the importance of each source is likely to vary, depending on the dominant control on river nutrient loads in a specific region.

Table 9 Model sensitivities to the changes in inputs and components examined based on Pearson correlation coefficients (r) and Nash–Sutcliffe model efficiency coefficient (NSE) for the log-transformed measurement-based vs. simulated loads across major world rivers for the year 1990 (range for the years 1990–2000 in parentheses).

The parameter sensitivity tests show relatively insensitive responses of river nutrient loads to the changes in the rate coefficients of decay processes that break down complex organic compounds into simpler organic or inorganic compounds. The changes in the rate coefficients for hydrolysis, nitrification, and denitrification, however, have relatively large impacts on river nutrient loads. This is especially true for the freshwater N cycle, which includes an additional loss pathway to the atmosphere via denitrification unlike the freshwater P cycle. The role of freshwater denitrification in global river N loads, however, has not been explicitly investigated by previous global watershed modeling studies. Our sensitivity results imply a prominent role of freshwater denitrification in determining the amount of N loss to the atmosphere vs. N exports to the coastal ocean at both global (Table 8) and regional (Fig. 8) scales.

Algae dynamics play a significant role in determining the relative composition of inorganic vs. organic nutrients in freshwaters. Decreasing the algal mortality rate constant by half enhances algal uptake, decreasing DIN (the sum of NO${}_{\mathrm{3}}^{-}$ and NH${}_{\mathrm{4}}^{+}$) and IP (the sum of PO${}_{\mathrm{4}}^{\mathrm{3}-}$ and PIP) by −11 % and −33 %, respectively, while it increases ON (the sum of DON and PON) and OP (the sum of DOP and POP) by 24 % and 33 %, respectively. Similarly, doubling the maximum photosynthesis rate or chlorophyll-a-specific initial slope of the photosynthesis–light curve enhances algal uptake, decreasing DIN by −6 or −7 % and IP by −21 % or −25 %, respectively, while it increases ON by 18 % or 21 % and OP by 23 % and 28 %, respectively. The opposite holds for the parameter changes that reduce algal uptake. An analysis of the model sensitivity simulations, wherein the dynamic contributors to light extinction (i.e., ISS, POM, and Chl in Eq. 23) were removed, further suggests that proper light limitation of algal growth is also important for skillful estimates of freshwater inorganic vs. organic nutrients. Removing the dynamic light shading component leads to a ∼ 26 % and ∼ 35 % overestimation of ON and OP loads and an underestimation of DIN and IP loads by ∼ 10 % and ∼ 32 %, respectively. Inorganic nutrient levels are suppressed by invigorated algal populations and more nutrients end up in organic forms. Algal controls thus offer an effective means of calibrating the mix of inorganic and organic constituents. We note uncertainty associated with the fractions that partition nutrient fluxes from algae mortality to different nutrient forms, which appears to have a modest effect on organic nutrient loads (Table 8).

Finally, additional uncertainty tests have shown the relatively insensitive responses of river loads to a broad range in (1) the fractions that divide externally prescribed TN and TP inputs into different N and P forms (see Sect. 3.1, Table S10), (2) the fractions that partition fluxes from complex organic nutrient decomposition to simpler organic vs. inorganic nutrients, and (3) the temperature correction factor values that account for the temperature effect on freshwater biogeochemical reactions.

4.5 Discussion on uncertainties and future work

The inputs and transport of solids and nutrients through the terrestrial–freshwater system and transformations within it are governed by complex and interlinked physical, chemical, and biological processes. The understanding of these processes varies greatly, as does the degree of their inevitable simplifications within LM3-FANSY. We have highlighted the numerous uncertainties and simplifications in the model description and result presentation. Here, we discuss the prominent uncertainties that will be prioritized in future work.

There are several significant uncertainties in modeling the soil erosion and associated fluxes of solids and PON to rivers and the coastal ocean. Our global river PON load estimate (7, 7–8 TgN yr⁻¹) is lower than that of Global NEWS 2 (13.5 TgN yr⁻¹, Table 5), but confidence in both estimates is low, without explicit evaluations due to relatively limited direct measurements of PON. However, both simulated global river SS loads (10, 10–11 Pg yr⁻¹) and global litter and soil N storage (86 PgN) in LM3-FANSY are at the lower bounds of previous estimates (11–27 Pg yr⁻¹ for SS loads, Table 5; 95 (70–820) PgN for N storage, Post et al., 1985). While river TKN loads, which include PON, agree fairly well with the measurement-based estimates across various rivers (Fig. 3), it seems probable that our global estimate is on the low end. Several factors may have contributed to this. The current relatively coarse-resolution of global implementation herein inevitably “glosses over” areas of peak soil erosion. The single vertical layer formulation of LM3 omits any potential interactions between the vertical distribution of soil erosion and that of litter and soil N storage. An implementation of LM3-FANSY at higher resolution capturing a larger number of rivers will allow an expansive evaluation against a larger number of observations and facilitate a better assessment of these uncertainties.

The challenges of modeling particulates continue once they have entered the freshwater system. The Rouse-number-dependent transport criterion from Pelletier (2012) was adapted to simulate the deposition–resuspension fluxes between the suspended matter (i.e., ISS, PON, POP, and PIP) and benthic sediments (i.e., Sed, SedN, and SedP). The criterion was designed to primarily simulate suspended loads, typically accounting for >80 % of total (i.e., suspended and bed) loads from most large (>∼100 km²) river basins (Pelletier, 2012; Turowski et al., 2010), without explicitly modeling benthic sediments. We acknowledge that our simplified benthic sediment component resulting from adapting the Pelletier approach drives uncertainties in modeling the suspended matter and benthic sediments, including important diagenesis, other biogeochemical transformations, and physical processes (e.g., mineralization, denitrification, and net long-term organic burial) that occur within the benthic sediments. An implementation of more sophisticated benthic sediment dynamics with improved observation-based constraints (Chapra et al., 2008; Di Toro, 2001) and bed-load transport processes is thus subject to critical future work.

Uncertainties associated with sediment dynamics and bed-load transport are compounded by the relatively simple representation of lakes and the exclusion of anthropogenic hydraulic controls like damming, irrigation, and diversion that affect many rivers. For model evaluation, if available, we used the natural water discharges of GEMS-GLORI when calculating loads and yields from the GEMS-GLORI concentrations (see Sect. 3.3). Large dams or reservoirs, however, have been shown to impound solids and nutrients to substantially decrease their loadings to rivers (Vorosmarty et al., 2003). Thus, despite the relatively low global river SS loads in this first implementation of LM3-FANSY, the lack of such sediment trapping may have induced overestimations of solid and nutrient loads from river basins including large dams or reservoirs. As a representative example, the Colorado River Basin is known for nearly complete trapping of solids due to large reservoir construction and flow diversion (Vorosmarty et al., 2003). LM3-FANSY does not capture such extreme trapping. As a result, the Colorado River SS load simulated by LM3-FANSY (99 232 kt yr⁻¹) is more consistent with the corresponding load calculated by using the “natural” water discharge of GEMS-GLORI (120 010 kt yr⁻¹) than with the load calculated by using the “actual” water discharge (649 kt yr⁻¹). Although use of the actual water discharges is found to not significantly alter the cross-watershed evaluations (Fig. S6), such anthropogenic hydraulic controls are expected to further increase in the future (Seitzinger et al., 2010). It will thus be important to consider the effects of such controls in future work.

There is also significant room for further model development and improvement. An improved representation of lakes (e.g., vertical layering) is necessary to better resolve algal processes and associated transformations between inorganic and organic nutrient phases. Modeling large lakes with the ocean component of GFDL's Earth system model (Adcroft et al., 2019) is one of our priority developments, particularly given the importance of algae as a control on the relative proportions of inorganic vs. organic nutrients in freshwaters. An initial configuration for the US Great Lakes is currently under development.

There is also a need to pursue advances to provide a more comprehensive and consistent approach of modeling the coupled N, C, and P cycles across the terrestrial and freshwater continuum of LM3-FANSY. Expansion of LM3 to include terrestrial P dynamics will be targeted to improve estimates of litter and soil P storage and fluxes to streams and rivers, generating mechanistically consistent estimates of N:P ratios of nutrient loads reaching coastal systems. Priority enhancements will also include integration of freshwater C, alkalinity, and silicon dynamics with the current solid, algae, and nutrient dynamics of FANSY to simulate the impacts of river inputs on coastal C budgets and acidification and better understand a role of silicate limitation in the development and persistence of many HABs worldwide.

The current version of LM3-FANSY simulates denitrification emissions from freshwaters but does not include processes that explicitly separate N₂O and N₂ emissions from the freshwater denitrification. As freshwater N₂O emissions have been recognized as an increasingly important greenhouse gas source (Yao et al., 2020), it will be important to differentiate N₂ and N₂O emission processes in future work.

Although LM3-FANSY is capable of producing river solids and nutrients in various forms and units, some disagreements between the modeled and measurement-based estimates remain. Many observational studies have noted the uncertainties associated with measurement methods, location, and frequency that likely contribute to these disagreements. While in this new model development study, we have particularly focused on evaluation of LM3-FANSY's capability to capture broad-scale differences in yields, loads, and concentrations across a globally distributed set of rivers, we have also attempted to include a meaningful analysis of interannual variability and trends to bolster this assessment. It is difficult, however, to comprehensively analyze geographic differences, temporal variability, trends, underlying mechanisms, and the model's external (e.g., climate forcings) vs. internal (e.g., “land N memory effects”, Lee et al., 2021) forcing effects in the same paper because data sources are disparate and a substantial fraction of rivers we searched for were too small to be well resolved by our relatively coarse-resolution global framework. We thus focused on the largest and most comprehensively observed rivers supplying climate-scale (multi-decadal) information to provide a meaningful (albeit limited) complement to the spatial patterns. Additional time-varying constraints are thus needed to build additional confidence in projected changes. A commitment to fostering long-term, frequent solid and nutrient sampling and further scrutiny of the impacts of using different load estimation methods on time series data are both essential to better constrain the model. Finally, all of these model improvement efforts will be facilitated by implementing LM3-FANSY at higher resolution to allow comparison against extensive measurements from smaller rivers across the world. Resolution alone, however, cannot address model deficiencies, but simultaneous improvements in the other areas described above are needed to improve the model.

5 Conclusion

Our comparisons of process-based LM3-FANSY outputs with measurement-based estimates across major world rivers demonstrate skillful simulations for most riverine constituents despite being restricted to a universal parameter set – the same parameters for all the basins (i.e., without tuning of each basin). LM3-FANSY represents a significant step forward in terms of capacity to model coupled algae, SS, N, and P dynamics in freshwaters at a process-based, spatially explicit, global scale. Although this study is focused on model development and descriptions of the coupled freshwater SS, N, and P cycles, the capability of LM3 to simulate changes in vegetation and soil C–N storage in response to many terrestrial dynamics under subannual to centurial historical climate and land use changes (M. Lee et al., 2016, 2019, 2021) allows applications of LM3-FANSY for studies of temporal (subannual to multi-year) variability and long-term trends in global and regional water pollution. Therefore, LM3-FANSY v1.0 can serve as a baseline for studies aimed at understanding the effects of terrestrial perturbations on coastal eutrophication. The mechanistic modeling framework of LM3-FANSY is also well suited to make future projections by use of a new generation of future socioeconomic and climate scenarios over centuries, though we acknowledge that further work is needed to fully resolve underlying mechanisms.