Global integrated modeling framework of riverine dissolved inorganic nitrogen with seasonal variation

: Understanding patterns and seasonal variations of exces‐ sive nutrients in surface water from anthropogenic activi‐ ties is important for pollution control. In this study, we developed an integrated biogeochemical modeling frame‐ work for nitrogen exchanges among the atmosphere, terres‐ trial, and aquatic ecosystems. A land surface model, a ter‐ restrial nitrogen cycle model, and a riverine hydrodynamics model incorporated with a river temperature model were consolidated and driven by multiple nitrogen sources related to anthropogenic activities. We estimated the global nitrogen loading and transporting in global rivers, with con‐ sideration of seasonal variations, and the validation demon‐ strates the reliability of the proposed model. The total dis‐ solved inorganic nitrogen (DIN) flow rate is accumulated following rivers and it has high total DIN loads even in regions with low population density but large basin area, such as those at high latitudes. This study successfully improves estimation of nitrogen loading on global scale with consideration of seasonal variation. Our results show consistent trends with the observed data of DIN concentra‐ tions in global rivers, where all above variables are greatly affected by seasonal variations. The results also reflect the monthly-variant nitrogen inputs help produce closer DIN concentration estimates to observations, which will possi‐ bly stress the need for further study on seasonal variability of anthropogenic emissions.


INTRODUCTION
For many decades, one of the most popular topics in biogeochemical cycle and nutrient transport has been the cycle of nitrogen, which is a requisite and highly demanded element for living organisms on Earth. However, since 1970, the amount of biologically available nutrients entering the terrestrial ecosystem has more than doubled in comparison to preindustrial era levels (Galloway et al., 2008). Excess Correspondence to: Yizhou Huang, Department of Multidisciplinary Science, Graduate School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro-ku, Tokyo 153-8902, Japan. E-mail: hyz@rainbow.iis.u-tokyo.ac.jp nitrogen from anthropogenic activities has interfered with the biogeochemical nitrogen cycle of natural ecosystems, resulting in ozone layer depletion, surface water eutrophication, soil acidification, nitrate pollution of water bodies, and even impacts to human health (Ravishankara et al., 2009;Duan et al., 2013;Fowler et al., 2013;Steffen et al., 2015). Specifically, water quality in terms of nitrate (NO 3 -) leaching from point and non-point sources is a relevant environmental issue in the world (Galloway et al., 1995;Galloway and Cowling, 2002;Mena-Rivera et al., 2017). Understanding the connection between anthropogenic nitrogen discharge and changes to nature, and quantifying them in water bodies, are critical for humans in order to secure natural sustainability under further social development (Schlesinger, 2009;Steffen et al., 2015;Best, 2019).
To address the above-mentioned problems, numerical models have recently been widely used to estimate transport of certain chemical substances in rivers and streams Seitzinger et al., 2005;He et al., 2011;Beusen et al., 2015;Liu et al., 2019). At the very beginning, Meybeck (1982) started the questions about carbon, nitrogen, and phosphorus running through world rivers. Then the research on nitrogen cycle, primarily focused on river basin scale, has been brought to attention Dumont et al., 2005;Seitzinger et al., 2005). The research group of Global Nutrient Export from Watersheds (Global NEWS) established multiform and future prediction models of riverine nutrient export (Seitzinger et al., 2010). However, the methods which treated entire river basins as a basic unit failed to illustrate the spatial variations and leaching and transport of the nitrogen within the basin (Howarth et al., 2002;Bouwman et al., 2005;Harrison et al., 2005). For improving studies in this field, He et al. (2011) used an integrated biogeochemical modeling framework with detailed information on nitrogen leaching or nitrate concentration in individually distributed grids. Beusen et al. (2015) also incorporated Integrated Model to Assess the Global Environment-Global Nutrient Model (IMAGE-GNM) with a distributed hydrological model to simulate global nitrogen and phosphorus transportation. These models and integrated studies have contributed greatly to understanding the spatial distribution of nitrogen within a river basin and have linked calibrated model in order. However, the hydrological processes are still missing in water bodies. To overcome these shortages, Nevison et al. (2016) and Liu et al. (2019) used Community Land Model (CLM) with the River Transport Model (RTM) built inside, along with a coupler of the Community Earth System Model (CESM) to simulate riverine nitrogen transport. However, CLM likely underestimates leaching and runoff of all forms of nitrogen because it includes temperate-zone managed crops only (Liu et al., 2019), where manure plays an important role in soil nitrogen budget . In addition to the hydrological processes, thermal process can also play a role in nitrogen transport in rivers. Although considered in certain global model simulations (Liu et al., 2019), nitrogen transport with thermal process was rarely presented with monthly validations or reliable monthly reproductions and predictions.
The aim of this study is to develop an integrated biogeochemical modeling framework, including materials exchange among the atmosphere, terrestrial and aquatic ecosystems, and to investigate the nitrogen load transport in global rivers, with consideration of seasonal variations. This framework was built to explore the behavior of the anthropogenic released nitrate by incorporating several models. Various point and non-point nitrogen sources such as synthetic fertilizer, manure production, atmospheric deposition and municipal wastewater were introduced as loading dataset. The seasonal variation which impacts the monthly nitrate transport mechanism in rivers was also discussed for a more comprehensive understanding of nitrogen fate and transport.

