Interpreting the n_load parameter

Science, Research & Data
1

Hello INVEST community! :seedling:

I’m looking for clarification about the n_load parameter.

I’m a little confused and unsure if there is contradictory information in the NDR documentation. I interpreted n_load as the N input associated with each land use type. For example, for agricultural land, N input would be the amount of fertilizer, manure, and atmospheric deposition. Locally, we have agricultural N use data, and I have used N inputs as the n_load without applying a correction factor to account for on-pixel nutrient retention. With this approach, and with modeling proportions of subsurface N, I’ve gotten a pretty good fit when calibrating against empirical N export measurements (slope: 0.92, intercept: 8,800 kg). I based this approach/understanding on the below excerpt from page 131 of the INVEST documentation.

"Note 2: Load values may be expressed either as the amount of nutrient applied (e.g. fertilizer, livestock waste, atmospheric deposition); or as “extensive” measures of contaminants, which are empirical values representing the contribution of a parcel to the nutrient budget (e.g. nutrient export running off urban areas, crops, etc.) In the latter case, the load should be corrected for the nutrient retention from downstream pixels of the same LULC. For example, if the measured (or empirically derived) export value for forest is 3 kg.ha-1.yr-1 and the retention efficiency is 0.8, users should enter 15(kg.ha-1.yr-1) in the n_load column of the biophysical table; the model will calculate the nutrient running off the forest pixel (n_export) as 15(1-0.8) = 3 kg.ha-1.yr-1."*

However, I just read this excerpt from the bottom of pg. 126-127, which makes me believe that I should apply a correction factor to the n_load values I have used.

“Loads are the sources of nutrients associated with each pixel of the landscape. Consistent with the export coefficient literature (California Regional Water Quality Control Board Central Coast Region, 2013; Reckhow et al., 1980), load values for each LULC are derived from empirical measures of nutrient export (e.g., nutrient export running off urban areas, crops, etc.). If information is available on the amount of nutrient applied (e.g. fertilizer, livestock waste, atmospheric deposition), it is possible to use it by estimating the on-pixel nutrient use (and apply this correction factor to obtain the load parameters).”

Which is the correct approach? ​Reading similar posts from other users, I’m inclined to say the second excerpt is accurate, but I wanted to double check.

If I followed the wrong approach, is there any validity in what I’ve done, or will I need to reassign N loads and recalibrate the model?

1 Like
2

Hi @atershy -

Yes, I agree that those paragraphs give contradictory information. I’ve only used the NDR model a few times, and reading these made me doubt my interpretation. So I’m trying to get someone on our science team to verify which is correct, and I’ll update the user guide appropriately and get back to you. Apologies for the delay.

~ Stacie

3

Hi again @atershy -

Well, your question led to a lot of consternation here at NatCap. The colleagues I’ve talked with all have used applied values of fertilizer (and/or atmospheric deposition etc) as the n_load parameter. We have not adjusted these using the retention coefficient.

But the contradiction in the User Guide text makes us wonder if something different is actually going on in the code that we weren’t aware of. The person who created the model has moved on from NatCap, and we are in the process of looking at the code and figuring things out. Apologies for the confusion, but given that we’re all juggling multiple projects right now, it might take a while before we have a definitive answer.

In the meantime, I suggest using fertilizer application without correction, since that’s been standard practice here. It’s good to hear that this approach gives a pretty good fit with calibration. We’ll update this thread with our findings.

~ Stacie

4 Likes
4

Hi Stacie,
Have there been any updates to the information here?
Thanks,
Dave

5

