本發(fā)明所屬技術(shù)領(lǐng)域為“海水入侵內(nèi)陸邊界控制方法”領(lǐng)域。具體講,涉及模擬海平面上升情況下海水入侵內(nèi)陸邊界的處理方法。

背景技術(shù)：

世界上大部分的人口分布在沿海地區(qū)，近年來,受地殼運(yùn)動、自然壓密、氣候變暖引起的海平面上升、大量抽用地下水等自然及人為因素的影響，致使沿海各地區(qū)已不同程度地遭到海水入侵,嚴(yán)重制約沿海地區(qū)的社會經(jīng)濟(jì)發(fā)展，引起了世界各沿海國家的充分重視。最近，由于海平面上升引起的海水入侵問題的研究也受到越來越多的關(guān)注，其中，內(nèi)陸邊界條件對模擬海平面上升引起的海水入侵有至關(guān)重要的影響。因此，開展內(nèi)陸邊界條件對模擬海平面上升引起的海水入侵影響的深入研究，尋求一種合適的內(nèi)陸邊界條件的處理方法，具有重要的科學(xué)價值和社會意義。

雖然流量控制系統(tǒng)邊界和水頭控制系統(tǒng)邊界在數(shù)值模擬中是最常應(yīng)用的兩種邊界條件，但事實上，以上兩種邊界都有一定的局限性，因為隨著海平面的上升，真實的模擬區(qū)域內(nèi)陸邊界既不是水頭恒定，也不是流量恒定，而是一個間接的過程，因此有人提出了第三種邊界條件——依靠水頭的流量邊界條件，這種邊界條件在一些特地場地也進(jìn)行了海平面上升對海水入侵的影響評估，而且也有人提出了解析解，但是卻沒有基于假定的理想的概念模型的一般性的研究，更缺乏基于理想的數(shù)值模擬在預(yù)測海平面上升對海水入侵的三種不同內(nèi)陸邊界條件的對比并得出合適的模擬海平面上升情況下海水入侵內(nèi)陸邊界的處理方法。

技術(shù)實現(xiàn)要素：

為克服現(xiàn)有技術(shù)的不足，本發(fā)明旨在在深入分析地下多相流運(yùn)移機(jī)理的基礎(chǔ)上，建立地下水氣-液二相流及溶質(zhì)運(yùn)移模型，來模擬分析不同內(nèi)陸邊界條件對沿海地區(qū)海平面上升引起的海水入侵的效果。本發(fā)明采用的技術(shù)方案是，模擬海平面上升情況下海水入侵內(nèi)陸邊界的處理方法，步驟如下：

步驟一：建立地下水氣-液二相流及溶質(zhì)運(yùn)移模型，包括基本控制方程及輔助方程，假定系統(tǒng)一直處于恒溫狀態(tài)。具體如下：

模型的基本控制方程為：

式中，M^κ表示κ組分的累積質(zhì)量密度，組分是指純水w、參考鹽水b和空氣a，F(xiàn)^κ為κ組分的流量密度，包括平流流量密度和擴(kuò)散流量密度q^κ為組分κ的源匯項；

(1)累積質(zhì)量密度等于β相中組分κ的質(zhì)量之和，其表達(dá)式為：

式中，β表示流體相，流體相為氣相或液相，液相為純水和參考鹽水的混合物，表示孔隙率，S_β為β相的飽和度，ρ_β為β相的密度，p_atm為κ組分占β相的質(zhì)量百分?jǐn)?shù)；

(2)平流流量p_g的表達(dá)式為：

式中，F(xiàn)_β為相的平流流量，符合達(dá)西定律，其多相流形式的表達(dá)式為：

其中，k為固有滲透系數(shù)，k_rβ為β相的相對滲透系數(shù)，μ_β為β相的粘滯性系數(shù)，為梯度符號，p_β為β相的孔隙壓力，g為重力加速度矢量；

(3)擴(kuò)散流量的表達(dá)式為：

式中，為κ組分在β相中的分子擴(kuò)散系數(shù)；τ₀τ_β(S_β)為迂曲度，τ₀是多孔介質(zhì)特性參數(shù)，τ_β＝τ_β(S_β)是飽和度的函數(shù)；

分子擴(kuò)散系數(shù)隨壓力p和溫度T變化的計算式為：

式中，θ為溫度系數(shù)；

非飽和地下水流及熱流傳輸TOUGH的第二個版本TOUGH2中提供了三種迂曲度模型：

τ₀τ_β＝τ₀k_rβ(S_β) 相對滲透率模型

米林頓夸克模型

τ₀τ_β＝τ₀S_β 常量擴(kuò)散 (7)

