说明书考虑Lisse现象的浅层地下水位预测方法
技术领域
本发明涉及预测浅层地下水位波动技术领域；具体讲，涉及考虑Lisse现象的浅层地下水 位预测方法。
技术背景
地下水作为水资源的重要组成部分，以水质较好、分布广泛以及便于开发等优点，成为 理想的供水水源之一，所以采用观测井对地下水位进行实时监测也越来越重要。然而对于降 雨条件下，地下水位变化较大较快的情形，观测井的监测数据可能存在偏差。因此，采用数 学模型模拟地下水位，更能精细地反映出水位的波动情况，且避免了不必要的人力和物力。
在传统意义上，地下水位与土壤中饱和区的顶部相一致，然而近年来，随着对饱和-非饱 和渗流研究的不断深入，越来越多的研究发现，地下水位以上存在饱和区，地下水位以下存 在非饱和区，因此，根据饱和度来定义地下水位是不合理的。本文采用孔隙水压力等于零的 线(面)，即浸润线，来代表地下水位，这也是土壤中水分流入观测井的必要条件。
Lisse现象即为，在强降雨条件下，观测井中水位快速上升，上升高度为湿润峰深度的若 干倍，且水位衰退须持续几小时至几天不等，而此时地下水并未获得真正的补给。其内在原 因为，雨水封闭了土体表面，湿润峰以下的气体受到压缩，导致非饱和区中孔隙气压力增加， 由此井中水位受到压迫而上升。因此，不能准确地识别出Lisse现象，可能会过量估计地下水 量，进而影响地下水的开采工作。
发明内容
为克服现有技术的不足，模拟分析降雨条件下，考虑Lisse现象的地下水位的波动情况。 为此，本发明采取的技术方案是，考虑Lisse现象的浅层地下水位预测方法，包括如下步骤：
步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一 直处于恒温状态，具体如下：
模型的基本控制方程为：
∂ M κ ∂ t = - div F κ + q κ - - - ( 1 ) ]]>
式中，Mκ表示包括空气a和水w的κ组分的累积质量密度；Fκ为κ组分的平流流量； qκ为κ组分的源汇项，t表示时间；
(1)累积质量密度Mκ表示β相中各组分κ的质量之和，其表达式为：
式中，β表示液相l或气相g，其中液相包括液态孔隙水及溶解的空气，气相包括干燥空 气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；Sβ为β相的饱和度；ρβ为β相 的密度；
(2)平流流量Fκ的表达式为：
F κ = Σ β X β κ F β - - - ( 3 ) ]]>
式中，Fβ为β相的平流流量，遵循达西定律，其表达式为：
F β = - K ρ β k rβ ( S l ) μ β ( ▿ p β - ρ β g ) - - - ( 4 ) ]]>
其中，K为固有渗透率；krβ为β相的相对渗透率，与液相饱和度Sl有关；μβ为β相的 粘滞性系数；pβ为β相的压力；为β相的压力梯度；g为重力加速度矢量；
毛细压力pc为液相压力pl与气相压力pg之间的差值：
pc＝pl-pg (5)
模型的本构关系包括毛细压力-饱和度关系与相对渗透率-饱和度关系，其中关于毛细压 力-饱和度关系曲线采用van Genuchten模型简称VG模型表示：
pc＝-p0[(S*)-1/λ-1]1-λ(-pmax≤pc≤0) (6)
式中，p0为进气值；λ为模型拟合参数，与土体均匀程度有关；S*为有效饱和度，表示 为S*＝(Sl-Slr)/(Sls-Slr)，Sl为液相饱和度，Slr为残余液相饱和度，Sls为饱和液相饱和度；
关于相对渗透率-饱和度的关系，其中液相相对渗透率krl采用van Genuchten-Mualem模 型简称VG-M模型表征：
k rl = S * [ 1 - ( 1 - ( S * ) 1 / λ ) λ ] 2 ( S l < S ls ) 1 ( S l ≥ S ls ) - - - ( 7 ) ]]>
气相相对渗透率krg采用科里(Corey)提出的表达式：
k rg = 1 - k rl ( S gr = 0 ) ( 1 - S ^ ) 2 ( 1 - S ^ 2 ) ( S gr > 0 ) - - - ( 8 ) ]]>
式中，Sgr为残余气相饱和度；
步骤二：模型求解：以TOUGH2/EOS3为工具，将变量分为主要变量和次要变量，其中， 次要变量可通过主要变量求得。空间上采用积分形式的有限差分方法(IFDM)进行离散，时 间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉斐逊 (Newton-Raphson)迭代方法，最后得到大型稀疏系数矩阵的线性方程组；
步骤三：模型边界条件确定：边界条件包括狄利克雷(Dirichlet)边界条件和黎曼 (Neumann)边界条件，其数学处理方法如下：
(1)Dirichlet边界条件
将边界条件单元的体积设为1×1050m3。包括空气边界和已知水头边界：
①对于空气边界，仅有气相，其主要变量为pg，和T，其中pg为气相压力，为 空气占气相的质量百分数，取值为0.9999，T为系统温度；
②对于已知水头边界，仅有液相状态，其主要变量为pl，和T，其中pl为液相压力， 为液相中空气所占的质量百分数，取值为1.0×10-10，T为系统温度；
(2)Neumann边界条件
Neumann边界条件描述的是系统与外界的流量交换情况，流入为正，流出为负；
①对于降雨入渗条件，可通过mr(t)＝ρwAeQw(t)，设置一个适当大小的源汇项来实现， 其中，源汇项mr表示单位时间内通过土体表面法向的水分质量，流入为正，流出为负；Ae为 土体的有效面积；Qw为降雨强度，ρw：水的密度；
②对于不透水边界或不入流边界，看作一类特殊的Neumann边界条件，边界上的流量设 为零；
步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件， 模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：
(1)渗流场变化：降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力 (pg-patm)/ρwg增加；接近地表处的孔隙水压力(pl-patm)/ρwg大于零，形成暂态饱和区； 在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，浸润线即pl＝0的线或 面上升，降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰下 方的孔隙水压力也恢复到初始值，当湿润峰到达初始地下水位时，初始地下水位上方形成暂 态饱和区，浸润线再次上升；
(2)地下水位变化：降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结 束时达到最大值，而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse 现象引起的；当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下 水获得真正的补给，因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下 水开采工作中，没能正确地识别Lisse现象，可能会造成地下水的过量开采。
步骤二进一步具体为：
(1)空间上采用积分有限差分法(IFDM)进行离散
首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的 积分形式进行空间离散，对于任意单元n，单元体积为Vn，边界面为Γn，单元的质量平衡方 程的积分形式如下：
d dt ∫ V n M κ d V n = ∫ Γ n F κ · nd Γ n + ∫ V n q κ d V n - - - ( 9 ) ]]>
式中，n为表面单元dΓn的单位法向量，方向为指向控制单元体内为正；
引入体积平均值：
∫ V n M κ d V n = V n M n κ - - - ( 10 ) ]]>
∫ V n q κ d V n = q n κ V n - - - ( 11 ) ]]>
式中，为Mκ，qκ在单元体积Vn上的平均值；
表面单元Γn上的面积分可近似为其所包含的各个表面Anm上面积分的平均值之和：
∫ Γ n F κ · nd Γ n = Σ m A nm F nm κ - - - ( 12 ) ]]>
式中，m为与单元n相邻的所有单元；Anm为单元n和m相邻的交界面；为Fκ在面Anm上内法线方向的平均值；
将式(10)、(11)和(12)代入到式(9)中，得到一组关于时间的一阶微分方程组
d M n w dt = 1 V n Σ m A nm F nm κ + q n κ - - - ( 13 ) ]]>
(2)时间上采用一阶向后差分方法进行离散
采用一阶向后差分方法，对式(13)进行时间上的离散，得到任一单元的全隐式非线性 方程组：
R n κ , k + 1 = M n κ , k + 1 - M n κ , k - Δt V n { Σ m A nm F nm κ , k + 1 + V n q n κ , k + 1 } - - - ( 14 ) ]]>
式中，引入了组分κ的余量上标k和k+1分别表示两相邻的时间步长指标；△t为 时间步长；
(3)Newton-Raphson迭代方法
对式(14)运用Newton-Raphson迭代方法进行线性化，引入迭代指标p，对式中的余量 在迭代步p+1处进行泰勒级数展开，只保留一阶项，推导出包含2×NEL(计算域 内单元数)个方程的线性方程组，并且以两迭代步的增量xi,p+1-xi,p为未知量；最后得到大型 稀疏系数矩阵的线性方程组：
- Σ i ∂ R n κ , k + 1 ∂ x i | p ( x i , p + 1 - x i , p ) = R n κ , k + 1 ( x i , p ) - - - ( 15 ) . ]]>
与已有技术相比，本发明的技术特点与效果：
(1)建立非饱和土的水-气二相流模型，模拟了降雨条件下，考虑Lisse现象的浅层地 下水位的波动情况。
(2)本发明避免了野外监测的高费用以及结果偏差问题，模拟结果合理可信，避免了由 于Lisse现象造成的地下水过量开采的后果。
附图说明
图1土体几何形状。
图2(a)毛细压力-饱和度关系曲线；(b)相对渗透率-饱和度关系曲线。
图3降雨过程中不同时刻的土体剖面(a)液相饱和度；(b)毛细压力；
(c)孔隙气压力；(d)孔隙水压力的分布情况。
图4降雨结束后不同时刻的土体剖面(a)液相饱和度；(b)毛细压力；
(c)孔隙气压力；(d)孔隙水压力的分布情况。
图5(a)B点压力水头和气体流速；(b)C点压力水头和液相饱和度；(c)D点压力水 头随时间的变化情况。
图6Lisse现象原理图。
图7地下水位随时间的变化情况
图8本发明的模拟方法流程图
具体实施方式
本发明旨在基于非饱和土中多孔介质的多相流理论，在同时考虑水相和气相流动规律及 其相互作用的基础上，建立了非饱和土的水-气二相流模型，来模拟分析降雨条件下，考虑 Lisse现象的地下水位的波动情况。
本发明采用的技术方案是：
步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一 直处于恒温状态。具体如下：
模型的基本控制方程为：
∂ M κ ∂ t = - div F κ + q κ - - - ( 1 ) ]]>
式中，Mκ表示κ组分(空气a和水w)的累积质量密度；Fκ为κ组分的平流流量；qκ为 κ组分的源汇项。
(1)累积质量密度Mκ表示β相中各组分κ的质量之和，其表达式为：
式中，β表示液相(l)或气相(g)，其中液相包括液态孔隙水及溶解的空气，气相包括 干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；Sβ为β相的饱和度；ρβ为β相的密度。
(2)平流流量Fκ的表达式为：
F κ = Σ β X β κ F β - - - ( 3 ) ]]>
式中，Fβ为β相的平流流量，遵循达西定律，其表达式为：
F β = - K ρ β k rβ ( S l ) μ β ( ▿ p β - ρ β g ) - - - ( 4 ) ]]>
其中，K为固有渗透率；krβ为β相的相对渗透率，与液相饱和度Sl有关；μβ为β相的 粘滞性系数；pβ为β相的压力；为β相的压力梯度；g为重力加速度矢量。
毛细压力pc为液相压力pl与气相压力pg之间的差值：
pc＝pl-pg (5)
模型的本构关系包括毛细压力-饱和度关系与相对渗透率-饱和度关系。其中关于毛细压 力-饱和度关系曲线采用van Genuchten模型(VG模型)表示：
pc＝-p0[(S*)-1/λ-1]1-λ(-pmax≤pc≤0) (6)
式中，p0为进气值；λ为模型拟合参数，与土体均匀程度有关；S*为有效饱和度，可表 示为S*＝(Sl-Slr)/(Sls-Slr)，Sl为液相饱和度，Slr为残余液相饱和度，Sls为饱和液相饱和度。
关于相对渗透率-饱和度的关系，其中液相相对渗透率krl采用van Genuchten-Mualem模 型(VG-M模型)表征：
k rl = S * [ 1 - ( 1 - ( S * ) 1 / λ ) λ ] 2 ( S l < S ls ) 1 ( S l ≥ S ls ) - - - ( 7 ) ]]>
气相相对渗透率krg采用Corey(科里)提出的表达式：
k rg = 1 - k rl ( S gr = 0 ) ( 1 - S ^ ) 2 ( 1 - S ^ 2 ) ( S gr > 0 ) - - - ( 8 ) ]]>
式中，Sgr为残余气相饱和度。
步骤二：模型求解：以TOUGH2/EOS3为工具，将变量分为主要变量和次要变量，其中， 次要变量可通过主要变量求得。空间上采用积分形式的有限差分方法(IFDM)进行离散，时 间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用Newton-Raphson(牛顿- 拉斐逊)迭代方法，最后得到大型稀疏系数矩阵的线性方程组；具体如下：
(1)空间上采用积分有限差分法(IFDM)进行离散
首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的 积分形式进行空间离散。对于任意单元n，单元体积为Vn，边界面为Γn，单元的质量平衡方 程的积分形式如下：
d dt ∫ V n M κ d V n = ∫ Γ n F κ · nd Γ n + ∫ V n q κ d V n - - - ( 9 ) ]]>
式中，n为表面单元dΓn的单位法向量，方向为指向控制单元体内为正。
引入体积平均值：
∫ V n M κ d V n = V n M n κ - - - ( 10 ) ]]>
∫ V n q κ d V n = q n κ V n - - - ( 11 ) ]]>
式中，为Mκ，qκ在单元体积Vn上的平均值。
表面单元Γn上的面积分可近似为其所包含的各个表面Anm上面积分的平均值之和：
∫ Γ n F κ · nd Γ n = Σ m A nm F nm κ - - - ( 12 ) ]]>
式中，m为与单元n相邻的所有单元；Anm为单元n和m相邻的交界面；为Fκ在面Anm上内法线方向的平均值。
将式(10)、(11)和(12)代入到式(9)中，得到一组关于时间的一阶微分方程组
d M n w dt = 1 V n Σ m A nm F nm κ + q n κ - - - ( 13 ) ]]>
(2)时间上采用一阶向后差分方法进行离散
采用一阶向后差分方法，对式(13)进行时间上的离散，得到任一单元的全隐式非线性 方程组：
R n κ , k + 1 = M n κ , k + 1 - M n κ , k - Δt V n { Σ m A nm F nm κ , k + 1 + V n q n κ , k + 1 } - - - ( 14 ) ]]>
式中，引入了组分κ的余量上标k和k+1分别表示两相邻的时间步长指标；△t为 时间步长；式子右端的平流流量项和源汇项均采用新的时间步长值，这种全隐式方法提高了 模型求解的数值稳定性。
(3)Newton-Raphson迭代方法
对式(14)运用Newton-Raphson迭代方法进行线性化。引入迭代指标p，对式中的余量 在迭代步p+1处进行泰勒级数展开，只保留一阶项，推导出包含2×NEL(计算域 内单元数)个方程的线性方程组，并且以两迭代步的增量xi,p+1-xi,p为未知量；最后得到大型 稀疏系数矩阵的线性方程组：
- Σ i ∂ R n κ , k + 1 ∂ x i | p ( x i , p + 1 - x i , p ) = R n κ , k + 1 ( x i , p ) - - - ( 15 ) ]]>
步骤三：模型边界条件确定：边界条件包括Dirichlet(狄利克雷)边界条件和Neumann (黎曼)边界条件，其数学处理方法如下：
(1)Dirichlet边界条件
在Dirichlet边界条件的数学处理方法中，一般将边界条件单元的体积看作非常大，即认 为边界条件单元与相邻土体单元的流量交换，不影响边界条件单元的状态，因此可将边界条 件单元的体积设为1×1050m3。包括空气边界和已知水头边界：
①对于空气边界，仅有气相，其主要变量为pg，和T，其中pg为气相压力，为 空气占气相的质量百分数，取值为0.9999，T为系统温度；
②对于已知水头边界，仅有液相状态，其主要变量为pl，和T，其中pl为液相压力， 为液相中空气所占的质量百分数，取值为1.0×10-10，T为系统温度。
(2)Neumann边界条件
Neumann边界条件描述的是系统与外界的流量交换情况，流入为正，流出为负，可以是 常量，也可以随时间变化。
①对于降雨入渗条件，可通过mr(t)＝ρwAeQw(t)，设置一个适当大小的源汇项来实现， 其中，源汇项mr表示单位时间内通过土体表面法向的水分质量，流入为正，流出为负；Ae为 土体的有效面积；Qw为降雨强度。
②对于不透水边界或不入流边界，可看作一类特殊的Neumann边界条件，边界上的流量 设为零。
步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件， 模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：
(1)渗流场变化：降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力 (pg-patm)/ρwg增加；接近地表处的孔隙水压力(pl-patm)/ρwg大于零，形成暂态饱和区； 在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，浸润线(pl＝0的线或 面)上升。降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰 下方的孔隙水压力也恢复到初始值。当湿润峰到达初始地下水位时，初始地下水位上方形成 暂态饱和区，浸润线再次上升。
(2)地下水位变化：降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结 束时达到最大值。而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse 现象引起的。当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下 水获得真正的补给。因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下 水开采工作中，没能正确地识别Lisse现象，可能会造成地下水的过量开采。
下面结合附图，针对模拟降雨条件下，考虑Lisse现象的浅层地下水位的波动情况进行具 体说明，其步骤如下：
(1)建立数学模型：模型选取均质各向同性的土体作为研究对象，模拟范围长100m， 宽100m，厚度为2m，其剖面图见图1。地下水位位于地表以下0.6m，观测点B、C、D分别 位于地表以下0.35m，0.58m和0.60m。土体材料特性参数取值见表1。毛细压力-饱和度、相 对渗透率-饱和度关系曲线分别采用VG模型和VG-M模型，如图2所示。假定含水层系统处 于恒温状态，T＝25℃。
表1土体材料特性参数取值
模型的基本控制方程为：
∂ M κ ∂ t = - div F κ + q κ - - - ( 1 ) ]]>
式中，Mκ表示κ组分(空气a和水w)的累积质量密度；Fκ为κ组分的平流流量；qκ为 κ组分的源汇项。
(1)累积质量密度Mκ表示β相中各组分κ的质量之和，其表达式为：
式中，β表示液相(l)或气相(g)，其中液相包括液态孔隙水及溶解的空气，气相包括 干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；Sβ为β相的饱和度；ρβ为β相的密度。
(2)平流流量Fκ的表达式为：
F κ = Σ β X β κ F β - - - ( 3 ) ]]>
式中，Fβ为β相的平流流量，遵循达西定律，其表达式为：
F β = - K ρ β k rβ ( S l ) μ β ( ▿ p β - ρ β g ) - - - ( 4 ) ]]>
其中，K为固有渗透率；krβ为β相的相对渗透率，与液相饱和度Sl有关；μβ为β相的 粘滞性系数；pβ为β相的压力；为β相的压力梯度；g为重力加速度矢量。
毛细压力pc为液相压力pl与气相压力pg之间的差值：
pc＝pl-pg (5)
模型的本构关系包括毛细压力-饱和度关系与相对渗透率饱和度关系。其中关于毛细压力 -饱和度关系曲线采用van Genuchten模型(VG模型)表示：
pc＝-p0[(S*)-1/λ-1]1-λ(-pmax≤pc≤0) (6)
式中，p0为进气值；λ为模型拟合参数，与土体均匀程度有关；S*为有效饱和度，可表 示为S*＝(Sl-Slr)/(Sls-Slr)，Sl为液相饱和度，Slr为残余液相饱和度，Sls为饱和液相饱和度。
关于相对渗透率-饱和度的关系，其中液相相对渗透率krl采用van Genuchten-Mualem模 型(VG-M模型)表征：
k rl = S * [ 1 - ( 1 - ( S * ) 1 / λ ) λ ] 2 ( S l < S ls ) 1 ( S l ≥ S ls ) - - - ( 7 ) ]]>
气相相对渗透率krg采用Corey(科里)提出的表达式：
k rg = 1 - k rl ( S gr = 0 ) ( 1 - S ^ ) 2 ( 1 - S ^ 2 ) ( S gr > 0 ) - - - ( 8 ) ]]>
式中，Sgr为残余气相饱和度。
(2)模型求解：为了保证模型计算结果的精确度，竖直方向上，模型网格尺寸在饱和区 为0.05m，非饱和区为0.01m；水平方向上，网格尺寸为1m。以TOUGH2/EOS3为工具，空 间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方 法进行离散，见式(9)；模型的线性化采用Newton-Raphson(牛顿-拉斐逊)迭代方法，最后 得到大型稀疏系数矩阵的线性方程组，见式(10)：
R n κ , k + 1 = M n κ , k + 1 - M n κ , k - Δt V n { Σ m A nm F nm κ , k + 1 + V n q n κ , k + 1 } - - - ( 9 ) ]]>
式中，引入了组分κ的余量分别为Mκ，qκ在单元体积Vn上的平均值； m为与单元n相邻的所有单元；Anm为单元n和m相邻的交界面；为Fκ在面Anm上内法线 方向的平均值；上标k和k+1分别表示两相邻的时间步长指标；△t为时间步长。
引入迭代指标p，对式(9)中的余量在迭代步p+1处进行泰勒级数展开，只保 留一阶项，推导出包含2×NEL(计算域内单元数)个方程的线性方程组，并且以两迭代步 的增量xi,p+1-xi,p为未知量；最后得到大型稀疏系数矩阵的线性方程组：
- Σ i ∂ R n κ , k + 1 ∂ x i | p ( x i , p + 1 - x i , p ) = R n κ , k + 1 ( x i , p ) - - - ( 10 ) ]]>
(3)模型边界条件确定：首先计算稳定渗流情况，作为降雨入渗的初始条件。在稳定渗 流情况下，其边界条件为：地下水位以下的侧面边界为已知水头边界，即 pl＝patm+ρwg(1.8-Z)，T＝25℃；地表的边界条件为空气边界，即 pg＝patm，T＝25℃；区域底部及其它边界为不透水边界。模型的初始条件设 为土体水相饱和，主要变量pl＝1.013×105pa，T＝25℃，在此定解条件下， 运行模型直至达到毛细压力-重力平衡状态，即为降雨模型的初始条件。
以稳定渗流状态作为降雨入渗模拟的初始条件，研究区域的侧面与底部为不透水边界， 地表为空气边界。区域顶部根据源汇项公式mr(t)＝ρwAeQw(t)，施加45mm/h的雨强，持续 时间为1h。
(4)根据非饱和土的水-气二相流模型，分析降雨条件下，非饱和区中的渗流场及考虑 Lisse现象的浅层地下水位的波动情况：
由于研究土体为砂土，因此其毛细管厚度很小，可以忽略不计。暂态饱和区是一种由非 饱和状态过渡到完全饱和状态的区域，此时毛细压力仍然存在；孔隙水压力((pl-patm)/ρwg) 和孔隙气压力((pg-patm)/ρwg)均为正值，且孔隙水压力小于孔隙气压力。图3为降雨过 程中不同时刻的液相饱和度、毛细压力和孔隙压力的分布情况。降雨开始后，湿润峰下移， 并在降雨结束时达到0.19m。毛细压力的变化趋势与饱和度变化相似，在湿润区负毛细压力 减小，湿润区下方仍保持初始值。由于入渗水流下移，湿润峰下方的空气受到压缩，因此非 饱和区的孔隙气压力增加，并在降雨结束时达到最大值，0.13m。接近地表处的孔隙水压力显 著增加，且大于零，说明此处形成了暂态饱和区；在湿润峰下方，由于孔隙气压力的增加， 孔隙水压力偏离了初始值，且偏离值的大小等于相应的孔隙气压力的变化。
图4为降雨结束后不同时刻的液相饱和度、毛细压力和孔隙压力的分布情况。降雨结束 之后，湿润峰继续下移，在6h至10h之间到达初始地下水位。在湿润锋到达初始地下水位之 前，由于非饱和区中的气体可以溢出地表，因此初始地下水位以上的孔隙气压力减小为零， 且湿润区下方的孔隙水压力恢复至初始值；而由于蒸发作用，接近地表处的孔隙水压力小于 零，暂态饱和区恢复到初始的非饱和状态。当湿润峰到达初始地下水位时，其上方的孔隙气 压力开始逐渐增加；孔隙水压力开始偏离初始值；同时由于封闭气泡，毛细压力仍然存在， 即在地下水位上方，形成暂态饱和区。
图5(a)为观测点B处孔隙压力、毛细压力和气体流速随时间的变化情况。在降雨过程 中，孔隙气压力一开始迅速增大，之后由于地下水的阻碍缓慢增加，并在降雨结束时达到最 大值0.13m。孔隙水压力的变化趋势与孔隙气压力基本一致。由于湿润峰未到达B点，因此 其毛细压力和饱和度均保持不变。降雨一结束，孔隙气压力和孔隙水压力均迅速减小，分别 接近于零和初始值。之后孔隙水压力与毛细压力的变化趋势一致。在降雨开始后的4小时， B点的负毛细压力开始减小，说明此时湿润峰到达该处，饱和度增加。由B点竖直方向的气 体流速随时间的变化情况来看，其中气体向上为正，向下为负，降雨一开始，由于雨水下渗， B点气体受到压缩，向下进入土体内部，导致孔隙气压力迅速增加；之后由于受到地下水的 限制，气体流速快速减小，因此孔隙气压力缓慢增加；当降雨结束时，气流反转，气体溢出 地表，造成B点孔隙气压力减小为零；当湿润峰到达B点时，气流向上，这是因为气体受到 水流驱替。在整个入渗过程中，B点的孔隙水压力均为负值，说明B点一直处于非饱和状态。
图5(b)观测点C处孔隙压力、毛细压力和液相饱和度随时间的变化情况。降雨过程中， C点孔隙气压力与B点变化情况相似。其孔隙水压力由初始的负值逐渐增加，并在降雨结束 时达到最大值0.106m，说明C点水分获得气流提供的能量，克服了毛管力，能够自由流动并 传递水压力。而C点的饱和度和毛细压力在降雨过程中基本保持不变。当降雨结束时，C点 孔隙气压力和孔隙水压力均迅速减小，分别接近于零和初始值。在7h时，C点饱和度增加， 说明湿润锋到达该点，其孔隙水压力和孔隙气压力也开始逐渐增加，但孔隙水压力始终小于 孔隙气压力，毛细压力存在。在18h时，孔隙水压力变为正值，浸润线上升。当湿润锋到达 C点之后，且其孔隙水压力大于零时，C点处于暂态饱和状态。
图5(c)观测点D处孔隙压力和毛细压力随时间的变化情况。由于D点位于初始地下水 位处，因此其初始的孔隙气压力和孔隙水压力均为零，且在整个入渗过程中，毛细压力保持 为零，孔隙水压力等于孔隙气压力。降雨过程中，D点孔隙气压力和孔隙水压力增加，且在 降雨结束时达到最大值，0.13m。降雨结束之后，孔隙气压力和孔隙水压力均恢复到零。在 7h时，由于湿润锋到达该点，两者开始逐渐增加。在整个入渗过程中，D点始终处于饱和状 态，孔隙水压力的增加有两点原因：降雨过程中的气流推动和获得真正补给。
图6为Lisse现象的原理解释，地下水位以下的土体和观测井可看作连通器。观测点D 位于土体中的地下水位处，点E位于井中，且与D点相平。根据连通器原理，D点的液相压 力等于E点的液相压力即：
p l D = p l E - - - ( 11 ) ]]>
在降雨开始之前，井中的水位应与地下水位相平，D点和E点的液相压力为：
p l D = p l E = p atm - - - ( 12 ) ]]>
在降雨过程中，由于初始地下水位和湿润峰之间的空气受到压缩，因此D点的气相压力 增加(如图5c)。对于任一时刻，D点增加的气相压力为水分提供了额外的能量，导致浸 润线上升。此时，D点的水相压力等于其气相压力：
p l D = p g E = p atm + Δ p g D - - - ( 13 ) ]]>
根据连通器原理，土体中的孔隙水将会流入井中，直至井中水位与浸润线相平。此时， 假设井中水位上升高度为△H，则E点的水相压力为：
p l E = p atm + ρ w gΔH - - - ( 14 ) ]]>
由方程(11)，水位上升高度△H可表示为：
ΔH = Δ p g D / ρ w g - - - ( 15 ) ]]>
此时，湿润峰并未到达初始地下水位，水位上升则表示发生了Lisse现象。而且，Lisse 现象是由地下水位处的孔隙气压力控制的。当湿润峰到达初始地下水位时，地下水获得真正 的补给，浸润线再次上升。
图7为浸润线上升高度随时间的变化情况，其中初始地下水位为基准线。浸润线的变化 趋势和初始地下水位处的孔隙气压力变化趋势一致：降雨过程中，首先快速增加，之后缓慢 增加，并在降雨结束时达到最大值0.13m。而此时湿润峰并未到达初始地下水位处，说明水 位的快速上升是由Lisse现象引起的。当降雨结束后，浸润线快速回到初始地下水位，之后在 7h处逐渐上升，说明此时地下水获得真正的补给。因此，当用浸润线代表地下水位时，水位 波动包括两部分原因：Lisse现象和真正补给。
综上所述，对本实施例总结如下：
(1)针对降雨条件下，考虑Lisse现象的浅层地下水位的波动，采用孔隙水压力等于零 的线(面)来代表地下水位，并利用数值模拟的方法研究了水位的变化情况。首先利用非饱 和土的水-气二相流模型，模拟了稳定渗流情况，作为降雨入渗的初始条件，之后分析了降雨 条件下，非饱和土中渗流场及水位的变化情况。
(2)对非饱和区中渗流场的分析表明，降雨过程中，非饱和区的气体受到湿润峰下移的 压缩，孔隙气压力增加；接近地表处的孔隙水压力大于零，形成暂态饱和区；在湿润峰下方， 由于孔隙气压力的增加，孔隙水压力偏离了初始值。降雨结束之后，非饱和区中的气体可以 自由溢出，孔隙气压力减小为零，湿润峰下方的孔隙水压力也恢复到初始值。当湿润峰到达 初始地下水位时，初始地下水位上方形成暂态饱和区，浸润线上升。
(3)降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结束时达到最大值。 而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的。当降雨 结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下水获得真正的补给。 因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下水开采工作中，没能 正确地识别Lisse现象，可能造成地下水过量开采的后果。然而，由于土体材料特性、降雨强 度和地形特点的不同，Lisse现象的形式多种多样，水位回退也不尽相同，因此需要相应开展 很多的研究工作。