Custom function for rain energy

[claudiacagn]

Hello, for an application of the "rfactor" in a pilot area in Italy, we are interested in using the definitions of the rain energy per unit by: 1) Brown and Foster (1987) 2) McGregor et al. (1987) instead of the Flanders/Belgium-calibrated equation. I have defined the custom functions below, which should replace the built-in function "rain_energy_per_unit_depth". Is it possibile to integrate the custom functions in a new release of the library? Thank you in advance, Claudia

def rain_energy_per_unit_depth_brown_and_foster(rain): """Calculate rain energy per unit depth according to Brown and Foster.

Parameters
----------
rain : numpy.ndarray
    Rain (mm)

Returns
-------
energy : float
    Energy per unit depth.

Notes
-----
The rain energy per unit depth :math:`e_r` (:math:`\\text{J}.\\text{mm}^{-1}.
\\text{m}^{-2}`) is defined by [1] [2]:

.. math::

    e_r = 29*(1-0.72*exp(-0.05*i_r)

with

 - :math:`i_r` the rain intensity for every 10-min
   increment (mm :math:`\\text{h}^{-1}` ).

References
----------
.. [1] Brown, L.C., Foster, G.R., 1987. 
Storm erosivity using idealized intensity distributions. 
Transactions of the ASAE 30, 0379–0386. https://doi.org/10.13031/2013.31957.

.. [2] Renard, K.G., Foster, G.R., Weesies, G.A., McCool, D.K., Yoder, D.C., 1997, 
Predicting soil erosion by water: a guide to conservation planning with the revised universal soil loss equation (RUSLE), 
Agriculture Handbook. U.S. Department of Agriculture, Washington. 
https://www.ars.usda.gov/ARSUserFiles/64080530/RUSLE/AH_703.pdf

"""
rain_energy = 29 * (1 - 0.72 * np.exp(-0.05*rain))
return rain_energy.sum()

def rain_energy_per_unit_depth_mcgregor(rain): """Calculate rain energy per unit depth according to McGregor.

Parameters
----------
rain : numpy.ndarray
    Rain (mm)

Returns
-------
energy : float
    Energy per unit depth.

Notes
-----
The rain energy per unit depth :math:`e_r` (:math:`\\text{J}.\\text{mm}^{-1}.
\\text{m}^{-2}`) is defined by [1]:

.. math::

    e_r = 29*(1-0.72*exp(-0.08*i_r)

with

 - :math:`i_r` the rain intensity for every 10-min
   increment (mm :math:`\\text{h}^{-1}` ).

References
----------
.. [1] McGregor, K.C., Bingner, R.L., Bowie, A.J. and Foster, G.R., 1995. Erosivity index values for 
      northern Mississippi. Transactions of the ASAE, 38(4), pp.1039-1047.

"""
rain_energy = 29 * (1 - 0.72 * np.exp(-0.08*rain))
return rain_energy.sum()

[Sachagobeyn]

Hi @claudiacagn:

This looks quit well, after your submission we realized we need to do a refactor of the naming and the units before we can proceed to include this. See https://github.com/watem-sedem/rfactor/issues/82

In addition, you will need to define a unit test, so to make sure we can automatically validate your results.

[claudiacagn]

Hi @Sachagobeyn, here attached a test case test_case_accumoli.zip

Do you have a possible timeline for the enhancement and modification? Thank you! Claudia

[Sachagobeyn]

@claudiacagn : see https://github.com/watem-sedem/rfactor/pull/87 for functionality custom energy function, I didn't have time to check units in energy functions. I'll leave this open until https://github.com/watem-sedem/rfactor/issues/81 is fixed

[claudiacagn]

Hi @Sachagobeyn, I tested the new functionality in JupyterLab and it seems to work. Thank you! I'll wait also for the download version. Bests, Claudia

[Sachagobeyn]

Dear @claudiacagn ,

Could you confirm Brown and Foster + McGregor implementation (same units as Verstraeten):

