土壤溶质异常输运的时间分数阶对流扩散模型

A Time-Fractional Advection-Diffusion Model for Anomalous Diffusion of Solute in Soil

摘要: 土壤系统经常表现出复杂的性质并导致溶质迁移的异常扩散。基于Skaggs等人的模型，本文研究开发蒸腾和根系吸水条件下的时间分数阶对流扩散方程(FADE)模型，以模拟根区的异常扩散并进行解析求解。模拟表明，时间分数阶对流扩散模型与整数阶对流扩散模型的数值结果在表面土壤附近出现偏差，随后随着时间的推移逐渐向下移动，偏差随深度逐渐扩大，较小的α对应较高的浓度曲线，说明土壤中溶质储层较强，导致溶质运移速度较慢，即存在亚扩散。

Abstract: Soil systems often exhibit complex properties and lead to abnormal diffusion of solute transport. Based on the model of Skaggs et al., this paper develops a time fractional advection-diffusion equation (FADE) model under transpiration and root water absorption conditions to simulate abnormal diffusion in the root zone and solves it analytically. The simulation shows that the numerical results of the time fractional advection-diffusion model and the integer advection-diffusion model deviate near the surface soil, and then gradually move downward with time. The deviation gradually expands with depth, and the smaller one corresponds to a higher concentration curve, indicating that the solute reservoir in the soil is strong, resulting in a slower solute migration rate, that is, there is sub-diffusion.

1. 引言

人们越来越关注土壤中化学品、重金属和肥料的合理使用。根区是植被吸收土壤–植物–大气连续体(SPAC)资源的主要区域[1]。了解和量化根区溶质运输对于预测植物养分吸收、肥料运输至关重要，从而有助于资源利用效率、可持续农业和环境保护。

单根或根系统对水和养分的吸收是从两个角度进行建模的。在单根建模中，将单根视为均匀的圆柱形结构，水和养分的吸收集中在根表面，并在根表面形成浓度梯度，导致水和养分向根扩散，将根对水分和养分的吸收包含在方程的汇项[2] [3]。在根系建模中，水由于重力、蒸腾作用、蒸发、降雨和灌溉而向下移动，水流由理查兹方程描述，水流驱动的养分输送由另一个对流扩散方程描述，其中汇项代表根系吸收的养分[4] [5]。由于理查兹方程的非线性使得计算变得复杂，所以一些根系模型也将含水量视为常数，但根系垂直吸水量是根密度的函数，因此养分吸收方程变为变系数方程，例如Skaggs等人的模型。

土壤系统通常表现出异质性，导致溶质的非菲克或异常传输[6]。例如，微孔、粘土矿物和其他土壤基质区域往往会吸附溶质颗粒，延迟化学移动并导致溶质迁移的亚扩散非菲克输送；但由根和互连裂缝组成的大孔隙可以为水提供非平衡优先路径，水流增强导致溶质的超扩散非菲克输送[7] [8]。异常扩散通常不满足菲克定律。因此，用遵循菲克定律的经典对流扩散模型来有效地描述根区溶质的非菲克或异常传输是不合理的。鉴于此，非局域传输模型(主要包括分数阶对流扩散模型)是解释非Fickian传输过程的可行模型[9]。分数阶对流扩散模型主要有三种类型：时间分数阶对流扩散、空间分数阶对流扩散和时空分数阶对流扩散，分别可以捕获异常的亚扩散、超扩散及其混合[10]。Skaggs等人开发了一种描述根区溶质迁移的对流扩散模型，但无法捕捉溶质迁移中的非菲克现象[11]。本研究尝试修改Skaggs等人的模型来模拟根区根系溶质的异常输运。

2. 模型和解析方法

2.1. 溶质运输模型

Skaggs等人[11]提出了一个具有可变系数和有限域内的混合边界条件的对流扩散模型。在该模型中土壤中的含水量近似为常数，根系吸水引起附近的水梯度(水压)，从而改变扩散系数和对流系数。

$R\frac{\partial c}{\partial t}=\frac{\partial }{\partial z}\left[D\left(z\right)\frac{\partial c}{\partial z}\right]-\left[v_0-v_{1}B\left(z\right)\right]\frac{\partial c}{\partial z}+\left[\left(1-\gamma \right)v_{1}b\left(z\right)\right]c,$ (1)

$v_{0}c\left(z,t\right)-D_0\frac{\partial c}{\partial z}=v_0{c}_0,z=0,$ (2)

$\frac{\partial c}{\partial z}=0,z=L,$ (3)

