潜在蒸散发计算(33个PET模型)

背景

目前已有100多个PET模型,但不同模型对PET估算及干旱风险评估的影响尚不明确。本研究选择在中国地区广泛应用的FAO56 Penman-Monteith (FAO56)、Priestley-Taylor (PT)、Thornthwaite (Th)、Hargreaves (Har)等模型(Li et al., 2022; Liu et al., 2022)以及适用于中国相似气候条件的Kharrufa (Kh)、Abtew (Ab)和Jensen-haise (JH)等模型(Shirmohammadi-Aliakbarkhani and Saberali, 2020; Zhao et al., 2021),共33个模型进行比较。其中,基于温度的模型7个,基于辐射的模型12个,基于质量传递的模型9个以及综合模型5个。以上模型在日尺度上估算PET(Th模型在月尺度),并通过累加至月尺度和年尺度。表1和补充S1详细介绍了这些模型的公式和输入,输入变量参考了FAO提供的计算方法(Singer et al., 2021)。

author = "Weiqi Liu(西北院刘伟琦)"
copyright = "Copyright (C) 2025 Weiqi Liu"
license = "NIEER"
version = "2025.07"
Reference paper == "Impact of the potential evapotranspiration models on drought monitoring"(共享不易,欢迎引用)

githubhttps://github.com/DroughtMonitor/Calculation-of-potential-evapotranspiration/tree/main

代码

变量准备

import aquacropeto as peto
import numpy as np
import xarray as xa
import matplotlib.pyplot as plt
import math
import calendar
import geopandas as geo
import salem
from sklearn.linear_model import LinearRegression
from sklearn.metrics import r2_score
from mpl_toolkits.basemap import Basemap
import scipy.stats as stats
import xskillscore as xs
import scipy.ndimage
from scipy import stats
import calendar
import pandas as pd


def avp_from_tdew(tdew: float) -> float:
    """
    Estimate actual vapour pressure (ea) from dew point temperature.

    Parameters
    ----------
    tdew : float
        Dew point temperature [°C]

    Returns
    -------
    float
        Actual vapour pressure [kPa]
    """
    return 0.6108 * math.exp((17.27 * tdew) / (tdew + 237.3))


def svp_from_t(t: float) -> float:
    """
    Estimate saturation vapour pressure (es) from air temperature.

    Parameters
    ----------
    t : float
        Air temperature [°C]

    Returns
    -------
    float
        Saturation vapour pressure [kPa]
    """
    return 0.6108 * math.exp((17.27 * t) / (t + 237.3))


def delta_svp(t: float) -> float:
    """
    Estimate the slope of the saturation vapour pressure curve at a given temperature.

    Parameters
    ----------
    t : float
        Air temperature [°C]

    Returns
    -------
    float
        Slope of saturation vapour pressure curve [kPa °C⁻¹]
    """
    es = svp_from_t(t)
    return (4098 * es) / ((t + 237.3) ** 2)


def psy_const(atmos_pres: float) -> float:
    """
    Calculate the psychrometric constant.

    Parameters
    ----------
    atmos_pres : float
        Atmospheric pressure [kPa]

    Returns
    -------
    float
        Psychrometric constant [kPa °C⁻¹]
    """
    return 0.000665 * atmos_pres


def atm_pressure(altitude: float) -> float:
    """
    Estimate atmospheric pressure from altitude.

    Parameters
    ----------
    altitude : float
        Elevation above sea level [m]

    Returns
    -------
    float
        Atmospheric pressure [kPa]
    """
    tmp = (293.0 - (0.0065 * altitude)) / 293.0
    return (tmp ** 5.26) * 101.3


def calc_lambda(tmean: float) -> float:
    """
    Calculate latent heat of vaporization (lambda).

    Parameters
    ----------
    tmean : float
        Mean daily air temperature [°C]

    Returns
    -------
    float
        Latent heat of vaporization [MJ kg⁻¹]
    """
    return 2.501 - 0.002361 * tmean


