陆地与大气之间的物质与能量的交换，在全球气候变化与大气环流中有显著作用，并且也是各类生态研究中的重要研究内容^[1]。地气之间能量交换的一个重要约束条件为地表热量平衡。地表热量平衡公式为：

$ {R_n} = H + {\rm{LE}} + {G_0} $

(1)

式中，R_n为净辐射通量，W·m^–2，是总收入辐射通量与总支出辐射通量的差值；H为感热通量，W·m^–2，是物体在加热或冷却过程中，温度升高或降低而不改变其原有相态所需吸收或放出的热量；LE为潜热通量，W·m^–2，是物质发生相变，且温度不发生变化时吸收或放出的热量；G₀为地表热通量，W·m^–2，反映了地面与下层土壤之间的热交换，是引起地表温度变化的热通量^[2]。现有相关文献对地表热量平衡的研究多以净辐射通量、感热通量和潜热通量为主。虽然地表热通量与其他三项通量相比数值很小，但是它是瞬时能量平衡的重要贡献者^[3]。地表热通量受到地表覆盖类型、土壤含水量状况以及太阳辐射等因素的影响而变化，不同情况下地表热通量在地表热平衡中所占比例不同^[4-5]。地表热通量作为地表能量平衡的重要组成之一，是众多学科领域的一项重要参数。

通常通过测量土壤某一深度处的土壤热通量数值，可得到地表热通量。土壤热通量的准确测量在农业、气象、生态、地理与土木工程等领域均有重要的研究意义^[6-7]。目前，获取土壤热通量的主要方法有量热法、梯度法、组合法和热通量板法。其中，热通量板在测量过程中存在接触热阻造成的误差：当将热通量板放入某一介质中或表面时，由于二者材料不同，即热特性上的差异，使得热通量板的存在增大或减少该介质本身的热阻，导致土壤热通量的测量值与真实值产生偏差^[8]。针对这一问题，Philip^[9]在1961年基于稳态热传导理论，提出了修正参数的概念。

通过分析常用土壤热通量测定方法的准确度，以及测量过程中可能导致误差的主要因素，可以帮助规范和改进测量技术，从而提高土壤热通量的测量精度，推进农、林、业、气象、生态等学科相关研究的顺利进行。本文以热通量板法为研究对象，模拟研究热通量板测量精度的影响因素，并且检验Philip在1961年提出的修正公式的适用性。

1 材料与方法 1.1 热通量板的相关原理

Carslaw和Jaeger在1959年推导出不同定解的条件下，一维坐标系下不同位置z（m）处温度T（K）和热通量G_c（W·m^–2）随时间变化的解析解。对于已知热特性参数的半无限长固体介质，初始温度为0 K，在端点处施加恒定热流F₀（W·m^–2），即定解的条件设为：

$ {G_c}\left( {0, t} \right) = {F_0}, T\left( {z, 0} \right) = 0 $

(2)

这种情况下，不同位置处温度和热通量随时间的变化分别为：

$ T\left( {z{\rm{, }}t} \right) = \frac{{{\rm{2}}{F_{\rm{0}}}}}{\lambda }\left[ {{{\left( {\frac{{\alpha t}}{\pi }} \right)}^{\frac{1}{2}}}{{\rm{e}}^{ - \frac{{{z^{\rm{2}}}}}{{{\rm{4}}\alpha t}}}} - \frac{z}{{\rm{2}}}erfc\frac{z}{{2\sqrt {\alpha t} }}} \right] $

(3)

$ {G_c}\left( {z, t} \right) = {F_0}erfc\frac{z}{{2\sqrt {\alpha t} }} $

(4)

式中，α为热扩散系数，m²· s^–1；erfc（x）为余误差函数。在室内实验设计中，可以通过满足定解条件式（1）~式（2）来实现，应用式（1）~式（4）作为热通量的理论参考值。

在实际测量中，通常将热通量板水平埋设在地表以下2.5~10 cm范围内^[10]。当有热流穿过热通量板，板内的填充材料充当热阻，使板的上下表面产生温度差，创造一个热表面和一个冷表面，板内上下表面的热敏电阻感受温度，通过热电偶将温度差转换为电压差输出，通过下面公式计算出热通量：

$ {G_p} = \frac{{{V_p}}}{{{E_{sen}}}} $

(5)

式中，G_p为热通量，W·m^–2；V_p为输出电压，V；E_sen为热通量板参数，mV·m²·W^–1。

由于实际应用中的热通量板和土壤热特性不一致，导致热流在热通量板边缘有扭曲变形^[11]，进而导致在应用式（1）~式（5）时，实际热流和该公式给出的预测值有较大的误差。1961年，Philip基于稳态热传导理论，提出了修正参数的概念：

$ f = \frac{\varepsilon }{{{\rm{1 + }}\left( {{\rm{1}} - \varepsilon } \right)H}} $

(6)

式中，f为修正参数，ε为土壤热导率与热通量板热导率的比值，H取决于热通量板的形状和尺寸（直径或长宽、厚度等）。

