检索项 检索词 到
  土壤学报  2023, Vol. 60 Issue (3): 599-609  DOI: 10.11766/trxb202210220581
0

引用本文  

张建彬, 高志球, 童兵, 等. 土壤温度预报方程研究进展. 土壤学报, 2023, 60(3): 599-609.
ZHANG Jianbin, GAO Zhiqiu, TONG Bing, et al. Progress of Soil Temperature Prediction Equation. Acta Pedologica Sinica, 2023, 60(3): 599-609.

基金项目

国家自然科学基金项目（41875013）资助

通讯作者Corresponding author

高志球，E-mail：zgao@nuist.edu.cn

作者简介

张建彬（1983—），男，新疆米东人，博士，高级工程师，数值预报模式及暴雨模拟。E-mail：20211103015@nuist.edu.cn
Contents              Abstract              Full text              Figures/Tables              PDF
土壤温度预报方程研究进展
张建彬1, 高志球1,2, 童兵3, 王琳琳2    
1. 南京信息工程大学大气物理学院, 南京 210044;
2. 中国科学院大气物理研究所大气边界层物理与大气化学国家重点实验室, 北京 100029;
3. 中国气象科学研究院灾害天气国家重点实验室, 北京 100081
收稿日期：2022-10-22；收到修改日期：2022-11-28；优先数字出版日期(www.cnki.net)：2023-02-01
基金项目：国家自然科学基金项目（41875013）资助
作者简介：张建彬（1983—），男，新疆米东人，博士，高级工程师，数值预报模式及暴雨模拟。E-mail：20211103015@nuist.edu.cn.
通讯作者Corresponding author：高志球，E-mail：zgao@nuist.edu.cn.
摘要：土壤温度（尤其是地表温度）是陆地和大气之间相互作用中关键的物理量，在地球系统中扮演了十分重要的角色。土壤温度预报技术一直是陆面模式、数值天气预报和气候预测中核心科学问题。本文系统回顾了土壤温度预报方程的研究进展，从经典的热传导方程到考虑了土壤水分垂直运动物理过程的热传导-对流方程，从用单一正弦波逼近到用傅里叶级数逼近地表温度日变化，从假设对流参数无日变化为常数到考虑其日变化，着重概述了土壤热传导-对流方程的创建、改进及求解。最后，本文对热传导-对流方程在地表能量平衡、土壤水分垂直运动、水通量和地震、冻土热传输研究中的应用进行了回顾。同时指出，全相态的土壤水和植物根系对热传导-对流方程的影响是土壤温度预报方程未来的研究方向。
关键词：土壤温度    热传导    热传导-对流    预报    
Progress of Soil Temperature Prediction Equation
ZHANG Jianbin1, GAO Zhiqiu1,2, TONG Bing3, WANG Linlin2    
1. School of Atmospheric Physics, Nanjing University of Information Engineering, Nanjing 210044, China;
2. State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China;
3. State Key Laboratory of Severe Weather, Chinese Academy of Meteorological Sciences, Beijing 100081, China
Foundation item: Supported by the National Natural Science Foundation of China(No. 41875013)
Abstract: Soil temperature (especially surface temperature) is a key physical quantity in the interaction between land and atmosphere, and plays a very important role in the earth system. Soil temperature prediction technology has always been the core scientific problem in land surface model, numerical weather prediction and climate prediction. This paper systematically reviews the research progress of soil temperature prediction equation, from the classical heat conduction equation to the heat conduction convection equation that takes into account the physical process of vertical movement of soil moisture, from the single sine wave approximation to the Fourier series approximation of the daily change of surface temperature, from the assumption that the diurnal change of convection parameters is constant to the consideration of its diurnal change, and emphatically summarizes the creation, improvement and solution of the soil heat conduction convection equation. Finally, this paper reviews the application of heat conduction convection equation in the study of surface energy balance, vertical movement of soil moisture, water flux, earthquake and frozen soil heat transfer. At the same time, it is pointed out that the influences of soil water phases and plant roots on the heat conduction-convection equation is warranted for the future research of soil temperature prediction equation.
Key words: Soil temperature    Heat conduction    Heat conduction-convection    Prediction    

土壤温度不仅在地表能量分配中扮演着重要角色，而且影响着土壤物理和化学以及植物生长等过程。因此，准确监测或预报土壤温度尤其是土壤表面温度一直是大气、土壤、水文和生物过程数值模型研发的重要环节。自从2003年，根据理论模型[1]或半经验模型[2]，Gao等[3]在半无限介质一维导热方程的解析解[4]的基础之上，首次提出了土壤热传导-对流方程，并通过Laplace变换推到出其数值解和利用谐波方法获得其解析解以来。围绕这两个土壤温度方程的分析和应用研究持续不断，并不断取得研究进展（包括：土壤热扩散率的确定、土壤温度预报、近地层能量平衡、土壤水垂直迁移和地震研究等）[5-19]。本文将对这些研究进展进行系统总结。

土壤温度预报直接影响陆面过程模式的表现。在陆面过程模型中，“强迫-恢复”方法是核心方程，而地表温度是最关键的诊断量[20-22]。最近50年来，随着观测数据的不断积累，土壤温度预报的研究持续不断，并不断取得新进展。在浅层地热应用、植物生长以及地表能量平衡中，通过土壤的热传递非常重要。本文主要聚焦土壤热传导-对流方程的创建、改进和求解，及其在不同研究中的应用进行评述。

1 土壤温度预报方程研究概述 1.1 土壤热传导

在前人基础上[4，23-24]，Gao等[3]推导了土壤热传导方程。由于分子传导控制着传输过程，他们用均匀体中的傅立叶热传导定律来表示土壤中z深度处的地下热通量Qg（W m–2），即，通量取决于土壤温度梯度，如下所示：

$ {Q_{\text{g}}} = - \lambda \partial T/\partial z, $ (1)