$c\left(z,t\right)=c_1,t=0,$ (4)

其中，$b\left(z\right)$ 是归一化根系吸水密度，它是关于深度的指数函数，$B\left(z\right)$ 是累积吸水分布，$D\left(z\right)$ 是扩散系数，

$b\left(z\right)=\left(1/\delta \right){\text{e}}^{-z/\delta },B\left(z\right)={\displaystyle {\int }_0^{z}b\left(\epsilon \right)\text{d}\epsilon },D\left(z\right)={\alpha }_L\left(v_0-v_{1}B\left(z\right)\right),$ (5)

其他参数和变量参见[11]。引入无量纲变量，

$C=c/c_0,T=v_{0}t/RL,Z=z/L,$ (6)

得到模型(1)~(4)的无量纲形式，

$\frac{\partial C}{\partial T}=\frac{1}{Pe}\frac{\partial C}{\partial T}-V_1\left(Z\right)\frac{\partial C}{\partial Z}+V_2\left(Z\right)\frac{\partial C}{\partial Z}+Q\left(Z\right)C,$ (7)

$C\left(Z,T\right)-p_0\frac{\partial C}{\partial Z}=1,Z=0,$ (8)

$\frac{\partial C}{\partial Z}=0,Z=1,$ (9)

$C\left(Z,T\right)=c_1/c_0,T=0,$ (10)

其中，$Pe=v_0{\alpha }_L/D\left(Z\right)$ 是Péclet数，$V_1\left(Z\right)={\alpha }_L{v}_1/v_0$ ，$V_2\left(Z\right)=\left[v_0-v_{1}B\left(z\right)\right]v_0$ ，$Q\left(Z\right)=L\left[\left(1-\gamma \right)v_{1}b\left(z\right)\right]/v_0$ ，$p_0=D_0/v_{0}L$ 。

考虑到土壤异质性对溶质运输的影响，我们将在接下来的小节中基于模型(1)-(4)建立时间分数阶对流扩散模型。

2.2. 时间分数阶对流扩散模型

在正常扩散(即布朗运动)中，均方位移(MSD)与时间呈线性关系，

$\langle x^2\left(t\right)\rangle \sim t.$

在异常扩散中，MSD是时间的幂函数，

$\langle x^2\left(t\right)\rangle \sim t^{\varphi },t\ne 1.$

我们引用时间分数阶导数，并将方程(8)重写为分数阶版本，用于捕捉亚扩散，

$\frac{{\partial }^{\beta }C}{\partial T^{\beta }}=\frac{1}{Pe}\frac{{\partial }^{2}C}{\partial T^2}-V_1\left(Z\right)\frac{\partial C}{\partial T}+V_2\left(Z\right)\frac{\partial C}{\partial Z}+Q\left(Z\right)C,$ (11)

其中，我们采用Caputo分数阶时间导数，而不是Riemann-Liouville导数，因为后者会产生负解，而这在溶质迁移中是不允许的[12]。Caputo分数阶时间导数由下式表示：

$\frac{{\partial }^{\alpha }C}{\partial T^{\alpha }}=\frac{1}{\Gamma \left(1-\alpha \right)}{\displaystyle {\int }_0^T\frac{\partial C\left(Z,\tau \right)}{\partial \tau }\left(T-\tau \right)\text{d}\tau },$ (12)

其中，$\alpha$ 是时间分数阶数($0\le \alpha \le 1$ )，${\partial }^{\alpha }C/\partial T^{\alpha }$ 是Caputo分数阶导数，$\Gamma (·)$ 是伽马函数。

因为时间分数阶模型的扩散跳跃仍然遵循经典的Fickian扩散，因此时间分数阶模型的边界条件和初始条件与条件(8)~(10)相同。在此，我们不在复写边界条件。

目前，大多数解析方法处理常数系数时间分数阶方程[13] [14]，而不是变系数时间分数阶方程，例如(1)和(7)。我们将尝试用拉普拉斯变换和数值逆拉普拉斯变换近似求解模型(11)、(8)~(10)，然后通过比较模型的数值和解析模拟结果，证明模型(11)、(8)~(10)的数值格式具有较高的精度。

解析方法

模型(9)、(6)~(8)的拉普拉斯变换为

$\left(s\overline{C}-C\left(Z,0\right)\right)s^{\alpha -1}=\frac{1}{Pe}\frac{{\partial }^2\overline{C}}{\partial Z^2}-V_1\left(Z\right)\frac{\partial \overline{C}}{\partial Z}+V_2\left(Z\right)\frac{\partial \overline{C}}{\partial Z}+Q\left(Z\right)\overline{C},$ (13)