对于边长为L，厚度为T的方形薄板：

$ H = 1 - 1.70\frac{T}{L} $

(7)

对于直径为D、厚度为T的圆形薄板：

$ H = 1 - 1.92\frac{T}{D} $

(8)

修正后的热通量G_Phi为：

$ {G_{{\rm{Phi}}}} = {G_p} \cdot f $

(9)

Philip修正公式是基于稳态热传导条件，即土壤内各点温度不随时间的变化而变化。而在实际测量当中，环境温度多为瞬态而非稳态，此时Philip修正公式是否还适用，有待验证。

1.2 计算机模拟设计

本研究应用COMSOL版本中的传热模块进行模拟。稳态分析指在模拟所给定的条件下，达到稳定不变（温度、热量流动等状态不变）时对模型的温度分布、热流分布等的分析；瞬态分析指在模拟所给定的条件和指定时间内对温度分布、热流分布等随时间变化的分析。本文的模拟条件不考虑热对流的存在，所以采用传热模块中的传导部分，分别建立模型，进行稳态分析和瞬态分析。本文选取两种土壤分别在三种含水量下进行计算机模拟，所用参数来自参见文献[12]，具体参数值见表 1。

本研究采用HFP01型热通量板为模拟对象。在土壤中测量热通量时，将周围土壤以及热通量板看成一个圆柱体集合体，从而简化为3D柱状模型。其中，土柱的直径设定为40 cm，厚度为10 cm；热通量板的直径为8 cm，厚度为0.5 cm。将热通量板水平放置，且埋藏深度分别设定为1 cm，2 cm和5 cm。由于圆柱体为中心对称图形，所以可将维度从3D简化为2D轴对称。如图 1所示，由于轴对称特性只需模拟对称轴右侧区域，以便减小计算量。土柱上边界为恒定热流边界，恒定热通量设定为150 W· m^–2[13]；两侧边界为绝热边界；下边界为狄式边界，恒定温度为20℃。模拟土壤的子域参数设定如表 1所示；热通量板的子域参数设定依据文献[14]，热导率为0.8 W·m^–1· K^–1，热容为1 300 J·kg^–1·K^–1，密度为1 700 kg·m^–3，初始温度均为20℃。在瞬态模拟研究中，时间步长设定为0︰0.5︰1 800，即从t=0 s时刻运行至t=1 800 s时刻，共计30 min，以0.5 s为计算间隔。稳态模拟是模型温度及热流达到稳定时的状态，所以无需设定时间步长。

本文通过瞬态模拟，分析土壤热通量板测量误差及影响因素。土壤热通量测量值为HFP01板中直径为3 cm圆形中心区域的热通量平均值（图 1中的①）；为了消除热通量板边界的热流扭曲造成的误差，本研究将土壤热通量的理论值设定为埋藏深度处接近板右侧边界的热通量平均值（图 1中的②），通过比较测量值与理论值的差值，分析土壤含水量与装置埋藏深度对测量误差的影响，并且对模型的温度场与热流分布进行了稳态模拟，分析测量误差的来源。

本文通过稳态模拟分析两种土壤在不同含水量的情况下，在三个埋藏深度处测量值、Philip修正值与理论值相对误差，比较不同条件对Philip修正效果的影响。由于目前的热通量板中央的传感器尺寸设置不同，本文还模拟了当热通量板内的传感器半径分别为1 cm、2 cm、3 cm和4 cm时对热通量测量误差的影响。

此外，本文将瞬态模拟与稳态模拟相结合，利用稳态模拟中计算出的Philip修正参数，计算出瞬态模拟的热通量的瞬态修正值，并与稳态模拟中的理论值进行对比，从而考察Philip修正公式在瞬态条件下的适用性。

1.3 室内实验设计

室内实验采用取自官厅水库的砂土，并设置3个水分梯度：0.00 g·kg^–1，20.0 g·kg^–1和150 g·kg^–1。分别对三个含水量的容重、热导率、体积热容量等物理参数进行测量，各参数均重复6次，结果如表 2所示。

实验在有机玻璃容器内进行，容器的内径为12 cm，高度为10 cm，厚度为0.5 cm，中间深度5 cm处安置加热膜，加热膜的直径为12 cm，厚度为0.27 cm，加热膜上方安置热通量板装置，热通量板水平安放在加热膜2 cm下方，如图 1整体区域所示。热通量板装置与数据采集仪（型号CR1000，Campbell Scientific公司，Logan，UT）连接，采样频率设定为1Hz。对于不同的土壤含水量，均设置6次重复，取平均值作为测量值，再应用Philip（1961）修正公式进行校正得到修正值。通过满足定解的条件式（1）~式（2），应用式（1）~式（4）作为理论参考值。本文获得的相对误差均通过计算测量值与理论值的偏差得到。此外，为了表述方便，文中出现的符号如表 3所示。

2 结果与讨论 2.1 热通量板法测量精度的主要影响因素

