Description of the tool

GREENeR is an open-source R package for assessing annual time series of nutrient loads in a river network and at the basin outlets, and contributions of nutrient sources to these loads. The package provides tools and methods to apply the model Geospatial Regression Equation for European Nutrient losses (GREEN; (B. Grizzetti et al. 2005; B. Grizzetti, Bouraoui, and De Marsily 2008; Bruna Grizzetti, Bouraoui, and Aloe 2012; Bruna Grizzetti et al. 2021)) to an area of interest. A brief description of the model, including sources and parameters, can be found at the end of this document.

The package includes functions to perform graphical summaries of model inputs, calibrate model parameters, run sensitivity analysis checks, and visualize model inputs and outputs through graphs and maps of total loads and contributions by source. The package works for both total nitrogen (TN) and total phosphorus (TP). It allows the analysis of different scenarios of nutrient input in the river basin or region of study. The package is parallel-capable to alleviate the computational burden in large basins.

An example case: the Lay river Basin

GREENeR functionalities are illustrated for the Lay river basin, in France. Please note that the basin data are extracted from a pan-European dataset (Bruna Grizzetti et al. 2021) for the purpose of showing the package tools and the analysis steps for a generic region. Since the dataset may not reflect local sources correctly, implications of results in terms of nutrient fluxes may be incorrect and are not evaluated nor commented herein.

The package includes two nutrient type scenarios (one for TN and one for TP), with disaggregated information on nutrient inputs for all catchments that compose the basin (spatial resolution of the Catchment Characterisation Model; (De Jager and Vogt 2007)). Both input (TN and TP) time series span from 1990 to 2018. The dataset includes also spatial information in vectorial format to enable mapping of tool results.

Loading and description of the scenario

# load the GREENeR package
library(GREENeR)

GreeneR needs information on the topology of the catchment and nutrient sources. All information is contained in two data.frames. Note that data.frame structures and column names must be respected when creating a new dataset, otherwise the functions will fail.

The topology is contained in the data.frame “catch_data_TN” (or catch_data_TP; in the example dataset, the Lay river Basin topology is presented).

# load the topological sequence of catchments and complementary info
data(catch_data_TN)
head(catch_data_TN)
#>     HydroID To_catch Shreve LakeFrRet NrmLengthKm
#> 1    361076   361195      1         0    0.010270
#> 28   361521   362009      4         0    0.059295
#> 55   362311   362426      2         0    0.011319
#> 84   361187   361195      2         0    0.011948
#> 109  361383   361521      1         0    0.005731
#> 136  362149   362175      1         0    0.002374

The fields are defined as follow:

• HydroID (integer): unique identifier of the catchment

• To_catch (integer): unique identifier of the catchment to which the catchment goes. Note that for the outlet To_catch== -1

• Shreve (integer positive): this indicates the Shreve order of the topological sequence in the stream network. GREEN resolves sequentially all loads for increasing numbers of Shreve order (headwaters, i.e. catchments that do not receive water, are Shreve == 1, catchments that receive waters only from Shreve 1 are Shreve == 2, etc till the outlet)

• LakeFrRet (fraction, 0-1): lake retention fraction (Lret in eqs 1, 6, 7). Note, this is different for TN or TP module

• NrmLengthKm (double): normalized length of catchment reach (eq 5)

The nutrient sources in each catchment per year are contained in the data.frame “annual_data_TN” for the TN scenario (“annual_data_TP” for the TP scenario).

# load the sources of TN for each year and catchment
data(annual_data_TN)
head(annual_data_TN)
#>   BasinID YearValue HydroID NextDownID   Atm    Min    Man   Fix  Soil    Sd
#> 1  291994      1990  361076     361195 9.608 28.111 21.379 2.909 2.108 0.514
#> 2  291994      1991  361076     361195 9.608 28.295 20.318 2.888 2.108 0.589
#> 3  291994      1992  361076     361195 9.608 26.482 21.514 2.917 2.108 0.665
#> 4  291994      1993  361076     361195 9.608 24.060 19.861 2.868 2.108 0.741
#> 5  291994      1994  361076     361195 9.608 23.725 20.199 2.890 2.108 0.816
#> 6  291994      1995  361076     361195 9.608 24.271 21.772 2.830 2.108 0.892
#>      Ps YearlyMass ForestFraction InvNrmRain
#> 1 0.314         NA      0.1700787   0.081113
#> 2 0.251         NA      0.1700787   0.075192
#> 3 0.188         NA      0.1700787   0.069081
#> 4 0.126         NA      0.1700787   0.066110
#> 5 0.063         NA      0.1700787   0.047314
#> 6 0.000         NA      0.1700787   0.061390

