彭曼公式推导与Python实现：从能量平衡到蒸散发计算

简介彭曼公式完整推导文档面向水文学、气象学及相关专业学生与研究者旨在系统梳理开放水面蒸发量估算核心公式的数学逻辑与物理背景。文档从能量平衡原理出发逐步推导KL-HElvrw方程详细说明短波净辐射、长波净辐射、对流热交换及蒸发潜热等项的含义并引入饱和水蒸气压斜率Δ来规避水温数据缺失问题推导过程严谨完整。资源为单个docx文件大小仅24KB内容精炼、公式排版清晰适合对照学习或教学备课使用。目前已有3736人学习得到学习者认可。通过学习该文档读者可深入理解Penman方程的来源与适用前提掌握HKHva(Ts-Ta)到Δ代换的关键技巧并熟悉Magnus经验公式等辅助计算方法有助于后期直接应用或进一步研究改进模型。1. 彭曼公式推导从能量平衡到可计算的蒸散发量在灌溉预报和作物生长模拟里参考蒸散量是避不开的输入。很多人直接调用现成库可一旦遇到小时尺度、传感器高度或高海拔气压公式里的每个权重就变得说不清。彭曼公式能同时吸收净辐射和空气干燥度对蒸发的影响而完整推导过程的价值在于让你看清Δ、γ、风函数是怎么挤进同一个式子的。下面从能量平衡和水汽扩散两个基本方程出发逐步消去表面温度得到经典彭曼公式然后用Python实现并给出验证方法。适合做水文模拟、农业气象和气候影响评估的工程师。2. 先摆证据能量平衡方程与空气动力学方程的物理基础彭曼公式的推导只需要两个方程一个是地表能量平衡另一个是湍流扩散下的水汽通量方程。它们单独看都简单合起来却会引入一个表面温度Ts这是推导中最容易绕晕的地方。下面先分别建立方程再统一变量和单位。2.1 能量平衡净辐射到底分给了谁蒸发面接收的净辐射Rn去向只有三个加热空气的显热通量H蒸发耗掉的潜热通量λE以及向下传入土壤或水体的热通量G。于是Rn H λE GE是蒸发速率单位kg/(m²·s)λ是水的蒸发潜热单位J/kg因此λE单位是W/m²。仪器能直接测到Rn和G但H和λE很难分开测这正是彭曼公式要解决的联立能量方程和水汽方程把λE单独解出来。日尺度计算中G通常可以取0因为白天的加热和夜间的冷却在24小时里基本抵消。但推导时不急着扔掉它否则后面线性化会少一项。小时尺度或裸土条件下G可能达到Rn的20%到30%那时必须单独估算。常见做法是把它建模为土壤热容和表层温度差的乘积或者在裸土场景用实测热流板数据但彭曼公式推导本身不关心G的取值只要求它作为已知输入参与计算。2.2 空气动力学方程水汽和热量都靠湍流输送水面或湿润冠层持续蒸发出水汽这些水汽必须由湍流带入大气否则蒸发面附近会迅速饱和。空气动力学方法用一个阻力ra来描述这个过程的“难易程度”ra越大水汽越难离开蒸发面。潜热通量写成λE (ρCp / γ) * (es(Ts) - ea) / ra感热通量是同一套湍流机制在输送热量所以也使用同一个raH ρCp * (Ts - Ta) / ra两个式子都依赖表面温度Ts。这里es(Ts)是表面温度对应的饱和水汽压ea是参考高度的实际水汽压。由于蒸发面湿润水汽压可以直接取饱和值。对实际作物冠层来说气孔会额外施加一个表面阻力这个点放到第5章再讨论经典彭曼推导先假设表面完全湿润。两个方程联立后Ts是未知量。综合气象站一般只观测参考高度气温Ta表面温度不直接可得。所以推导的主要工作就是找到一种方式把Ts从方程里消去这就是第3章线性化要处理的问题。符号含义典型单位Rn净辐射通量密度W/m²日尺度常用MJ/(m²·day)G土壤或水体热通量W/m² 或 MJ/(m²·day)H显热通量W/m²λE潜热通量W/m²ρ空气密度kg/m³Cp定压比热J/(kg·°C)γ干湿表常数kPa/°Cra空气动力学阻力s/mTs蒸发面温度°CTa参考高度气温°Ces(Ts)表面温度下的饱和水汽压kPaea参考高度实际水汽压kPa2.3 从单位体系看为什么需要干湿表常数干湿表常数γ是连接温度梯度和水汽压梯度的桥梁定义为γ Cp * P / (0.622 * λ)标准大气压101.3 kPa、气温20°C附近λ取2.45 MJ/kgγ约为0.066 kPa/°C。温度变化对γ的影响不超过5%所以很多实现直接把它固定为常数只有高海拔站点才需要用实际气压重算。下面一段Python会贯穿后续所有实现先单独写出来import math def saturation_vapor_pressure(T): # T: 气温单位°C # 返回饱和水汽压单位kPa return 0.6108 * math.exp(17.27 * T / (T 237.3)) def slope_sat_vapor_pressure(T): # 饱和水汽压曲线在T处的斜率单位kPa/°C es saturation_vapor_pressure(T) return 4098 * es / (T 237.3) ** 2 def psychrometric_constant(P101.3): # P: 气压kPa缺省为标准大气压 # 返回干湿表常数单位kPa/°C lam 2.45 # MJ/kg cp 1.013e-3 # MJ/(kg·°C) return cp * P / (0.622 * lam)Tetens公式在-40到50°C范围内精度足够气象站平均气温基本不会超出。slope函数直接用解析导数而不是数值差分避免了步长选择的问题。psychrometric_constant中Cp和λ都使用MJ单位所以输出自动是kPa/°C。后续所有实现都会从这三个函数里取数。3. 核心推导线性化饱和水汽压并消去表面温度第2章的方程里有一个非线性项es(Ts)它是Ts的指数函数。如果不做任何近似两个方程虽然也能数值求解但需要迭代还会面临初始值选取的问题。彭曼的巧妙之处是在气温Ta附近对饱和水汽压曲线做一阶展开把一个隐式问题变成显式公式。下面按步骤展开。3.1 为什么要对饱和水汽压做线性化在Ta处对es(T)做一阶泰勒展开es(Ts) ≈ es(Ta) Δ(Ts - Ta)Δ是饱和水汽压曲线在Ta处的斜率也就是第2章代码里的slope_sat_vapor_pressure。它的物理含义是气温每升高1°C饱和水汽压增加多少kPa在0到40°C之间从0.044增加到0.37左右不是常数。线性化有一个前提就是|Ts - Ta|不能太大。对水面和湿润冠层来说表面温度与气温通常相差不超过10°C这一项导致的误差相对于饱和水汽压差D可以接受。但对干燥裸土或极端干旱冠层表面温度可能比气温高20°C以上经典彭曼公式的偏差就会显现。那种场景要用更完整的显式阻力模型而不是怪公式本身。3.2 感热通量改写为辐射项的函数先从感热通量方程解出Ts - TaH ρCp * (Ts - Ta) / ra代入能量平衡 Rn H λE GTs - Ta (Rn - G - λE) * ra / (ρCp)这个结果可以直观理解净辐射减去潜热和土壤热通量后剩余能量会用来把蒸发面加热到高于气温的温度。如果蒸发非常旺盛λE很大Ts - Ta就会变小对应蒸发冷却现象。后面要用这个关系式替换水汽方程里的Ts。3.3 把两个方程装进一个式子把线性化后的es(Ts)代入潜热通量方程λE (ρCp / γ) * [D Δ(Ts - Ta)] / ra其中D es(Ta) - ea是气温下的饱和水汽压与实际水汽压之差也就是常说的饱和水汽压差。再把3.2的Ts - Ta表达式代入λE (ρCp / γ) * [D Δ * (Rn - G - λE) * ra / (ρCp)] / ra展开后把含λE的项都移到等号左侧λE * (1 Δ / γ) Δ / γ * (Rn - G) ρCp / (γ * ra) * D两边同乘γ/(γ Δ)得到λE Δ/(Δ γ) * (Rn - G) γ/(Δ γ) * f(u) * D这里f(u) ρCp / (γ * ra)是一个速度函数实际计算时用风速的经验函数替代。式子两边再除以蒸发潜热λ就得到以水层深度表示的开阔水面蒸发量E Δ/(Δ γ) * (Rn - G)/λ γ/(Δ γ) * f(u) * D这就是完整推导的终点。可以看到公式里没有任何剩余的表面温度未知量只需要输入气温、湿度、风速和净辐射。3.4 对公式中每一项的物理意义右侧第一项是辐射项权重为Δ/(Δγ)。气温越高Δ越大辐射项的权重越高说明在温暖气候下蒸发主要由能量决定。第二项是空气动力学项权重为γ/(Δγ)它乘上风函数f(u)和饱和水汽压差D代表干燥大气和风联合驱动下的蒸发能力。在低温环境比如10°C时Δ约0.082 kPa/°CΔ/(Δγ)只有0.55左右空气动力学项的权重上升到0.45。在35°C时Δ约0.35 kPa/°C辐射项权重升到0.84。这个权重变化解释了为什么热带湿润地区只靠辐射估算蒸发也能得到不错结果而温带和干旱区必须把风与湿度纳入。下面一个表整理不同情形下各项的反应情形辐射项变化空气动力学项变化对E的总体影响气温升高权重增大Rn可能增大D增大权重略降增大风速增大基本不变f(u)增大增大湿度增大基本不变D减小减小净辐射增强直接增大间接影响小增大这套对照关系在调试代码时很有用。如果某个输入变化后E的反应方向与表里相反一定是在符号或单位上出了错。4. 用Python把彭曼公式算出来并核对单位推导的终点是一个显式公式但能不能用对取决于单位。很多工程计算错误不是公式写成而是净辐射用了W/m²、风速用了10m高度、湿度用了露点而不是相对湿度。下面给出一套可执行的Python实现并用手算检查每个中间量。4.1 输入变量与单位约定本节使用经典彭曼公式的水层深度写法E Δ/(Δγ) * (Rn - G)/λ γ/(Δγ) * f(u2) * (es - ea)Rn和G统一为MJ/(m²·day)λ取2.45 MJ/kg因此第一项计算出来是kg/m²/day数值上近似等于mm/day。es、ea和D单位都为kPa。风函数f(u2)在经典写法中单位是mm/day/kPa输入风速为2m高度处平均风速u2单位m/s。风函数没有统一标准不同文献给出的系数差异很大。我这里演示用的是Penman风函数的一种常用简化形式f(u2) 0.26 * (1 0.536 * u2)请注意这个式子适用于中纬度湿润半湿润地区换到强风区或高原后应该用本地率定系数。如果你要的是参考作物ET0建议直接用第5章的FAO-56公式而不是抠经典风函数。4.2 完整实现代码前面第2章的三个函数继续使用这里再补上风函数和逐日计算函数import math def saturation_vapor_pressure(T): return 0.6108 * math.exp(17.27 * T / (T 237.3)) def slope_sat_vapor_pressure(T): es saturation_vapor_pressure(T) return 4098 * es / (T 237.3) ** 2 def psychrometric_constant(P101.3): lam 2.45 cp 1.013e-3 return cp * P / (0.622 * lam) def wind_function(u2): # u2: 2m高度风速m/s # 返回风函数单位mm/day/kPa return 0.26 * (1 0.536 * u2) def penman_day(Ta, RH, u2, Rn, G0.0, P101.3): # Ta: 平均气温°C # RH: 平均相对湿度0~100 # u2: 2m高度风速m/s # Rn: 净辐射MJ/m2/day # G: 土壤热通量MJ/m2/day日尺度取0 es saturation_vapor_pressure(Ta) ea es * RH / 100.0 D es - ea delta slope_sat_vapor_pressure(Ta) gamma psychrometric_constant(P) lam 2.45 radiation_term (delta / (delta gamma)) * (Rn - G) / lam aerodynamic_term (gamma / (delta gamma)) * wind_function(u2) * D return radiation_term aerodynamic_term代码逻辑按公式顺序展开没有任何隐藏换算。参数说明D是实际水汽压项的核心它把相对湿度转成kPa避免从露点温度反算时多一层误差。日尺度上G默认0只有在小时尺度或裸土场景才需要显式传入G。气压P只影响γ海拔越高γ越小空气动力学项权重降低但影响通常不超过10%。4.3 用一组气象数据核对计算结果取夏季晴天典型数据Ta25°CRH60%u22m/sRn25MJ/(m²·day)G0P101.3kPa。运行et penman_day(Ta25, RH60, u22, Rn25) print(et)输出约7.7mm/day。中间量可以手动核对es3.17kPaea1.90kPaD1.27kPaΔ0.188kPa/°Cγ0.066kPa/°Cf(u2)0.537mm/day/kPa。辐射项约7.53mm/day空气动力学项约0.18mm/day加和就是最终结果。从这个例子可以看出湿润夏季里辐射项占绝对主导风在湿度差很小时贡献有限。放一个中间量表方便观察变量数值单位es3.17kPaea1.90kPaD1.27kPaΔ0.188kPa/°Cγ0.066kPa/°Cf(u2)0.537mm/day/kPa辐射项7.53mm/day空气动力学项0.18mm/dayE07.71mm/day手动核算时的四舍五入会导致尾数有0.05以内的误差这属于正常现象。如果你算出的辐射项小了一半大概率是Rn用了W/m²而没除以0.0864换算到MJ/(m²·day)。4.4 常见单位坑最容易出错的地方有三个。第一个是净辐射气象站通常给出太阳总辐射Rs需要先换算成净辐射Rn相当于(1-反照率)乘以总辐射再扣除净长波辐射不能直接拿总辐射代入。第二个是风速高度公式只认2m高度如果传感器在10m需要按对数风廓线换算强风条件下高度误差会被放大。第三个是相对湿度直接用RH和气温计算ea即可不要用露点温度反查水汽压两者在温差大的日子会引入系统偏差。5. 从经典彭曼到FAO-56 Penman-Monteith的扩展与验证经典彭曼假设蒸发面完全湿润。真实的作物冠层并非如此气孔会控制水汽散失干旱时蒸腾还会进一步下降。Monteith在经典推导基础上加入表面阻力rs并把空气动力学阻力ra显式留在公式里得到Penman-Monteith公式。FAO-56推荐的参考作物蒸散量就是它的实用版本单位更统一更适合作为标准。5.1 引入表面阻力后的方程结构湿润表面的潜热通量用的是es(Ts)覆盖整个冠层表面。引入表面阻力rs后水汽路径变成两段串联从叶片气孔到冠层表面是rs从冠层表面到参考高度是ra。于是潜热通量表达式变为λE (ρCp / γ) * (es(Ts) - ea) / (ra rs)感热通量仍走纯粹的空气路径H ρCp*(Ts - Ta)/ra。把这两个式子连同能量平衡方程联立在线性化es(Ts)后可以推出λE (Δ(Rn - G) ρCp * D / ra) / (Δ γ*(1 rs/ra))当rs0时分母变成Δ γ分子和经典彭曼的阻力形式严格一致。所以Penman-Monteith是经典公式的自然推广经典彭曼只是它在中性湿润表面下的特例。rs越大分母越大蒸发受到的气孔抑制越强这对应干旱胁迫场景。5.2 FAO-56公式与经典彭曼的对照FAO-56把参考作物定义为一种假想的、高度0.12m、表面阻力70s/m、反照率0.23的草地。把这个设定代入Penman-Monteith并用24小时平均气温简化热传导项得到ET0 [0.408Δ(Rn-G) γ * (900/(Ta273)) * u2 * D] / [Δ γ*(1 0.34*u2)]0.408是蒸发潜热λ的倒数因为FAO-56直接输出mm/day。900/(Ta273)来自空气动力学阻力ra的经验近似。还有一个细节值得注意在FAO-56的标准化设定里空气动力学阻力约等于208/u2表面阻力rs取70s/m于是rs/ra 70*u2/208 ≈ 0.336u2这也是分母中0.34u2的来源。0.34u2已经内嵌了参考作物的气孔特性所以不要在FAO-56公式外再乘一个风函数那样会重复计入空气动力学阻力。两个公式的适用边界有必要说清对比点经典彭曼FAO-56 PM表面阻力无rs0隐含70s/m适用对象开阔水面或充分湿润下垫面参考作物ET0时间尺度日、月、小时均可用设计为24小时风函数经验式需率定由0.34u2内置单位风险风函数版本多易错常数统一风险低日常做农田灌溉推荐用FAO-56但研究水面蒸发和地表能量闭合时回到经典阻力表达式会更灵活不会被内置的参考作物假设限制。5.3 用Python对比同一组气象数据扩展实现FAO-56代码非常简洁def fao56_et0(Ta, RH, u2, Rn, G0.0, P101.3): es saturation_vapor_pressure(Ta) ea es * RH / 100.0 D es - ea delta slope_sat_vapor_pressure(Ta) gamma psychrometric_constant(P) numerator 0.408 * delta * (Rn - G) gamma * (900 / (Ta 273)) * u2 * D denominator delta gamma * (1 0.34 * u2) return numerator / denominator用4.3同一组输入跑一遍结果约8.1mm/day和经典彭曼的7.7mm/day相比高5%左右。差异主要来自空气动力学项的处理经典水面公式的风函数是经验系数而FAO-56把参考作物气孔阻力内嵌到分母里。两者在常规气象条件下差异通常在10%到20%以内如果差异超过30%先检查输入单位而不是公式结构。这段代码也很适合用来做回归测试把同一组数据分别喂给penman_day和fao56_et0并在断言里锁住偏差范围防止后续改动破坏公式。6. 彭曼公式推导结果的自检灵敏度分析与量级校验公式实现了怎么确认没写错。我用三个方法做自检。第一个是符号检查对E分别求风速和净辐射的偏导数。风速和净辐射增加蒸散发都应增大除非饱和水汽压差D为负。第二个是极端场景检查把Rn和G设为0公式只剩空气动力学项此时E可正可负但不应出现巨大的正值或负值。第三个是量级检查在25°C、相对湿度60%、风速2m/s、净辐射20MJ/(m²·day)时E应在5到8mm/day超过这个范围就查单位。下面这段代码用数值差分检验灵敏度放在单元测试里很实用h 1e-6 dE_du (penman_day(25, 60, 2 h, 25) - penman_day(25, 60, 2, 25)) / h dE_dRn (penman_day(25, 60, 2, 25 h) - penman_day(25, 60, 2, 25)) / h print(dE_du, dE_dRn)以第4章的风函数实现∂E/∂u2大概在0.05到0.2mm/day每m/s∂E/∂Rn大概在0.25到0.4mm/day每MJ/(m²·day)。如果你的风速导数是0或者负风函数或符号必然有错如果辐射导数超过0.5通常是Rn或λ的单位不匹配。这个方法同样适用于FAO-56和自定义的小时尺度版本。把这两个导数的合理范围写进测试断言后续任何人改动单位换算或替换风函数只要导数跑偏就会立刻报错。这样彭曼公式的推导和实现才算真正闭环。本文还有配套的精品资源点击获取