$\overline{C}\left(0,s\right)-p_0\frac{\partial \overline{C}}{\partial Z}=\frac{1}{s},Z=0,$ (14)

$\frac{\partial \overline{C}}{\partial Z}=0,Z=1,$ (15)

其中，s，$\overline{C}\left(Z,s\right)$ 分别是时间t和$C\left(Z,s\right)$ 的Laplace变换。模型(13)~(15)的解为

$\overline{C}\left(Z,s\right)=\frac{s^{\alpha -1}{c}_1/c_0}{s^{\alpha }-Q\left(Z\right)}+{\text{e}}^{H_1\left(Z,s\right)Z}A{}_1\left(s\right)+{\text{e}}^{H_2\left(Z,s\right)Z}{A}_2\left(s\right),$ (16)

其中，

$H_1\left(Z,s\right)=\frac{Pe\left(V_1\left(Z\right)+V_2\left(Z\right)\right)-\sqrt{Pe\left(-4Q\left(Z\right)+4s^{\alpha }+Pe{\left(V_1\left(Z\right)+V_2\left(Z\right)\right)}^2\right)}}{2},$ (17)

$H_2\left(Z,s\right)=\frac{Pe\left(V_1\left(Z\right)+V_2\left(Z\right)\right)+\sqrt{Pe\left(-4Q\left(Z\right)+4s^{\alpha }+Pe{\left(V_1\left(Z\right)+V_2\left(Z\right)\right)}^2\right)}}{2},$ (18)

$A_1\left(s\right)=\frac{{\text{e}}^{H_2\left(Z,s\right)}\left(Q\left(Z\right)\left(s+\left(-1+c_1/c_0\right)s^{\alpha }\right)\right)H_2\left(Z,s\right)}{s\left({\text{e}}^{H_2\left(Z,s\right)}{H}_2\left(Z,s\right)\left(-1+p_0{H}_1\left(Z,s\right)\right)+{\text{e}}^{H_1\left(Z,s\right)}\left(H_1-H_1{H}_2{p}_0\right)\right)\left(s^{\alpha }-Q\left(Z\right)\right)},$ (19)

$A_2\left(s\right)=\frac{{\text{e}}^{H_2\left(Z,s\right)}\left(Q\left(Z\right)\left(s+\left(-1+c_1/c_0\right)s^{\alpha }\right)\right)H_2\left(Z,s\right)}{s\left({\text{e}}^{H_1\left(Z,s\right)}{H}_1\left(Z,s\right)\left(-1+p_0{H}_2\left(Z,s\right)\right)+{\text{e}}^{H_2\left(Z,s\right)}\left(H_2-H_1{H}_2{p}_0\right)\right)\left(s^{\alpha }-Q\left(Z\right)\right)}.$ (20)

式(16)结构复杂，难以直接利用留数定理和复积分对其取逆拉普拉斯变换。我们打算利用Zakian逆拉普拉斯变换来求得其解。从而得到式(16)的解析解为

$C\left(Z,T\right)=\frac{2}{T}{\displaystyle \sum _{m=1}^n\mathrm{Re}\left\{K_m\left(\frac{{\left(a_m/T\right)}^{\alpha -1}{c}_1/c_0}{{\left(a_m/T\right)}^{\alpha }-Q\left(Z\right)}+{\text{e}}^{H_{1}Z}{A}_1\left(a_m/T\right)+{\text{e}}^{H_{2}Z}{A}_2\left(a_m/T\right)\right)\right\},}$ (21)

其中，$K_m$ 和$a_m$ 的参数可以在文献[13]中找到。

一些工作已经证明拉普拉斯变换和Zakian逆拉普拉斯变换可以为对流扩散方程问题提供精确的结果[15] [16]。因此，当经典和分数阶对流扩散方程问题的数值模拟结果足够近似时，可以认为数值格式具有较高的精度，这种近似是由分数阶而不是差分格式引起的。

3. 数值模拟

在数值方法中，我们对一阶空间导数项应用向后差分，对二阶空间导数项应用二阶中心差分，时间分数阶导数项应用L-1格式来近似，模型(11)、(8)~(10)的数值格式请参考附录。附录中只要满足在时间和空间步长足够小的情况下，其数值解会一致收敛于精确解。模型参数值来自于文献[11]，例如，$R=1$ ，$v_0=1\text{cm}\cdot {\text{d}}^{-1}$ ，$\gamma =0.1$ ，$v_1=0.75\text{cm}\cdot {\text{d}}^{-1}$ ，${\alpha }_L=10\text{ cm}$ ，$\delta =10\text{ cm}$ ，$L=\text{120 cm}$ 。模型的计算是建立在无量纲的基础上，但为了实际方便，我们图中仍以有量纲形式显示。