步驟二：模型求解：以TOUGH2的水、鹽水和空氣的混合模塊EOS7為工具，空間上采用積分形式的有限差分方法IFDM進(jìn)行離散，時間上采用一階向后差分的全隱式方法進(jìn)行離散，模型的線性化采用牛頓-拉弗森Newton-Raphson迭代方法，最后得到大型稀疏系數(shù)矩陣的線性方程組；

步驟三：模型邊界條件確定：模型計算中的邊界條件包括狄利克雷Dirichlet邊界條件和諾伊曼Neumann邊界條件兩種，其數(shù)學(xué)處理方法如下：

(1)Dirichlet邊界條件

Dirichlet邊界條件上，邊界條件單元的主要變量在計算過程中保持不變：

①對于空氣邊界，其邊界條件單元的相態(tài)為僅有氣相狀態(tài),變量為p_g,X_b,T，對于與大氣接觸的邊界上p_g＝p_atm，X_b＝0.0；T為氣溫；

②對于已知水頭邊界，包括地下淡水水頭邊界和海水水頭邊界，其邊界條件單元的相態(tài)均為僅有液相狀態(tài),變量為p_g，X_b，T，由于液相飽和狀態(tài)下毛細(xì)壓力為零，因此有p_g＝p_atm+ρ_lgh，其中ρ_l為淡水或海水的密度，h為邊界上的水深；地下淡水水頭邊界上X_b等于零，海水邊界上的X_b是根據(jù)參考鹽水的特性以及海水的鹽度計算出來；

(2)Neumann邊界條件

Neumann邊界條件描述的是系統(tǒng)與外界的流量交換情況，邊界條件單元的單位流量對應(yīng)式中的源匯項，流入為正：

①對于不透水邊界，是一類特殊的Neumann邊界條件，邊界上的流量為零，積分形式的有限差分法的處理很簡單，就是不設(shè)置與之相鄰的邊界條件單元，則與不透水面單元的流量交換為零；

②對于降雨入滲條件的處理方法，則設(shè)置一個適當(dāng)大小的源匯項來實現(xiàn)；

(3)依靠水頭的流量邊界條件

設(shè)區(qū)域外的已知水位又稱參考水頭為H_ref，水流流入或者流出研究區(qū)域的流量大小與研究區(qū)內(nèi)的水頭H和參考水頭H_ref之間的水頭差成比例，并且用如下線性方程式表示：

q＝-C(H-H_ref) (15)

其中q是從內(nèi)陸側(cè)流入研究區(qū)的單寬淡水補(bǔ)給量，C是導(dǎo)水系數(shù)。

而研究區(qū)域淡水側(cè)的等效淡水水頭H_f跟淡水補(bǔ)給量q_f在承壓含水層和非承壓含水層中分別有如下關(guān)系：

承壓含水層：

非承壓含水層：

其中K是含水層的滲透系數(shù)，α＝ρ_f/(ρ_s-ρf)是相對密度，ρ_f和ρ_s分別是淡水密度和海水密度，L和B分別為研究區(qū)域含水層的長度和厚度，H_s是距離含水層底部的平均海水深度。

把公式(15)分別帶入公式(16)和(17)中，并且讓H＝H_f和q＝q_f，可得出如下關(guān)系式：

承壓含水層：

非承壓含水層：

模型求解具體步驟如下：

(1)空間上采用積分有限差分法IFDM進(jìn)行離散

首先將計算域離散成子單元，其性質(zhì)由形心點代表，分別對各個單元的質(zhì)量平衡方程的積分形式進(jìn)行空間離散，對于任意單元n，單元體積為V_n，邊界面為Γ_n，單元的質(zhì)量平衡方程的積分形式如下：

式中，n為表面單元dΓ_n的單位法向量，指向控制單元體內(nèi)為正；

引入體積平均值，有

式中，為M^κ，q^κ在V_n上的平均值。

Γ_n上的面積分可近似為其所包含的各個表面A_nm的面積分的平均值之和，有

式中，m為與單元n相鄰的所有單元，A_nm是單元n和m相鄰的交界面，是F^κ在面A_nm上的內(nèi)法線方向的平均值；

將式(9)、(10)和(11)代入到式(8)中，得到一組關(guān)于時間的一階微分方程組

(2)時間上采用一階向后差分方法進(jìn)行離散

對式(12)的時間微分采用一階向后差分方法，得到任意單元的全隱式非線性方程組，見式(13)：

式中，引入了組分κ＝w,g,a的余量△t為時間步長，上標(biāo)k和k+1表示兩相鄰的時間步長指標(biāo)；右端的流量項和源匯項均采用新的時間步長值；