两种土壤在三种含水量情况下，均为埋藏深度越深测量准确度越低。如图 2（a）所示，含水量为0.00 g·kg^–1时，土壤的热导率小于热通量板的热导率，使得测量热通量大于土壤实际热通量，相对误差为正；含水量为200 g·kg^–1和400 g·kg^–1时，相对误差为负。两种情况的相对误差绝对值均随着埋藏深度的增大而增大。这是由于在计算机模拟中，表面被设定为恒定热源，埋深越深则热通量板的边界效应越明显，热流扭曲越大，从而偏离理论值越远。

两种土壤在三种埋藏深度的情况下，含水量影响大致相同，干燥土壤的测量误差较大；具有一定含水量的土壤的测量误差较小。如图 2（b）所示，当热通量板埋藏深度为5 cm时，相对误差由小到大依次排列为：含水量为200 g·kg^–1的黏土、含水量为400 g·kg^–1的黏土、含水量为400 g·kg^–1的砂土。

结果表明，土壤热通量测量的相对误差大小取决于土壤热导率与热通量板热导率的差异程度，热通量板的热导率越接近土壤的热导率，测量结果的准确度越高。

2.2 计算机模拟热通量板的误差来源

热通量板放入土壤后，在热通量板的边界会产生温度与热流的不均匀分布。当土壤热导率高于热通量板热导率时，在热通量板与土壤交界处，土壤中的温度梯度小于热通量板内温度梯度，所以热通量板的放置影响了土壤中温度的分布。其热流分布如图 3（a）所示，土壤的热通量高而流经热通量板的热通量低。以热通量板的径向作为横坐标，土壤热通量为纵坐标，如图 4（a）所示，在交界处，土壤热通量对应有峰值出现，而流经热通量板的热通量则出现低谷。当土壤热导率低于热通量板热导率时，以埋藏深度为1 cm的含水量为0.00 g·kg^–1的砂土为例，如图 3（b）和图 4（b）所示，变化规律相反：在交界处，土壤热通量对应有低谷出现，而流经热通量板的热通量出现峰值。

上述两种情况显示，由于热通量板和土壤热导率的差异，导致在热通量板边界处，土壤热通量与热通量板中热通量的不同，即在相同的埋深处，热通量板边界处产生了热通量的突变值。

2.3 Philip修正公式的效果

稳态模拟结果表明（如表 4所示），传感器半径越小，Philip修正效果越好；修正前后相对误差均随深度增加而增大；Philip修正公式在具有一定含水量的土壤中应用性较好，而应用于干燥土壤时较差。瞬态模拟结果表明，Philip修正公式在干燥土壤中应用性较差，在具有一定含水量的土壤中应用性较好（如图 5（a）所示）；在埋深2 cm和5 cm处的应用性好于埋深1 cm（如图 5（b）所示）。

对三种含水量的砂土进行室内实验，应用式（1）~式（9）的结果作为理论值。以Carslaw理论值G_E-C作为横坐标，以测量值G_E-P和G_E-Phi作为纵坐标，对数据进行过原点的直线拟合分析，并以1︰1线为参考，分析Philip的修正效果，如表 5所示。在土壤含水量为20.0 g·kg^–1时，HFP01的测量准确度较高，所以修正公式的修正效果并不明显；在土壤含水量为0.00 g·kg^–1和150 g·kg^–1时，虽然修正后的拟合直线斜率均较未修正的拟合直线斜率更加接近于1，但是校正后的拟合直线斜率与1仍存在偏离，说明修正值与理论值之间仍存在偏差。

计算机模拟与室内实验结果表明，由于Philip的稳态假设处理，修正后仍存在误差。因为Philip的解是热传导方程在稳态条件下给出的，并不能应用于瞬态问题，特别是当土壤温度未达到稳定分布时。所以，在土壤温度随时间波动的情况下，Philip的解会偏离瞬态热传导方程给出的解，并导致测量误差。

3 结论

本文采用有限元仿真软件COMSOL进行计算机模拟，直观呈现出热通量板测量准确性的影响因素。通过计算机模拟与室内实验结合，考察了Philip修正公式在瞬态条件下的适用性。结果表明：和文献所述不同^[15]，热通量板测量结果的准确度并非随着埋藏深度增加而提高；本实验结果为深度越深准确度越低；土壤热导率是热通量板测量准确度的主要影响因素，热通量板的热导率越接近土壤的热导率，测量结果的准确度越高。计算机模拟证实了前人的猜测^[16-17]：由于热通量板与土壤热导率的差异，使得两种材料中温度分布不同，进而由于温度梯度的差异导致在热通量板的边界产生热流的不均匀分布。Philip修正公式在干燥的土壤情况下，修正准确性较差；在具有一定含水量的土壤情况下，修正准确性较好。Philip修正理论有局限性，是因为在修正公式推导过程中，为简化运算做了很多近似，并且在土壤温度随时间波动的情况下，Philip的解会偏离瞬态热传导方程给出的解，导致了校正值与理论值间仍存在差异的现象，表明此公式有待优化。