def wind_speed_2m(ws: float, z: float) -> float:
    """
    Convert wind speed measured at height z to wind speed at 2 meters above ground level.

    Parameters
    ----------
    ws : float
        Wind speed measured at height z [m s⁻¹]
    z : float
        Height of wind measurement above the ground [m]

    Returns
    -------
    float
        Wind speed adjusted to 2 meters height [m s⁻¹]
    """
    return ws * (4.87 / math.log((67.8 * z) - 5.42))


def clip_zeros(arr: np.ndarray) -> np.ndarray:
    """
    Replace negative values with 0.

    Parameters
    ----------
    arr : np.ndarray
        Input array (e.g., from Pandas Series or xarray DataArray)

    Returns
    -------
    np.ndarray
        Array with negative values clipped to 0
    """
    return np.where(arr < 0, 0, arr)

综合模型

import numpy as np

def fao56_penman_monteith(net_rad: float, t: float, ws: float,
                           svp: float, avp: float, delta_svp: float,
                           psy: float, shf: float = 0.0) -> float:
    """
    FAO-56 Penman-Monteith equation for reference evapotranspiration (ETo).

    Parameters
    ----------
    net_rad : float
        Net radiation at crop surface [MJ m⁻² day⁻¹]
    t : float
        Air temperature at 2 m height [°C]
    ws : float
        Wind speed at 2 m height [m s⁻¹]
    svp : float
        Saturation vapour pressure [kPa]
    avp : float
        Actual vapour pressure [kPa]
    delta_svp : float
        Slope of saturation vapour pressure curve [kPa °C⁻¹]
    psy : float
        Psychrometric constant [kPa °C⁻¹]
    shf : float, optional
        Soil heat flux [MJ m⁻² day⁻¹], by default 0.0

    Returns
    -------
    float
        Reference evapotranspiration [mm day⁻¹]
    """
    term1 = 0.408 * (net_rad - shf) * delta_svp / (delta_svp + psy * (1 + 0.34 * ws))
    term2 = 900 * ws * (svp - avp) * psy / ((t + 273) * (delta_svp + psy * (1 + 0.34 * ws)))
    pet = term1 + term2
    return clip_zeros(pet)


def penman(lambd: float, net_rad: float, wind: float, delta_svp: float,
           svp: float, avp: float, psy: float,
           aw: float = 1, bw: float = 0.537, ku: float = 6.43,
           shf: float = 0.0) -> float:
    """
    Standard Penman equation for potential evapotranspiration.

    Returns PET [mm day⁻¹]
    """
    fu = ku * (aw + bw * wind)
    denominator = lambd * (delta_svp + psy)
    term1 = delta_svp * (net_rad - shf) / denominator
    term2 = psy * (svp - avp) * fu / denominator
    pet = term1 + term2
    return clip_zeros(pet)


def FAO24_penman(lambd: float, net_rad: float, wind: float, delta_svp: float,
                 svp: float, avp: float, psy: float,
                 aw: float = 1, bw: float = 0.864, ku: float = 2.7,
                 shf: float = 0.0) -> float:
    """
    FAO-24 version of Penman equation.

    Returns PET [mm day⁻¹]
    """
    fu = ku * (aw + bw * wind)
    denominator = lambd * (delta_svp + psy)
    term1 = delta_svp * (net_rad - shf) / denominator
    term2 = psy * (svp - avp) * fu / denominator
    pet = term1 + term2
    return clip_zeros(pet)


def FAO_ppp17_Penman(lambd: float, net_rad: float, wind: float,
                     delta_svp: float, svp: float, avp: float,
                     psy: float, Tmax: float, Tmin: float,
                     aw: float = 1, ku: float = 6.43, shf: float = 0.0) -> float:
    """
    FAO PPP17 version of Penman with temperature-dependent wind coefficient.

    Returns PET [mm day⁻¹]
    """
    delta_T = Tmax - Tmin
    if delta_T < 12:
        bw = 0.54
    else:
        bw = 0.54 + 0.35 * (delta_T - 12) / 4

    fu = ku * (aw + bw * wind)
    denominator = lambd * (delta_svp + psy)
    term1 = delta_svp * (net_rad - shf) / denominator
    term2 = psy * (svp - avp) * fu / denominator
    pet = term1 + term2
    return clip_zeros(pet)