(3)Newton-Raphson迭代方法

對式運(yùn)用Newton-Raphson迭代方法進(jìn)行線性化，引入迭代指標(biāo)p，對式(13)中的余量在迭代步p+1處進(jìn)行泰勒級數(shù)展開，只保留一階項，得到包含3×NEL個方程的線性方程組，NEL為計算域內(nèi)單元數(shù)，并且以兩迭代步的增量為未知量；最后得到大型稀疏系數(shù)矩陣的線性方程組，如式(14)：

還包括模型驗證步驟：運(yùn)用地下水氣-液二相流及溶質(zhì)運(yùn)移模型，模擬分析承壓含水層受咸、淡水共同作用的靜力平衡情況，并與已有文獻(xiàn)的計算結(jié)果進(jìn)行對比驗證；

還包括如下分析步驟：根據(jù)地下水氣-液二相流及溶質(zhì)運(yùn)移模型，以咸-淡水靜力平衡情況作為初始條件，采用依靠水頭的流量邊界，分析模擬海平面上升情況下海水入侵的效果：

(1)入侵程度變化：內(nèi)陸邊界采用依靠水頭的流量邊界時，隨著海平面的上升過程中，用相對鹽度等值線表征的海水入侵范圍始終介于流量恒定和水頭恒定中間，具體如下：當(dāng)內(nèi)陸邊界為流量恒定時，隨海平面的上升，海水入侵程度基本與初始情況保持不變；而當(dāng)內(nèi)陸邊界為水頭恒定時，隨著海平面的上升，海水入侵范圍明顯增大；當(dāng)內(nèi)陸邊界為依靠水頭的流量邊界時，隨著海平面的上升，海水入侵范圍有所增大，但入侵程度明顯比水頭恒定情況小；

(2)研究區(qū)域左右兩側(cè)等效水頭差和左側(cè)補(bǔ)給量的變化：內(nèi)陸邊界采用依靠水頭的流量邊界時，隨著海平面的上升，研究區(qū)域左右兩側(cè)等效水頭差和左側(cè)補(bǔ)給量都在逐漸減小，但是內(nèi)陸邊界為流量恒定時，以上兩項都基本保持不變，而水頭恒定時以上兩項減少的幅度要大于依靠水頭的流量邊界。

本發(fā)明的特點及有益效果是：

(1)建立了地下水氣-液二相流及溶質(zhì)運(yùn)移模型，可以用來模擬不同內(nèi)陸邊界條件對，由海平面上升引起的海水入侵的影響。

(2)本發(fā)明采用了基于假定的理想的概念模型來模擬海平面上升下海水入侵的內(nèi)陸邊界處理方法，這種依靠水頭的流量邊界，更加符合實際，因為它的補(bǔ)給會隨著海平面的上升而逐漸減小，它是由參考水頭(遠(yuǎn)處參考的某一水面)與海平面之間的高度差以及含水層的性質(zhì)共同決定的。內(nèi)陸邊界條件對海水入侵有至關(guān)重要的影響，因此提出一種更加合理的模擬海平面上升下海水入侵的內(nèi)陸邊界處理方法，對于合理預(yù)測沿海地區(qū)海水入侵的位置和移動情況，有效防治地下水污染，保護(hù)有限的地下淡水資源具有重大的現(xiàn)實意義，而且對于水資源管理者更好的評估和管理海平面上升對海水入侵潛在的長期影響具有一定的指導(dǎo)意義。

附圖說明：

圖1含水層的一般概念模型：(a)承壓含水層，(b)非承壓含水層。

圖2 Henry’s問題的研究范圍與邊界條件。

圖3穩(wěn)態(tài)下相對鹽度為0.5的等值線位置及其與其他文獻(xiàn)結(jié)果的對比圖。

圖4流量恒定(FC)和水頭恒定(HC)下，達(dá)到穩(wěn)定時相對鹽度為0.5的等值線分布。