# load the sources of TP for each year and catchment
data(annual_data_TP)
head(annual_data_TP)
#>   BasinID YearValue HydroID NextDownID    Bg   Min   Man    Sd    Ps YearlyMass
#> 1  291994      1990  361187     361195 0.082 8.577 3.108 0.330 0.315         NA
#> 2  291994      1991  361187     361195 0.082 8.300 2.819 0.319 0.300         NA
#> 3  291994      1992  361187     361195 0.082 7.688 3.393 0.307 0.285         NA
#> 4  291994      1993  361187     361195 0.082 6.614 3.660 0.296 0.270         NA
#> 5  291994      1994  361187     361195 0.082 6.119 4.508 0.284 0.256         NA
#> 6  291994      1995  361187     361195 0.082 6.025 4.676 0.273 0.241         NA
#>   ForestFraction InvNrmRain
#> 1              0   0.081765
#> 2              0   0.075672
#> 3              0   0.068371
#> 4              0   0.066009
#> 5              0   0.047398
#> 6              0   0.061488

In the annual data frames some fields are in common:

• BasinID (integer): unique identifier of the Basin. This may be useful when the region of interest encompasses more than one basin; all catchments belonging to a same Basin share the same BasinID.
• YearValue (integer): the year for which data are defined
• HydroID (integer positive): unique catchment identifier (matches the ID in catch_data)
• NextDownID (integer positive): corresponds to To_catch
• ForestFraction: Non-agricultural land cover in the catchment (fraction, FF in eqs 6 and 7)
• InvNrmRain: inverse of normalized rainfall (eq 4)
• YearlyNMass: observed annual total nitrogen load (TN ton/yr) from monitoring station data where available (NA when no data is present for that catchment and year). These data allow calibration of model parameter; thus sufficient observed loads are essential for the calibration process.

The other fields specify nutrient sources, and differs in the two nutrient scenarios. In the case of TN, fields are:

• Atm: Annual nitrogen deposition from atmosphere (ton/yr)
• Min: Annual amount of nitrogen from mineral fertilisers (ton/yr)
• Man: Annual amount of nitrogen in manure fertilisers (ton/yr)
• Fix: Annual amount of nitrogen fixation by leguminous crops and fodder (ton/yr)
• Soil: Annual amount of nitrogen fixation by bacteria in soils (ton/yr)
• Sd: Nitrogen input from scattered dwellings (ton/yr)
• Ps: Nitrogen input from point sources (ton/yr)

In the case of TP fields are:

• Bg: annual amount of phosphorus background losses (ton/yr; eq 7)
• Min: annual amount of phosphorus mineral fertilisers (ton/yr)
• Man: annual amount of phosphorus in manure fertilisers (ton/yr)
• Sd: phosphorus input from scattered dwellings (ton/yr)
• Ps: phosphorus input from point sources (ton/yr)

For some functionalities of the GREENeR package, the geographical information of the catchments is required in the form of an sf (simple feature) object (the sf package, Pebesma & Bivand, implements the Simple Features open standard)

# load the geographical information of the basin (require for some functionalities)
data(sh_file)
head(sh_file)
#> Simple feature collection with 6 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 3440600 ymin: 2645700 xmax: 3476900 ymax: 2691300
#> Projected CRS: ETRS89-extended / LAEA Europe
#>   OBJECTID HydroID NextDownID BasinID JunctionID AreaSqkm DrainAreaS Shreve
#> 1    71204  368447         -1  291994   20837332    43.68    1971.10     95
#> 2   350329  362484     362498  291994   33061724     0.58       0.58      1
#> 3   350337  368467     368429  291994   33061732    21.37      21.37      1
#> 4   350345  362497     362454  291994   33061740    25.52      25.52      1
#> 5   350346  362499     362498  291994   33061741     0.03     158.18      8
#> 6   350347  362498     362537  291994   33061742     1.44     160.20      9
#>     FecID   CsbID SuppID  Csb0ID Shape_STAr Shape_STLe
#> 1 6031419 7000265  71204 7100074   43680000      59800
#> 2 6031445 7000265 350329 7100074     580000       3600
#> 3 6031420 7000265 350337 7100074   21370000      32600
#> 4 6031450 7000265 350345 7100074   25520000      36400
#> 5 6031445 7000265 350346 7100074      30000        800
#> 6 6031445 7000265 350347 7100074    1440000       8000
#>                         geometry
#> 1 MULTIPOLYGON (((3447500 264...
#> 2 MULTIPOLYGON (((3476600 268...
#> 3 MULTIPOLYGON (((3440600 265...
#> 4 MULTIPOLYGON (((3458000 268...
#> 5 MULTIPOLYGON (((3476700 268...
#> 6 MULTIPOLYGON (((3476000 269...