其中，负号表示热流沿着温度梯度$\partial T/\partial z$向下，$\lambda $为土壤热分子传导率（W m–1 K–1），T为土壤温度（K），z为土壤的垂直坐标（正向下）。又基于Garratt[25]的研究，忽略厚度为z的薄层土壤中的水平热传导，假设土壤温度随时间的变化取决于流入和流出该层的热通量的差异。他们通过对浅层土壤热量守恒的研究，发现了表层土壤热通量Q0与土壤小深度热通量Qg之间的关系，并用热力学第二定律得出了简单的预测方程：

$ {C_{\text{g}}}\partial T/\partial t = - \partial {Q_{\text{g}}}/\partial z, $ (2)

其中，左侧项表示单位体积土壤中随时间变化的能量变化率；右边项是由于土壤中的分子传导而产生的单位体积能量的净输入；t为以秒为单位的时间。

在此基础上，假设研究区域内的土壤具有（ⅰ）各向同性和均质性，以及（ⅱ）含水量随深度不变或其变化对Cg和$\lambda $的影响可忽略不计。此外，假设能量交换仅发生在垂直方向上且方程（2）可简化为固体柱内热传导方程，由此推导出了土壤热传导方程，即：

$ \partial T/\partial t = k{\partial ^2}T/\partial {z^2} $ (3)

其中，$k = \lambda /{C_{\text{g}}}$（m2 s–1）为土壤热扩散率。

1.2 土壤热对流方程

热量可以通过传导、对流和辐射在土壤中传递[26-27]；然而，大多数土壤温度变化发生在地表附近的浅层内[26]，辐射可以忽略。因此，土壤温度和含水量的分布是研究土壤热性质的关键变量。为此，在土壤热传导方程基础上，Gao等[3]用Qv表示由水的垂直运动引起的热通量，并假设液体渗透速率为W（m3 s–1 m–2，正向上），将Qv表示为：

$ {Q_{\text{v}}} = {C_{\text{w}}}W\theta \Delta T, $ (4)

其中，Cw为水的体积热容量（4.2 MJ W–3 K–1），θ为土壤体积含水量（m3 m–3），而ΔT为土壤中水温的垂直差（℃）。又根据热力学第二定律可以得到：

$ {C_{\text{g}}}\partial T/\partial t = - \partial {Q_{\text{v}}}/\partial z, $ (5)

进而可以得到土壤热对流方程，

$ \partial T/\partial t = W\partial T/\partial z, $ (6)
1.3 土壤热传导-对流方程的创建和求解

假设土壤中热传导和对流是各自独立的物理过程，Gao等[3]合并热传导方程（3）和对流方程（6），得到以下土壤热传导-对流方程：

$ \frac{{\partial T}}{{\partial t}} = k\frac{{{\partial ^2}T}}{{\partial {z^2}}} + W\frac{{\partial T}}{{\partial z}}, $ (7)

方程（7）左边为土壤温度随时间的变化；$ {C}_{\text{g}} $（${C_{\text{g}}} = \rho \times c$，$\rho $为土壤密度（kg m–3））为土壤体积热容（J m–3 K–1），c为土壤比热（J kg–1 K–1），z为土壤的垂直坐标（正向下）。总体对流参数$W = - {C_{\text{w}}}/{C_{\text{g}}}\omega \theta $，$w$为液态水通量密度（m3 s–1 m–2）。

浅层土壤中k会随土壤深度的增加而变化。Gao等[7]进一步考虑了土壤热扩散率的垂直不均匀性，推导得到$W = \partial k/\partial z - {C_{\text{w}}}/{C_{\text{g}}}\omega \theta $，包括土壤总体热扩散率垂直梯度（$\partial k/\partial z$）和土壤中液态水通量密度（$ - {C_{\text{w}}}/{C_{\text{g}}}\omega \theta $）两部分。

土壤热传导-对流方程的求解可以用谐波法和拉普拉斯变换法。具体过程如下：

1）谐波法

基于Hillel[28]的假设，Gao等[3]将地表温度的日变化视为一个纯正弦函数：${\left. {T\left( t \right)} \right|_{z = 0}} = {T_0} + $ $A\sin \omega t, (t > 0)$其中T0是一个常数，A为土壤表面温度变化的振幅（℃），ω为地球旋转的角速度（$2\pi /p = 7.292 \times {10^{ - 5}}$rad s–1）。在此边界条件下，土壤热传导-对流方程可改写为：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial T}}{{\partial t}} = k\frac{{{\partial ^2}T}}{{\partial {z^2}}} + W\frac{{\partial T}}{{\partial z}}(t > 0, z > 0)} \\ {{{\left. T \right|}_{z = 0}} = {T_0} + A\sin \omega t\left( {t \geqslant 0} \right)} \end{array}} \right., $ (8)

我们通过一系列谐波运算得到方程（8）的解析解为：

$ \begin{gathered} T\left( {z, t} \right) = {T_0} + A\exp \left[ {\left( { - \frac{W}{{2k}} - \frac{{\sqrt 2 }}{{4k}}\sqrt {{W^2} + \sqrt {{W^4} + 16{k^2}{\omega ^2}} } } \right)z} \right] \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\sin \left[ {\omega t - z\frac{{\sqrt 2 \omega }}{{\sqrt {{W^2} + \sqrt {{W^4} + 16{k^2}{\omega ^2}} } }}} \right], \hfill \\ \end{gathered} $ (9)

方程（9）适用于上边界条件为随时间正弦变化的半无限介质。土壤热扩散系数k和对流参数W可通过测量z1和z2深度两个不同土层的土壤温度获得，公式为：