Hi Stacie et al.,
I’d like to follow up on this thread by better understanding how nutrient loads for a given pixel are calculated. From my understanding of the documentation and the above discussion,

  1. the load for a given land use type can be calculated from the sum of the inputs (say, fert + manure + deposition).
  2. The modified load to a given pixel is the load x RPI to account for higher transport potential with higher precipitation (eq. 42 in the on-line documentation).
  3. NDR then estimates how much of that modified load reaches the stream by calculating the retention efficiency of that pixel and all downslope pixels on the flowpath to the stream. NDRo,i = 1 - eff’i (eq. 47). The NDRo,i is then modified by Borselli’s k and the topographic index IC, to account for hydrological connectivity, to get the NDRi for a given pixel (eq. 46.
  4. So, a given pixel’s export to the stream = LOADi x NDRi (Equation in Fig. 5).

What I’ve had difficulty understanding is that I don’t see anywhere that the LOADi to a given pixel is augmented by nutrients carried to it from upslope pixels. That is, say you have a forest pixel that is part of a riparian buffer, with an eff_N of 0.8, and a LOAD of 5 kg N/ha/yr from deposition. NDR would estimate that pixel’s export to be 5 * (1-0.8) = 1 kg N/ha/yr (for simplicity, here I’m not accounting for conversions of pixel sizes to ha, IC, k, and critical lengths longer than a pixel size, etc.). But say that forest buffer pixel has a crop pixel just upslope of it. The crop pixel has a LOAD of 100 kg N/ha/yr and eff_N = 0.5.

  1. NDR calculates the N export of the crop pixel as 100 (1 - eff_Nforest) = 20 kg N/ha/yr (from eq. 48, where the eff,down,i > eff,i), and the export from the forest pixel as 1 kg N/ha/yr, for a total of 21 kg N/ha/yr.

  2. But with the cropland only retaining 50% of the applied N, that means that 50 kg N/ha/yr are leaving the crop pixel and going to the forest pixel (assuming that all the N leaves in runoff and none in gaseous form). So, the true load to the forest pixel is 50 from runoff + 5 from deposition = 55 kg N/ha/yr. Assuming, as NDR does, that the forest pixel has a constant eff_N no matter how much N comes in (more on that later), export = 55 * (1 - 0.8) = 11 kg N/ha/yr. So, of the N load to the crop pixel, the total making it to the stream is LOADcrop*(1-EFFcrop)*(1-EFFforest) = 100 * 0.5 * 0.2 = 10 kg N/ha/yr, rather than 20.

  3. And, if we want to know the ecosystem service of N retention by the forest buffer, NDR tells us that the buffer has just absorbed 4 kg N/ha/yr [LOADforest * (1 - EFFforest)], but really it has absorbed 44 kg N/ha/yr.

  4. This gets more complicated if one assumes multiple crop pixels upslope of the single pixel of forest buffer. In that case, the first crop pixel would retain 50 kg N/ha/yr and export 50 to the downslope crop pixel, which would now have a load of 150 kg N/ha/yr, of which it would retain 75 and pass on 75 to the forest, whose load now becomes 80, of which it would retain 64 and pass on 16 to the stream. But NDR would calculate the exported N as 100 * (1 - 0.8) from the first crop pixel + 100*(1-0.8) from the second crop pixel + 5 * (1 - 0.8) from the forest pixel = 41 kg N/ha/yr, with the forest still only retaining 4, rather than 64.

  5. These discrepancies are important for our application, since we’re trying to estimate the ecosystem service provided by the forest buffer and test different buffer scenarios. Can we just not use NDR for this application?
    a. It has led to some counter-intuitive results, as @atershy pointed out in a recent post. We added more forested riparian buffer to predominantly agricultural watersheds, but ended up with less total N retention by forested buffers.
    b. How could this be - even with NDR’s formulation? More land area of forest should lead to more total N retention by forest, no matter what, correct?

  6. As an aside, it seems contrary to the idea of N saturation to allow a forest or cropland to have a constant eff_N, no matter what the load to it is. So, we’re trying to figure out if there is a way to decrement eff_N for a pixel if LOADi for that pixel exceeds a certain threshold (e.g., whatever the maximum absolute retention capacity of that LULC type is, in kg N/ha/yr). But maybe this can’t work in NDR, if NDR never augments loads to a pixel based on upslope runoff. Thoughts?

Sorry for the long post! Just trying to work through the nuts and bolts to better interpret our results and make sure the model is capable of doing what we’re asking of it.
Thanks,
Dave

1 Like
6

Hey @dhooper,

Thanks for taking the time and care to work through the NDR logic. It’ll take us a little bit to parse and get something satisfying back, but we will be taking a look. In the meantime, one of our software engineers logged the steps the code is taking as it pertains to load_n, back when this topic first was posted and I thought it might be helpful to share that in this thread.

…here’s what I’m seeing relating to the load_n parameter as it’s traced through the model. I’ll just list out the files that it affects below as a list of pixel-stack raster calculator operations:

load_n[lucode] * cell_area_ha => load_n.tif
load_n.tif * runoff_proxy_index.tif => modified_load_n.tif
modified_load_n.tif * (1 - proportion_subsurface_n[lucode]) => surface_load.tif
surface_load_n.tif * ndr_n.tif => n_surface_export.tif ( ndr_n.tif is not in any way affected by load_n)
n_surface_export.tif + n_subsurface_export.tif => n_total_export.tif

And since the subsurface calculations derive from modified_load_n.tif, which uses load_n, those are calculated as:

modified_load_n.tif * proportion_subsurface_n[lucode] => subsurface_load_n.tif
subsurface_load_n.tif * sub_ndr_n.tif => n_subsurface_export.tif

In short, I do not see any place in the code where load_n or its derivatives are used in a way that is different from what is defined in the user’s guide.

More soon!

7

Hi @dcdenu4 ,
Thanks for your quick reply and for looking further into this. The code you posted helps confirm what I understood from the documentation. In thinking through this more yesterday, and trying to explain it to my colleagues, here’s what I wrote to them:

What became clearer to me is how NDR splits its calculations into two parts, Loads and Retention efficiencies (whose inverse is nutrient delivery ratio):

  1. Nutrient load and modified load: the modified load only depends on the inputs (LOAD) from the biophysical table and how easily those get transported as represented by the runoff proxy index, which reflects a given pixel’s precipitation relative to average precip in the entire raster. That is, modified load of pixel i = LOADi x RPIi. It does not include any inputs from upslope pixels.

a… Conceptually, I find this approach a bit confusing for a couple of reasons:

i. the way I think about it, higher runoff doesn’t really change the LOAD of nutrients to any given pixel, so much as potentially change the efficiency of capture of that load (which, to my mind, would be in the second set of calculations on retention efficiencies; see below).

ii. I would think that “modified load” would reflect direct inputs (fert + manure + deposition + N-fixation) to a pixel based on its LULC type, as in the biophysical table, PLUS any N transported to that pixel by runoff from upslope pixels. But that apparently does not occur. The brief reply showing the code in response to my post confirms this.

  1. The block of calculations on Surface NDR (including eq. 46 & 47) relates to modification of the efficiency of N capture by a given pixel of a given LULC type. The documentation describes an iterative process for assigning eff’i, the maximum retention of N added to a given pixel by that pixel AND its downstream flow path. This concept, and how it is encapsulated in equation 48 of the documentation is what gave rise to my Excel calculations and post. That is, the topographic parameters (ICo, ICi, and k in eq. 47) don’t modify the nutrient load at all. They just modify the proportion of nutrients lost by a pixel to the stream via the downstream flow path (NDRi, in eq. 46).

a. A key feature of 48 is that if the retention efficiency of a downstream pixel (eff’downi) is greater than that of the pixel in question (eff_N from the biophysical table, = effLULCi in Eq. 48), the model seems to assume that the pixel in question then has that same retention efficiency – that is, all nutrients added to pixeli (say a crop pixel) will eventually be absorbed by downstream pixelj (say a buffer pixel) at the efficiency of the buffer pixelj. It does this rather than having pixeli absorb part of the nutrients added to it and then having the leftovers runoff and augment inputs to pixelj. This leads to the discrepancies in retention that I pointed out in my initial post to the conversation board.

  1. Ultimately, the model then calculates the nutrients from a given pixel that make it to the stream as “modified load” x NDRi (see the equation for pixel export on Fig. 5, and the code in the response to my post).

So, to my mind, this makes for a confusing conceptualization of watershed nutrient retention. It also seems to make the calculations of N retention by riparian buffers that we’re trying to do very difficult or impossible, because the only retention recorded for those pixels comes from the direct load to that LULC type from the biophysical table, even though the model attributes even more retention to them by applying their eff’i to the upslope crop pixel, as in eq. 48.

I hope that additional explanation helps to clarify some of the points that I found difficult to understand in the model implementation. In case it helps, I’ve also linked an Excel file working through the calculations in eq. 48 and 49 that helped me better understand how those approaches influenced the estimates of eff’i for a given pixel. If I’m wrong about any of this, please let me know!
Thanks!
Dave

8

Hi NDR folks,
I’m wondering if you have any thoughts on the issues raised above.

Thanks,
Dave

9

Hey @dhooper,

Thanks for nudging this. We’ve reach out to @Perrine, one of the model authors, and are waiting to hear back. Sorry for the delay!

Doug

10

Dear @dhooper (and colleagues),

Lots of interesting thoughts in this thread! This reminds me of some rich conversations we had about the best approach for the InVEST nutrient model.

A few points to respond to @dhooper’s last message. Point 4 is the most important one (with respect to your intended application):

  1. Modified loads: They were introduced in the model to reflect the fact that, all else being equal, an area with higher runoff would wash off more nutrients than an area with less runoff. This is because there is less opportunity for vegetation to take up nutrients (as discussed in Endreny et al. 2003: https://onlinelibrary.wiley.com/doi/10.1111/j.1752-1688.2003.tb01569.x)

Note: Modified loads do not reflect upslope pollutants. This is based on the conceptual model where we estimate the contribution of each pixel to the stream (inspired by the work of Borselli et al. for sediment transport, cf. InVEST’s sediment model). With this framework, if we added the load from upstream pixel j to the load of pixel i, we would double count load_j (which would be represented in pixel j AND pixel i)

  1. Load values (in response to @atershy): I confirm that we’ve used amount of nutrients applied (cf. also Redhead et al. 2018: https://www.sciencedirect.com.remotexs.ntu.edu.sg/science/article/pii/S0048969717320909)

  2. Calcs on surface NDR: That is correct, they don’t modify the loads, but rather the NDRsurf value. NDRsurf, in turn, is what modifies load_surf (as per equation 54)

  3. Retention of riparian buffers: In response to point 7 in your previous post (" And, if we want to know the ecosystem service of N retention by the forest buffer, NDR tells us that the buffer has just absorbed 4 kg N/ha/yr [LOADforest * (1 - EFFforest)], but really it has absorbed 44 kg N/ha/yr.")

→ The model does not estimate the retention of riparian buffers (or other landscape element). This is indeed an important limitation. Computationally, it would take another algorithm to calculate the contribution of each pixel to the retention (ie. NDR factor) of all upstream pixels flowing through that particular pixel. We’ve discussed that with the software team back in the days but I don’t think anyone had time to work on this since then.

→ This limitation is discussed in the User’s guide in the section “Evaluating Nutrient Retention Services” (copied below). I believe such an approach could be used to assess the retention of particular landscape elements, but I agree that it’s less ideal than having a model output representing the retention by a given pixel.

The NDR model does not directly quantify the amount of nutrient retained on the landscape. However, if you have scenarios that are being compared with current conditions, the nutrient retention service may be estimated by taking the difference in nutrient export between the scenario and current conditions. This quantifies the difference in nutrient reaching a stream, based on the changes in land cover/climate/etc present in the scenario, which provides a way of evaluating impacts to downstream uses such as drinking water.

To calculate per pixel nitrogen retention services within a single scenario, we recommend subtracting n_total_export.tif from the modified_load_n.tif result located in the intermediate output folder. Similarly, per pixel phosphorus retention services can be calculated by subtracting p_surface_export.tif from modified_load_p.tif . Use the .gpkg output to quantifty watershed scale nutrient retention services by subtracting the n_total_export result from (n_surface_load + n_subsurface_load ) for nitrogen and p_surface_export from p_surface_load for phosphorus.

I hope this helps!

5 Likes
11

Hello @atershy @dhooper and @Perrine -

NatCap science, software and analysts had a convening about NDR yesterday. The main thing that came out of it is the realization that we should be using export coefficients, not applied nutrients, for the n_load and p_load parameters.

This is because the model calculations currently do not take into account retention on the pixel itself, they only include retention downslope of each pixel. But in reality, if we apply nutrient to a pixel of cropland, the crop will use/retain some of the nutrient, and what remains will move downslope. Since this is not currently included in the model, we should be correcting for it when creating our n_load and p_load parameters. If you’ve read the previous posts, yes, that means that even much of NatCap’s own use of the model has used incorrect loading values, since we often used applied nutrient without correction.

So we’re doing two things to improve this situation:

  1. I have updated the User Guide to be correct and consistent with its guidance for how the model currently works, saying this:

Data sources may provide loading values as either the amount of applied nutrient (e.g. fertilizer, livestock waste, atmospheric deposition); or as “extensive” measures of contaminants, which are empirical values representing the contribution of a parcel to the nutrient budget (e.g. nutrient export running off urban areas, crops, etc.) In the case of having applied nutrient values, they should be corrected for the nutrient retention provided by the pixel itself, using the application rate and retention efficiency value (eff_n or eff_p) for that land cover type. For example, if the nitrogen application rate for an agricultural LULC class is 10 kg/ha/year, and the retention efficiency is 0.4, you should enter a value of 6.0 into the n_load column of the biophysical table. If you have “extensive”/nutrient export values, then you may use them directly in the biophysical table without correction.

  1. We will be updating the NDR model to allow users to enter either applied nutrient or export coefficient values for each LULC class in the biophysical table, with additional columns specifying whether each value represents applied nutrient or export coefficient. If the value is applied nutrient, the model will calculate (application rate x retention efficiency) as described above. If the value is export coefficient, then the model will use the value as is, without correction. The User Guide will also be updated.

I’m not sure how long the NDR update will take, but I will try to remember to post about it here when it is released.

Apologies for the confusion, and providing incorrect guidance. We’ve all learned from this one!

~ Stacie

4 Likes
12

@swolny I just edited your post of the UG excerpt with the UG edits you made yesterday, for anyone reading this thread in the future.

2 Likes