A8 - 能量平衡、光合作用与叶温

同样的光照下，叶片会因空气温度、湿度、风速和蒸腾条件不同而达到不同温度。本例把一片叶子放在顶部开放的腔体中，从辐射平衡和固定气孔阻力开始，再验证光合作用、动态气孔与叶温的耦合。四部分共用同一套几何场景，分别与解析值或独立方程求解结果比较。

场景与输入

拖动可旋转，也可单独打开三维场景。腔体四壁与地面温度固定，叶片温度由能量平衡求解；热壁条件将东墙单独加热至 340 K。

条件	设置
腔体	4 m × 4 m，高 3 m，顶部开放，无周期边界
壁面和地面	短波反射率 0，长波发射率 1，温度为环境温度
叶片	1 m × 1 m，中心 (2, 2, 1) m，由两个三角形组成；每个三角形单独求温度
光学属性	反射率 0.2，透射率 0；长波不透射
太阳光	550 nm，带宽 1 nm；表中 E 为法向光谱辐照度，单位 W m⁻² nm⁻¹
天空	无大气、无短波天空光；下行长波为 σTₑ⁴
初始叶温	280 K，另用 260 K 与 330 K 检查初值影响
长波采样	每三角形每面 32、128、512；热壁为 512、2048、8192、32768

下表每行是一组条件，三个后端使用相同输入。除第 4 行外，叶片水平、太阳天顶角为 0°。第 4 行叶片向东抬升 30°，太阳位于东侧、天顶角 30°，入射余弦为 cos 60°。

序号	条件	E	环境温度 / K	发射率	辐射面数	气孔	其他设置
1	同温平衡	0	300	0.95	2	关闭	—
2	低温弱光	100	280	0.95	2	关闭	—
3	低发射率强光	400	285	0.70	2	关闭	—
4	倾斜叶片	400	300	0.95	2	关闭	入射角 60°
5	单面交换	200	295	0.90	1	关闭	—
6	显热交换	400	300	0.95	2	关闭	气温 25°C，湿度 50%
7	下表面蒸腾	400	300	0.95	2	下表面	气温 25°C，湿度 50%
8	高湿度	400	300	0.95	2	下表面	湿度 90%
9	双面蒸腾	400	300	0.95	2	两面	湿度 50%
10	风速控制	400	300	0.95	2	下表面	风速 4 m/s，叶宽 0.05 m，按叶宽与风速计算阻力
11	局部热壁	400	300	0.95	2	下表面	东墙 340 K

第 1–5 组的气温设为各自的解析辐射平衡温度，使解处的显热恰好为零。第 6–11 组气温均为 25°C，除第 8 组外相对湿度均为 50%。指定的空气交换阻力为 100 s/m，气孔阻力为 150 s/m；第 10 组由叶宽和风速计算空气交换阻力。叶片土壤热通量 G 为零。

理论参考

第一步：辐射平衡点

设叶片反射率为 r、辐射交换面数为 n、发射率为 ε，环境温度为 Tₑ。吸收的短波和长波分别为

\[ Q=(1-r)E\Delta\lambda\cos i,\qquad L_{abs}=n\varepsilon\sigma T_e^4. \]

关闭蒸腾并令气温等于辐射平衡温度，解处 H = LE = G = 0，因而

\[ T_r=\left(T_e^4+\frac{Q}{n\varepsilon\sigma}\right)^{1/4}. \]

既比较叶温，也分别比较吸收短波、吸收长波、发射长波和显热。第一组还检查无太阳光时叶片能否回到 300 K 的同温平衡。

第二步：显热与蒸腾

给定气孔阻力，独立求解

\[ Q+L_{abs}-n\varepsilon\sigma T^4-H(T)-LE(T)=0, \]

\[ H=C_H(T-T_a),\qquad LE=C_E[e_s(T)-e_a], \]

\[ e_s(T)=610.78\exp\left[\frac{17.27(T-273.15)}{T-273.15+237.3}\right], \qquad e_a=\frac{RH}{100}e_s(T_a). \]

温度用 K，水汽压用 Pa。指定空气交换阻力 rₐ 时，C_H = ρcₚ/rₐ，该阻力表示整个表面的显热交换；按叶宽计算时，rₐ = 200√(d/max(u, 0.1))，C_H = nρcₚ/rₐ。潜热系数为 C_E = nₛρcₚ/[γ(rₐ+rₛ)]，nₛ 是有气孔的面数。