$ \left\{ {\begin{array}{*{20}{l}} {k = \frac{{{{\left( {{z_1} - {z_2}} \right)}^2}\omega \ln \left( {\frac{{{A_2}}}{{{A_1}}}} \right)}}{{\left( {{\Phi _2} - {\Phi _1}} \right)\left[ {{{\left( {{\Phi _2} - {\Phi _1}} \right)}^2} + {{\left( {\ln \left( {\frac{{{A_2}}}{{{A_1}}}} \right)} \right)}^2}} \right]}}} \\ {W = \frac{{\omega \left( {{z_2} - {z_1}} \right)}}{{{\Phi _2} - {\Phi _1}}}\left[ {\frac{{2{{\ln }^2}\left( {\frac{{{A_2}}}{{{A_1}}}} \right)}}{{{{\left( {{\Phi _2} - {\Phi _1}} \right)}^2} + {{\left( {\ln \left( {\frac{{{A_2}}}{{{A_1}}}} \right)} \right)}^2}}} - 1} \right]} \end{array}} \right. $ (10)

其中A1和A2，Φ1和Φ2分别为z1和z2深处土壤温度的振幅（K）和初始相位（rad）。

2）拉普拉斯变换法

微气象和水文研究中，传统的假设是初始土壤温度垂直剖面是恒定的，而不是随深度线性变化。然而，在现实中，平均土壤温度可能会随着土壤深度的增加而有所不同。Gao等[3]设定初始条件${\left. T \right|_{t = 0}} = {T_0} - \gamma z, \left( {z \geqslant 0} \right)$，其中T0为表面平均温度（K），γ（$\gamma = \mathop \sum \limits_{i = 1}^n {T_{0i}} - \mathop \sum \limits_{i = 1}^n \frac{{{T_i}}}{z}$，且γ=0是理想状态）为土壤温度衰减率（K m–1），z为土壤深度（m）。在此基础上，他们将边界条件设置为与谐波法相同，得到：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial T}}{{\partial t}} = k\frac{{{\partial ^2}T}}{{\partial {z^2}}} + W\frac{{\partial T}}{{\partial z}}(t > 0, z > 0)} \\ {{{\left. T \right|}_{t = 0}} = {T_0} - \gamma z, \left( {z \geqslant 0} \right)} \\ {{{\left. T \right|}_{z = 0}} = {T_0} + A\sin \omega t\left( {t \geqslant 0} \right)} \end{array}} \right., $ (11)

通过设置T=T*+T0-γz-γWt，应用拉普拉斯变换方案，并设${T^{**}} = {T^*}{e^{\frac{W}{{2k}}z + \frac{{{W^2}}}{{4k}}t}}$得到：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial {T^{**}}}}{{\partial t}} = k\frac{{{\partial ^2}{T^{**}}}}{{\partial {z^2}}}} \\ {t = 0, {T^{**}} = 0} \\ {z = 0, {T^{**}} = \left( {A\sin \omega t + \gamma Wt} \right){e^{\frac{{{W^2}}}{{4k}}t}}} \end{array}} \right., $ (12)

然后对T**进行拉普拉斯变换，并设$\bar T = \mathcal{L}\left[ {{T^{**}}} \right]$。

在此基础上，他们通过一系列推导、运算得到：

$ T = {T_0} - \gamma z - \gamma Wt + {T^*}, $ (13)
$ \begin{gathered} T\left( {z, t} \right) = {T_0} - \gamma z - \gamma Wt + \frac{{z{e^{\frac{W}{{2k}}z\frac{{{W^2}}}{{4k}}t}}}}{{2\sqrt {k\pi } }} \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\mathop \smallint \limits_0 \frac{{Asin(\omega \tau + \Phi + \gamma W\tau ){e^{\frac{{{W^2}}}{{4k}}\tau }}}}{{{{\left( {t - \tau } \right)}^{3/2}}}}{e^{\frac{{{z^2}}}{{4k\left( {t - \tau } \right)}}}}d\tau . \hfill \\ \end{gathered} $ (14)
1.4 土壤热传导-对流方程的改进和求解

考虑到W取日平均值不能体现土壤水分垂直迁移的日变化，Wang等[29]将W的日变化考虑到热传导-对流方程中，在Shao等[30]的基础上，用正弦函数比较土壤上边界条件，建立了以下方程：

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial T}}{{\partial t}} = k\frac{{{\partial ^2}T}}{{\partial {z^2}}} + \left( {a + b\sin (\omega t)} \right)\frac{{\partial T}}{{\partial z}}} \\ {T\left( {z, 0} \right) = f\left( z \right)} \\ {T\left( {\infty , t} \right) = {T_1}} \\ {T\left( {0, t} \right) = {T_0} + \mathop \sum \limits_i^n {A_i}\sin \left( {i\omega t + {\phi _i}} \right), \left( {t \geqslant 0} \right)} \end{array}} \right., $ (15)

令${T^*} = T\left( {z, t} \right) - {T_1}$可得到同性的边界条件，经此变换上述方程、边界条件和初始条件变为，

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial {T^*}}}{{\partial t}} = k\frac{{{\partial ^2}{T^*}}}{{\partial {z^2}}} + \left( {a + b\sin \omega t} \right)\frac{{\partial {T^*}}}{{\partial z}}} \\ {{T^*}\left( {0, t} \right) = \left( {{T_0} - {T_1}} \right) + \mathop \sum \limits_{i = 1}^n {A_i}\sin (i\omega t + {\Phi _i}), i = 1, 2, ..., n} \\ {{T^*}\left( {\infty , t} \right) = 0} \\ {{T^*}\left( {z, 0} \right) = f\left( z \right) - {T_1}} \end{array}} \right., $ (16)

经过不同推演（方程（17）-方程（35）的推导请见附录），要得到$V\left( {p, 0} \right)$，需要确定$f\left( Z \right)$，即初始时刻土壤温度廓线。如此，就要根据实际情况而定，因为不同区域或者质地土壤温度廓线肯定有所不同。

最后根据傅立叶逆变换得到，$V\left( {Z, t} \right) = \frac{2}{\pi }$ $\mathop \smallint \limits_0^\infty V\left( {p, t} \right)\sin (pZ)dp$，所以

