重力观测通常会受到降雨、蒸发、地下水渗透和流动等影响，尤其是局部地区地下水变化的干扰^[1-4]。因此，为防止地下水重力效应对地球物理学或动力学信号的掩盖，有必要对精密重力观测进行地下水改正。许多学者对地下水重力效应进行模拟^[5-8]，其中Kazama等^[7]收集日本北部Isawa扇形地区水文、气象资料，利用一维地下水渗透方程，模拟该地区超导台站局部区域土壤含水率的时空分布及其导致的重力变化，取得与超导重力观测吻合较好的结果。但在其使用有限差分法解算渗透方程时采用显式差分格式，虽然计算简单，但受到收敛条件的限制，使得时间和空间步长的选取不自由。有限差分法除显式差分外，还有隐式差分和中心差分两种无条件收敛的差分格式^[9]。本文在Kazama等^[7]的基础上，探讨不同差分格式对地下水模拟结果以及相应重力改正的影响。

1 一维地下水渗透方程

假设土壤性质均匀，且不考虑水平方向上的地下水流动，根据质量守恒原理和达西定律，有一维非饱和层的非线性渗透方程，即Richards方程^[9-10]：

$ \begin{array}{*{20}{c}} {\frac{{\partial \theta }}{{\partial t}} = \frac{\partial }{{\partial z}}\left[ {D\left( \theta \right)\frac{{\partial \theta }}{{\partial z}} + K\left( \theta \right)} \right] = }\\ {\frac{\partial }{{\partial z}}\left[ {D\left( \theta \right)\frac{{\partial \theta }}{{\partial z}}} \right] + \frac{{\partial K\left( \theta \right)}}{{\partial z}}} \end{array} $

(1)

式中，θ为土壤含水率，t为时间，z为垂直方向坐标。K(θ)和D(θ)分别为垂向渗透系数和扩散系数，均是含水率的函数^[11-12]，且与土壤性质参数有关。土壤性质参数包括垂向饱和渗透系数K_S、饱和扩散系数D_S、常数项a、b以及最大、最小含水率θ_max和θ_min，取值见表 1^[7]。

渗透方程的求解需要加入一定的边界条件和初始条件。这里给定的2个边界条件分别是地表处的渗透速率和地下水面处的含水率，初始条件设定为定水头稳定流条件下的稳态解^[7]。

2 有限差分解算

按深度把土壤均匀分为若干层，用θ_i^j、v_i^j分别表示垂向第i层某时刻t_j的含水率和渗透速率，利用有限差分法可将式(1)差分化。

2.1 显式差分

Kazama等^[7]的显式差分用时间上的向前差分公式代替$ \frac{{\partial \theta }}{{\partial t}} $，用空间上的一阶中心差分公式代替$ \frac{\partial }{{\partial z}}\left[{D\left( \theta \right)\frac{{\partial \theta }}{{\partial z}}} \right] $和$ \frac{{\partial K\left( \theta \right)}}{{\partial z}} $，舍去无穷小量后有：

$ \begin{array}{*{20}{c}} {\frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} = D_{i + 1/2}^j\frac{{\theta _{i + 1}^j - \theta _i^j}}{{\left[ {\Delta z} \right]^2}} + D_{i - 1/2}^j\frac{{\theta _{i - 1}^j - \theta _i^j}}{{{{\left[ {\Delta z} \right]}^2}}} + }\\ {\frac{{K_{i + 1/2}^j - K_{i + 1/2}^j}}{{\Delta z}}} \end{array} $

(2)

式中，θ_i^j+1是垂向第i层下一个时刻t_j+1的含水率，可以由该时刻t_j的含水率分布直接求得。式(2)的截断误差为O(Δt)+O(Δz)。

根据求解的地下水分布，利用牛顿引力定律和平面半无限空间近似计算地下水重力效应^{[7, 13]}：

$ {g_{\rm{w}}} = 2{\rm{ \mathsf{ π} }}{\rho _{\rm{w}}}G\Delta z \cdot \sum\limits_{i = 1}^N {{\theta _i}} $

(3)

式中，g_w为地下水重力效应，Δz是水平含水层厚度，G为万有引力常数，ρ_w为水的密度。

式(2)中系数D_i+1/2^j、D_i－1/2^j、K_i+1/2^j、K_i－1/2^j都不是整数层的参数值，涉及到层间参数的取值问题。层间参数取值的加权公式有多种^[14]。