def kimberly_penman(ws: float, net_rad: float, jj: int, avp: float,
                    svp: float, lambd: float, delta_svp: float,
                    psy: float, g: float = 0.0, ku: float = 2.62) -> float:
    """
    Kimberly Penman equation.

    Parameters
    ----------
    jj : int
        Day of year (DOY)
    
    Returns PET [mm day⁻¹]
    """
    # Empirical seasonal wind function
    w = ((0.4 + 1.4 * np.exp(-((jj - 173) / 58) ** 2)) +
         (0.605 + 0.345 * np.exp(-((jj - 243) / 80) ** 2))) * ws

    denominator = lambd * (delta_svp + psy)
    term1 = delta_svp * (net_rad - g) / denominator
    term2 = ku * psy * (svp - avp) * w / denominator
    pet = term1 + term2
    return clip_zeros(pet)

基于温度模型

import numpy as np
import calendar

# Default number of days in each month
_MONTHDAYS = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
_LEAP_MONTHDAYS = [31, 29, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]


def thornthwaite(monthly_t, monthly_mean_dlh, year=None):
    """
    Estimate monthly PET using Thornthwaite (1948) method.

    Parameters
    ----------
    monthly_t : list or np.ndarray
        Mean daily temperature for each month [°C]
    monthly_mean_dlh : list or np.ndarray
        Mean daily daylight hours for each month [hours]
    year : int, optional
        Year (used for leap year check), default assumes non-leap year

    Returns
    -------
    list
        Monthly potential evapotranspiration [mm/month]
    """
    if len(monthly_t) != 12:
        raise ValueError(f'monthly_t should have 12 values, got {len(monthly_t)}')
    if len(monthly_mean_dlh) != 12:
        raise ValueError(f'monthly_mean_dlh should have 12 values, got {len(monthly_mean_dlh)}')

    month_days = _LEAP_MONTHDAYS if (year and calendar.isleap(year)) else _MONTHDAYS

    adj_t = [max(0, t) for t in monthly_t]  # Negative temps set to 0

    I = sum((t / 5.0) ** 1.514 for t in adj_t if t > 0)
    a = (6.75e-7 * I**3) - (7.71e-5 * I**2) + (1.792e-2 * I) + 0.49239

    pet = []
    for Ta, L, N in zip(adj_t, monthly_mean_dlh, month_days):
        pet_month = 1.6 * (L / 12.0) * (N / 30.0) * ((10.0 * Ta / I) ** a) * 10.0
        pet.append(pet_month)
    return pet


def blaney_criddle(tmean, py, a=-1.55, b=0.96):
    """
    Estimate PET using Blaney-Criddle method.

    Parameters
    ----------
    tmean : float or np.ndarray
        Mean daily air temperature [°C]
    py : float or np.ndarray
        Relative daylength ratio (fraction of max annual)
    a : float
        Calibration coefficient
    b : float
        Calibration coefficient

    Returns
    -------
    float or np.ndarray
        PET [mm/day]
    """
    pet = a + b * (py * (0.457 * tmean + 8.128))
    return clip_zeros(pet)


def kharrufa(tmean, py):
    """
    Estimate PET using Kharrufa method.

    Parameters
    ----------
    tmean : float or np.ndarray
        Mean daily temperature [°C]
    py : float or np.ndarray
        Relative daylength ratio

    Returns
    -------
    float or np.ndarray
        PET [mm/day]
    """
    pet = 0.34 * py * np.power(np.maximum(tmean, 0), 1.3)
    return clip_zeros(pet)