圖5在承壓含水層中三種不同內(nèi)陸邊界在不同海平面上升情況下，達(dá)到穩(wěn)態(tài)時相對鹽度為0.5的等值線分布圖：(a)流量恒定系統(tǒng)邊界(FC systems),(b)依靠水頭的流量邊界(GHB,和(c)水頭恒定系統(tǒng)邊界(HC systems)。

圖6承壓和非承壓含水層中海平面上升前穩(wěn)態(tài)下相對鹽度為0.5的等值線分布。

圖7在非承壓含水層中三種不同內(nèi)陸邊界在不同海平面上升情況下，達(dá)到穩(wěn)態(tài)時相對鹽度為0.5的等值線分布圖：(a)流量恒定系統(tǒng)邊界(FC systems),(b)依靠水頭的流量邊界(GHB,和(c)水頭恒定系統(tǒng)邊界(HC systems)。

圖8本發(fā)明的模擬方法流程圖。

具體實施方式

本發(fā)明步驟如下：

步驟一：建立地下水氣-液二相流及溶質(zhì)運(yùn)移模型，包括基本控制方程及輔助方程，假定系統(tǒng)一直處于恒溫狀態(tài)。具體如下：

模型的基本控制方程為：

式中，M^κ表示κ組分(純水w、參考鹽水b和空氣a)的累積質(zhì)量密度，F(xiàn)^κ為κ組分的流量密度，包括平流流量密度和擴(kuò)散流量密度q^κ為組分κ的源匯項。

(1)累積質(zhì)量密度等于β相中組分κ的質(zhì)量之和，其表達(dá)式為：

式中，β表示流體相(氣相或液相)，液相為純水和參考鹽水的混合物，表示孔隙率，S_β為β相的飽和度，ρ_β為β相的密度，p_atm為κ組分占β相的質(zhì)量百分?jǐn)?shù)。

(2)平流流量p_g的表達(dá)式為：

式中，F(xiàn)_β為相的平流流量，符合達(dá)西定律，其多相流形式的表達(dá)式為：

其中，k為固有滲透系數(shù)，k_rβ為β相的相對滲透系數(shù)，μ_β為β相的粘滯性系數(shù)，p_β為β相的孔隙壓力，g為重力加速度矢量。

(3)擴(kuò)散流量的表達(dá)式為：

式中，為κ組分在β相中的分子擴(kuò)散系數(shù)；τ₀τ_β(S_β)為迂曲度，τ₀是多孔介質(zhì)特性參數(shù)，τ_β＝τ_β(S_β)是飽和度的函數(shù)。

分子擴(kuò)散系數(shù)隨壓力p和溫度T變化的計算式為：

式中，θ為溫度系數(shù)，通?？扇≈禐?.8。

TOUGH2中提供了三種迂曲度模型：分子擴(kuò)散系數(shù)τ₀一般取1.0

τ₀τ_β＝τ₀k_rβ(S_β) (相對滲透率模型)

(米林頓夸克模型)

τ₀τ_β＝τ₀S_β (常量擴(kuò)散) (7)

步驟二：模型求解：以非飽和地下水流及熱流傳輸TOUGH的第二個版本TOUGH2的水、鹽水和空氣的混合模塊EOS7為工具，空間上采用積分形式的有限差分方法(IFDM)進(jìn)行離散，時間上采用一階向后差分的全隱式方法進(jìn)行離散，模型的線性化采用牛頓-拉弗森Newton-Raphson迭代方法，最后得到大型稀疏系數(shù)矩陣的線性方程組；具體如下：

(1)空間上采用積分有限差分法(IFDM)進(jìn)行離散

首先將計算域離散成子單元，其性質(zhì)由形心點代表，分別對各個單元的質(zhì)量平衡方程的積分形式進(jìn)行空間離散。對于任意單元n，單元體積為V_n，邊界面為Γ_n，單元的質(zhì)量平衡方程的積分形式如下：

式中，n為表面單元dΓ_n的單位法向量，指向控制單元體內(nèi)為正。

引入適當(dāng)?shù)捏w積平均值，有

式中，為M^κ，q^κ在V_n上的平均值。

Γ_n上的面積分可近似為其所包含的各個表面A_nm的面積分的平均值之和，有

式中，m為與單元n相鄰的所有單元，A_nm是單元n和m相鄰的交界面，是F^κ在面A_nm上的內(nèi)法線方向的平均值。

將式(9)、(10)和(11)代入到式(8)中，得到一組關(guān)于時間的一階微分方程組

(2)時間上采用一階向后差分方法進(jìn)行離散

對式(12)的時間微分采用一階向后差分方法，得到任意單元的全隱式非線性方程組，見式(13)：

式中，引入了組分κ＝w,g,a的余量△t為時間步長，上標(biāo)k和k+1表示兩相鄰的時間步長指標(biāo)；右端的流量項和源匯項均采用新的時間步長值，這種全隱式方法提高了模型求解的數(shù)值穩(wěn)定性。

(3)Newton-Raphson迭代方法

對式運(yùn)用Newton-Raphson迭代方法進(jìn)行線性化。引入迭代指標(biāo)p，對式(13)中的余量在迭代步p+1處進(jìn)行泰勒級數(shù)展開，只保留一階項，得到包含3×NEL(計算域內(nèi)單元數(shù))個方程的線性方程組，并且以兩迭代步的增量為未知量；最后得到大型稀疏系數(shù)矩陣的線性方程組，如式(14)：

