背景
目前已有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"(共享不易,欢迎引用)
github:https://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)