Overview of integrated model framework
In this study we proposed an integrated model framework for nitrogen loads from different sources, and the transporting of nitrogen under effects of water temperature in global rivers (Figure 1). A land surface model − MATSIRO (Takata et al., 2003), a terrestrial nitrogen cycle model (TNCM) based on Lin et al. (2000) and Suga et al. (2005), and a riverine hydrodynamic model -CaMa-Flood (Yamazaki et al., 2011) incorporated with a river temperature model − HEAT-LINK (Tokuda et al., 2019), were consolidated for achieving such evaluations. For input datasets (Table I), the integrated model framework first produced global soil temperature, soil moisture, runoff, and runoff temperature from MATSIRO, which is driven by atmosphere forcing. The nitrogen cycle and the riverine hydrodynamics were built based on MATSIRO outputs. Nitrogen input (Text S1) introduced in this study consists of sources from synthesis fertilizer application (Nishina et al., 2017), manure production (Zhang et al., 2017), external atmospheric deposition (Tian et al., 2018), urban wastes flow (Morée et al., 2013), and aquaculture pollutions . The nitrate leaching from various land types was calculated by TNCM. The riverine hydrodynamic model was aimed to carry nitrogen from anthropogenic and natural sources through riverine network. Each grid cell in the riverine network acquired water containing dissolved inorganic nitrogen (DIN) from grid cells on upstream, as well as the DIN diffused from nearby land and point sources pollution within the grid cell. Meanwhile, the river temperature model controlled freezing and thawing status and affected material transport throughout a year. Lastly, annual total nitrogen amount at river outlets and monthly