1) 算术平均值：

$ D_{i \pm 1/2}^j = \frac{{D_i^j + D_{i \pm 1}^j}}{2},K_{i \pm 1/2}^j = \frac{{K_i^j + K_{i \pm 1}^j}}{2} $

(4)

2) 几何平均值：

$ D_{i \pm 1/2}^j = \sqrt {D_i^j \times D_{i \pm 1}^j} ,K_{i \pm 1/2}^j = \sqrt {K_i^j \times K_{i \pm 1}^j} $

(5)

3) 调和平均值：

$ D_{i \pm 1/2}^j = 2\frac{{D_i^j \times D_{i \pm 1}^j}}{{D_i^j + D_{i \pm 1}^j}},K_{i \pm 1/2}^j = 2\frac{{K_i^j \times K_{i \pm 1}^j}}{{K_i^j + K_{i \pm 1}^j}} $

(6)

4) 上游值：

$ D_{i \pm 1/2}^j = \left\{ {\begin{array}{*{20}{c}} {D_{i \pm 1}^j,D_{i \pm 1}^j \ge \theta _i^j}\\ {D_i^j,\theta _i^j > \theta _{i \pm 1}^j} \end{array}} \right.,K_{i \pm 1/2}^j = \left\{ {\begin{array}{*{20}{c}} {K_{i \pm 1}^j,\theta _{i \pm 1}^j \ge \theta _i^j}\\ {K_i^j,\theta _i^j > \theta _{i \pm 1}^j} \end{array}} \right. $

(7)

为比较不同加权公式的模拟效果，采用相同的差分格式、相同的模型参数、相同的初始条件和不同的层间加权公式，分别对渗透问题进行模拟计算，获得计算区域(地表至地下10 m)的土壤含水率分布及其重力效应。

图 1为不同层间参数加权公式对地下50 cm处含水率模拟结果的影响^[7]。可以发现，Kazama等^[7]的模拟结果与采用几何平均值法的结果差异为零，两者具有等价性；与上游值法差异最大，达到0.01；与调和平均值法差异次之；与算数平均值法差异更小，小于0.001。

图 2为不同层间参数加权公式对地下水重力效应模拟的影响^[7]，表 2给出其标准差。可以发现，图 2与图 1有着类似的特征，Kazama等^[7]的模拟结果与几何平均值法差异为零，而与上游值法差异最大，最大达0.15 μGal，标准差为0.064 μGal；与调和平均值法差异次之，标准差为0.024 μGal；与算数平均值法差异更小，小于0.05 μGal，标准差为0.015 μGal。相对于该台站变化幅度达8 μGal的地下水重力效应^[7]来说，不同加权公式对地下水重力效应的影响在1.9%以内。

2.2 隐式差分和中心差分

隐式差分用时间上的向后差分公式代替$ \frac{{\partial \theta }}{{\partial t}} $，舍去无穷小量后有：

$ \begin{array}{*{20}{c}} {\frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} = D_{i + 1/2}^{j + 1}\frac{{\theta _{i + 1}^{j + 1} - \theta _i^{j + 1}}}{{{{\left[ {\Delta z} \right]}^2}}} + D_{i - 1/2}^{j + 1}\frac{{\theta _{i - 1}^{j + 1} - \theta _i^{j + 1}}}{{{{\left[ {\Delta z} \right]}^2}}} + }\\ {\frac{{K_{i + 1/2}^{j + 1} - K_{i - 1/2}^{j + 1}}}{{\Delta z}}} \end{array} $

(8)

式中，θ_i－1^j+1、θ_i^j+1和θ_i+1^j+1分别是垂向第i-1、i和i+1层下一个时刻t_j+1的含水率，不可以由该式直接从时刻t_j的含水率分布求得，而要通过联立形成三对角方程组，用追赶法进行求解。式(8)的截断误差也为O(Δt)+O(Δz)。

中心差分用时间上的中心差分公式代替$ \frac{{\partial \theta }}{{\partial t}} $，舍去无穷小量后有：

$ \begin{array}{*{20}{c}} {\frac{{\theta _i^{j + 1} - \theta _i^j}}{{\Delta t}} = D_{i + 1/2}^{j + 1/2}\frac{{\theta _{i + 1}^{j + 1/2} - \theta _i^{j + 1/2}}}{{{{\left[ {\Delta z} \right]}^2}}} + }\\ {D_{i - 1/2}^{j + 1/2}\frac{{\theta _{i - 1}^{j + 1/2} - \theta _i^{j + 1/2}}}{{{{\left[ {\Delta z} \right]}^2}}} + \frac{{K_{i + 1/2}^{j + 1/2} - K_{i - 1/2}^{j + 1/2}}}{{\Delta z}}} \end{array} $

(9)

式中，θ_i－1^j+1/2、θ_i^j+1/2和θ_i+1^j+1/2分别是垂向第i-1、i和i+1层在t_j+1/2时刻的含水率，用它们在t_j+1和t_j时刻含水率的平均值代替，对D_i+1/2^j+1/2、D_i－1/2^j+1/2、K_i+1/2^j+1/2和K_i－1/2^j+1/2采用同样的处理。该式也无法直接求解，需通过联立形成三对角方程组，用追赶法进行求解。式(9)的截断误差为O([Δt]²)+O(Δz)，与式(8)具有相同的形式，两者的区别在于三对角方程的系数和常数项的不同。中心差分格式的截断误差比显式和隐式差分格式的要小，因此具有更高的精度。

由式(8)和(9)可知，隐式差分与中心差分公式中的扩散系数和渗透系数都是待求含水率的函数，因此这两方程是非线性的，这就牵涉到方程的线性化问题。常用的线性化方法有显式法、预测-校正法和迭代法^[9]。

以隐式差分法为例，对3种线性化方法进行说明。显式法是用上一时刻的参数值D_i^j(或K_i^j)代替方程中的D_i^j+1(或K_i^j+1)。预测-校正法是先对未知数θ_i^j+1进行预测，然后在其基础上得到校正值。这里先用显式差分法得到t_j+1时刻的预测值θ_i^j+1，根据土壤参数与含水率的关系曲线得到D_i^j+1(或K_i^j+1)的值，然后按隐式差分解算三对角方程组得到t_j+1时刻的校正值θ_i^j+1。迭代法是先取t_j时刻的D_i^j(或K_i^j)作为D_i^j+1(或K_i^j+1)的估计值，然后按照隐式差分解算三对角方程组，得到t_j+1时刻含水率θ_i^j+1的第一次迭代值θ_i^j+1(1)，根据θ_i^j+1(1)及土壤参数与含水率的关系曲线得到新的D_i^j+1(或K_i^j+1)值，再次按照隐式差分解算三对角方程组，得到第二次迭代值θ_i^j+1(2)。重复迭代过程，直到各层前后两次迭代值之差小于所规定的允许误差e为止，即

$ \max \left| {\theta _i^{j + 1\left( m \right)} - \theta _i^{j + 1\left( {m - 1} \right)}} \right| \le e $

(10)

为比较不同差分形式和线性化方法的结果，采用相同的层间加权公式、模型参数、初始条件和不同的差分形式或线性化方法，分别对Kazama等^[7]中的渗透问题进行模拟计算。表 3列出7种差分形式和线性化方法的组合，其中Kazama等^[7]使用的是组合形式1(显式差分法)。图 3为不同差分格式和线性化方法组合形式对地下水重力效应模拟的影响，表 4给出了其标准差。从图 3看，各曲线随着时间增加慢慢分散，差异越来越大，说明各方法之间的差异存在一定的累积效应。在计算的639 d里，与组合形式1(显式差分法)的差异中，最大的是组合形式4(隐式迭代法)，其最大差异值为0.12 μGal；而标准差最大的是组合形式7(中心迭代法)，其差异的标准差为0.022 μGal。相对于该台站变化幅度达8 μGal的地下水重力效应而言，不同计算方法对地下水重力效应的影响在1.5%以内。

3 结语

本文比较了不同差分格式对一维地下水渗透模拟结果及其相应地下水重力效应的影响。在日本北部Isawa扇形地区超导台站，不同层间参数加权公式最大能够引起约0.15 μGal的重力效应差异，影响在1.9%以内；不同差分格式和线性化方法组合最大能够引起约0.12 μGal的重力效应差异，影响在1.5%以内。该结果可为一维地下水模拟及重力效应改正算法的选取提供参考。