步驟三：模型邊界條件確定：模型計算中的邊界條件包括狄利克雷Dirichlet邊界條件和諾伊曼Neumann邊界條件兩種，其數(shù)學(xué)處理方法如下：

(1)Dirichlet邊界條件

Dirichlet邊界條件上，邊界條件單元的主要變量在計算過程中保持不變，為此，設(shè)定邊界條件單元的體積非常大，如1×10⁴⁰m³。當(dāng)邊界單元的體積相對于土體單元很大時，與土體單元的流量交換將不會改變邊界單元的主要變量值。

①對于空氣邊界，其邊界條件單元的相態(tài)為僅有氣相狀態(tài),主要變量為p_g,X_b,T，對于與大氣接觸的邊界上p_g＝p_atm，X_b＝0.0；T為氣溫。

②對于已知水頭邊界，包括地下淡水水頭邊界和海水水頭邊界，其邊界條件單元的相態(tài)均為僅有液相狀態(tài),主要變量為p_g、X_b、T，由于液相飽和狀態(tài)下毛細(xì)壓力為零，因此有p_g＝p_atm+ρ_lgh，其中ρ_l為淡水或海水的密度，h為邊界上的水深；地下淡水水頭邊界上X_b等于零，海水邊界上的X_b可根據(jù)參考鹽水的特性以及海水的鹽度計算出來。

(2)Neumann邊界條件

Neumann邊界條件描述的是系統(tǒng)與外界的流量交換情況，邊界條件單元的單位流量對應(yīng)式中的源匯項，流入為正，可以是常量，也可以隨時間變化。

①對于不透水邊界，是一類特殊的Neumann邊界條件，邊界上的流量為零，積分形式的有限差分法的處理很簡單，就是不設(shè)置與之相鄰的邊界條件單元，則與不透水面單元的流量交換為零。

②對于降雨入滲條件的處理方法，則可設(shè)置一個適當(dāng)大小的源匯項來實現(xiàn)。

(3)依靠水頭的流量邊界條件

依靠水頭的流量邊界條件受研究區(qū)域外的已知水位的較大地表水體(如濕地、河流、水庫等)影響，區(qū)域外的已知水位又稱參考水頭H_ref。水流流入或者流出研究區(qū)域的流量大小與研究區(qū)內(nèi)的水頭H和參考水頭H_ref之間的水頭差成比例，并且可以用如下線性方程式表示：

q＝-C(H-H_ref) (15)

其中q是從內(nèi)陸側(cè)流入研究區(qū)的單寬淡水補(bǔ)給量，C是導(dǎo)水系數(shù)。

而研究區(qū)域淡水側(cè)的等效淡水水頭H_f跟淡水補(bǔ)給量q_f在承壓含水層和非承壓含水層中分別有如下關(guān)系：

承壓含水層：

非承壓含水層：

其中K是含水層的滲透系數(shù)，α＝ρ_f/(ρ_s-ρf)是相對密度，ρ_f和ρ_s分別是淡水密度和海水密度，L和B分別為研究區(qū)域含水層的長度和厚度，H_s是距離含水層底部的平均海水深度。

把公式(15)分別帶入公式(16)和(17)中，并且讓H＝H_f和q＝q_f，可得出如下關(guān)系式：

承壓含水層：

非承壓含水層：

步驟四：模型驗證：運(yùn)用地下水氣-液二相流及溶質(zhì)運(yùn)移模型，模擬分析承壓含水層受咸、淡水共同作用的靜力平衡情況，并與已有文獻(xiàn)的計算結(jié)果進(jìn)行對比驗證；

步驟五：根據(jù)地下水氣-液二相流及溶質(zhì)運(yùn)移模型，以咸-淡水靜力平衡情況作為初始條件，采用依靠水頭的流量邊界，分析模擬海平面上升情況下海水入侵的效果：

(1)入侵程度變化：內(nèi)陸邊界采用依靠水頭的流量邊界時，隨著海平面的上升過程中，用相對鹽度等值線表征的海水入侵范圍始終介于流量恒定和水頭恒定中間。具體如下：當(dāng)內(nèi)陸邊界為流量恒定時，隨海平面的上升，海水入侵程度基本與初始情況保持不變；而當(dāng)內(nèi)陸邊界為水頭恒定時，隨著海平面的上升，海水入侵范圍明顯增大；當(dāng)內(nèi)陸邊界為依靠水頭的流量邊界時，隨著海平面的上升，海水入侵范圍有所增大，但入侵程度明顯比水頭恒定情況小。