Exploring scenario inputs

The library includes functions to examine the nutrient sources in a Basin, and explore how these are applied, either globally, or distributed over time or in space.

A first summary is provided by the input_plot() function, which - when option “B” is given, shows the mean annual nutrient loads per source for the whole period of analysis.

# the barplot for the Lay TN and TP Scenarios
input_plot(annual_data_TN, sh_file, "Lay", "B")
input_plot(annual_data_TP, sh_file, "Lay", "B")

Barplot for the Lay TN scenario.

Barplot for the Lay TP scenario.

By selecting option “D”, input_plot() returns density plots of nutrient loads per catchment. These are useful to see the distribution of nutrient source inputs in the Basin.

# the density plots for the Lay TN and TP Scenarios
input_plot(annual_data_TN, sh_file, "Lay", "D")
input_plot(annual_data_TP, sh_file, "Lay", "D")

Density plot for the Lay TN scenario.

Density plot for the Lay TP scenario.

The function input_Tserie() allows to examine the temporal evolution of the inputs. There are three options: displaying cumulative total amounts (“gr1”), total amounts (“gr2”), or total amounts per unit area (“gr3”). The latter graph is particularly useful when catchment sizes are when the range of catchment areas is large to compare nutrient application intensity between catchments.

In the Lay river example, these graphs show, among other things, how the amount of P from mineral fertilization (Min) has decreased over time, whereas the amount of P from manure has increased.

# the time serie plot type 1 (areas)
input_Tserie(catch_data_TN, annual_data_TN, sh_file, "Lay", "gr1")
input_Tserie(catch_data_TP, annual_data_TP, sh_file, "Lay", "gr1")

Time series plot type 1 (areas) for TN and TP.

# the time serie plot type 2 (lines)
input_Tserie(catch_data_TN, annual_data_TN, sh_file, "Lay", "gr2")
input_Tserie(catch_data_TP, annual_data_TP, sh_file, "Lay", "gr2")

Time series plot type 2 (lines) for TN and TP.

Time series plot type 3 (areas by year and km2) for TN and TP.

There is a fourth alternative for the function input_Tseries(), “gr4”. This option compares the levels of nutrient inputs in three zones of the basin (upper, middle and lower part). The upper zone includes the catchments up to the median (50%) shreve value, the middle zone includes the catchments from 50% to 75% of the shreve values, and the lower part includes the catchments from 75% to 100%. At the top of each graph, the corresponding area of each part is indicated. In the example below, it can be seen that differences in N inputs between the upper and lower part of the basin are noticeable over the whole period.

# the time serie plot type 4 (by km2 and Shreve)
input_Tserie(catch_data_TN, annual_data_TN, sh_file, "Lay", "gr4")
input_Tserie(catch_data_TP, annual_data_TP, sh_file, "Lay", "gr4")

Time series plot type 4, thatcompares the levels of nutrient inputs in three zones of the basin.

The function input_maps() shows nutrient inputs distribution in the Basin. The function generates a map for each nutrient source, plus a map (in yellow-green scale) of the total diffuse nutrient sources, which can be contrasted to the point sources map. The first option (“gr1”) shows total inputs (kt/year).

# The title of the plot
mapTitle <- "Lay Basin"

# the Input Load Map by source type 1 (kt/year)
input_maps(catch_data_TN, annual_data_TN, sh_file, mapTitle, "gr1", legend_position = 3)

Input map by source type 1 for TN scenario

For the TP scenario:

# the Input Load Map by source type 1 (kt/year)
input_maps(catch_data_TP, annual_data_TP, sh_file, mapTitle, "gr1", legend_position = 3)

Input map by source type 1 for TP scenario

Option “gr2” shows inputs per area (in kt/y/km2):

# the Input Load Map by source type 2 (kt/y/km2)
input_maps(catch_data_TN, annual_data_TN, sh_file, mapTitle, "gr2", legend_position = 3)

Input map by source type 2 for TN scenario

# the Input Load Map by source type 2 (kt/y/km2)
input_maps(catch_data_TP, annual_data_TP, sh_file, mapTitle, "gr2", legend_position = 3)

