Linking global terrestrial and ocean biogeochemistry with processbased, coupled freshwater algae-nutrient-solid dynamics in LM3-FANSY v1.0 by Lee et al.

This paper presents the development of a coupled river (and lake) biogeochemical model for N and P on the global scale, based on the coupling of LM3, the freshwater module of FANSY. The P input and some of the N sources are taken from another data source (IMAGE-GNM).

The authors claim that the new model FANSY is incorporated in the Land Model (LM3). This coupling results in a more process-based representation of biochemistry model FANSY-LM3 which links the terrestrial and ocean biochemistry. The purpose of the LM3 model: “It captures processes including changes in vegetation functioning, plant-soil nitrogen cycling, and streamflow dynamics. The land model sheds light on the effects of these processes on atmospheric physics and chemistry and, conversely, the effects of changing climate and CO2 concentration on the land.” (from https://www.gfdl.noaa.gov/land-model/). It is clearly a goal of this coupled system (FANSY-LM3) to capture changes and streamflow dynamics. The FANSY model has indeed a lot of dynamics incorporated (and are described here), so this model is well suited to couple with the LM3 model. However, the results that are presented here are not dynamic, but are on one single time period. How can the reader discover whether FANSY is capable of producing reasonable results when it is used in combination with the LM3 model? If the purpose of the model only is to simulate the chosen time period, then I would not take this model, but the Global NEWS model. The advantage of the Global NEWS model is that it is published and well cited and the results are very transparent! The current presentation of the result does not convince me to start using a dynamic process-based model. Also the sensitivity analyses are focused on the Pearson correlation coefficient of that same time period on a limited set of rivers. This approach does make it very hard to understand how the model is reacting on different inputs. It is not clear to me, why these 64 basins are chosen and whether they are a good representation of the global river basins. Also the use of different time periods to do this comparison is creating all kinds of uncertainties. I would suggest to do a proper sensitivity analyses (not a one factor analysis but change all parameters at once) on the important model output parameters to show and understand the mechanistic behavior of the model. The current validation of the model is not good enough and is not convincing. Additional validation on long term data sets is needed to show that the dynamic patterns are well simulated.

I would advise to reject the current manuscript or major revisions under the condition that there will be a better sensitivity analyses including a more extensive validation over a longer time period.

 

Major comments

1.    It is very important to know if the calibration data and validation overlap. However, details about calibration data are absent in this paper.

2.    There should be an overview of the loads that flow into the river, the removal loads to the different sectors (sediments, air, …) and the export to the coastal waters for N and P.

3.    The consistency of the P input compared to the N input is a concern. As long as this is not clear, I would suggest to avoid to calculate and present the N:P ratio.

4.    In line 438-440 it is clearly states that FANSY does not simulate the processes in reservoirs. This should be discussed including what this means when this model is used in “unprecedented scenarios”.

5.    It is not appropriate to use the measurement data unevenly collected during for 1970s-1990s to validate model results for 1982-2010, especially because there are a lot of published global measurement data available for validation for the period 1982-2010. This introduced high temporal uncertainty may be larger than the model uncertainty itself. Substantial work is needed to re-do the model validation to eliminate this uncertainty.

Specific comments

Line 22-23: How you do evaluate which are “better global estimates of nutrient inputs to rivers”? As far as I know, there is no way to validate nutrient inputs to rivers, and it’s hard to say which is a better estimate.

Line 32-38: This sentence is too long, and extends for 7 lines. Please re-organize it into several short sentences to increase readability.

Line 40: “a comprehensive model” is not clear. In terms of what aspect?

Line 48 and 52: Term use: NO3- instead of “NO3”, same for NH4 and PO4, same for everywhere in the text.

Line 43-45: “In particular, inorganic nutrients, which are characterized by higher bioavailability than organic forms (Sipler and Bronk, 2004), have been recognized as critical drivers of algal blooms (including non-HABs) and hypoxic events (Kemp et al., 2005).” It is “excessive inorganic nutrients” rather than simply “inorganic nutrients” that may drive algal blooms and hypoxic blooms.

Line 45:”Meanwhile” should be “Besides”

Line 51-53: “Furthermore, nutrient and algae dynamics are strongly linked with solid dynamics through phosphate (PO4) sorption/desorption interactions with solid particles (McGechan and Lewis, 2002) and algae growth reduction due to light shading by suspended solids (SS) (Dio Toro, 1978) ”: There are more processes that link nutrient and algae dynamics with solid dynamics, and what you list here are only two of them.

Line 55-57: “Projecting global freshwater biogeochemistry changes requires process-based models that are robust under unprecedented conditions expected in the next century.” This paper does not involve any future projections, but “unprecedented conditions expected in the next century” and similar statements have appeared many times in this paper. This paper does not show that this model can handle “unprecedented” situations. Therefore, I would advise to remove this kind of sentences.

Line 59: “Prior applications of process-based freshwater biogeochemistry models, such as …, have generally been limited to small watersheds”. This is not true. Many process-based models have been developed to simulate freshwater biogeochemistry in inland waters more than “small watersheds”, such as DLEM (Tian et al., 2020; Yao et al., 2020; Bian et al., 2022), LOAC (Akbarzadeh et al., 2019), INCA(Whitehead et al., 1998a; Whitehead et al., 1998b; Wade et al., 2002), IMAGE-DGNM (Vilmin et al., 2020; Vilmin et al., 2022), for large river basins and global scale. A better understanding and summary of processed freshwater biogeochemistry models are needed.

Line 69-70: Indeed IMAGE-GNM does not differentiate the different forms, but therefore the IMAGE-DGNM is developed (Vilmin et al., 2020; Vilmin et al., 2022).

Line 73-75: “IMAGE-GNM takes a mass balance approach to calculate soil nutrient budgets, which at times rests on simple scaling without potential dynamical feedbacks (e.g., an estimation of litter from floodplains to rivers as 50% of total net primary production (Beusen et al., 2015)).” Firstly, floodplain vegetation is not calculated as an item in the soil budget in IMAGE-GNM but as a direct source to surface water. Secondly, the use of 50% of net primary production (NPP) for wetlands and floodplains can be uncertain, but the area of floodplain and wetlands due to hydrology, as well as the simulated NPP are both spatially explicit and time-dependent. It is therefore not correct to say that calculation of nutrients from floodplain vegetation to rivers in IMAGE-GNM is just simplified without dynamic feedback. Not to mention that for all the other soil and water budget, IMAGE-GNM simulates them in a mass-balanced and dynamic way. 

Line 86-89: A sentence too long to read.

Line 90: “performance assessment against measurement-based global and regional estimates across world major rivers”. Ambiguous. The global result cannot be compared with “major river (s)” to assess model performance.

Line 119-123: These sentences are unclear: First sentences states that sewage and aquaculture are from prescribed sources. Last line claims all P input is not from LM3? Please make more clear.

Line 132-133: “Algae chlorophyll a (CHL), algae C (Ca), algae P (Pa), and algae dry matter (Da) are diagnosed from algal N (Na) and, in the case of CHL, nutrient and light conditions”, not clear. It’s better to separate it into two sentences and describe the two situations.

Line 135-138: “The 5 prognostic P variables include the same organic forms as for N (dissolved organic P (DOP), particulate organic P (POP), and sedimentary organic P (SedP)), but includes dissolved and particulate inorganic forms (PO4 and PIP, Sec. 2.2.4) rather than the oxidation state distinction as done for N. ” This sentence is too long and not so readable. Consider improving it. Perhaps add “variables” at the end of this sentence, after “N”.

Line 140: Ambiguous expression of “variable processes”: do you mean processes which are variable, or variable’s processes, or variables and processes? 

Line 4 in Table 1: From Line 128 in the text, I think you mean Sed is the abbreviation for “bottom sediment”. But in Table 1, it is not clear and differs from the text. Modification is needed to ensure the term use is consistent and clear.

Table 1 and text: This is a biogeochemistry model which includes quite some biogeochemical processes and different chemical forms. Therefore, it is not appropriate to use the common symbols of Na, Ca and Pa, as new abbreviations of the terms in this study. Please consider revising these symbols to ensure clarity and avoid confusion.

Figure 1: (1) The pathway of Algae dry matter to SS is quite strange. Do you mean that this part of organic matter will go to the pool of SS in water or finally Sed in sediment? Then this results in permanent additive pools and the pathway, which will only increase. This means that you may have partitioned the algae organic matter directly into one part which will be degraded, and the other refractory part which will not participate in biogeochemical reactions. The refractory part after partitioning will always increase the SS (water) and Sed pools (bottom sediment) pools, but will not add to the total PON and POP pools, which is quite strange. Or another possible assumption for this is that the dead algae are counted as SS which only participates in the physical processes, while in fact the dead algae already become NH4+ (see your pathway of another algae mortality), which should not add to the SS pool.

Besides, you have SedN and SedP specially for N- and P-contained sediment. The interactions seem somewhat unclear and not self-explanatory.

(2) algae growth and mortality are not the common definitions themselves (see Line 216-217), which calls for earlier explanations or notes to inform the readers from the very beginning. Otherwise this is too confusing. 

Line 156 – 157: Equation 1 and 2. The parameter Ri is used in both equations with the same sign. Is that correct? Is Ri not part of Fin or Fout? Can you explain this better?

Line 167: Which rainfall is meant here? Annual, monthly, weekly, daily? Clarify.

Line 170: C1 is calibrated. How and on what? Please elaborate on this.

Line 1 in Table 2: “Many reactions are reported at 20℃” is not appropriate for the use of the ideal temperature (Tref) of all reactions. Detailed citations of mechanism studies are needed. Also, the dynamic feedback, which has been advertised so many times, cannot be well examined with the same Tref for every biogeochemical process related to every element form driven by every organism group. Many mechanism studies have shown the obvious difference in Tref of nitrification, denitrification, mineralization, primary production, etc. Using the same Tref constraint for all biogeochemical processes related to all element forms driven by organism groups will yield a wrong estimation of the process rates, element and algae processing and the consequent nutrient concentrations. 

Line 6 in Table 2 and Line 192: “R” symbol is used before, as term reaction. It is confusing to use the same letter for two different terms. Change one of them.

Equation 8: Please refer to the original reference (Ludwig and Probst). This equation is very old and only based on 19 observations. Are there not more observations to improve this?

Line 216-217: the definitions of algae growth and mortality differ from other studies! This should be mentioned to inform readers, or at least tell them these processes are “generalized” or “net”. Otherwise, it is too confusing.

Line 273: it is hard to understand how equation #23 calculates the combined the net flux of a generalized mortality (non-predatory mortality + grazing + settling + excretion) as mentioned in Line 216-217. Contradictorily, according to Line 271-277, this equation does not reflect so-called “generalized mortality”, but grazing (death because of predation).

Line 291: “assuming C weight content …. 1500 kg m3”. Does this mean that you use a global number for bulk density and for the C content of the soil? I hope not…

Line 297-296: I wonder why PON cannot be calculated from TKN-NH4-DON for this calibration purpose.

Line 299-302: I wonder whether or not fPON can or should be used for compensating hydraulic controls like dams and reservoirs. I’m afraid this is not very appropriate for a mechanistic model. This factor is constant (Table 2) and from the calibration of measurements. But it is unclear how this calibration was performed. Line 295 claims that fPON  slows down the N movement, but that that not make any sense…. Also the “limited measurement-based PON estimates” makes it a bit a mystery.

Line 328: see comment for Line 273. I don’t find any rationality for using a part of the mortality calculation of predation-caused algae death as DON excretion to NH4+. Also, I don’t know how fDON used here is calibrated, which is very strange to be done because generally it is unknown how much DON excretion will account for in algal mortality due to predation.

Equations 35 – 38: What is CPN?

Line 367: Calculation for 1900 – 2010. But Beusen et al. (2015) has provided only data up to 2000. How did you calculate the last ten years?

Line 368: Confusing: “reported point and nonpoint N and P inputs to rivers”  Can you better describe what you used?

Line 370: 1990 calibrated (constant) from partitioning factors for instream processes may be not suitable for the period 1982-2010 considering the possible temporal changes. Calibration process is not clear here. 

Line 373-381: This atmospheric N deposition is calculated here, and can be different from Beusen et al. 2015. However, In Table 3 it is clearly written that atmospheric N deposition is from Beusen et al. 2015. So what are you using? Please clarify.

Table 3: So LM3-FANSY simulates soil N budget and soil budget and their delivery to rivers and lakes, but no P budget. This means that P is not coupled to the LM3. How consistent are the assumptions of P of Beusen et al. (2015) with the assumptions of LM3? Do they both have agricultural land on the same spots? This is a major concern of me, because one of the results is the N:P ratio. If the assumptions do not match, then the N:P ratio does not mean anything. Perhaps the N:P ratio is a bridge to far for this model at the moment…..  

Table 3: The numbers in Table 3 are not the same as in in Vilmin et al. (2018). The fractions for wastewater are dependent on the treatment level. Here it is assumed that whole the world has primary treatment, which is not true. The fractions of agricultural surface runoff and aquaculture are different than Vilmin et al. (2018). How is this determined?

Line 434-445: The use of observation data for this model is very strange. There are a lot of published global measurement data available for validation for the period 1982-2010, but this paper intentionally chooses the measurement data during the 1970s-1990s which Global NEWS used. This introduced high temporal uncertainty may be larger than the model uncertainty itself. It is not clear how the comparison is made. The timestep is 30 minutes. So how is the value that is used for calibration constructed? Discharge weighed? Mouth of the river or the location of the observation? Only the year 1990? Or average over yearly concentrations of 1970 – 1990 or 1982 – 2010? This is very unclear. The time period difference is very strange and is also a concern. A lot has happened with the N and P loads/concentrations in the period 1982 – 2010.

Line 444-445: “Cross-watershed …. and trends.” Here the purpose of the model is defined. I do not agree with this. If this is true, then you don’t need to spinup 11000 years, and do all the dynamic calculations. Just use the Global NEWS model. This model is far to complicated for this purpose and this purpose does not fit in the LM3 framework. 

Line 450-452: It is unclear whether it is because of the calibrated parameter of instream solid dynamics that improves the model performance, because the results without calibrated parameter are not shown.

Table SI1 – SI9: Nice to have all the observations together, but I miss the results from FANSY in these tables. Then the reader can also see, which rivers are doing fine or not.

Table 4:  Again, “R” is used as a new abbreviation for another term in the same paper.

Table 4 and 5: Not clear how this 15% is applied. Is the change applied only for 1990 and the calculation is only done for 1990? Constant over the year or not? Or total N or P and then some forms get more change than others? Please elaborate on this in the methods section.

Tables 4-5: hard to interpret. To me, this sensitivity does not make any sense. First of all because the output parameter is the Pearson correlation coefficient. I can not understand this, how this helps me to understand the model mechanics. Better to do this on output parameters of the model. Then also an explanation could be added on why this difference is big or not. Secondly, I think that a one factor analysis is very limited for this kind of models. 

Table 6:  Needs to mention the estimates from IMAGE-GNM are for 2000 instead of the 1990s as in other studies.

Table 6:  LM3-FANSY results needed to be added in this table for straightforward comparison.

Line 485-486: This is not true. global river chemical measurement data are not limited, but they are simply not used in this paper.

Line 495-498: Explanations are needed for the estimate ranges like 36.4-41.3. Not clear whether this is because of temporal change or uncertainty or something different.

Line 507: This model estimates a global export of 6.5-7.8 TgP. But if I correctly understand this article, you used the same riverine input as Beusen et al (2015) who estimated a global export of 4 TgP. So there is problem here. The global input of P is around 9 TgP (stated in Beusen et al. (2015). So the P loss in the river is around 30% whereas the P loss in Beusen et al is above 50%.  This is a huge difference and should be explained here. Also the N loss should be checked!

Line 521: should be “Arctic”.

Line 551: grammar error: “are the among the …”

Figures 6-8: The maps are for 1990 but the comparisons of concentrations, loads and yields are performed for 1982-2010. Why is there a difference in time between the two result analyses?

References

Akbarzadeh, Z., Maavara, T., Slowinski, S., Van Cappellen, P., 2019. Effects of damming on river nitrogen fluxes: a global analysis. Global Biogeochem. Cycles 33, 1339-1357.

Beusen, A.H.W., Bouwman, A.F., Van Beek, L.P.H., Mogollón, J.M., Middelburg, J.J., 2016. Global riverine N and P transport to ocean increased during the 20th century despite increased retention along the aquatic continuum. Biogeosciences 13, 2441-2451.

Beusen, A.H.W., Van Beek, L.P.H., Bouwman, A.F., Mogollón, J.M., Middelburg, J.J., 2015. Coupling global models for hydrology and nutrient loading to simulate nitrogen and phosphorus retention in surface water. Description of IMAGE-GNM and analysis of performance. Geosci. Model Dev. 8, 4045–4067.

Bian, Z., Pan, S., Wang, Z., Yao, Y., Xu, R., Shi, H., Kalin, L., Anderson, C., Justic, D., Lohrenz, S., Tian, H., 2022. A Century-Long Trajectory of Phosphorus Loading and Export From Mississippi River Basin to the Gulf of Mexico: Contributions of Multiple Environmental Changes. Global Biogeochem. Cycles 36, e2022GB007347.

Tian, H., Xu, R., Pan, S., Yao, Y., Bian, Z., Cai, W.-J., Hopkinson, C.S., Justic, D., Lohrenz, S., Lu, C., Ren, W., Yang, J., 2020. Long-term trajectory of nitrogen loading and delivery from Mississippi River Basin to the Gulf of Mexico. Global Biogeochem. Cycles 34, e2019GB006475.

Vilmin, L., Mogollón, J. M., Beusen, A. H. W., and Bouwman, A. F.: Forms and subannual variability of nitrogen and phosphorus loading to global river networks over the 20th century, Glob. Planet. Change, 163, 67–85, 2018.

Vilmin, L., Bouwman, A.F., Beusen, A.H.W., van Hoek, W.J., Mogollón, J.M., 2022. Past anthropogenic activities offset dissolved inorganic phosphorus retention in the Mississippi River basin. Biogeochemistry 161, 157-169.

Vilmin, L., Mogollón, J.M., Beusen, A.H.W., van Hoek, W.J., Liu, X., Middelburg, J.J., Bouwman, A.F., 2020. Modeling process‐based biogeochemical dynamics in surface fresh waters of large watersheds with the IMAGE‐DGNM framework. Journal of Advances in Modeling Earth Systems 12, 1-19.

Wade, A.J., Durand, P., Beaujouan, V., Wessel, W.W., Raat, K.J., Whitehead, P.G., Butterfield, D., Rankinen, K., Lepisto, A., 2002. A nitrogen model for European catchments: INCA, new model structure and equations. Hydrol. Earth Syst. Sci. 6, 559-582.

Whitehead, P.G., Wilson, E.J., Butterfield, D., 1998a. A semi-distributed Integrated Nitrogen model for multiple source assessment in Catchments (INCA): Part I - Model structure and process equations. Sci. Total Environ. 210-211, 547-558.

Whitehead, P.G., Wilson, E.J., Butterfield, D., Seed, K., 1998b. A semi-distributed integrated flow and nitrogen model for multiple source assessment in catchments (INCA): Part II - Application to large river basins in south Wales and eastern England. Sci. Total Environ. 210-211, 559-583.

Yao, Y., Tian, H., Shi, H., Pan, S., Xu, R., Pan, N., Canadell, J.G., 2020. Increased global nitrous oxide emissions from streams and rivers in the Anthropocene. Nat. Clim. Chang. 10, 138-142.

In this study, Lee et al. intend to develop a comprehensive freshwater model FANSY which can simulate the dynamics of algae, solid and nutrient dynamics in inland waters, and the fresh water model has been incorporated into the Land Model LM3. The topic is very important and should be attractive to a broad of readers. The paper is well organized and the methods and materials have been described adequately. This study involves tremendous of model development work. I appreciate the effects of the authors on developing such a model.

Nonetheless, I still have some concerns on the algorithms used in this study to simulate the transfer & transformation of sediment, N and P through the inland water systems. In specific:

1) To calculate the erosion-induced POC loss from soil to inland waters, why not calculate the POC loss rate based on the soil SOC concentration and the soil erosion rate, especially you are going to develop a process-based model? The empirical equation from Busen et al., 2005 might be too coarse to be used here for calculating the POC loss rate, and it cannot represent the effects of SOC dynamics on the riverine POC flux. In addition, many studies have proved the effectiveness of calculating POC loss rate based on the soil erosion rate and the SOC concentration, such as Tian et al., 2015 (Anthropogenic and climatic influences on carbon fluxes from eastern North America to the Atlantic Ocean: A process-based modeling study) and Zhang et al., 2022 (Estimating the lateral transfer of organic carbon through the European river network using a land surface model. Earth System Dynamics)

2) PON generally flows together with POC, as POC, PON and POP are all contained in the POM (particulate organic matter). Almost all particulate matter loss from land to rivers is caused by soil erosion, rather than the water drainage. Thus, it does not make sense to calculate PON loss based on the soil water drainage. A more reasonable solution is that, loss of all particulate matters (e.g. sediment, POC, POC, POP) from land to rivers should be calculated based on the soil erosion rate and the concentrations of these matters in soil; loss of all dissolved matters (DOC, DON, DOP) should be calculated based on soil water drainage (i.e. the leaching processes). Tian et al., 2015 & Zhang et al., 2022 can be good references.