(2)研究區(qū)域左右兩側(cè)等效水頭差和左側(cè)補(bǔ)給量的變化：內(nèi)陸邊界采用依靠水頭的流量邊界時，隨著海平面的上升，研究區(qū)域左右兩側(cè)等效水頭差和左側(cè)補(bǔ)給量都在逐漸減小，但是內(nèi)陸邊界為流量恒定時，以上兩項都基本保持不變，而水頭恒定時以上兩項減少的幅度要大于依靠水頭的流量邊界。

下面結(jié)合附圖，針對模擬海平面上升下海水入侵的內(nèi)陸邊界處理方法進(jìn)行具體說明，其步驟如下：

(1)建立模擬海水入侵的數(shù)學(xué)模型：該數(shù)學(xué)模型描述的是承壓含水層未受污染的淡水受海水入侵的情況，即經(jīng)典的Henry’s問題。模型選取200米長、100米深、1米厚的含水層為研究對象，其上、下兩層不透水頂板、底板之間的含水層均質(zhì)、各向同性，內(nèi)陸側(cè)邊界上有恒定的淡水流入量(或者等效的淡水作用下的靜水壓力)，海水側(cè)邊界上分布著密度較大的海水作用下的靜水壓力，是一個理想化的準(zhǔn)二維模型。其研究范圍及邊界條件見圖2，模型相關(guān)參數(shù)取值見表1。假定含水層系統(tǒng)處于恒溫狀態(tài)(T＝25℃)。

模型的基本控制方程為：

式中，M^κ表示κ組分(純水w、參考鹽水b和空氣a)的累積質(zhì)量密度，F(xiàn)^κ為組分κ的流量密度，包括平流流量密度和擴(kuò)散流量密度q^κ為組分κ的源匯項。

①累積質(zhì)量密度等于β相中組分κ的質(zhì)量之和，其表達(dá)式為：

式中，β表示流體相(氣相或液相)，液相為純水和參考鹽水的混合物，表示孔隙率，S_β為β相的飽和度，ρ_β為β相的密度，p_atm為κ組分占β相的質(zhì)量百分?jǐn)?shù)。

②平流流量p_g的表達(dá)式為：

式中，F(xiàn)_β為相的平流流量，符合達(dá)西定律，其多相流形式的表達(dá)式為：

其中，k為固有滲透系數(shù)，k_rβ為β相的相對滲透系數(shù)，μ_β為β相的粘滯性系數(shù)，p_β為β相的孔隙壓力，g為重力加速度矢量。

③擴(kuò)散流量的表達(dá)式為：

式中，為κ組分在β相中的分子擴(kuò)散系數(shù)；τ₀τ_β(S_β)為迂曲度，τ₀是多孔介質(zhì)特性參數(shù)，τ_β＝τ_β(S_β)是飽和度的函數(shù)。

分子擴(kuò)散系數(shù)隨壓力p和溫度T變化的計算式為：

式中，θ為溫度系數(shù)，通?？扇≈禐?.8。

TOUGH2中提供了三種迂曲度模型：分子擴(kuò)散系數(shù)τ₀一般取1.0

τ₀τ_β＝τ₀k_rβ(S_β) (相對滲透率模型)

(米林頓夸克模型)

τ₀τ_β＝τ₀S_β (常量擴(kuò)散) (7)

(2)模型求解：對計算域進(jìn)行網(wǎng)格剖分，含水層被劃分為5m×5m×1m的8節(jié)點的均勻六面體單元，共800個單元，1722個節(jié)點。以TOUGH2/EOS7為工具，空間上采用積分形式的有限差分方法(IFDM)進(jìn)行離散，時間上采用一階向后差分的全隱式方法進(jìn)行離散，見式(13)；模型的線性化采用Newton-Raphson迭代方法，最后得到大型稀疏系數(shù)矩陣的線性方程組，見式(14)：

式中，引入了組分κ＝w,g,a的余量表示M^κ在控制體積V_n上的平均值；表示q^κ在控制體積V_n上的平均值；m為與單元n相鄰的所有單元；A_nm是單元n和m相鄰的交界面的面積；是F^κ在面A_nm上的內(nèi)法線方向的平均值；△t為時間步長；上標(biāo)k和k+1表示兩相鄰的時間步指標(biāo)。

引入迭代指標(biāo)p，對式(13)中余量在迭代步p+1處進(jìn)行泰勒級數(shù)展開，只保留一階項，得到包含3×NEL(計算域內(nèi)單元數(shù))個方程的線性方程組，并且以兩迭代步的增量為未知量：

