PET calculation with the Hargreaves Modificated method (Droogers and Allen, 2002)

[JavierDiezSierra]

Generic Issue

• xclim version: 0.39.0
• Python version: 3.8.13
• Operating System: Ubuntu 22.04.02

Description

In the context of the Copernicus Climate Atlas (https://atlas.climate.copernicus.eu/atlas) we (https://github.com/SantanderMetGroup) are calculating the SPEI with different PET (Potential Evapotranspiration) methods.

We are using xclim during the index computation workflow to ensure process traceability and also because we want to use tools that are aligned with the Copernicus roadmap.

The xclim.indices.potential_evapotranspiration() function includes several methods that we are testing, comparing the results with our in-house functions. Currently, we have validated two of the available methods ("thornthwaite48" and "hargreaves85"), and we would like to implement the modificated "hargreaves85" method (Droogers and Allen 2002; https://doi.org/10.1023/A:1015508322413), which includes precipitations as a proxy for humidity (when precipitation occurs, the air contains humidity, thus reducing PET).

If you are interested in including this new method, we can submit a pull request and provide you with the documentation. Bellow, you'll find our in-house function to calculate the modified method which uses precipitation. This function also implements the "hargreaves85" method because we noticed that the formulation used in xclim is not exactly the same as the one available in the original paper (Hargreaves and Zohrab A. Samani 1985). The xclim formulation (

xclim/xclim/indices/_conversion.py

Line 1285 in bb39cb1

def potential_evapotranspiration(

) includes a parameter (lv = 2.5), which differs from the formulation in Hargreaves et al (2002). The inverse of “lv” (1 / 2.5 = 0.4) coincides approximately with the constant 0.408 in Hargreaves et al. (2002)]. This constant is used to convert radiation to evaporation equivalents in mm (kg/MJ).

Many thanks for your help!

def potential_evapotranspiration(tasmin, tasmax, pr = None, method = 'hg85'):
    """
    Calculate the Potential Evapotranspiration (PET) index using either the HG (Hargreaves 1985) or the MHG (Droogers and Allen 2002) method.

    Parameters
    ----------
    tasmin : xr.Dataset
        Input dataset with daily or monthly 'tasmin'.
    tasmax : xr.Dataset
        Input dataset with daily or monthly 'tasmax'.
    pr : xr.Dataset
        Input dataset with daily or monthly 'pr' for the MHG method. Note that for the MHG method, 
        the suitable temporal scale is monthly. Results might differ when calculated at daily or monthly aggregations 
        (daily-to-monthly aggregated results are NOT identical to monthly ones)
    methods : string
        if 'hg85' then the HG method (Hargreaves 1985) is applied
        if 'mhg' then the MHG method (Droogers and Allen, 2002) is applied
    Returns
    -------
    xarray.Dataset
        Dataset containing the calculated Potential Evapotranspiration index.
    """

    # Extracting temperature variables from the dataset
    tasmin = convert_units_to(tasmin, "degC")
    tasmax = convert_units_to(tasmax, "degC")

    tas = (tasmin + tasmax) / 2        
    tr = tasmax - tasmin
    tr = tr.where(tr>0, 0)
    ra = compute_extraterrestrial_solar_radiation(tasmin, "MJ m-2 d-1")
    
    ra = ra * 0.408 # Is used to convert the radiation to evaporation equivalents in mm (kg/MJ)

    if method == 'mhg':
        # Hargreaves Modificado 2002; Droogers and Allen, 2002 (precipitation in mm/month)
        pr = convert_units_to(pr, "mm/month")
        ab = tr - 0.0123 * pr
        pet_index = (
            0.0013 * ra * (tas + 17.0) * ab ** 0.76
        )
        pet_index = xr.where(np.isnan(ab**0.76), 0, pet_index)

    elif method == 'hg85':
        # Hargreaves et al. (1985) and Hargreaves (1994)
        pet_index = 0.0023 * ra * (tas + 17.8) * tr ** 0.5

    pet_index = pet_index.clip(0)
    pet_index = pet_index.to_dataset(name="pet")
    pet_index['pet'].attrs['units'] = 'mm'
        
    return pet_index

Code of Conduct

• I agree to follow this project's Code of Conduct

[aulemahal]

Hi @JavierDiezSierra,
Yes we are interested! I am quite happy to hear that xclim is being part of such a project.

We will happily welcome a PR to add this new method and the other fix to the PET computation. The transition from your function to xclim's implementation looks relatively easy to do. Thanks!

[Zeitsperre]

Hi @JavierDiezSierra, it's always great to see more and more research groups making use of xclim!

For the modified Hargreaves constant, would it make more sense to propose that it be called via droogersallen02/da02?

If I understand correctly, the hg85 calculation you propose here is more precise than the one currently implemented. Would it also make sense to modify our existing method?

We can help with getting the code and documentation consistent once the PR is opened. Thanks again!

[JavierDiezSierra]

@Zeitsperre and @aulemahal thank you very much for responding so quickly and for showing interest in including this new method.

@Zeitsperre, I think the codes (droogersallen02/da02) you propose for the new method are perfect. This method includes precipitation and is more accurate than the hg85 for arid climates. However, I think it's interesting to have both available in xclim in case precipitation data is not available. I will make a pull request throughout next week and help you document the method (input units, temporal aggregation, etc).

@aulemahal, I agree that the transition from my function to xclim's implementation is relatively easy to do ;-)

Many many thanks!

[JavierDiezSierra]

Hi @aulemahal,

I have just created the pull request: JavierDiezSierra:da02