def rain_energy_verstraeten(rain):
    """Calculate rain energy per unit depth according to Salles/Verstraeten.

    Parameters
    ----------
    rain : numpy.ndarray
        Rain (mm)

    Returns
    -------
    energy : float
        Energy per unit depth.

    Notes
    -----
    The rain energy per unit depth :math:`e_r` (:math:`\\text{MJ}.\\text{mm}^{-1}.
    \\text{ha}^{-1}`) for an application for Flanders/Belgium is defined
    by [1]_ , [2]_ and [3]_:

    .. math::

        e_r = 0.1112i_r^{0.31}

    with

     - :math:`i_r` the rain intensity for every 10-min
       increment (mm :math:`\\text{h}^{-1}` ).

    References
    ----------
    .. [1] Salles, C., Poesen, J., Pissart, A., 1999, Rain erosivity indices and drop
        size distribution for central Belgium. Presented at the General Assembly of
        the European Geophysical Society, The Hague, The Netherlands, p. 280.
    .. [2] Salles, C., Poesen, J., Sempere-Torres, D., 2002. Kinetic energy of rain and
        its functional relationship with intensity. Journal of Hydrology 257, 256–270.
    .. [3] Verstraeten, G., Poesen, J., Demarée, G., Salles, C., 2006, Long-term
        (105 years) variability in rain erosivity as derived from 10-min rainfall
        depth data for Ukkel (Brussels, Belgium): Implications for assessing soil
        erosion rates. Journal Geophysysical Research, 111, D22109.
    """
    rain_energy = 0.1112 * ((rain * 6.0) ** 0.31) * rain
    return rain_energy.sum()

def rain_energy_brown_and_foster(rain):
    """Calculate rain energy per unit depth according to Brown and Foster.

    Parameters
    ----------
    rain : numpy.ndarray
        Rain (mm)

    Returns
    -------
    energy : float
        Energy per unit depth.

    Notes
    -----
    The rain energy per unit depth :math:`e_r` (:math:`\\text{MJ}.\\text{mm}^{-1}.
    \\text{ha}^{-1}`) is defined by [4] and [5]:

    .. math::

        e_r = 0.29*(1-0.72*exp(-0.05*i_r)

    with

     - :math:`i_r` the rain intensity for every 10-min
       increment (mm :math:`\\text{h}^{-1}` ).

    References
    ----------
    .. [4] Brown, L.C., Foster, G.R., 1987. Storm erosivity using idealized intensity
        distributions. Transactions of the ASAE 30, 0379–0386.
        https://doi.org/10.13031/2013.31957.

    .. [5] Renard, K.G., Foster, G.R., Weesies, G.A., McCool, D.K., Yoder, D.C., 1997,
        Predicting soil erosion by water: a guide to conservation planning with the
        revised universal soil loss equation (RUSLE),  Agriculture Handbook. U.S.
        Department of Agriculture, Washington.
        https://www.ars.usda.gov/ARSUserFiles/64080530/RUSLE/AH_703.pdf
    """
    rain_energy = 0.29 * (1 - 0.72 * np.exp(-0.05 * rain))
    return rain_energy.sum()

def rain_energy_mcgregor(rain):
    """Calculate rain energy per unit depth according to McGregor.

    Parameters
    ----------
    rain : numpy.ndarray
        Rain (mm)

    Returns
    -------
    energy : float
        Energy per unit depth.

    Notes
    -----
    The rain energy per unit depth :math:`e_r` (:math:`\\text{MJ}.\\text{mm}^{-1}.
    \\text{ha}^{-1}`) is defined by [6]:

    .. math::

        e_r = 0.29*(1-0.72*exp(-0.08*i_r)

    with

     - :math:`i_r` the rain intensity for every 10-min
       increment (mm :math:`\\text{h}^{-1}` ).

    References
    ----------
    .. [6] McGregor, K.C., Bingner, R.L., Bowie, A.J. and Foster, G.R., 1995.
        Erosivity index values for northern Mississippi. Transactions of the ASAE,
         38(4), pp.1039-1047. 10.13031/2013.27921

    """
    rain_energy = 0.29 * (1 - 0.72 * np.exp(-0.08 * rain))
    return rain_energy.sum()

It is now implemented in https://github.com/watem-sedem/rfactor/pull/94

