邊文英

(河北省地礦局國土資源勘查中心，石家莊 050000)

0 引言

全球有約三分之一的陸地位于干旱與半干旱地區(qū)，而且世界上近一半國家都在不同程度上受到干旱缺水問題的困擾。我國也不例外，尤其是西北地區(qū)干旱缺水問題更為突出，并且研究區(qū)蒸發(fā)十分強(qiáng)烈，降水稀少[1-4]。因此為了更好的研究在此條件下干旱區(qū)地下水的動態(tài)變化，本文在認(rèn)真研究該區(qū)域水文地質(zhì)條件特征的基礎(chǔ)上，提出通過解析法計(jì)算傍河地下水在蒸發(fā)條件下典型剖面的潛水面分布情況，為研究此類地區(qū)水文地質(zhì)條件提供了新思路。

本方法通過概化邊界條件、建立數(shù)學(xué)模型，分別按照定水頭和定流量兩種補(bǔ)給方式，計(jì)算干旱區(qū)傍河地下水蒸發(fā)型一維穩(wěn)定流潛水面的分布情況，此方法能夠通過較少的參數(shù)，獲取剖面地下水潛水面的分布情況，并且精度較高，對于具有相似水文地質(zhì)條件地區(qū)可以直接利用本次計(jì)算成果，從而能夠大大節(jié)省時間、減少工作量，為研究局部剖面地下水潛水面分布提供依據(jù)。

1 數(shù)學(xué)模型及其求解方法

本文以選擇一個傍河地下水蒸發(fā)型穩(wěn)定流典型剖面為例進(jìn)行介紹，典型剖面圖見圖1?，F(xiàn)根據(jù)剖面的具體情況采用解析法對傍河地下水蒸發(fā)型一維穩(wěn)定流進(jìn)行分析。

1.1 數(shù)學(xué)模型

根據(jù)研究區(qū)的含水層結(jié)構(gòu)、地下水系統(tǒng)的劃分以及當(dāng)?shù)氐湫偷臍夂蛱卣鳎_定研究區(qū)的數(shù)學(xué)模型為源匯項(xiàng)只包括潛水蒸發(fā)的一維穩(wěn)定流的潛水模型。

模型的水流控制方程為：

(1)

圖1 典型剖面圖Figure 1 Typical section

模型的邊界條件：

h(x=0)=h0

(2)

(3)

式中：h為潛水位(m)；x為水平坐標(biāo)；t為時間(d)；K為x方向的滲透系數(shù)(m/d)；Eg為潛水面的蒸發(fā)量(m/d)；h0典型剖面所傍河的水位高程；L為地下水位穩(wěn)定時典型剖面距離河道的距離。

傍河地下水蒸發(fā)型穩(wěn)定流解析解的具體求解方法如下。

1.2 定水頭補(bǔ)給潛水面分布

1.2.1 線性蒸發(fā)定水頭補(bǔ)給潛水面分布

根據(jù)潛水蒸發(fā)公式的不同選擇分為線性蒸發(fā)公式和非線性蒸發(fā)公式兩種方式求其解析解。線性蒸發(fā)根據(jù)定水頭補(bǔ)給的數(shù)學(xué)模型，通過整理變形的不同，又分為兩種。下面首先介紹線性蒸發(fā)公式的第一種變形求解方法。

(1)線性蒸發(fā)公式第一種變形

線性蒸發(fā)公式：

(4)

式中Dmax為潛水蒸發(fā)為零時的地下水極限埋深，E0為水面蒸發(fā)強(qiáng)度，h為潛水含水層厚度，z為潛水含水層底板到地面的高度。

定水頭補(bǔ)給-數(shù)學(xué)模型：

(5)

h(x=0)=h0

(6)

(7)

將(4)代入(5)：

(8)

(9)

(10)

設(shè)y=h2，則(8)轉(zhuǎn)化成(11)

(11)

y″-By=A

(12)

方程(12)為二階常系數(shù)非齊次線性微分方程，因此若想求(12)的通解，歸結(jié)為求對應(yīng)的齊次方程y″-By=0的通解和非齊次方程(12)本身的一個特解。

齊次方程y″-By=0的通解：

為常數(shù))

(13)

非齊次方程(12)的特解：

(14)

則方程(12)的通解為：

(15)

即：

(16)

根據(jù)控制方程(6)、(7)可求得常數(shù)C1、C2的值：

(17)

(18)

將(17)和(18)值代入(16)，整理得(19)：

(19)

再將(9)和(10)代入(19)得最終定水頭補(bǔ)給線性蒸發(fā)的潛水面形狀方程:

(20)