Figure 1. Numerical and analytical solutions of the solute concentration with depth calculated by models (11), (8)~(10) at t = 1 day

图1. 在时间t = 1天模型(11)、(8)~(10)计算得出的溶质浓度与深度关系的数值和解析解曲线

(a) t = 5 (b) t = 50

Figure 2. The profiles of solute concentration with depth calculated by model (11), (8)~(10) at different simulated time

图2. 在不同模拟时间下模型(11)、(8)~(10)计算得出的溶质浓度与深度关系的曲线

(a) z = 3 (b) z = 15

Figure 3. The curves of solute concentration versus time at different depths calculated by models (11), (8)~(10), with a simulation time of 50 days

图3. 在模拟时间为50天模型(11)、(8)~(10)计算得出的不同深度下溶质浓度与时间关系的曲线

模型(11)、(8)~(10)的解析解(21)与数值解除图1中坡度较大的根面附近外，其余情况一致，此例说明数值格式(A2)~(A5)具有较高的精度。时间分数阶对流扩散模型(11)、(8)~(10)与整数阶对流扩散模型的数值结果出现了偏差：首先出现在表面土壤附近，随后随着时间的推移逐渐向下移动(图2)，偏差随深度逐渐扩大(图3)，较小的$\alpha$ 对应较高的浓度曲线，说明土壤中溶质储层较强，导致溶质运移速度较慢，即存在亚扩散，但并不明显(图2、图3)。

4. 结论

本文研究扩展了Skaggs等人的工作，用于模拟多孔和土壤中水通量向下的根系吸收养分的情况。我们建立了一个时间分数阶对流扩散模型，用Caputo导数来描述亚扩散方向。本文利用拉普拉斯变换和数值逆拉普拉斯变换求解时间分数阶对流扩散模型，结果表明，时间分数阶对流扩散模型与整数阶对流扩散模型的结果在表面土壤附近出现偏差，且偏差随着时间的推移逐渐向下移动，但在短期模拟时间内亚扩散现象不明显。本项研究模型以含水量均匀为假设条件进行建模的，这使得模型在运用于实际农田时非常有限，在未来工作中，我们将对这一问题进行深入研究。

附 录

在本附录中，我们将给出模型(9)、(6)~(8)的数值格式。将空间域和时间域均分成网格点，$Z_i=i\Delta Z\left(i=1,2,\cdots ,N\right)$ 为空间网格点，$\Delta Z$ 是空间步长。$T^n=n\Delta t\left(n=0,1,\cdots \right)$ 为时间网格点，$\Delta t$ 是时间步长。我们用$C_i^n$ 表示网格点上$C\left(Z_i,T^n\right)$ 的近似值。

时间分数阶导数采用L-1格式近似[17]，

$\frac{{\partial }^{\alpha }C}{\partial T^{\alpha }}=\frac{\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \sum _{j=0}^n\left[{\left(j+1\right)}^{1-\alpha }-j^{1-\alpha }\right]}\left[C_i^{n+1-j}-C_i^{n-j}\right]=\frac{\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \sum _{j=0}^n{b}_j^{\alpha }}\left[C_i^{n+1-j}-C_i^{n-j}\right],$ (A1)

其中，$b_j^{\alpha }={\left(j+1\right)}^{1-\alpha }-j^{1-\alpha }$ 。

从而，模型(11)、(8)~(10)的数值格式为

$\frac{\Delta t^{-\alpha }}{\Gamma \left(2-\alpha \right)}{\displaystyle \sum _{k=1}^n{b}_k^{\alpha }\left[C_i^{n+1-k}-C_i^{n-k}\right]}=\frac{1}{Pe}\frac{C_{i+1}^n-2C_i^n+C_{i-1}^n}{\Delta Z^2}-V_1\left(Z_i\right)\frac{C_i^n-C_{i-1}^n}{\Delta Z}+V_2\left(Z_i\right)\frac{C_{i+1}^n-C_i^n}{\Delta Z}+Q\left(Z_i\right)C_i^n,$ (A2)

$C_1^n-p_0\frac{C_2^n-C_1^n}{\Delta Z}=1,$ (A3)

$\frac{C_N^n-C_{N-1}^n}{\Delta Z}=0,$ (A4)

$C_i^0=c_1/c_0.$ (A5)