def hamon(tmean, dl):
    """
    Estimate PET using Hamon method.

    Parameters
    ----------
    tmean : float or np.ndarray
        Mean daily temperature [°C]
    dl : float or np.ndarray
        Daylight hours [h]

    Returns
    -------
    float or np.ndarray
        PET [mm/day]
    """
    pet = np.power(dl / 12, 2) * np.exp(tmean / 16)
    return clip_zeros(pet)


def linacre(tmean, elevation, lat, tdew):
    """
    Estimate PET using Linacre method.

    Parameters
    ----------
    tmean : float or np.ndarray
        Mean temperature [°C]
    elevation : float
        Site elevation [m]
    lat : float
        Latitude [°]
    tdew : float or np.ndarray
        Dew point temperature [°C]

    Returns
    -------
    float or np.ndarray
        PET [mm/day]
    """
    th = tmean + 0.006 * elevation
    pet = (500 * th / (100 - lat) + 15 * (tmean - tdew)) / (80 - tmean)
    return clip_zeros(pet)


def romanenko(tmean, rh):
    """
    Estimate PET using Romanenko method.

    Parameters
    ----------
    tmean : float or np.ndarray
        Mean daily temperature [°C]
    rh : float or np.ndarray
        Relative humidity [%]

    Returns
    -------
    float or np.ndarray
        PET [mm/day]
    """
    pet = 4.5 * ((1 + tmean / 25) ** 2) * (1 - rh / 100)
    return clip_zeros(pet)


def schendel(tmean, rh):
    """
    Estimate PET using Schendel method.

    Parameters
    ----------
    tmean : float or np.ndarray
        Mean daily temperature [°C]
    rh : float or np.ndarray
        Relative humidity [%]

    Returns
    -------
    float or np.ndarray
        PET [mm/day]
    """
    pet = 16 * tmean / rh
    return clip_zeros(pet)

基于辐射模型

import numpy as np

def priestley_taylor(rn, dlt, gamma, lambd, alpha=1.26, g=0):
    """Priestley-Taylor method for PET estimation.

    Parameters
    ----------
    rn : float or array
        Net radiation [MJ m-2 d-1]
    dlt : float or array
        Slope of saturation vapour pressure curve [kPa °C⁻¹]
    gamma : float or array
        Psychrometric constant [kPa °C⁻¹]
    lambd : float or array
        Latent heat of vaporization [MJ kg⁻¹]
    alpha : float, optional
        Priestley-Taylor coefficient [-]
    g : float or array, optional
        Soil heat flux [MJ m-2 d-1]

    Returns
    -------
    pet : float or array
        Potential evapotranspiration [mm d⁻¹]
    """
    pet = (alpha * dlt * (rn - g)) / (lambd * (dlt + gamma))
    return clip_zeros(pet)

def makkink(rs, dlt, gamma, lambd, alpha=0.61, b=-0.012):
    """Makkink method for PET estimation."""
    pet = alpha * dlt * rs / ((dlt + gamma) * lambd) + b
    return clip_zeros(pet)

def jensen_haise(tmean, rs, lambd, cr=0.025, tx=-3):
    """Jensen-Haise method for PET estimation."""
    pet = rs / lambd * cr * (tmean - tx)
    return clip_zeros(pet)

def turc(tmean, rs, rh, k=0.013):
    """Turc method for PET estimation (converted radiation units)."""
    rs = rs * 23.885  # MJ/m2/day -> cal/cm2/day
    c = np.where(rh >= 50, 1, 1 + (50 - rh) / 70)
    pet = k * c * tmean * (rs + 50) / (tmean + 15)
    pet = np.where(tmean <= 0, 0, pet)
    return clip_zeros(pet)

def mcguinness_bordne(tmean, ra, lambd, k=0.0147):
    """McGuinness-Bordne method for PET estimation."""
    pet = k * ra * (tmean + 5) / lambd
    return clip_zeros(pet)