(3)模型邊界條件確定：假定系統(tǒng)處于恒溫狀態(tài)T＝25℃在整個模擬過程中保持不變。在初始條件下，計算域均為液相飽和狀態(tài)，毛細(xì)壓力可以忽略，相對滲透系數(shù)取值1.0，主要變量為p_g、X_b、T。因此有p_g＝p_atm+ρ_lgh，其中ρ_l為淡水或海水的密度，h為邊界上的水深。對于淡水側(cè)(左側(cè))，作用淡水源匯項，將Q_in＝V_indρ_w(d為含水層深度100m，ρ_w為淡水密度)等于7.639×10^-3kg/s，均勻分布在左側(cè)邊界上，主要變量p_g＝p_atm+ρ_wg(100-z)，X_b＝0.0，和T＝25℃；對于海水側(cè)(右側(cè))，受海水作用，主要變量p_g＝p_atm+ρ_sg(100-z)，X_b＝1.0，和T＝25℃；其中，p_atm為大氣壓力，等于1.013×10⁵Pa。初次計算的初始條件為：所有單元液相飽和，壓力均為大氣壓力與靜水(淡水)壓力之和即p_g＝p_atm+ρ_wg(100-z)，X_b＝0.0，

以上實際上屬于流量恒定的情況，為了使初始穩(wěn)態(tài)條件相同，當(dāng)水頭恒定時，由公式(16)計算可知，在承壓含水層中恒定的等效水頭H_f＝102.57米(q_f＝0.66m²/d且H_s＝100m)。因此其內(nèi)陸邊界條件設(shè)置為P_l＝P_atm+ρ_wg(102.57-z),X_b＝1.0 and然而依靠水頭的流量的邊界條件的設(shè)置中，假定參考水頭H_ref＝105m，因此由公式(15)知，C＝0.2716m/s(q_f＝0.66m²/d且H_f＝102.57m)，且其初始條件的設(shè)置與流量恒定完全相同。

(4)模型驗證：

在靜力平衡情況下，相對(海水)鹽度X_b＝0.5的等值線位置及其與已有文獻(xiàn)計算結(jié)果的對比如圖3。由圖可見，本模型的計算結(jié)果與其他文獻(xiàn)的結(jié)果基本一致，從而驗證了計算模型模擬海水入侵問題的有效性。

(5)根據(jù)地下水氣-液二相流及溶質(zhì)運(yùn)移模型，分析不同內(nèi)陸邊界條件對模擬海平面上升的海水入侵的效果：

以咸-淡水靜力平衡情況作為初始條件，由圖4可知三種不同內(nèi)陸邊界條件的海水入侵情況是等效的(圖4沒有標(biāo)注依靠水頭的流量(GHB)邊界情況，因為初始情況下，GHB邊界與流量邊界完全相同)。

為了比較不同內(nèi)陸邊界條件對海平面上升的海水入侵的影響，右側(cè)海水邊界的海水位分別設(shè)置為100.25,100.50,100.75 and 101.0m(對應(yīng)于海平面分別上升0.25,0.50,0.75 and 1.0m)。對于流量恒定和水頭恒定的內(nèi)陸邊界條件，左側(cè)(內(nèi)陸邊界)補(bǔ)給流量q_f＝0.66m²/d和等效淡水水頭H_f＝102.57m分別保持不變。而對于水頭依靠的流量內(nèi)陸邊界條件，由公式(19)可知，補(bǔ)給流量q_f隨著海水位H_s的增大而減小，具體的q_f分別為0.61,0.57,0.52,和0.48m²/d對應(yīng)于H_s分別為100.25,100.50,100.75和101.0m。對于水頭依靠的流量內(nèi)陸邊界條件的設(shè)置跟流量恒定條件的設(shè)置類似，除了隨海水位變化的q_f。

圖5比較了在承壓含水層中三種不同內(nèi)陸邊界在不同海平面上升情況下，達(dá)到穩(wěn)態(tài)時相對鹽度為0.5的等值線分布。結(jié)果表明，內(nèi)陸邊界條件對模擬海平面上升引起的海水入侵有至關(guān)重要的影響。水頭恒定(HC)的內(nèi)陸邊界條件對海平面上升引起的海水入侵最為敏感、脆弱，因為隨著海平面的上升，左右兩側(cè)的水力梯度逐漸減小；而流量恒定(FC)的內(nèi)陸邊界條件對模擬海平面上升引起的海水入侵卻幾乎不變，因為海平面上升時左側(cè)淡水水頭同等上升，水力梯度幾乎保持不變；而依靠水頭的流量邊界(GHB)條件介于以上兩種邊界之間，因為它的淡水補(bǔ)給會隨著海平面的上升而逐漸減小，它是由參考水頭(遠(yuǎn)處參考的某一水面)與海平面之間的水頭差以及含水層的性質(zhì)共同決定的。