根據(jù)所求出的定水頭補(bǔ)給線性蒸發(fā)的潛水面形狀方程，在理想狀態(tài)下得到其潛水面形狀，見圖2。假設(shè)h0=50m，z=53m，k=20m/d，Dmax=5m，E0=0.005 479m/d，L=1 000m。

圖2 定水頭補(bǔ)給潛水面形狀圖Figure 2 Phreatic water table form under constantwater head recharge

(2) 線性蒸發(fā)公式第二種變形

線性蒸發(fā)公式：

(21)

定水頭補(bǔ)給—數(shù)學(xué)模型：

(22)

h|x=0

(23)

(24)

將(21)代入(22)最終得(25)：

(25)

(26)

(27)

h″-Bh=A

(28)

方程(28)為二階常系數(shù)非齊次線性微分方程，因此若想求它的通解，先求對應(yīng)的齊次方程h″-Bh=0的通解和非齊次方程(28)本身的一個特解。

齊次方程h″-Bh=A的通解：

為常數(shù))

(29)

非齊次方程(28)特解：

(30)

則方程(25)的通解為：

(31)

根據(jù)控制方程(23)、(24)得C1、C2的值：

(32)

(33)

將(32)和(33)代入(31)，整理得：

(34)

再將(26)和(27)代入(34)得最終定水頭補(bǔ)給的線性蒸發(fā)的潛水面形狀方程:

(35)

根據(jù)所求出的定水頭補(bǔ)給線性蒸發(fā)的潛水面形狀方程，在理想狀態(tài)下得到其潛水面形狀，見圖3。假設(shè)h0=50m，z=53m，k=20m/d，Dmax=5m，E0=0.005 479m/d，L=1 000m。

圖3 定水頭補(bǔ)給潛水面形狀圖Figure 3 Phreatic water table form under constantwater head recharge

1.2.2 非線性蒸發(fā)定水頭補(bǔ)給潛水面分布

由于非線性蒸發(fā)的公式比較復(fù)雜，直接求其解析解存在一定難度，因此本文采用龍格-庫塔方法(數(shù)學(xué)手冊，2005)[5]。下面簡單介紹龍格-庫塔方法的計(jì)算公式和步驟。

(1)龍格-庫塔方法

使用龍格-庫塔方法計(jì)算，需要兩個前提，即存在含多個未知函數(shù)的微分方程組：

(36)

(37)

初始條件：y(x0)=y0,z(x0)=z0

(38)

龍格-庫塔方法具體計(jì)算公式:

(39)

其中：

先計(jì)算k1、k2、k3、k4、l1、l2、l3、l4，再計(jì)算yn+1和zn+1。

(2)非線性蒸發(fā)潛水蒸發(fā)面的求解

目前有關(guān)潛水蒸發(fā)的經(jīng)驗(yàn)公式較多,比較典型的主要有:阿維里揚(yáng)諾夫的潛水穩(wěn)定蒸發(fā)拋物線型公式(1958),葉水庭等人(1977)提出的潛水蒸發(fā)指數(shù)型公式;隨著研究的深入,無作物條件下計(jì)算潛水蒸發(fā)的阿維里揚(yáng)諾夫公式、指數(shù)型公式已被用于有作物條件下潛水蒸發(fā)的計(jì)算,特別是阿維里揚(yáng)諾夫經(jīng)驗(yàn)公式的應(yīng)用更為廣泛。但是由于阿維里揚(yáng)諾夫公式當(dāng)大氣蒸發(fā)力較大時會過高估計(jì)潛水蒸發(fā)量，而指數(shù)型公式可以彌補(bǔ)這一問題[6-12]。因此本次非線性潛水蒸發(fā)選取葉水庭公式(1977)，葉水庭公式如下：

Eg=E0e-αD=E0e-α(z-h)

(40)

式中，α為衰減系數(shù)，通過實(shí)測資料確定;D為潛水埋深,E0為水面蒸發(fā)，z為地面高程，h為潛水面高程。

(41)

h(x=0)=h0

(42)

(43)

將(40)代入(41)得：

(44)

將(44)做無量綱化處理最終得(45)：

(45)

其中：xD=x/L,hD=h/h0,zD=z/h0

(46)

無量綱數(shù)學(xué)模型-定水頭補(bǔ)給：

(47)

hD(xD=0)=1

(48)

(49)

將(47)分解為2個常微分方程，構(gòu)成方程組

(50)

(51)

在excel表格中從xD=0開始計(jì)算,知道hD(0)和y(0)。先假設(shè)一個y(0)=C(常數(shù))。計(jì)算一系列hD(x)和y(x)的數(shù)值，直到獲得y(xD=1)=0，可能不等于零。調(diào)整y(0)的猜測值，直到y(tǒng)(x=1)=0或無限接近于0為止。計(jì)算結(jié)果見圖4。