采用 σ = 5.670374419 × 10⁻⁸ W m⁻² K⁻⁴、ρ = 1.225 kg/m³、cₚ = 1004 J kg⁻¹ K⁻¹、γ = 66.5 Pa/K。参考脚本用十进制高精度运算和二分法求根，不调用 LESS 的能量平衡求解器。前两部分采用固定气孔阻力，后两部分进一步验证光合作用控制的气孔响应。

热壁与三维长波交换

第 11 组按三角形面积平均的视角因子计算热壁输入：

\[ L_{abs,j}=2\varepsilon\sigma T_e^4+ \varepsilon\sigma F_{j,east}(340^4-T_e^4). \]

每个接收点到矩形壁面的视角因子由多边形投影立体角计算，再对叶片三角形做 Gauss 积分。两个三角形到东墙的双面视角因子分别为 0.180504129346 和 0.135629828566，因此叶温也略有不同。32 阶与 64 阶积分的最大变化小于 10⁻¹²。视角因子的定义见 Sandia 的说明。

验证结果

后端	等温环境叶温误差 / K	热壁叶温误差 / K	能量闭合残差 / W m⁻²	结果
optix	0.000221	0.008672	0.003326	通过
vulkan	0.000221	0.008641	0.002930	通过
embree	0.000221	0.008672	0.003342	通过

33 组最终结果全部通过。表中误差取每个后端全部对应条件、两个叶片三角形的最大值；闭合残差独立按 |S_abs + L_abs − L_emit − H − LE − G| 计算，单位均按单面叶面积归一化。

等温环境要求叶温误差不超过 0.005 K、吸收短波误差不超过 0.005 W/m²、吸收长波误差不超过 0.02 W/m²，其余通量误差不超过 0.05 W/m²。热壁条件允许叶温误差 0.1 K、长波和其他通量误差 3 W/m²。所有条件的能量闭合残差均要求不超过 0.02 W/m²。

采样与迭代收敛

等温环境提供几乎无采样波动的对照；热壁条件则检验长波采样。图中用同一个种子逐级增加采样，最终一级相对前一级的叶温变化不超过 0.05 K。最终预算另用三个种子重复，叶温标准差不超过 0.00832 K；三个后端的均值与波动幅度一致。

低发射率强光和热壁条件另将温度/能量迭代容差从 0.01 K / 0.5 W/m² 收紧到 0.001 K / 0.05 W/m²，再收紧到 0.0001 K / 0.01 W/m²。最后两级叶温变化不超过 0.000214 K，平均通量变化不超过 0.002708 W/m²。从 260 K 与 330 K 开始求解，相对最终结果的叶温差均小于 0.005 K。

光合作用与动态气孔

第三步：固定叶温

叶片仍为反射率 0.2、不透射的平面，使用总吸收光子驱动 C₃ 光合作用。固定叶温后，分别改变光照、CO₂、湿度和温度，检验光合速率与气孔导度。吸收的 PAR 光子通量由 Qλ/(hcNₐ) × 10⁶ 换算，单位为 μmol 光子 m⁻² s⁻¹。

条件	E / W m⁻² nm⁻¹	叶温 / K	CO₂ / μmol mol⁻¹	空气湿度	方程形式
暗条件	0	298.15	400	50%	默认平滑共限制、显式 Jmax
弱光	25	298.15	400	50%	默认平滑共限制
中等光照	100	298.15	400	50%	默认平滑共限制、显式 Jmax
强光	400	298.15	400	50%	默认平滑共限制、显式 Jmax
高 CO₂	400	298.15	800	50%	默认平滑共限制
干燥空气	400	298.15	400	30%	默认平滑共限制
高叶温	400	308.15	400	50%	默认平滑共限制

气温为 25°C。Vcmax25 = 60 μmol m⁻² s⁻¹，Rd25 = 0.9 μmol m⁻² s⁻¹，Ball–Berry 斜率 m = 9、截距 b = 0.01 mol H₂O m⁻² s⁻¹；显式 Jmax 形式取 Jmax25 = 120 μmol m⁻² s⁻¹。气压 101325 Pa、氧浓度 209 mmol mol⁻¹。默认平滑共限制对应 formulation="scope"，另一种为 formulation="jmax"。