為了研究非承壓含水層的模擬內(nèi)陸邊界對海平面上升的海水入侵的情況，首先把Henry問題擴(kuò)展到非承壓區(qū)域，保持相同的邊界條件和參數(shù)，但是上層邊界為潛水自由水面且增加了5米，網(wǎng)格離散化范圍為0.25到1米不等。因此非承壓含水層研究區(qū)域的長為200m，厚度為105米。為了更好的研究海平面上升0到1米的對海水入侵的情況，網(wǎng)格在100米上下被細(xì)化。

非承壓含水層邊界條件的設(shè)置類似于承壓含水層，對于流量恒定的內(nèi)陸邊界條件，左側(cè)單寬補(bǔ)給q_f＝0.66m²/d保持不變。對于水頭恒定的內(nèi)陸邊界條件，等效淡水水頭由公式(17)可知，H_f＝102.54m(q_f＝0.66m²/d且H_s＝100m)，而對于依靠水頭的流量的邊界條件，由公式(15)知，C＝0.2681m/s(q_f＝0.66m²/d，H_ref＝105且H_f＝102.57m)。當(dāng)右側(cè)海水邊界的海水位分別設(shè)置為100.25,100.50,100.75 and 101.0m(對應(yīng)于海平面分別上升0.25,0.50,0.75 and 1.0m)時，對于流量恒定和水頭恒定的內(nèi)陸邊界條件，左側(cè)(內(nèi)陸邊界)補(bǔ)給流量q_f＝0.66m²/d和等效淡水水頭H_f＝102.54m分別保持不變。而對于水頭依靠的流量內(nèi)陸邊界條件，由公式(21)可知，補(bǔ)給流量q_f隨著海水位H_s的增大而減小，具體的q_f分別為0.62,0.57,0.53,和0.48m²/d對應(yīng)于H_s分別為100.25,100.50,100.75和101.0m。

由圖6承壓和非承壓含水層中海平面上升前穩(wěn)態(tài)下相對鹽度為0.5的等值線分布可知，這四條曲線基本重合，表明海水入侵過程不論在承壓還是非承壓含水層中，其初始穩(wěn)態(tài)條件下相一致。對于海平面上升情況下，不同內(nèi)陸邊界對海水入侵的影與承壓含水層情況基本一致(圖7)。

綜上所述，對本實施例總結(jié)如下：

(1)基于標(biāo)準(zhǔn)的理想的概念模型，采用數(shù)值模擬的方法研究了不同內(nèi)陸邊界條件對海平面上升引起的海水入侵效果。首先利用地下水氣-液二相流及溶質(zhì)運(yùn)移模型模擬分析了經(jīng)典的Henry問題的咸-淡水靜力平衡情況，通過與以前文獻(xiàn)結(jié)果的對比，驗證了所采用模型的有效性。然后基于咸-淡水靜力平衡情況，模擬分析了不同內(nèi)陸邊界條件對海平面上升引起的海水入侵效果。

(2)對模擬結(jié)果的分析表明，內(nèi)陸邊界條件對模擬海平面上升引起的海水入侵有至關(guān)重要的影響。水頭恒定(HC)的內(nèi)陸邊界條件對海平面上升引起的海水入侵最為敏感、脆弱，因為隨著海平面的上升，左右兩側(cè)的水力梯度逐漸減??；而流量恒定(FC)的內(nèi)陸邊界條件對模擬海平面上升引起的海水入侵卻幾乎不變，因為海平面上升時左側(cè)淡水水頭同等上升，水力梯度幾乎保持不變；而依靠水頭的流量邊界(GHB)條件介于以上兩種邊界之間且更符合實際。因為它的淡水補(bǔ)給會隨著海平面的上升而逐漸減小，它是由參考水頭(遠(yuǎn)處參考的某一水面)與海平面之間的水頭差以及含水層的性質(zhì)共同決定的。

(3)依靠水頭的流量內(nèi)陸邊界條件是一種更加合理的模擬海平面上升下海水入侵的內(nèi)陸邊界處理方法，這對于合理預(yù)測沿海地區(qū)海水入侵的位置和移動情況，有效防治地下水污染，保護(hù)有限的地下淡水資源具有重大的現(xiàn)實意義，而且對于水資源管理者更好的評估和管理海平面上升對海水入侵潛在的長期影響具有一定的指導(dǎo)意義。