def abtew(tmax, rs, lambd, k=0.01786):
    """Abtew method for PET estimation."""
    pet = k * rs * tmax / lambd
    return clip_zeros(pet)

def doorenbos_pruitt(rs, delta, rh, wind, psy, lambd):
    """Doorenbos-Pruitt method for PET estimation."""
    b = -0.3
    a = (1.066 - 0.0013 * rh - 0.0002 * rh * wind +
         0.045 * wind - 0.0000315 * rh ** 2 - 0.0011 * wind ** 2)
    pet = a * (delta / (delta + psy)) * rs / lambd + b
    return clip_zeros(pet)

def irmak(rs, tmean):
    """Irmak method for PET estimation."""
    pet = -0.611 + 0.149 * rs + 0.079 * tmean
    return clip_zeros(pet)

def oudin(tmean, ra, lambd, k1=100, k2=5):
    """Oudin method for PET estimation."""
    pet = ra * (tmean + k2) / lambd / k1
    pet = np.where((tmean + k2) > 0, pet, 0)
    return clip_zeros(pet)

def hargreaves(tmin, tmax, tmean, et_rad):
    """Hargreaves method for PET estimation."""
    pet = 0.0023 * (tmean + 17.8) * np.sqrt(tmax - tmin) * 0.408 * et_rad
    return clip_zeros(pet)

def droogers_allen(tmin, tmax, tmean, et_rad):
    """Droogers-Allen variation of Hargreaves method."""
    pet = 0.0025 * (tmean + 16.8) * np.sqrt(tmax - tmin) * 0.408 * et_rad
    return clip_zeros(pet)

def allen(tmin, tmax, tmean, et_rad):
    """Allen variation of Hargreaves method."""
    pet = 0.003 * (tmean + 20) * ((tmax - tmin) ** 0.4) * 0.408 * et_rad
    return clip_zeros(pet)

def dorji(tmin, tmax, tmean, et_rad):
    """Dorji method for PET estimation."""
    pet = 0.002 * (tmean + 33.9) * ((tmax - tmin) ** 0.296) * 0.408 * et_rad
    return clip_zeros(pet)

基于质量传递模型

def Dalton(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Dalton equation.

    Parameters
    ----------
    wind : float or array-like
        Wind speed at 2 m height [m s-1].
    es : float or array-like
        Saturation vapour pressure [kPa].
    ea : float or array-like
        Actual vapour pressure [kPa].

    Returns
    -------
    pet : float or array-like
        Reference evapotranspiration [mm day-1].
    """
    pet = (0.3648 + 0.07223 * wind) * (es - ea)
    return clip_zeros(pet)


def Trabert(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Trabert equation.
    """
    pet = 0.3075 * (wind ** 0.5) * (es - ea)
    return clip_zeros(pet)


def Meyer(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Meyer equation.
    """
    pet = (0.375 + 0.05026 * wind) * (es - ea)
    return clip_zeros(pet)


def Rohwer(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Rohwer equation.
    """
    pet = 0.44 * (1 + 0.27 * wind) * (es - ea)
    return clip_zeros(pet)


def Penman_mass(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Penman mass transfer method.
    """
    pet = 0.35 * (1 + 0.98 / (100 * wind)) * (es - ea)
    return clip_zeros(pet)


def Albrecht(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Albrecht equation.
    """
    pet = (0.1005 + 0.297 * wind) * (es - ea)
    return clip_zeros(pet)


def Brockamp_Wenner(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Brockamp-Wenner equation.
    """
    pet = 0.543 * (wind ** 0.456) * (es - ea)
    return clip_zeros(pet)


def WMO(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the WMO equation.
    """
    pet = (0.1298 + 0.0934 * wind) * (es - ea)
    return clip_zeros(pet)


def Mahringer(wind, es, ea):
    """
    Estimate reference evapotranspiration (ETo) using the Mahringer equation.
    """
    pet = 0.15072 * ((3.6 * wind) ** 0.5) * (es - ea)
    return clip_zeros(pet)

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值