Incorporation of TNCM and HEAT-LINK
The TNCM has been developed to take account of the mass balance of nitrogen in both organic soil and vegetation of the ecosystem. The terrestrial nitrogen cycle partially includes biological processes that rely on a different set of environmental factors. The internal parameters of TNCM, such as storage and fluxes, are the same as the model structure of the nitrogen cycle model by Lin et al. (2000Lin et al. ( , 2001, and two alternate versions revised by Suga et al. (2005) and He et al. (2009). In addition, this study also introduced the nitrogen source of manure production and atmospheric deposition ( Figure S1). The model contains four variables for nitrogen: nitrogen in vegetation (N veg , unit: Megagram nitrogen km -2 , hereafter Mg N km -2 ), organic nitrogen matter in soil (N som , Mg N km -2 ), ammonium (N amm , Mg N km -2 ), and nitrate (N nit , Mg N km -2 ). The nitrogen balance for each process is shown below: where, bfix is flux of nitrogen fixation as in nitrogen, uptk is flux of nitrogen uptake by plant, litf is flux of litter-fall from leaf, trunk, and root as in nitrogen, mnrl is flux of soil organic matter mineralization as in nitrogen, manu is the amount of manure production from livestock, depoamm (in+ex) is flux of internal and external nitrogen deposition as in ammonium, vola is flux of ammonia volatilization, ntrf is flux of nitrification, deponit (in+ex) is flux of internal and external nitrogen deposition as in nitrate, dnit is flux of nitrification and gas emission process, lech is flux of nitrate leaching, fertamm is the amount of ammonium in fertilizer, and fertnit is the amount of nitrate in fertilizer. All of these fluxes are in units of Mg N km -2 d -1 . hvst is the ratio of harvested crops. For natural ecosystem, all external nitrogen inputs such as depo (ex) , manu, and fert were set to zero. The function of coupled riverine hydrodynamic model (Tokuda et al., 2019) in this integrated model is to create a network of river channels and to simulate riverine transport of the nutrients, as well as to simulate seasonal variation under temperature change for global rivers in a more realistic manner. For instance, the freeze-thaw cycle of a river can directly impact river discharge, and thus nitrogen transport. The governing equations in HEAT-LINK follow the one-dimensional mass conservation law, the Saint-Venant equation for momentum conservation, and the conservation law of energy. More detail of the models and data used in this study are specified in the Supplemental Material.

Model setting
Similar research recently performed validations of DIN load in 1995 (He et al., 2011;Liu et al., 2019;Tian et al., 2020). Therefore, the initial database set up in this study was based on the year of 1995. An experimental group of constant nitrogen input was added for discussion. An additional run for 2001-2010 was also implemented for 10-year simulations which can show recent temporal changes of the DIN. The spin-up period of the entire integrated model was set as 20 years. The initial condition for the river water temperature was the air temperature on the first day of 1995, and CaMa-Flood assumes the sea-level elevation for the initial condition.

Data analysis
In addition to R 2 , the determination of efficiency, to evaluate model performance we also evaluated root mean Where O is the mean of observed data, O i is observed data, S i is the simulated data, and n is the number of observed data. We consider pRMSE under 50% to be acceptable in perspective of the global scale of the integrated model (Beusen et al., 2015).

Model validation with annual DIN load
We validated the integrated model in terms of the annual DIN load and yield (load per area) for the year of 1995. Figure 2 shows the comparisons between simulated results and observations for 30 major rivers in different continents (10 in Asia, 7 in Europe, 8 in North America, 3 in South America, 1 in Africa, and 1 in Australia; also see Figure 3). The bar graphs are sorted in terms of the basin size from the largest to the smallest in Figure 2a and 2b. Results show that the basin size is one of the factors for DIN load, as higher DIN load is observed for larger rivers. However, higher average DIN yield can also lead to larger DIN load at river outlets, even though the basin size is smaller (e.g. Pearl River and Rhine River). Low DIN loads are observed in Murray and Rio Grande due to the low yields there, however the low load in Tornionjoki is probably due to its small size.
The scatter plots in Figure 2c and 2d indicate a linear relationship between the simulated and observed DIN load and yield, respectively. The regression lines between simulated and observed DIN load and yield are very close to the 1:1 line. The coefficient of determination R 2 is 0.86 and 0.96 for the two variables, respectively. The pRMSE for DIN load was 20.4% and 13.5% for the DIN yield. From here it can be said the integrated model reproduced annual

Result of global spatial distribution
The global spatial distribution of total DIN flow rates in 1995 is illustrated in Figure 3. The DIN flow rate is generally high for rivers with large annual discharge including Amazon River, Mississippi River and Yangtze River. Besides the large river basins selected for validation, rivers in Western Europe, Northern China and South Asia also have high values in annual DIN flow rates. Thanks to the integration of riverine models, the total DIN flow rate is accumulative along the rivers, from upstream to the downstream. Additionally, the major rivers in higher latitude, where seasonal variations play a larger role than regions in lower latitude, contribute much higher DIN flow rates compared with nearby rivers.

River discharge and temperature
The river temperature directly affects the freeze-thaw cycle and thus leads to impact on the seasonality of river discharges, and further affects the concentration and transport of nutrients. Here we selected 4 major rivers (i.e. Ob, Mackenzie, Lena, and Yukon) with sufficient observed nutrient data to evaluate the model's ability in reproducing the temperature and discharge. Figure S2 compares the estimated monthly river temperature with the observations. In these high-latitude regions, rivers are covered with ice on the surface and the water temperature decreased to zero degrees during winter periods. Despite the cold season, the model can simulate the river temperature and show the general seasonal variations with the observed water temperature fairly well.
Discharge simulation is more difficult than temperature simulation because the bias in discharge is propagated from precipitation, runoff generation and river routing. It has been observed that the discharge simulation in the high latitudes remains challenging in many studies. We also compared the monthly river discharge with observations and the simulation at the selected rivers ( Figure S3). However, there is some shift between observed and simulated results. Also, the discharge is overestimated in the month when river thawing begins. This is because the current method of river ice processes and snow melting in the model is comparatively simple. When runoff increases at the beginning of the ice breaking season, the ice is pushed up and crashes. The floating ice flows downstream, and temporary ice jams are formed when some of the floating ice forms a cluster, especially in narrow sections. The model ignores such processes and instead considers the removal of surface ice in mechanisms related only to heat budget (Tokuda et al., 2019).
The coefficient of determination R 2 has been tested for discharge and temperature. The estimated river temperature works very well on this model with 0.98, 0.94, 0.95 and 0.99 for Ob, Mackenzie, Lena, and Yukon rivers, respectively. The estimated river discharge also works well for Ob and Lean rivers, with 0.87 and 0.64, respectively. The model does not predict well for Mackenzie and Yukon, with R 2 of 0.17 and 0.28, respectively. Despite this, they still show a suitable seasonal trend to DIN simulations.

Monthly average DIN concentration
We further assessed monthly average river DIN concentration against observations in selected rivers (Figure 4). The average concentration has an inversed seasonal variation compared to the river discharge, as the concentration usually reaches the highest level during the freezing season and has a lower level after ice breaking. Ob River had an acceptable result, except for a time shift when the lowest concentration occurs. Beside slight overestimation of concentration, Lena River had a one-month shift of the lowest concentration. The reason may be the uncertainty of river discharge estimation, which is mentioned above. Mackenzie River had a fairly close match in DIN concentration. Yukon River had an overall underestimation in DIN transport, but still had a shape that corresponds to the seasonal variation. Overall, the pRMSE for DIN concentration in rivers mentioned above are 27.0%, 9.2%, 47.3%, and 24.0%, respectively. Extra uncertainties and shortcomings of the model will be discussed in the Discussion section. From the perspective of monthly time series simulation, the model reproduced a reasonable seasonal variation of monthly average DIN concentration in major rivers.
We also demonstrate the impact of temporal resolutions  Figure 4. Except the Mackenzie River, the monthly variation had changed a lot compared to the timevariant N inputs (red line). The pRMSE for the simulation of the four rivers are 20.5%, 28.9%, 56.6%, and 27.0%, respectively. Although the pRMSE does not show a great increase, the monthly-variant N inputs produce closer DIN concentration estimates to the observations, comparing to the constant nitrogen inputs. The differences indicate that the time-variant inputs will to some degree benefit the DIN estimations. The improvement should work together with the introduced temperature model as the impact of temperature is different in nitrogen concentration and transport. For assessing the model performance for more recent conditions, we ran the simulation for a further 10 years (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010) for Mackenzie and Yukon rivers (due to data availability, Figure 5). In general, the simulated DIN concentration is in the similar degree of the observations and the pRMSE are 48.8% and 46.6% for Mackenzie and Yukon. The average simulated DIN concentration is around 0.03 mg L -1 in summer and 0.10 mg L -1 in winter, showing the importance of seasonal variations and the confidence of long-term simulations of this model.

DISCUSSION
There are still several aspects to be discussed about uncertainties of the integrated model framework, including the datasets and the integrated model.
Regarding the nitrogen input, several uncertainties could be pointed out. GDP, population, and crop area are used to estimate the global fertilizer distribution. However, even in the existing activity data, the national census has inherent spatial distribution (Winiwarter and Muik, 2010). For the global manure data, the livestock distribution was generated by using one-phase static GLIMS (Zhang et al., 2017) without variations in time. Moreover, the manure, municipal waste and aquaculture pollution datasets are in yearly resolution, which makes precise seasonal variation study difficult. Hence, adjustments of all inputs with spatial distribution at the subnational level and temporal variations at monthly resolution would be perfect.
As we know, the river hydrodynamic model governs nutrient transport in rivers, and the uncertainties that occur  in-stream will also create uncertainties in estimating material transport. First, the uncertainty of runoff generated from the land surface model can lead to biases in river discharge resulting in very different DIN concentration results. Second, human activities such as irrigation abstraction and dam regulation can also change the DIN concentration in rivers. Moreover, the river temperature will affect chemical changes in rivers (Duan et al., 2016) but it is not considered. Nutrient reaction with other components including other nutrients, chemicals, microbes, and other pollutants is also yet to be included.

CONCLUSIONS
It is very important to understand in-stream dynamics with nitrogen transport for future related studies. In this study, a fully coupled biogeochemical model with land surface model and river dynamic model was proposed and has been proven with data consistency and efficiency in terms of the temperature, river discharge, DIN concentration, and nitrogen load. This study successfully improves the estimation of nitrogen loading on a global scale with consideration of seasonal variation. Additionally, the results of monthly-variant nitrogen input suggests the need for further study on the seasonal variability of the anthropogenic emissions. This model framework provides a chance of incorporating future prediction information, in order to understand the consequences of climate change and anthropogenic activities to global nitrogen cycle over the long term.

ACKNOWLEDGMENTS
The authors are grateful for a MEXT scholarship by The Ministry of Education, Culture, Sports, Science, and Technology of JAPAN and acknowledge the support from the Japan Society for the Promotion of Science [KAKENHI:16H06291] and the Environment Research and Technology Development Fund [JPMEERF20202005] of the Environmental Restoration and Conservation Agency of Japan. We also gratefully acknowledge the United Nations Environment Programme Global Environment Monitoring System (GEMS) for providing us with monthly observed data of river water temperature and nutrients concentration, and the Global Runoff Data Centre (GRDC) for the daily observed discharge data.

SUPPLEMENTS
Text S1. Detail of the models and data used Figure S1. Flow chart for the terrestrial nitrogen cycle model, orange boxes represent added parts of fluxes (Based on Suga et al., 2005) Figure S2. Comparison of monthly simulated (red line) and observed (black dots) river temperature (°C) for selected rivers at specific stations Figure S3. Comparison of monthly simulated (red line) and observed (black line) river discharge (m 3 s -1 ) for selected rivers at specific stations