$ U\left( {z, t} \right) = \frac{2}{\pi }\mathop \smallint \limits_0^\infty V\left( {p, t} \right)\sin (p\left( {x - {p_1}\left( t \right)} \right)dp, $ (36)

所以，${T^*} = U\left( {z, t} \right)\exp \left( { - \frac{{{a_1}^2t}}{{4k}} - \frac{{{a_1}z}}{{2k}}} \right)$。最后$T\left( {z, t} \right) = {T_1} + {T^*}\left( {z, t} \right)$。由于方程（36）是显式，因此最后的结果$T\left( {z, t} \right)$必是显式而非隐式。

表 1汇总了土壤温度的不同预报方程，及其边界条件和解析解。

表 1 土壤温度预报方程及其边界条件、初始条件和解析解 Table 1 Soil temperature prediction equation and its boundary conditions, initial conditions and analytical solutions
2 土壤热传导-对流方程的应用

土壤热传导-对流方程被应用于近地表能量平衡、土壤水运动趋势、土壤水通量和地震研究中（表 2），下面将分别阐述各个应用：

表 2 土壤热传导-对流方程的应用 Table 2 Application of soil heat conduction-convection equation
2.1 在近地表能量平衡研究中的应用

Gao[5]在$\frac{{\partial T}}{{\partial t}} = k\frac{{{\partial ^2}T}}{{\partial {z^2}}} + W\frac{{\partial T}}{{\partial z}}$的基础上提出了考虑热传导-对流过程的地表土壤热通量计算方法：

$ {G_0} = {G_z} + {C_g}k\frac{{\partial T}}{{\partial z}} + {C_w}W\nabla T $ (37)

利用1998年GAME/Tibet试验那曲站数据，计算了土壤热扩散率和对流系数，进而利用以上方程计算土壤热通量，并将计算的土壤热通量带入地表能量平衡方程中，发现能量平衡闭合度由0.7提高到0.9。

利用土壤热传导-对流方程的解，Gao等[9]对湿润土壤和干土的土壤表面热通量和温度的日变化的相位进行了理论表征和比较，发现自然土壤的表面热通量和温度的相位差在0至pi/4之间；近地表能量闭合度对土壤表面热通量振幅（A4）和净辐射振幅（A1）的比值十分敏感，高的A4/A1比值导致近地表能量闭合度差。

2.2 在土壤水运动趋势研究中的应用

基于土壤热传导-对流方程，Gao等[10]确定了塔克拉玛干沙漠站0~0.2 m土壤层10个月的热扩散率和对流参数；并通过对比对流参数（$W = \partial k/\partial z - {C_w}/{C_g}\omega \theta $）及其两组分：热扩散率垂直梯度（$\partial k/\partial z$）和液态水通量密度（$ - {C_w}/{C_g}\omega \theta $），首次发现热传导-对流方法可以捕捉月尺度上沙漠土壤水的运动趋势。

2.3 在土壤水通量/地震研究中的应用

基于我国川西断裂带基岩温度测量网络采集的基岩温度数据，Liu等[16]利用振幅和位相法计算了该地区的土壤热扩散系数，定量研究了土壤温度日振幅衰减（${A_1}/{A_2}$）和相位后移${\varnothing _1} - {\varnothing _2}$与土壤水分运动的关系，进而认定相位后移（${\varnothing _1} - {\varnothing _2}$）和振幅衰减（${A_1}/{A_2}$）之间的差异可以作为水运动或对流换热的示踪剂。通过观测网获取的数据分析，他们还发现（${\varnothing _1} - {\varnothing _2}$）和（${A_1}/{A_2}$）的空间变化特征与该地区的温泉分布相一致，从而为青藏高原东缘断裂带的热液活动提供了新的证据。由于热对流过程和热传导过程一起影响着土壤温度周期振荡的相位和振幅，因此可以从多深度温度的周期变化信号中导出地下水通量的时间变化，从而能够识别地震引起的流体流动变化。Liu等[17]利用安装在青藏高原东缘鲜水河断裂带的六个不同位置的0至20 m之间不同深度的高灵敏度温度传感器测量的五年基岩温度时间序列，借助动态谐波回归技术提取了年振荡的时变振幅和相位，估计了附近两次大地震后地下水通量变化。从而揭示了两次地震的不同水文响应。这是首次使用温度-时间序列量化受地震影响的地下水通量连续变化。更深入的探索，可望在热流变化中找寻地震前兆信号。

2.4 在冻土热传输研究中的应用

上述研究未考虑冻土的温度预报方程。冻土的热属性与非冻土差异显著，冰的热扩散率（1.1 × 10–6 m2 s–1）远大于液态水的热扩散率（0.14×10–6 m2 s–1）[34]。基于青藏高原东部的两个站点的观测，Wang等[35]评估了冻融期通过添加土壤冰蓄热量和土壤冰相变潜热对地表热通量估算的改进，发现冰的大体积热容量和土壤冰相变潜热对土壤热通量的估算有特定的影响。最近，Tong等[36]利用青藏高原三个站点两层土壤温度3年的观测数据，基于热传导-对流方程确定了土壤热扩散率，发现非冻土热扩散率与土壤湿度关系显著，然而冻土时二者关系发生变化。至今尚未有很好的观测冻土热属性的方法。而冻土热属性是研究冻融交替的土壤水热传输的必要参数，也会通过影响土壤温度和含水量间接影响冻土温室气体的排放。随着全球变暖的加剧，高纬度或者高原上的永久冻土在转变为季节性冻土或者非冻土，对下垫面和大气之间的热量和物质（水汽和温室气体等）交换有显著的影响。未来需要在土壤温度预报方程中加入冻土的影响，结合具体物理过程和观测，考虑冻土时土壤热属性的时空变化，以提升土壤温度预报方程在冻土时的表现。

3 土壤温度预报方程的研究不足与展望

相对经典的仅仅考虑了热传导的土壤热传导方程，土壤热传导-对流方程考虑了更合理的土壤热传输过程。然而，即使是物理过程和边界条件等最完善的土壤热传导-对流方程[29]，依然未完全考虑全相态（液态、固态、气态）土壤水对热传输的影响。此外，土壤中热传输方程一般是基于裸土创建的，植物根系的影响鲜有考虑。然而，除了少数下垫面，绝大多数的土壤中都含有植物根系。因此，为进一步提升土壤温度预报方程表现，下一步需要考虑不同相态土壤水和植物根系对土壤热传输的定量影响。

4 总结

综上所述，土壤温度预报方程从最初只考虑热传导过程到热传导-对流耦合方程的创建和求解，边界条件从一个正弦波演变为傅里叶级数，土壤介质被视为均一到考虑土壤热属性随深度的垂直变化，土壤温度预报方程不断发展，越来越接近真实土壤的热传输过程。接着，我们阐述了热传导-对流的土壤温度预报方程的在近地表能量平衡、土壤水垂直迁移、地震和冻土研究中的应用。同时指出，土壤热传导-对流方程下一步还需要考虑全相态的土壤水和植物根系的影响，进而将更加完备的土壤预报方程耦合入天气和气候模式，旨在提升模式表现。

附录：土壤热传导-对流方程的求解过程

为了去掉方程中$a\frac{{\partial {T^*}}}{{\partial z}}$项可做如下代换，令${T^*} = U\left( {z, t} \right)\exp \left( {mt + nz} \right)$，则

$ \frac{{\partial {T^*}}}{{\partial t}} = \left( {\frac{{\partial U}}{{\partial t}} + Um} \right)\exp \left( {mt + nz} \right), $
$ \frac{{\partial {T^*}}}{{\partial z}} = \left( {\frac{{\partial U}}{{\partial z}} + Un} \right)\exp \left( {mt + nz} \right), $
$ \frac{{{\partial ^2}{T^*}}}{{\partial {z^2}}} = \left[ {\frac{{{\partial ^2}U}}{{\partial {z^2}}} + n\frac{{\partial U}}{{\partial z}} + n\left( {\frac{{\partial U}}{{\partial t}} + Un} \right)} \right]\exp \left( {mt + nz} \right), $

将以上三方程代入$\frac{{\partial {T^*}}}{{\partial t}} = k\frac{{{\partial ^2}{T^*}}}{{\partial {z^2}}} + $ $\left( {a + b\sin \omega t} \right)\frac{{\partial {T^*}}}{{\partial z}}$，

同时方程左右两边同除$\exp \left( {mt + nz} \right)$，得

$ \begin{gathered} \frac{{\partial U}}{{\partial t}} + Um = k\left( {\frac{{{\partial ^2}U}}{{\partial {z^2}}} + 2n\frac{{\partial U}}{{\partial z}} + U{n^2}} \right) + \left( {a + b\sin \omega t} \right) \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {\frac{{\partial U}}{{\partial z}} + Un} \right), \hfill \\ \end{gathered} $

