[发明专利]考虑Lisse现象的浅层地下水位预测方法

权利要求书

1.一种考虑Lisse现象的浅层地下水位预测方法，其特征是，包括如下步骤：

步骤一：建立非饱和土的水-气二相流模型，包括基本控制方程及本构关系，假定系统一直处于恒温状态，具体如下：

模型的基本控制方程为：

∂Mκ∂t=-divFκ+qκ---(1)]]>

式中，M^κ表示包括空气a和水w的κ组分的累积质量密度；F^κ为κ组分的平流流量；q^κ为κ组分的源汇项，t表示时间；

(1)累积质量密度M^κ表示β相中各组分κ的质量之和，其表达式为：

式中，β表示液相l或气相g，其中液相包括液态孔隙水及溶解的空气，气相包括干燥空气和水蒸气；为κ组分占β相的质量百分数；φ为孔隙率；S_β为β相的饱和度；ρ_β为β相的密度；

(2)平流流量F^κ的表达式为：

Fκ=ΣβXβκFβ---(3)]]>

式中，F_β为β相的平流流量，遵循达西定律，其表达式为：

Fβ=-Kρβkrβ(Sl)μβ(▿pβ-ρβg)---(4)]]>

其中，K为固有渗透率；k_rβ为β相的相对渗透率，与液相饱和度S_l有关；μ_β为β相的粘滞性系数；p_β为β相的压力；为β相的压力梯度；g为重力加速度矢量；

毛细压力p_c为液相压力p_l与气相压力p_g之间的差值：

p_c＝p_l-p_g    (5)

模型的本构关系包括毛细压力-饱和度关系与相对渗透率-饱和度关系，其中关于毛细压力-饱和度关系曲线采用van Genuchten模型简称VG模型表示：

pc=-p0[(S*)-1/λ-1]1-λ,(-pmax≤pc≤0)---(6)]]>

式中，p₀为进气值；λ为模型拟合参数，与土体均匀程度有关；S^*为有效饱和度，表示为S^*＝(S_l-S_lr)/(S_ls-S_lr)，S_l为液相饱和度，S_lr为残余液相饱和度，S_ls为饱和液相饱和度；

关于相对渗透率-饱和度的关系，其中液相相对渗透率k_rl采用van Genuchten-Mualem模型简称VG-M模型表征：

krl=S*[1-(1-(S*)1/λ)λ]2(Sl<Sls)1(Sl≥Sls)---(7)]]>

气相相对渗透率k_rg采用科里(Corey)提出的表达式：

krg=1-krl(Sgr=0)(1-S^)2(1-S^2)(Sgr>0)---(8)]]>

式中，S^=(Sl-Slr)/(1-Slr-Sgr),]]>S_gr为残余气相饱和度；

步骤二：模型求解：以TOUGH2/EOS3为工具，将变量分为主要变量和次要变量，其中，次要变量可通过主要变量求得，空间上采用积分形式的有限差分方法(IFDM)进行离散，时间上采用一阶向后差分的全隐式方法进行离散，模型的线性化采用牛顿-拉斐逊(Newton-Raphson)迭代方法，最后得到大型稀疏系数矩阵的线性方程组；

步骤三：模型边界条件确定：边界条件包括狄利克雷(Dirichlet)边界条件和黎曼(Neumann)边界条件，其数学处理方法如下：

(1)Dirichlet边界条件

将边界条件单元的体积设为1×10⁵⁰m³，包括空气边界和已知水头边界：

①对于空气边界，仅有气相，其主要变量为p_g，和T，其中p_g为气相压力，为空气占气相的质量百分数，T为系统温度；

②对于已知水头边界，仅有液相状态，其主要变量为p_l，和T，其中p_l为液相压力，为液相中空气所占的质量百分数，T为系统温度；

(2)Neumann边界条件

Neumann边界条件描述的是系统与外界的流量交换情况，流入为正，流出为负；

①对于降雨入渗条件，通过m_r(t)＝ρ_wA_eQ_w(t)，设置一个适当大小的源汇项来实现，其中，源汇项m_r表示单位时间内通过土体表面法向的水分质量，流入为正，流出为负；A_e为土体的有效面积；Q_w为降雨强度，ρ_w为水的密度；