Input map by source type 2 for TP scenario

Model calibration

The quality of model calibration depends on the quality and number of available observations. Before diving into calibration it is thus useful to start by exploring observations.

# Lay Basin scenario calibration for TN
# subset of data rows with reference values to be used in the calibration
TN_ref_values <- annual_data_TN[!is.na(annual_data_TN$YearlyMass),]
TN_ref_values
#>      BasinID YearValue HydroID NextDownID    Atm     Min     Man    Fix   Soil
#> 400   291994      2003  361226     361187 12.021  38.799  27.080  4.763  2.972
#> 401   291994      2004  361226     361187 11.832  39.125  27.297  4.526  2.960
#> 402   291994      2005  361226     361187 11.472  37.502  27.631  4.356  2.960
#> 969   291994      2007  363217     363876 25.208  84.777  74.632  9.129  7.660
#> 970   291994      2008  363217     363876 22.580  84.131  77.208  9.018  7.660
#> 971   291994      2009  363217     363876 23.018  71.262  86.496 10.419  7.660
#> 1793  291994      2007  362500     363555 34.932 108.107  95.171 11.641  9.768
#> 1794  291994      2008  362500     363555 30.246 107.284  98.455 11.500  9.768
#> 1795  291994      2009  362500     363555 31.590  90.874 110.299 13.286  9.768
#> 3016  291994      2007  363878     363855 10.129  33.955  29.892  3.656  3.068
#> 3017  291994      2008  363878     363855  8.744  33.696  30.923  3.612  3.068
#> 3018  291994      2009  363878     363855  9.047  28.542  34.644  4.173  3.068
#> 3633  291994      2007  365569     365715 66.515 258.580 227.638 27.844 23.364
#> 3634  291994      2008  365569     365715 59.035 256.612 235.494 27.506 23.364
#> 3635  291994      2009  365569     365715 59.319 217.360 263.824 31.779 23.364
#> 4532  291994      2003  366683     367700 35.532 191.333 133.540 23.489 14.656
#> 4533  291994      2004  366683     367700 36.139 193.511 135.009 22.385 14.640
#> 4534  291994      2005  366683     367700 32.287 185.481 136.662 21.546 14.640
#> 4535  291994      2006  366683     367700 37.763 161.308 144.156 21.797 14.640
#> 4536  291994      2007  366683     367700 37.096 162.028 142.639 17.447 14.640
#> 4780  291994      2008  366683     367700 32.642 160.794 147.562 17.236 14.640
#> 4781  291994      2009  366683     367700 32.023 136.199 165.314 19.913 14.640
#>         Sd    Ps YearlyMass ForestFraction InvNrmRain
#> 400  1.496 1.111     17.724    0.152901786   0.062910
#> 401  1.430 1.118     12.667    0.152901786   0.070694
#> 402  1.364 1.126      8.762    0.152901786   0.090762
#> 969  1.815 0.744    832.388    0.099490976   0.059976
#> 970  1.963 0.683    934.193    0.099490976   0.050057
#> 971  2.111 0.621    778.023    0.099490976   0.059081
#> 1793 4.305 4.011    351.973    0.062302006   0.060989
#> 1794 4.493 4.007    391.211    0.062302006   0.050472
#> 1795 4.682 4.003    283.156    0.062302006   0.059423
#> 3016 1.244 0.344    354.944    0.000000000   0.060897
#> 3017 1.095 0.505    418.703    0.000000000   0.050287
#> 3018 0.946 0.665    311.261    0.000000000   0.058547
#> 3633 8.045 4.348    235.181    0.019378028   0.059086
#> 3634 8.295 4.210    277.517    0.019378028   0.049728
#> 3635 8.546 4.072    190.167    0.019378028   0.056961
#> 4532 3.998 0.982   2685.982    0.001319609   0.067024
#> 4533 3.869 1.092   2476.694    0.001319609   0.073931
#> 4534 3.740 1.201   1555.908    0.001319609   0.088373
#> 4535 3.951 1.096   4353.326    0.001319609   0.055746
#> 4536 4.162 0.990   3665.792    0.001319609   0.064562
#> 4780 4.373 0.885   4288.808    0.001319609   0.054822
#> 4781 4.584 0.779   3264.154    0.001319609   0.064394

# number of reference data
nrow(TN_ref_values)
#> [1] 22

