基于Penman-Monteith公式的潜在蒸散发计算

1948年，彭曼将能量平衡与传质方法结合起来，根据日照、温度、湿度和风速的标准气候记录，导出了一个计算开放水面蒸发量的方程式。这个所谓的组合后，许多研究人员进一步发展了该方法，并通过引入阻力因子将其推广到了作物的面应用。对各种计算方法性能的分析表明，需要制定一个标准的计算潜在蒸散发 ETo 的方法。推荐使用粮农组织 Penman-Monteith 方法作为唯一标准方法。它是一种在广泛的位置和气候下极有可能正确预测 ETo 的方法，并且可以在数据缺失的情况下应用。

组合公式的 Penman-Monteith 形式为：
$$E{T}_0=\frac{0.408\Delta \ast (R_n-G)}{\Delta +\gamma \ast (1+0.34u_2)}+\frac{\frac{900}{T+273}+\gamma \ast u_2\ast (e_s-e_a)}{\Lambda +\gamma (1+0.34u_2)}$$
式中$$R_n$$为净辐射($$MJ\ast m^{-2}\ast d^{-1}$$):$$G$$为土壤热通量($$MJ\ast m^{-2}\ast d^{-1}$$),计算逐日蒸散发量时忽略：$$T$$为日平均温度($${}^{0}C$$)；$$u_2$$为2m高度处风速($$m\ast s^{-1}$$);$$es$$和$$ea$$分别为饱和水汽压与实际水汽压($$kPa$$);$$\Delta$$为饱和水汽压—温度曲线斜率($$kPa\ast {}^o{C}^{-1}$$)；$$\gamma$$为温度计常数($$kPa\ast {}^o{C}^{-1}$$)。采用上式可计算逐日$$E{T}_0$$。

先根据上式设计计算ET0的函数，其中Rn、u2、es、ea以及△需要再设计函数进行计算。

def pm_fao56(tmean, wind, rs=None, rn=None, g=0, tmax=None, tmin=None,
             rhmax=None, rhmin=None, rh=None, pressure=None, elevation=None,
             lat=None, n=None, nn=None, rso=None, a=1.35, b=-0.35, ea=None,
             albedo=0.23, kab=None, as1=0.25, bs1=0.5, clip_zero=True):
    """
    计算潜在蒸散量的函数

    输入参数：
    tmean:日平均气温 [°C]
    wind:日平均风速 [m/s]
    rs: 入射太阳辐射 [MJ/ ((m^2)*d)]
    rn:净辐射 [MJ m-2 d-1]
    g: 土壤热通量 [MJ m-2 d-1]
    tmax:日最高气温 [°C]
    tmin:日最低气温 [°C]
    rhmax: 日最大相对湿度 [%]
    rhmin: 最小日相对湿度 [%]
    rh: 日平均相对湿度 [%]
    pressure:大气压[kPa]
    elevation:观测点高程 [m]
    lat: 观测点纬度 [rad]
    n: 实际日照时间 [hour]
    nn: 最大可能日照时间或日照时数 [hour]
    rso: 晴空太阳辐射 [MJ m-2 day-1]
    a: 净长波辐射的经验系数 [无单位]
    b: 净长波辐射的经验系数 [无单位]
    ea: 实际水蒸气压 [kPa]
    albedo: 地表反照率 [无单位]
    kab:从as1，bs1推导出用于估算晴空辐射的系数  [degrees].
    as1: 回归常数，表示阴天时地外辐射到达地球的比例 (n = 0) [无单位]
    bs1: 地外辐射经验系数 [无单位]
    clip_zero: 如果为True，则将所有负值替换为0

    输出参数：潜在蒸散发 [mm/d].

    """
    if tmean is None:
        tmean = (tmax + tmin) / 2
    pressure = calc_press(elevation, pressure)
    gamma = calc_psy(pressure)
    dlt = calc_vpc(tmean)

    gamma1 = (gamma * (1 + 0.34 * wind))

    ea = calc_ea(tmean=tmean, tmax=tmax, tmin=tmin, rhmax=check_rh(rhmax),
                 rhmin=check_rh(rhmin), rh=check_rh(rh), ea=ea)
    es = calc_es(tmean=tmean, tmax=tmax, tmin=tmin)

    rn = calc_rad_net(tmean, rn, rs, check_lat(lat), n, nn, tmax, tmin, rhmax,
                      rhmin, rh, elevation, rso, a, b, ea, albedo, as1, bs1,
                      kab)

    den = dlt + gamma1
    num1 = (0.408 * dlt * (rn - g)) / den
    num2 = (gamma * (es - ea) * 900 * wind / (tmean + 273)) / den
    pet = num1 + num2
    pet = clip_zeros(pet, clip_zero)
    return pet

函数pm_fao56（)中calc_press()、calc_psy()、calc_vpc()、calc_ea()、 calc_es()以及calc_rad_net()是分别计算站点大气压强、湿度常数、空气温度下饱和蒸汽压曲线的斜率、实际水汽压、饱和水汽压以及净辐射的函数，clip_zeros()函数的作用是将负值设置为0。这些函数的设计不再赘述。

实例12：下面基于上述函数利用郑州市2016年的包含日平均温度、日最高温度、日最低温度以及太阳辐射等信息的数据进行郑州市2016年的潜在蒸散发估算。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from PemanFAO56 import PMcombination as PM
data = pd.read_csv("E:\ETProject3(方法分开)\PMFAO56_to_ET\data\data1\zz.txt", parse_dates=True, index_col="date")
et0_coagmet = data["et_asce0"]
meteo = pd.DataFrame({"tmean":data["tavg"],
                      "tmax":data["tmax"],
                      "tmin":data["tmin"],
                      "rhmax":data["rhmax"]*100,
                      "rhmin":data["rhmin"]*100,
                      "u2":data["windrun"]*1000/86400,
                      "rs":data["solar"]*86400/1000000})

tmean, tmax, tmin, rhmax, rhmin, wind, rs = [meteo[col] for col in meteo.columns]
lat = 34.37 * np.pi / 180

et0_fao56 = PM.pm_fao56(tmean, wind, rs=rs, elevation=1138, lat=lat,
                         tmax=tmax, tmin=tmin, rhmax=rhmax, rhmin=rhmin)

plt.figure()
matplotlib.rcParams['font.sans-serif'] = ['SimHei']
x=range(1,367)
y=et0_fao56
plt.plot(x,y,color='red')
plt.xlabel('Date day')
plt.ylabel('ET0 mm/d')
plt.title('ET0 by PM')
plt.show()

结果如下：