②对于不透水边界或不入流边界，看作一类特殊的Neumann边界条件，边界上的流量设为零；

步骤四：利用非饱和土的水-气二相流模型，以稳定渗流情况作为降雨入渗的初始条件，模拟分析降雨条件下，考虑Lisse现象的浅层地下水位的波动情况：

(1)渗流场变化：降雨过程中，非饱和区的气体受到湿润峰下移的压缩，孔隙气压力(p_g-p_atm)/ρ_wg增加；接近地表处的孔隙水压力(p_l-p_atm)/ρ_wg大于零，形成暂态饱和区；在湿润峰下方，由于孔隙气压力的增加，孔隙水压力偏离了初始值，浸润线(p_l＝0的线或面)上升，降雨结束之后，非饱和区中的气体可以自由溢出，孔隙气压力减小为零，湿润峰下方的孔隙水压力也恢复到初始值，当湿润峰到达初始地下水位时，初始地下水位上方形成暂态饱和区，浸润线再次上升；

(2)地下水位变化：降雨过程中，地下水位首先快速增加，之后缓慢增加，并在降雨结束时达到最大值，而此时湿润峰并未到达初始地下水位处，说明水位的快速上升是由Lisse现象引起的；当降雨结束后，浸润线快速回到初始地下水位，之后逐渐上升，说明此时地下水获得真正的补给，因此，地下水位的波动有两部分原因：Lisse现象和真正的补给，在地下水开采工作中，没能正确地识别Lisse现象，可能会造成地下水的过量开采。

2.如权利要求1所述的考虑Lisse现象的浅层地下水位预测方法，其特征是，步骤二进一步具体为：

步骤二进一步具体为：

(1)空间上采用积分有限差分法(IFDM)进行离散

首先将计算域离散成子单元，其性质由形心点代表，分别对各个单元的质量平衡方程的积分形式进行空间离散，对于任意单元n，单元体积为V_n，边界面为Γ_n，单元的质量平衡方程的积分形式如下：

ddt∫VnMκdVn=∫ΓnFκ·ndΓn+∫VnqκdVn---(9)]]>

式中，n为表面单元dΓ_n的单位法向量，方向为指向控制单元体内为正；

引入体积平均值：

∫VnMκdVn=VnMnκ---(10)]]>

∫VnqκdVn=qnκVn---(11)]]>

式中，M^κ，q^κ在单元体积V_n上的平均值；

表面单元Γ_n上的面积分可近似为其所包含的各个表面A_nm上面积分的平均值之和：

∫ΓnFκ·ndΓn=ΣmAnmFnmκ---(12)]]>

式中，m为与单元n相邻的所有单元；A_nm为单元n和m相邻的交界面；为F^κ在面A_nm上内法线方向的平均值；

将式(10)、(11)和(12)代入到式(9)中，得到一组关于时间的一阶微分方程组

dMnwdt=1VnΣmAnmFnmκ+qnκ---(13)]]>

(2)时间上采用一阶向后差分方法进行离散

采用一阶向后差分方法，对式(13)进行时间上的离散，得到任一单元的全隐式非线性方程组：

Rnκ,k+1=Mnκ,k+1-Mnκ,k-ΔtVn{ΣmAnmFnmκ,k+1+Vnqnκ,k+1}---(14)]]>

式中，引入了组分κ的余量上标k和k+1分别表示两相邻的时间步长指标；△t为时间步长；

(3)Newton-Raphson迭代方法

对式(14)运用Newton-Raphson迭代方法进行线性化，引入迭代指标p，对式中的余量在迭代步p+1处进行泰勒级数展开，只保留一阶项，推导出包含2×NEL(计算域内单元数)个方程的线性方程组，并且以两迭代步的增量x_i,p+1-x_i,p为未知量；最后得到大型稀疏系数矩阵的线性方程组：

-Σi∂Rnκ.k+1∂xi|p(xi,p+1-xi,p)=Rnκ,k+1(xi,p)---(15).]]>