[claudiacagn]

Hi @Sachagobeyn, the implementation should be: def rain_energy_brown_and_foster(rain): rain_energy = 0.29 (1 - 0.72 np.exp(-0.3rain)) rain return rain_energy.sum()

def rain_energy_mcgregor(rain): rain_energy = 0.29 (1 - 0.72 np.exp(-0.48rain)) rain return rain_energy.sum()

[SethCallewaert]

Hi @claudiacagn, Can you provide links to the corresponding articles for these formulae? From the papers I have read, the above-mentioned formulae are correct, except for mcgregor where I find the value of e=0.29(1-0.72 exp(-0.082 rain). see: https://www.researchgate.net/publication/317383994_Rainfall_erosivity_An_historical_review#pf6 or https://www.google.com/url?sa=t&source=web&rct=j&opi=89978449&url=https://esdac.jrc.ec.europa.eu/public_path/shared_folder/WorkshopSeoul/32.%2520Joon-Hak%2520Lee%2520-%2520A%2520study%2520of%2520Rainfall%2520Kinetic%2520Energy%2520Equations%2520based%2520on%2520Climate%2520Classification.pdf&ved=2ahUKEwj9op_iwLWGAxWZUaQEHTdrCwgQFnoECBUQAQ&usg=AOvVaw3a3Eb2mcqwT3i_id7v7uqm or https://www.ideals.illinois.edu/items/96627/bitstreams/310834/data.pdf

[claudiacagn]

Hi @SethCallewaert, thank you for the comment. You are right, the original manuscript by McGregor et al. (1995), which is freely available here (https://www.researchgate.net/publication/279892111_Erosivity_Index_Values_for_Northern_Mississippi) reports a coefficient of -0.082 (see page 1043). However, some implementation codes use the rounded version of -0.08. For the fully-extended implementation of mcgregor, the code should look like (if I'm not wrong):

def rain_energy_mcgregor(rain): rain_energy = 0.29 (1 - 0.72 np.exp(-0.492rain)) rain return rain_energy.sum()

Looking for your feedback,

Claudia

[Sachagobeyn]

Dear @claudiacagn

Can you clarify why you use 0.492 as coefficient?

Kind regards,

[claudiacagn]

For my understanding of the source code, "rain" in the energy formula should be intended as an intensity, and, in case of 10min-interval rainfall data, it should be multiplied by 6 to have rain per hour, so that 0.082*6 = 0.492. When doing so, I obtained expected values of the R-factor when testing with the custom function. Was it wrong? Bests, Claudia

[Sachagobeyn]

Dear @claudiacagn,

Very well spotted, I agree: I suggest we apply it in this manner,

def rain_energy_brown_and_foster(rain):
    ...
    rain_energy = 0.29 * (1 - 0.72 * np.exp(-0.05 * rain*6))
    return rain_energy.sum()

def rain_energy_mcgregor(rain):
    ...
    rain_energy = 0.29 * (1 - 0.72 * np.exp(-0.08 * rain*6))
    return rain_energy.sum()

[claudiacagn]

I agree. Also, I suggest to make explicit in the name that the energy functions refer to 10min-interval data, for example by changing the name into: rain_energy_brown_and_foster_10min

Not sure if it makes sense though....

[Sachagobeyn]

I'll add it in the docs. I would rather check for including it as a a parameter, see https://github.com/watem-sedem/rfactor/issues/97

[Sachagobeyn]

I agree. Also, I suggest to make explicit in the name that the energy functions refer to 10min-interval data, for example by changing the name into: rain_energy_brown_and_foster_10min

Not sure if it makes sense though....

Also not we added * rain for`both formula's:

rain_energy = 0.29 * (1 - 0.72 * np.exp(-0.05 * rain*6))* rain

since we are computing the total storm kinetic energy and not the rainfall energy per unit depth, see https://github.com/watem-sedem/rfactor/pull/94/commits/2ff5af3c3a9c2d7d176ea08f7b672810d3adfe31

[Sachagobeyn]

@SethCallewaert (and @claudiacagn): can you both review changes in https://github.com/watem-sedem/rfactor/pull/94