3) To my knowledge, riverine DON mostly originated from soil during the leaching processes. Decomposition of PON and algae mortality only contribute a small part of the total DON. However, the FANSY model in this study assumed all DON is originated from the decomposition of PON and algae mortality. This assumption does not make sense. I would suggest the authors to add the fluxes of DOC, DON and DOP from land to river during the leaching processes. Many studies have revealed that the fluxes of dissolved matters caused by leaching is much larger than the loss of particulate matters caused by soil erosion (Regnier et al., 2022; Nature, The land-to-ocean loops of the global carbon cycle)

4) I am surprised that the model does not consider the emission of N2O from inland waters. Both observation and simulation indicated that a larger part of the riverine N is loss to air in the form of N2O (Yao et al., Nature Climate Change, Increased global nitrous oxide emissions from streams and rivers in the Anthropocene)

5) Many of the forcing data used in this study were obtained from Beusen et at., 2015. It is better to provide more information about these data. Foe example, the spatial and temporal resolutions of these data, and the method used by Beusen et al., 2015 to produce these data.

6) Discussion on the uncertainties of this study is weak. Consider the simplification and ignorance of many processes that are closely related to the transfers of sediment, C, N & P. There should be more discussion on the potential uncertainties in the simulated results. For example, the uncertainties in simulated soil erosion rate, POC loss rate and the effects of dams and reservoirs.

7) Compared to existing models, a highlight point of the FANSY model in this study is that it explicitly simulated the biogeochemical cycles related to the algae. However, there is no evaluation of the simulated algae processes, and it is hard for the readers to determine whether the scheme used in the FANSY model is effective or reliable on simulating algae dynamics. Is it possible for the authors to collect some observation data to validate the simulated algae biomass, growth rate or mortality rate in this study.

Minor comments:

In Eq. 3, how is the slope extracted? Is the slope extracted from DEM at 1-degree resolution?

In Eq. 4, how the depth of river and lake is calculated? Do you have a map of area and lake as the forcing data of your model? If yes, is the area of river or lake fixed?