整理后

$ \begin{gathered} \frac{{\partial U}}{{\partial t}} = k\frac{{{\partial ^2}U}}{{\partial {z^2}}} + b\sin \omega t\frac{{\partial U}}{{\partial z}} + Unb\sin \omega t + \hfill \\ \;\;\;\;\;\;\;\left( {2nk + a} \right)\frac{{\partial U}}{{\partial z}} + U\left( {k{n^2} + an - m} \right), \hfill \\ \end{gathered} $ (17)

故$\left\{ {\begin{array}{*{20}{l}} {2nk + a = 0} \\ {k{n^2} + an - m = 0} \end{array}} \right.$

得：$n = - \frac{a}{{2k}}$，$m = - \frac{{{a^2}}}{{4k}}$。

故令${T^*} = U\left( {z, t} \right)\exp \left( { - \frac{{{a^2}t}}{{4k}} - \frac{{az}}{{2k}}} \right)$，则

$ \frac{{\partial {T^*}}}{{\partial t}} = \left( {\frac{{\partial U}}{{\partial t}} - U\frac{{{a^2}}}{{4k}}} \right)\exp \left( { - \frac{{{a^2}t}}{{4k}} - \frac{{az}}{{2k}}} \right), $ (18)
$ \frac{{\partial {T^*}}}{{\partial z}} = \left( {\frac{{\partial U}}{{\partial z}} - U\frac{a}{{2k}}} \right)\exp \left( { - \frac{{{a^2}t}}{{4k}} - \frac{{az}}{{2k}}} \right), $ (19)
$ \frac{{{\partial ^2}{T^*}}}{{\partial {z^2}}} = \left( {\frac{{{\partial ^2}U}}{{\partial {z^2}}} - \frac{a}{k}\frac{{\partial U}}{{\partial z}} + \frac{{{a^2}}}{{4{k^2}}}U} \right)\exp \left( { - \frac{{{a^2}t}}{{4k}} - \frac{{az}}{{2k}}} \right), $ (20)

因此，

$ \begin{gathered} \frac{{\partial U}}{{\partial t}} - U\frac{{{a^2}}}{{4k}} = k\left( {\frac{{{\partial ^2}U}}{{\partial {z^2}}} - \frac{a}{k}\frac{{\partial U}}{{\partial z}} + \frac{{{a^2}}}{{4{k^2}}}U} \right) + \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left( {a + b\sin \omega t} \right)\left( {\frac{{\partial U}}{{\partial z}} - U\frac{a}{{2k}}} \right), \hfill \\ \end{gathered} $

继而，

$ \begin{array}{l}\frac{\partial U}{\partial t}=k\left(\frac{{\partial }^{2}U}{\partial {z}^{2}}-\frac{a}{k}\frac{\partial U}{\partial z}+\frac{{a}^{2}}{4{k}^{2}}U\right)+a\frac{\partial U}{\partial z}-U\frac{{a}^{2}}{2k}+\\ b \ \mathrm{sin} \ (\omega t)\frac{\partial U}{\partial z}-b \ \mathrm{sin} \ (\omega t)\frac{a}{2k}U+U\frac{{a}^{2}}{4k}\\ =k\frac{{\partial }^{2}U}{\partial {z}^{2}}-\left(\frac{ab}{2k} \ \mathrm{sin} \ \omega t\right)U+b \ \mathrm{sin} \ \omega t\frac{\partial U}{\partial z},\end{array} $ (21)

则方程（16）变为，

$ \left\{ \begin{array}{l} \frac{\partial U}{\partial t}=k\frac{{\partial }^{2}U}{\partial {z}^{2}}-\left(\frac{ab}{2k} \ \mathrm{sin} \ \omega t\right)U\left(z, t\right)+b \ \mathrm{sin} \ \omega t\frac{\partial U}{\partial z}\\ U\left(0, t\right)=\mathrm{exp}\left(\frac{{a}^{2}t}{4k}\right)[\left({T}_{0}-{T}_{1}\right)+ \sum\limits_{i = 1}^n {A}_{i} \ \mathrm{sin} \ (i\omega t+{\Phi }_{i})], i=1, 2, \mathrm{...}, n\\ U\left(\infty , t\right)=0\\ U\left(z, 0\right)=\left[f\left(z\right)-{T}_{1}\right]\mathrm{exp}\left(\frac{az}{2k}\right)\end{array} \right. , $ (22)

为了进一步简化方程需去掉方程中$b\sin \omega t\frac{{\partial U}}{{\partial z}}$项，引进一长度以m为单位的参数

$ {p_1}\left( t \right) = \mathop \smallint \limits_0^t b\sin \omega tdt = \frac{b}{\omega }\left( {1 - \cos \omega t} \right), $ (23)

令$Z = z + {p_1}\left( t \right)$，得到$U\left( {z, t} \right) = U\left[ {Z - {p_1}\left( t \right), t} \right] = V\left( {Z, t} \right)$。

对U和V对t求偏导，可得$\frac{{\partial U}}{{\partial t}} = \frac{{\partial V}}{{\partial t}} + \frac{{\partial V}}{{\partial Z}}\frac{{\partial Z}}{{\partial t}}$由此可以得到

$ \frac{{\partial U}}{{\partial t}} = \frac{{\partial V}}{{\partial t}} + \frac{{\partial {p_1}\left( t \right)}}{{\partial t}}\frac{{\partial V}}{{\partial Z}} = \frac{{\partial V}}{{\partial t}} + b\sin \omega t\frac{{\partial V}}{{\partial Z}}, $ (24)

同样有，

$ \left\{ \begin{array}{l} \frac{\partial U}{\partial z}=\frac{\partial V}{\partial Z}\\ \frac{{\partial }^{2}U}{\partial {z}^{2}}=\frac{{\partial }^{2}V}{\partial {Z}^{2}}\end{array} \right., $ (25)

因此，结合方程（24）和方程（25），方程（22）变为，

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{\partial V}}{{\partial t}} = k\frac{{{\partial ^2}V}}{{\partial {Z^2}}} - \left( {\frac{{ab}}{{2k}}\sin \omega t} \right)V\left( {Z, t} \right)} \\ {V\left( {{p_1}\left( t \right), t} \right) = \exp \left( {\frac{{{a^2}t}}{{4k}}} \right)[\left( {{T_0} - {T_1}} \right) + \mathop \sum \limits_{i = 1}^n {A_i}\sin (i\omega t + {\Phi _i})], i = 1, 2, ..., n} \\ {V\left( {\infty , t} \right) = 0} \\ {V\left( {Z, 0} \right) = \left[ {f\left( Z \right) - {T_1}} \right]\exp \left( {\frac{{aZ}}{{2k}}} \right)} \end{array}} \right. , $ (26)