# distribution of the references by year and location
table(TN_ref_values$HydroID,TN_ref_values$YearValue)
#>         
#>          2003 2004 2005 2006 2007 2008 2009
#>   361226    1    1    1    0    0    0    0
#>   362500    0    0    0    0    1    1    1
#>   363217    0    0    0    0    1    1    1
#>   363878    0    0    0    0    1    1    1
#>   365569    0    0    0    0    1    1    1
#>   366683    1    1    1    1    1    1    1

In the example basin there are 22 observations (derived from publicly available dataset of European Environment Agency, EEA for 6 different stations from 2003 to 2009; note however that only one station has references for all years. (For the TP case there are 58 observations, collected in eight catchments from 1997 to 2012, with one complete time-series).

To run the calibration process for a scenario (function calib_green()), the following settings must be defined:

1. The expected range for each parameter. This range is defined by two vectors of three values, one for the lower limits and one for the upper limits of the three parameters. The values correspond to each of the parameters of the model in sequence: alpha_P (eq 4), alpha_L (eq 5), and sd_coeff (eqs 6 and 7).

2. The number of iterations to be performed during the calibration process. The higher the number of iterations, the more likely it is to achieve quality parameters, but the longer the computation time required. Although it depends on the catchment and the proposed intervals, it is recommended to run at least 200 iterations to have enough information to continue the calibration process.

3. The years to be included in the calibration process.

# Parameter for the calibration of the model
# the lower limits for all params (alpha_P, alpha_L, sd_coef)
low <- c(10, 0.000, 0.1) 
# the upper limits for all params (alpha_P, alpha_L, sd_coef)  
upp <- c(50, 0.08, 0.9)       

# number of iterations
n_iter <- 2000    

# years in which the model should be executed
years <- c(2003:2009)

In the example below, a calibration of TN in the Lay basin is performed with the parameters indicated above. The process is parallelised and will use all the cores of the computer. The computation time depends on the computer, the basin, and the number of iterations.

# execution of the calibration  
DF_calib_TN <- calib_green(catch_data_TN, annual_data_TN, n_iter, low, upp, years)

The calibration function applies a Latin Hypercube sampling scheme to the three parameters within the possible range (defined by lower and upper limits) and evaluates model performance (predictions against available observations) for the calibration period (specified in ‘years’) by calculating several “goodness-of-fit” metrics. The function returns a dataframe with parameters and goodness-of-fit scores that can be further analyzed.

The function applies the following goodness-of-fit metrics (Althoff and Rodrigues 2021):

• NSE: Nash-Sutcliffe efficiency
• rNSE: Relative Nash-Sutcliffe efficiency
• mNSE: Modified Nash-Sutcliffe efficiency
• cp: Coefficient of Persistence
• VE: Volumetric Efficiency
• KGE: Kling-Gupta efficiency
• d: Index of Agreement
• md: Modified index of agreement
• rd: Relative Index of Agreement
• r: Linear correlation coefficient
• R2: Coefficient of determination
• PBIAS: Percent bias
• MAE: mean absolute error
• RMSE: Root mean square error
• ME: Mean error
• MSE: Mean square error
• NRMSE: Normalized Root Mean Square Error

Exploring calibration results

Choosing the right goodness-of-fit metric to select the best set of parameters for a model largely depends on the overall study scope, the area of interest (e.g. high or low load, upper or lower catchment area), and the available observation dataset (size and quality). The library includes a series of functions to examine the calibration results and help select the most suitable set of parameters.

The calib_boxplot() function shows relationships between best parameter sets chosen according to one goodness-of-fit parameter (title of each boxplot) according to one goodness-of-fit metric (title of plot) in relation to others (x label). Only the six metrics that are used most frequently in hydrological calibration are included in this figure. In the lower panel, the figure shows as well the distribution of model parameters in the best parameter sets. The absolute best parameter set according to each goodness-of-fit metrics is marked as red dots in each boxplot.

# Generating the box plots  
rateBS <-  5  #  percent of parameters selected from the whole calibration set
calib_boxplot(DF_calib_TN, rateBS)

The calib_dot() function shows the distribution of parameters in relation to each other for a chosen goodness-of-fit metric. With this function the figure can be generated for all the previously described goodness-of-fit metrics.

# Generating the dot plots  
Gof_mes <- "NSE"
calib_dot(DF_calib_TN, Gof_mes)

In the example figure, the best results for alpha_P to achieve best NSE are in the range 0.25-0.30 (cyan color), while there are no clear patterns for the other two parameters. Also, no significant interaction between parameters can be appreciated. We can conclude that the most sensitive parameter in this case is alpha_P.

The scatter_plot() function shows all parameter realizations in the calibration dataset against a selected goodness-of-fit metric to visualize the influence of each parameter on model results. The result depends on the goodness-of-fit metric the user wants to consider. With this function the figure can be generated for all metrics generated by green_calib() function.

# Generating the scatter plots  
Gof_mes <- "NSE"
scatter_plot(DF_calib_TN, Gof_mes)

Selecting the best parameter sets

The best parameter set according to a selected goodness-of-fit metrics can be extracted with the select_params() function:

Gof_mes <- "NSE" # according the NSE goodness of fit metric
NSE_bestParams <- select_params(DF_calib_TN, Gof_mes)
NSE_bestParams
Param_NSE2 <- as.numeric(NSE_bestParams[2:4])

The function can be applied with different goodness-of-fit metric to obtain different potential model realizations:

Gof_mes <- "rNSE" 
rNSE_bestParams <- select_params(DF_calib_TN,Gof_mes)
rNSE_bestParams    
Param_rNSE2 <- as.numeric(rNSE_bestParams[2:4])

(note that actual values returned by the function could change at each calibration run)

Annual observed vs predicted values for a parameter set

Once a parameter set has been selected, it is useful to verify graphically the differences between the observations (ObsLoad) and the model predictions (PredictLoad). Included in the tool is the function simobs_annual_plot(), which allows a year-to-year comparison of observed and predicted values for a given parameter sets.

# annual scatter plot comparing observed vs modeled loads by year
label_plot <- "NSE best params for TN in the Lay" 
simobs_annual_plot(catch_data_TN, annual_data_TN, Param_NSE2[1], Param_NSE2[2], Param_NSE2[3], years, label_plot)

Comparing Observed vs Predicted for two sets of parameters

The best parameter set depends on the goodness-of-fit metric. It is advisable to pre-select alternative parameter sets (according to different GoF) and compare model results, before making the final selection of the parameter set. The compare_calib() function shows simultaneously all observations against modelled values obtained from two parameter sets.

setPlabels <- c("NSE", "rNSE")    
label_plot <- "Comparing two sets of parameters for Lay TN" 
compare_calib(catch_data_TN, annual_data_TN, Param_NSE2[1], Param_NSE2[2], Param_NSE2[3], Param_rNSE2[1], Param_rNSE2[2], Param_rNSE2[3], years, label_plot, setPlabels)

At the top left of the figure, the scores of six frequently used goodness-of-fit metrics for the two parameter sets (as identified in setPlabels) are reported.

Estimation of nutrient fluxes and source contribution in the basin

Once the most appropriate parameter set has been selected, it is possible to run the model to estimate nutrient loads across the basin.

The function region_nut_balance() runs the GREEN model with the selected parameter set, and returns the region nutrient mass balance of the nutrient fluxes for the whole simulation period. The results of this function can be visualized in a Sankey diagram with the function N4_sankey().

# Computing the nutrient balance  
years <- c(2003:2009)
Nut_bal_TN <- region_nut_balance(catch_data_TN, annual_data_TN, Param_NSE2[1], Param_NSE2[2], Param_NSE2[3], years)
# Plot the sankey plot with the result of the balance
sank <- N4_sankey(Nut_bal_TN)

Share of nutrients

The function green_shares() runs the GREEN model with the selected parameter set and returns for all catchments nutrient loads and contributions by each nutrient source in the simulation period. The results of the model can be examined by several functions.

The function nutrient_tserie(), with the option “gr1”, shows the temporal evolution of total load outputs for the upper, central and lower catchment area.

years <- c(2003:2009)
basin_sa <- green_shares(catch_data_TN, annual_data_TN, Param_NSE2[1], Param_NSE2[2], Param_NSE2[3], years)
    
# The title of the plot
plot_title <- "Time series Load Output for the Lay Basin"    
# Output Load Basin average time series (lines)
nutrient_tserie(basin_sa, sh_file, plot_title, "gr1")

The “gr2” option shows the temporal evolution of the total load at the basin outlet in the simulation period.

# Output Load Basin average time series (lines)
nutrient_tserie(basin_sa, sh_file, plot_title, "gr2")

The “gr3” option shows the temporal evolution of the total load at the basin outlet in the simulation period with contributions by sources.

# Total Load in the Basin Outlet by source apportionment time series (lines)
nutrient_tserie(basin_sa, sh_file, "Lay Basin", "gr3")

The function nutrient_maps() shows maps of nutrient loads in the basin. The option “gr1” shows loads by source (in logarithmic scale to improve the visualization as the differences in a basin may vary over several orders of magnitude).

# The title of the Map
map_title <- "Output Loads for the Lay Basin"
    
# Basin Output Load Maps by source
nutrient_maps(basin_sa, sh_file, map_title, "gr1", legend_position = 3)

The alternative “gr2” shows the total load at catchment outlets (in logarithmic scale to improve their visualisation, as the differences between sources and outlets are usually several orders of magnitude).

# Basin Output Total Load  Maps 
nutrient_maps(basin_sa, sh_file, map_title, "gr2", legend_position = 3)

The alternative “gr3” produces maps of total load by catchment area (in kt/km2/y) in the basin.

# Basin Output Specific Load by km2 Maps 
nutrient_maps(basin_sa, sh_file, map_title, "gr3", legend_position = 3)

GREEN model equations, sources and resources for input data

The Geospatial Regression Equation for European Nutrient losses (GREEN; B. Grizzetti et al. (2005); B. Grizzetti, Bouraoui, and De Marsily (2008); Bruna Grizzetti, Bouraoui, and Aloe (2012); Bruna Grizzetti et al. (2021)) is a conceptual model to assess nutrient loads (total nitrogen and total phosphorous) from a basin given diffuse and point emissions to the land and river network. The model comprises two major nutrient pathways: diffuse sources that undergo retention in the land phase (basin retention) before reaching the stream, and point sources that are directly discharged into surface waters. Once in the river network, all sources are reduced by the in-stream retention (river retention).

Diffuse sources include nutrient from mineral fertilisers, manure application, nitrogen crop and soil fixation, inputs from scattered dwellings (i.e. isolated houses and small agglomerations that are not connected to sewerage systems), and nitrogen atmospheric deposition (for nitrogen module) or background losses (for phosphorus module). Point sources consist of urban and industrial wastewater discharges. Basin retention is modelled as a decay function proportional to the inverse of the total annual precipitation in the catchment; river retention is estimated as a decay function proportional to the river length, considered as a proxy for water residence time. In addition, lake retention is simulated as a function of lakes residence time and average depth.

The basin is divided into spatial subunits (called catchments), which have a given area, a river reach, an in-let node, and an outlet node. The catchments are topologically connected from basin heads to the outlet in cascading sequence. The sequence of nutrient load accumulation through the stream network is defined by Shreve (1966) order. Nutrient input from the different sources, basin and river retention are simulated in each small catchment and routed through the river network. For each catchment i in the basin, the GREEN nutrient load \(L_i\) is estimated by the general equation:

\[\begin{equation} \label{eq:green_nutrient_load} L_i = (1 - Lret_i) * (DS_i * (1 - Bret_i) + PS_i + U_i) * (1 - Rret_i) \end{equation}\]

where:

• \(L_i\) = Nutrient load at the catchment outlet (out-node) (ton/yr)

Sources of nutrients are:

1. \(DS_i\) = Nutrient diffuse sources in the catchment (ton/yr)
2. \(PS_i\) = Nutrient point sources in the catchment (ton/yr)
3. \(U_i\) = Nutrient load from upstream catchments (ton/yr)

Sinks of nutrients are:

1. \(Lret_i\) = Lake retention (fraction). Lret can be defined according to Kronvang et al. (2004): \[\begin{equation} Lret_i = 1 - \frac{1}{[1 + (\frac{7.3}{z}) * RT]} \text{(for nitrogen)} \end{equation}\] \[\begin{equation} Lret_i = 1 - \frac{1}{[1 + (\frac{26}{z}) * RT]} \text{(for phosphorus)} \end{equation}\]

where:

• \(z\) average lake depth (m),
• \(RT\) = hydraulic residence time (yr).

At European scale, average lake depth and hydraulic residence time can be obtained from HydroLAKES database (https://www.hydrosheds.org/pages/hydrolakes, Messager et al. (2016)).

2. \(Bret_i\) = Basin retention (fraction): \[\begin{equation} \label{eq:basin_retention} Bret_i = 1 - \exp(-alpha_p * NrmInvRainy) \end{equation}\] where \(NrmInvRain\) is the inverse of annual precipitation (in mm) of year y, normalized by the maximum in the dataset: \(NrmInvRain = (1 / precipitationy ) / max(1 / precipitation))\)

3. \(Rret_i\) = River retention (fraction): \[\begin{equation} Rret_i = 1 - \exp(-alpha_l * NrmLengthKmi) \end{equation}\]

where \(NrmLengthKm\) is the length (in km) of the catchment reach, normalized by the maximum in the dataset: \(NrmLengthKm\) = (length of catchment \(i\) reach, in km) / \(max\)(Reach length in the basin).

Thus, the inverse of precipitation and the river length are scaled by maximum scaling (Frank and Todeschini 1994). Note that in this way basin retention varies from year to year according to annual precipitation, whereas river retention for a given catchment is constant. The equation is applied sequentially from the most upstream nodes to the outlet. The model parameters are:

1. the basin retention coefficient (\(alpha_p\)), which together with annual precipitation in each catchment regulates nutrient retention of diffuse sources \(\eqref{eq:basin_retention}\);
2. the river retention coefficient (\(alpha_l\)), which regulates the river retention per river km.
3. the fraction of domestic diffuse sources that reaches the stream network (\(sd_{coef}\)). Note that in Europe the amount ot nturients from scattered dwellings is a minor source of nutrients, therefore this parameter is not very sensitive, and a fixed value of \(0.50\) is recommended. In other regions of the world however this may not be the case, hence the model in GREENeR package defines this fraction as a calibration parameter.

\(\eqref{eq:green_nutrient_load}\) defines the general GREEN function. However, diffuse and point sources are defined differently for each nutrient module, i.e. nitrogen or phosphorus. More specifically in GREEN nitrogen model, for the total nitrogen load \(L_i\) of catchment \(i\) is estimated as:

\[\begin{equation} L_i = (1 - Lret_i) * ((MinN_i + ManN_i + FixN_i + SoilN_i + (1 - FF_i) * AtmN_i) * \\ (1 - Bret_i) + 0.38 * FF_i * AtmN_i + sd_{coef} * SdN_i + PsN_i + U_i) * (1 - Rret_i) \end{equation}\]

where:

• \(MinN_i\) = Annual amount of nitrogen from mineral fertilisers (ton/yr)
• \(ManN_i\) = Annual amount of nitrogen in manure fertilisers (ton/yr)
• \(FixN_i\) = Annual amount of nitrogen fixation by leguminous crops and fodder (ton/yr)
• \(SoilN_i\) = Annual amount of nitrogen fixation by bacteria in soils (ton/yr)
• \(AtmN_i\) = Annual nitrogen deposition from atmosphere (ton/yr)
• \(FF_i\) = Non-agricultural land cover in the catchment (fraction)
• \(SdN_i\) = Nitrogen input from scattered dwellings (ton/yr)
• \(PsN_i\) = Nitrogen input from point sources (ton/yr)
• \(U_i\) = Nitrogen load from upstream cacthments (ton/yr).

Note that nitrogen atmospheric deposition losses are split into two parts, i.e. inputs to agricultural land undergo the basin retention (which includes also the crop uptake), while in all other areas (FF fraction) they are reduced by a fixed rate, derived from the literature, before entering into the stream. Background losses for nitrogen are thus estimated as \(0.38 * FF * AtmN\). For an atmospheric deposition of \(10\) kgN/ha this corresponds to a background of \(3.8\) kgN/ha (in line with the values reported by Commission et al. (2004)).

In GREEN phosphorus model, the total annual phosphorus load \(L_i\) of catchment \(i\) is estimated by the equation: \[\begin{equation} L_i = (1 - Lret_i) * ((MinP_i + ManP_i + (1 - FF_i) * BgP_i) * (1 - Bret_i) + FF_i * BgP_i + \\ sd_{coef} * SdP_i + PsP_i + U_i) * (1 - Rret_i) \end{equation}\]

• \(MinP_i\) = annual amount of phosphorus mineral fertilisers (ton/yr)
• \(ManP_i\) = annual amount of phosphorus in manure fertilisers (ton/yr)
• \(FF_i\) = Non-agricultural land cover in the catchment (fraction)
• \(BgP_i\) = annual amount of phosphorus background losses (ton/yr)
• \(SdP_i\) = phosphorus input from scattered dwellings (ton/yr)
• \(PsP_i\) = phosphorus input from point sources (ton/yr)
• \(U_i\) = Phosphorus load from upstream catchments (ton/yr)

Note that phosphorous background losses are splited into two parts, with the inputs to agricultural land undergoing basin retention, while in all other areas (FF fraction) they are considered entering into the stream. Background losses for phosphorus are estimated at \(0.15\) kgP/ha (in line with the values reported by Commission et al. (2004)).

Assembling input data for running the GREEN model is time consuming. Data sources for application at European scale are described in Bruna Grizzetti et al. (2021). Domestic sources of nitrogen and phosphorus can also be found in Vigiak et al. (2020).