参考根据 Farquhar 等的 C₃ 光合模型及所选方程形式的温度响应，独立联立生化同化与气孔扩散：

\[ A_c=V_{cmax}\frac{C_i-\Gamma^*}{C_i+K_m},\qquad g_s=\max\left(b,b+m\frac{A_n h_s}{C_a}\right),\qquad A_n=\frac{g_s}{1.6}(C_a-C_i). \]

默认形式采用平滑共限制，将羧化速率、光限制速率和产物输出限制联立，计算净同化速率 Aₙ；hₛ 使用叶温下的表面相对湿度。其方程说明见生化模型文档。显式 Jmax 形式使用非直角双曲线计算电子传递速率，取羧化和电子传递两种限制的较小值，再减去暗呼吸；气孔关系使用环境湿度，并对扩散关系中的 Aₙ 取非负部分。暗条件的净同化为负的暗呼吸速率，气孔导度回到截距 b。

参考脚本在 CO₂ 浓度上直接二分求根，独立于 LESS 的原生迭代。完整速率、温度系数和两种形式均在参考与运行脚本中给出。

第四步：开放叶温耦合

采用强光条件，分别验证默认环境、高湿度 90%、CO₂ 800 μmol mol⁻¹、东墙 340 K，以及显式 Jmax 形式。叶温从 280 K 开始，气孔导度由光合作用计算，再换算为水汽阻力：

\[ r_s=\frac{P}{RTg_s},\qquad LE=\frac{\rho c_p}{\gamma}\frac{e_s(T)-e_a}{r_a+r_s}. \]

空气交换阻力仍为 100 s/m，只有叶片下表面有气孔。耦合形式采用 cₚ = 1005 J kg⁻¹ K⁻¹、γ = 66 Pa/K。参考在每个候选叶温下先求解光合作用与气孔，再对净辐射减去显热和潜热的方程二分求根。热壁仍使用前述两个三角形的独立视角因子。

验证分别检查净同化速率、气孔导度、叶温和通量误差；将 LESS 输出的 Aₙ、gₛ 代回生化与气孔方程检查碳交换闭合，再用输出的叶温和 gₛ 重新计算显热、潜热及能量残差。

后端	净同化误差 / μmol m⁻² s⁻¹	气孔导度误差 / mol m⁻² s⁻¹	叶温误差 / K	独立能量残差 / W m⁻²	结果
optix	0.004230	0.00007750	0.009891	0.003833	通过
vulkan	0.004228	0.00007752	0.009891	0.003867	通过
embree	0.004232	0.00007624	0.009891	0.003005	通过

新增 15 组条件在三个后端上共 45 组最终结果全部通过。净同化误差要求不超过 0.02 μmol CO₂ m⁻² s⁻¹，气孔导度误差不超过 0.0005 mol H₂O m⁻² s⁻¹，APAR 误差不超过 0.002 μmol 光子 m⁻² s⁻¹。均匀环境叶温误差限为 0.005 K，热壁为 0.1 K；相应通量误差限为 0.1 和 3 W/m²。输出通量的能量闭合限为 0.02 W/m²，按叶温与气孔导度独立重算的能量残差限为 0.05 W/m²。

耦合条件使用每三角形每面 32、128、512 个长波样本；热壁使用 512、2048、8192、32768、65536 个样本，最后两级叶温变化小于 0.05 K。默认形式、热壁和显式 Jmax 三组还逐级收紧温度、能量和生理迭代容差，最终取 0.0001 K、0.01 W/m² 和 0.0001 的相对生理变化限。收紧容差及从 260 K、330 K 开始求解，最终叶温变化不超过 0.000793 K，净同化变化不超过 0.000668 μmol m⁻² s⁻¹，气孔导度变化不超过 0.00000751 mol m⁻² s⁻¹。

复现

下载输入条件、墙面、水平叶片、倾斜叶片、主运行脚本与光合作用参考及运行脚本。保持 scripts/ 与 data/energy-balance/ 的相对位置，安装 LESS、NumPy 和 Matplotlib 后运行全部四部分：

python scripts/appendix_energy_balance.py

完整精度与收敛结果包含当前输入、独立参考、最终通量与叶温，以及采样和迭代检查。