上式可通过正弦傅立叶变化求解。因此对方程和初始边界条件做关于Z的正弦傅立叶变换，即

$ V\left( {p, t} \right) = \mathop \smallint \limits_0^\infty V\left( {Z, t} \right)\sin (pZ)dZ, $ (27)

其中，p为傅立叶变换参数。通过此变化，结果变为初边值的常微分方程，

$ \left\{ {\begin{array}{*{20}{l}} {\frac{{dV\left( {p, t} \right)}}{{dt}} = - \left[ {k{p^2} + \frac{{ab}}{{2k}}\sin (\omega t} \right)]V\left( {p, t} \right) + kp\exp (\frac{{{a^2}t}}{{4k}})\left[ {\left( {{T_0} - {T_1}} \right) + \mathop \sum \limits_{i = 1}^n {A_i}\sin (i\omega t + {\Phi _i})} \right]} \\ {V\left( {p, 0} \right) = \mathop \smallint \limits_0^\infty [f\left( Z \right) - {T_1}]\exp \left( {\frac{{aZ}}{{2k}}} \right)\sin \left( {pZ} \right)dZ} \end{array}} \right. $ (28)

为了求解方程（28），我们首先用分离变量法求解齐次方程，显式解析解表示为

$ \begin{array}{l}V=\mathrm{exp}(-k{p}^{2}t+{b}_{3} \ \mathrm{cos} \ \omega t)\left(M+J+c\right)=\\ \mathrm{exp}(-k{p}^{2}t+{b}_{3} \ \mathrm{cos} \ \omega t)\left({V}_{1}+{V}_{2}+{V}_{3}+{V}_{4}\right),\end{array} $ (29)

其中，

$ {V}_{1}={b}_{1}\mathrm{exp}({b}_{2}t)\sum\limits_{i = 1}^n \frac{{b}_{2}{A}_{i} \ \mathrm{sin} \ (i\omega t+{\Phi }_{i})-i\omega {A}_{i} \ \mathrm{cos} \ (i\omega t+{\Phi }_{i})}{{b}_{2}{}^{2}+{(i\omega )}^{2}}, $ (30)
$ {V}_{2}=-{b}_{1}{b}_{3}\sum\limits_{i = 1}^n \left[ \ \mathrm{cos} \ {\Phi }_{i}{V}_{2}{}'\left(i\right)+ \ \mathrm{sin} \ {\Phi }_{i}{V}_{2}{}''\left(i\right)\right], $ (31)
$ \begin{array}{l}{V}_{2}{}'\left(i\right)=\frac{1}{2}\frac{\mathrm{exp}({b}_{2}t)}{{b}_{2}{}^{2}+{(i+1)}^{2}{\omega }^{2}}\left\{{b}_{2}{A}_{i} \ \mathrm{sin} \ [\left(i+1\right)\omega t\left]-\left(i+1\right)\omega {A}_{i} \ \mathrm{cos} \ [\left(i+1\right)\omega t\right]\right\}\\ +\frac{1}{2}\frac{\mathrm{exp}({b}_{2}t)}{{b}_{2}{}^{2}+{(i-1)}^{2}{\omega }^{2}}\left\{{b}_{2}{A}_{i} \ \mathrm{sin} \ [\left(i-1\right)\omega t\left]-\left(i-1\right)\omega {A}_{i} \ \mathrm{cos} \ [\left(i-1\right)\omega t\right]\right\},\end{array} $ (32a)
$ \begin{array}{l}{V}_{2}{}''\left(i\right)=\frac{1}{2}\frac{\mathrm{exp}({b}_{2}t)}{{b}_{2}{}^{2}+{(i+1)}^{2}{\omega }^{2}}\left\{{b}_{2}{A}_{i} \ \mathrm{cos} \ [\left(i+1\right)\omega t\left]+\left(i+1\right)\omega {A}_{i} \ \mathrm{sin} \ [\left(i+1\right)\omega t\right]\right\}\\ +\frac{1}{2}\frac{\mathrm{exp}({b}_{2}t)}{{b}_{2}{}^{2}+{(i-1)}^{2}{\omega }^{2}}\left\{{b}_{2}{A}_{i} \ \mathrm{cos} \ [\left(i-1\right)\omega t\left]+\left(i-1\right)\omega {A}_{i} \ \mathrm{sin} \ [\left(i-1\right)\omega t\right]\right\},\end{array} $ (32b)
$ {V}_{3}=\frac{{b}_{4}}{{b}_{2}}\mathrm{exp}({b}_{2}t), $ (33)
$ {V}_{4}=-{b}_{4}{b}_{3}\mathrm{exp}({b}_{2}t)\frac{{b}_{2} \ \mathrm{cos} \ (\omega t)+\omega \ \mathrm{sin} \ (\omega t)}{{b}_{2}{}^{2}+{\omega }^{2}}。$ (34)

当$t = 0$时，

$ V\left( {p, 0} \right) = \exp (\frac{{{a_1}b}}{{2k\omega }})\left[ {{V_1}\left( {p, 0} \right) + {V_2}\left( {p, 0} \right) + {V_3}\left( {p, 0} \right) + {V_4}\left( {p, 0} \right) + c} \right] $

则，

$ c = \exp ( - \frac{{ab}}{{2k\omega }})V\left( {p, 0} \right) - \left[ {{V_1}\left( {p, 0} \right) + {V_2}\left( {p, 0} \right) + {V_3}\left( {p, 0} \right) + {V_4}\left( {p, 0} \right)} \right], $

其中，

$ V\left( {p, 0} \right) = \mathop \smallint \limits_0^\infty [f\left( Z \right) - {T_1}]\exp \left( {\frac{{aZ}}{{2k}}} \right)\sin \left( {pZ} \right)dZ。$ (35)
参考文献
[1]
De Vries D A. Thermal properties of soils//van Wijk W R. Physics of plant environment [M]. Amsterdam: North-Holland Publishing Company, 1963: 210-235. (0)
[2]
Johansen O. Thermal conductivity of soils [D]. Norwegian University of Science and Technology, 1975: 291. (0)
[3]
Gao Z, Fan X, Bian L. An analytical solution to one-dimensional thermal conduction-convection in soil[J]. Soil Science, 2003, 168: 99-107. DOI:10.1097/00010694-200302000-00004 (0)
[4]
Horton R, Wierenga P J, Nielsen D R. Evaluation of methods for determination apparent thermal diffusivity of soil near the surface[J]. Soil Science Society of American Journal, 1983, 47(1): 23-32. (0)
[5]
Gao Z. Determination of soil heat flux in a Tibetan short-grass prairie[J]. Boundary-Layer Meteorology, 2005, 114: 165-178. DOI:10.1007/s10546-004-8661-5 (0)
[6]
Gao Z, Chen T J, Hu Y. Impact of soil vertical water movement on the energy balance of different land surfaces[J]. International Journal of Biometeorology, 2007, 51: 565-573. DOI:10.1007/s00484-007-0095-6 (0)
[7]
Gao Z, Lenschow D H, Horton R, et al. Comparison of two soil temperature algorithms for a bare ground site on the Loess Plateau in China[J]. Journal of Geophysics Research: Atmospheric, 2008, 113: D18105. DOI:10.1029/2008JD010285 (0)
[8]
Gao Z, Horton R, Wang L, et al. An improved force-restore method for soil temperature prediction[J]. European Journal of Soil Science, 2008, 59: 972-981. DOI:10.1111/j.1365-2389.2008.01060.x (0)
[9]
Gao Z, Horton R, Liu H P. Impact of wave phase difference between soil surface heat flux and soil surface temperature on soil surface energy balance closure[J]. Journal of Geophysics Research: Atmospheric, 2010, 115: D16112. DOI:10.1029/2009JD013278 (0)
[10]
Gao Z, Tong B, Horton R, et al. Determination of desert soil apparent thermal diffusivity using a conduction‐convection algorithm[J]. Journal of Geophysical Research: Atmospheres, 2017, 122: 9569-9578. DOI:10.1002/2017JD027290 (0)
[11]
Wang L, Gao Z, Horton R. Comparison of six algorithms to determine the soil apparent thermal diffusivity at a site in the Loess Plateau of China[J]. Soil Science, 2010, 175(2): 51-60. DOI:10.1097/SS.0b013e3181cdda3f (0)
[12]
Huang F, Zhan W, Ju W, et al. Improved reconstruction of soil thermal field using two-depth measurements of soil temperature[J]. Journal of Hydrology, 2014, 519: 711-719. DOI:10.1016/j.jhydrol.2014.08.014 (0)
[13]
Liu Q Q, Wei D P, Sun Z T, et al. Review of related researches on calculation of shallow ground heat transfer (In Chinese)[J]. Progress in Geophysics, 2014, 29(6): 2510-2517. [刘迁迁, 魏东平, 孙振添, 等. 浅层地表热传输计算方法综述[J]. 地球物理学进展, 2014, 29(6): 2510-2517.] (0)
[14]
Falocchi M, Barontini S, Ranzi R. A parametrization of a steady periodic solution of the Fourier equation to model soil temperature dynamics[J]. Journal of Geophysical Research: Earth Surface, 2015. DOI:10.1002/2015jf003664 (0)
[15]
Gao Z, Russell E S, Missik J E C, et al. A novel approach to evaluate soil heat flux calculation: An analytical review of nine methods[J]. Journal of Geophysical Research: Atmospheres, 2017, 122: 6934-6949. DOI:10.1002/2017JD027160 (0)
[16]
Liu Q, Chen S, Jiang L, et al. Determining thermal diffusivity using near-surface periodic temperature variations and its implications for tracing groundwater movement at the eastern margin of the Tibetan Plateau[J]. Hydrological Processes, 2019, 33(8): 1276-1286. DOI:10.1002/hyp.13399 (0)
[17]
Liu Q, Chen S, Chen L, et al. Detection of groundwater flux changes in response to two large earthquakes using long-term bedrock temperature time series[J]. Journal of Hydrology, 2020, 590: 125245. DOI:10.1016/j.jhydrol.2020.125245 (0)
[18]
Tong B, Xu H, Horton R, et al. Determination of long-term soil apparent thermal diffusivity using near-surface soil temperature on the Tibetan Plateau[J]. Remote Sensing, 2022, 14(17): 4238. DOI:10.3390/rs14174238 (0)
[19]
Ioannidis T, Bakas N A. An analytical solution for the heat conduction–convection equation in non-homogeneous soil[J]. Boundary-Layer Meteorology. DOI:10.1007/s10546-022-00753-2 (0)
[20]
Deardorff J W. Efficient prediction of ground surface temperature and moisture inclusion of a layer of vegetation[J]. Journal of Geophysical Research, 1978, 83: 1889-1903. DOI:10.1029/JC083iC04p01889 (0)
[21]
Lin J D. On the force-restore method for prediction of ground surface temperature[J]. Journal of Geophysical Research, 1980, 85: 3251-3254. DOI:10.1029/JC085iC06p03251 (0)
[22]
Dickinson R E. The force-restore model for surface temperature and its generalization[J]. Journal of Climate, 1988, 1: 1086-1097. DOI:10.1175/1520-0442(1988)001<1086:TFMFST>2.0.CO;2 (0)
[23]
Carslaw H S, Jaeger J C. Conduction of heat in solids . 2nd ed[M]. UK: Oxford University Press, 1959: 510. (0)
[24]
Fan X G, Tang M C. A preliminary study on conductive and convective soil heat flux (In Chinese)[J]. Plateau Meteorology, 1994, 13(1): 14-19. [范新岗, 汤懋苍. 土壤传导-对流热通量计算的初步结果[J]. 高原气象, 1994, 13(1): 14-19.] (0)
[25]
Garratt J R. The atmospheric boundary layer [M]. Cambridge University Press, 1992: 117-119, 291-292. (0)
[26]
Rybach L, Muffler P. Geothermal systems: Principles and case histories . Hoboken: John Wiley and Sons, 1981. (0)
[27]
Stull R B. An introduction to Boundary Layer Meteorology [M]. Atmospheric Sciences Library, Springer Science & Business Media, 1988: 666. (0)
[28]
Hillel D. Introduction to soil physics . New York: Academic Press, 1982: 364. (0)
[29]
Wang L, Gao Z, Horton R, et al. An analytical solution to the one-dimensional heat conduction-convection equation in soil[J]. Soil Science Society of America Journal, 2012, 76(6): 1978-1986. DOI:10.2136/sssaj2012.0023N (0)
[30]
Shao M, Horton R, Jaynes D B. Analytical solution for one-dimensional heat conduction–convection equation[J]. Soil Science Society of America Journal, 1998, 62(1): 123-128. DOI:10.2136/sssaj1998.03615995006200010016x (0)
[31]
Jackson R D, Kirkham D. Method of measurement of the real thermal diffusivity of moist soil[J]. Soil Science Society of America Journal, 1958, 22(6): 479-482. DOI:10.2136/sssaj1958.03615995002200060001x (0)
[32]
van Wijk W R. Physics of plant environment . Amsterdam: North-Holland Publishing Co., 1963. (0)
[33]
Hu G, Zhao L, Wu X, et al. New Fourier-series-based analytical solution to the conduction–convection equation to calculate soil temperature, determine soil thermal properties, or estimate water flux[J]. International Journal of Heat and Mass Transfer, 2016, 95: 815-823. DOI:10.1016/j.ijheatmasstransfer.2015.11.078 (0)
[34]
Farouki O T. Thermal properties of soils, Special Report 81-1 [R]. Cold Regions Research and Engineering Lab Hanover, NH. 1981. (0)
[35]
Wang J, Luo S, Lv Z, et al. Improving ground heat flux estimation: Considering the effect of freeze/thaw process on the seasonally frozen ground[J]. Journal of Geophysical Research: Atmospheres, 2021, 126: e2021JD035445. (0)
[36]
Tong B, Gao Z, Horton R, et al. Soil apparent thermal diffusivity estimated by conduction and by conduction-convection heat transfer models[J]. Journal of Hydrometeorology, 2017, 18(1): 109-118. DOI:10.1175/JHM-D-16-0086.1 (0)