圖4 定水頭補(bǔ)給潛水面形狀圖Figure 4 Phreatic water table form under constantwater head recharge

1.3 定流量補(bǔ)給潛水面分布

定流量補(bǔ)給的潛水面分布的求解方法和定水頭求解方法和步驟基本相同，只是數(shù)學(xué)模型發(fā)生變化。

1.3.1 線性蒸發(fā)定流量補(bǔ)給潛水面分布

(1)線性蒸發(fā)公式第一種變形

線性蒸發(fā)公式：

(52)

式中Dmax為潛水蒸發(fā)為零時的地下水極限埋深，E0為水面蒸發(fā)強(qiáng)度，h為潛水層厚度，z為潛水層地板到地面的高度。

定流量補(bǔ)給-數(shù)學(xué)模型：

(53)

(54)

(55)

所求方程和定水頭求解方法一致，只是控制條件發(fā)生變化，即C1、C2變化，因此不再贅述。

即：

(56)

(57)

(58)

將(57)和(58)代入(56)中，整理得：

(59)

再將(9)和(10)代入(59)得最終定流量補(bǔ)給的線性蒸發(fā)的潛水面形狀方程：

(60)

根據(jù)所求出的定流量補(bǔ)給線性蒸發(fā)的潛水面形狀方程，在理想狀態(tài)下得到其潛水面形狀，見圖5。假設(shè)h0=50m，z=53m，k=20m/d，Dmax=5m，E0=0.005 479m/d，L=1 000m，q=1.33m3/d。

圖5 定流量補(bǔ)給潛水面形狀圖Figure 5 Phreatic water table form underconstant flow recharge

(2) 線性蒸發(fā)公式第二種變形

線性蒸發(fā)公式：

(61)

定流量補(bǔ)給—數(shù)學(xué)模型：

(62)

(63)

(64)

所求方程和定水頭一致，只是控制條件發(fā)生變化，即C1、C2變化。

即：

(65)

(66)

(67)

將(66)和(67)值代入(65)，整理得：

(68)

再將(26)和(27)代入(68)得最終定流量補(bǔ)給的線性蒸發(fā)的潛水面形狀方程:

(69)

根據(jù)所求出的定流量補(bǔ)給線性蒸發(fā)的潛水面形狀方程，在理想狀態(tài)下得到其潛水面形狀，見圖6。假設(shè)h0=50m，z=53m，k=20m/d，Dmax=5m，E0=0.005 479m/d，L=1 000m，q=1.63m3/d。

圖6 定流量補(bǔ)給潛水面形狀圖Figure 6 Phreatic water table form under constant flow recharge

1.3.2 非線性蒸發(fā)定流量補(bǔ)給潛水面分布

葉水庭公式：

Eg=E0e-αD=E0e-α(z-h)

(70)

式中，α為衰減系數(shù)，通過實(shí)測資料確定;D為潛水埋深,E0為水面蒸發(fā)，z為地面高程，h為潛水面高程。

數(shù)學(xué)模型：

(71)

(72)

(73)

將(70)代入(71)得：

(74)

將(74)做無量綱化處理得(75)：

(75)

(76)

其中：xD=x/L,hD=h/z

無量綱數(shù)學(xué)模型-定流量補(bǔ)給：

(77)

(78)

(79)

將(77)分解為2個常微分方程，構(gòu)成方程組

(80)

(81)

在excel表格中從xD=0開始計(jì)算,知道hD(0)和y(0)。先假設(shè)一個hD(0)=C(常數(shù))，y(0)=-q/K。計(jì)算一系列hD(x)和y(x)的數(shù)值，直到獲得y(xD=1)=0，可能不等于零。調(diào)整hD(0)的猜測值，直到y(tǒng)(xD=1)=0或無限接近于0為止。計(jì)算結(jié)果見圖7。

圖7 定流量補(bǔ)給潛水面形狀圖Figure 7 Phreatic water table form under constant flow recharge

2 結(jié)論

①通過對一維蒸發(fā)型穩(wěn)定流進(jìn)行分析和計(jì)算，得到了一維穩(wěn)定流在線性蒸發(fā)條件下的定水頭補(bǔ)給和定流量補(bǔ)給的解析解；

②通過龍格-庫塔方法，在excel中實(shí)現(xiàn)了一維穩(wěn)定流在非線性蒸發(fā)條件下的定水頭補(bǔ)給和定流量補(bǔ)給的潛水面形態(tài)；并且在計(jì)算的過程中分別對定水頭補(bǔ)給和定流量補(bǔ)給做了無量綱化處理，計(jì)算結(jié)果對傍河典型剖面具有一定適用性。