陽 玲，杜金月，王同科，趙志學(xué)，郝永紅

(1.天津師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院，天津 300387；2.天津師范大學(xué)天津市水資源與水環(huán)境重點(diǎn)實(shí)驗(yàn)室，天津 300387)

世界上有40%的人口居住在離海岸線100 km以內(nèi)的海濱城市，給近海海域帶來了大量污染，破壞了生態(tài)環(huán)境[1-3]。源于陸地的海洋污染有44%是來自于地下水海底排泄和地表徑流[1]。雖然地下水海底排泄的流量是地表河流入海徑流量的40%[2]，但是它的污染物濃度卻是地表徑流的10倍[3]。海底地下水排放被視為海洋污染最主要的方式之一[4-5]。預(yù)測海底地下水排放和輸送到海洋中的污染物對海洋生態(tài)和環(huán)境保護(hù)有著重要意義。

水文地質(zhì)參數(shù)識別又稱水文地質(zhì)逆問題求解，是建立地下水?dāng)?shù)值模型的關(guān)鍵步驟[6]。含水層的壓力傳導(dǎo)系數(shù)作為描述地下水溶質(zhì)運(yùn)移的關(guān)鍵參數(shù)被廣泛的研究[7-10]。傳統(tǒng)識別壓力傳導(dǎo)系數(shù)最廣泛的技術(shù)是抽水試驗(yàn)和振蕩試驗(yàn)[11-13]，但是這些方法有一個(gè)共同的缺點(diǎn)：需要耗費(fèi)大量的人力、財(cái)力和物力。為了更好地描述含水層的特性，有必要開發(fā)一種經(jīng)濟(jì)高效的方法來估計(jì)含水層的壓力傳導(dǎo)系數(shù)。

海洋和河流的潮汐可以看作自然的振蕩抽水試驗(yàn)，可通過分析潮汐和地下水響應(yīng)關(guān)系估計(jì)壓力傳導(dǎo)系數(shù)。Ferris估計(jì)了潮汐河含水層的壓力傳導(dǎo)系數(shù)，他使用了時(shí)間滯后和潮汐與地下水位相關(guān)的潮汐效率因素的方法[14]。Erskine估計(jì)了非承壓和高度滲透沿海含水層的壓力傳導(dǎo)系數(shù)，該含水層位于英國的一個(gè)核電站附近，他用了時(shí)間滯后和潮汐效率因素方法[15]。Zhou利用觀測井觀測地下水水位振幅和潮汐確定了T/S(給水度與導(dǎo)水系數(shù)之比),寫出了T/S的表達(dá)式[16]。在貝爾港沿海巖溶含水層的背景下，類似的壓力傳導(dǎo)系數(shù)的估計(jì)表明，水動(dòng)力環(huán)境包括基質(zhì)、裂隙和管道流動(dòng)系統(tǒng)。

然而，前人運(yùn)用地下水水位對潮汐的響應(yīng)識別含水層的參數(shù)多采用數(shù)值擬合的方法。在前人的研究基礎(chǔ)上，本文試圖獲得潮汐產(chǎn)生的壓力波在含水層中傳播的解析解，通過數(shù)值擬合估計(jì)含水層的壓力傳導(dǎo)系數(shù)。其能夠更加清楚地了解潮汐影響地下水的機(jī)理。

1 地下水模型

假設(shè)海濱區(qū)有一個(gè)水平、均質(zhì)、等厚的含水層，它位于兩弱透水層之間(圖1)。因?yàn)榈叵滤膲毫λ^高于含水層和海水的高度，所以地下水向海中排泄。

基于圖1的概念模型，建立笛卡爾直角坐標(biāo)系，x軸的原點(diǎn)位于內(nèi)陸，距海岸的水平距離為l m，且地下水水位是定水位H1。x軸是水平的，沿著含水層向海岸方向延伸，垂直方向?yàn)楹０毒€方向。在實(shí)際條件中左邊界為定水位邊界的可能性不大，更多的是流量邊界。假設(shè)有一條與海岸平行的運(yùn)河，保證左邊界為定水位邊界。該模型需假設(shè)有一條與海岸平行的運(yùn)河，且該運(yùn)河與海不相通。

圖1 濱海區(qū)含水層概念模型

承壓含水層中地下水運(yùn)動(dòng)的一維偏微分方程為：

(1)

式中：H——水頭/m；

x——水平坐標(biāo)/m；

t——時(shí)間/h；

T——含水層的滲透系數(shù)/(m·h-1)；

S——含水層儲水系數(shù)；

l——含水層的長度/m。

(2)

邊界條件為：

H|x=0=H1，H|x=l=H2+Asinωt，t>0，

(3)

式中：H1——內(nèi)陸含水層定水頭邊界/m；

H2——海水平均水位/m；

A——海水振蕩振幅/m；

ω——角速度/h-1。

2 地下水位對潮汐響應(yīng)的解析解

將式(2)～(3)轉(zhuǎn)化到復(fù)數(shù)域上[17-18],即將有關(guān)物理量從實(shí)平面擴(kuò)充到復(fù)平面[19]。變量x和t是相互獨(dú)立的。因此可以用分離變量法求解：

H(x,t)=u1(x,t)+u2(x,t)

(4)

(5)

將式(4)代入式(1)～(3)中：

(6)

u2(0,t)=0,u2(l,t)=Asinωtt>0

(7)

因?yàn)槭?7)關(guān)于時(shí)間是周期的，所以可以省略初始條件。把式 (4)～ (6)擴(kuò)張到復(fù)數(shù)域，假設(shè)u2(x,t)=lm(U(x,t)),代入式(6)～(7)，可以得到：

(8)

U(0,t)=0,U(l,t)=Aeiωt,t>0

(9)

(10)

(11)

最終求解結(jié)果：

(12)

式(12)可以簡化為：

(13)

其中：

(14)

(15)

3 識別含水層的壓力傳導(dǎo)系數(shù)

(16)

(17)

4 數(shù)值仿真

為了驗(yàn)證潮汐信號衰減率方法的正確性，了解該方法的精度，對式(13)中的參數(shù)賦值(表1)，進(jìn)行了數(shù)值仿真。賦值后由式(13)獲得的地下水水位對潮汐的響應(yīng)圖，見圖2。

表1 含水層參數(shù)

圖2 地下水對潮汐的響應(yīng)

為更加清晰地了解不同位置地下水響應(yīng)情況，圖3表達(dá)了x=30 m和x=80 m以及海水水位三種不同情況下地下水水位隨時(shí)間的變化。

圖3 x=30 m和x=80 m以及海水水位三種不同情況下地下水位隨時(shí)間的變化

為檢驗(yàn)潮汐衰減率r(x,t,ω)與余切函數(shù)cotωt之間的線性關(guān)系。當(dāng)x=80 m時(shí)地下水水位，r(x,t,ω)與cotωt呈線性相關(guān)(圖4)。

圖4 x=80 m時(shí)r(x,t,ω)與cotωt的線性關(guān)系

為了驗(yàn)證潮汐信號衰減率方法的正確性，模擬野外地下水水位的觀測過程，對x=80 m處的地下水?dāng)?shù)據(jù)進(jìn)行取樣，取樣時(shí)間步長為15 min，共獲得48個(gè)采樣點(diǎn)(圖5)。

由圖6可見，r(80,t,ω)與cotωt呈線性關(guān)系。由于采樣步長15 min，相對于樣本總體時(shí)間較小，圖6中部分點(diǎn)重合。

圖5 在x=80 m處，時(shí)間間隔為0.25 h,地下水對潮汐響應(yīng)的合成數(shù)據(jù)

圖6 潮汐信號衰減率r(x,t,ω)與cotωt的線性關(guān)系

圖7 用α和β估計(jì)初值

5 結(jié)論

(1)在濱海含水層中，地下水水位與海洋潮汐的交互作用類似于自然的振蕩抽水試驗(yàn)。當(dāng)漲潮時(shí)，相當(dāng)于給地下水注水，地下水水位上升。當(dāng)退潮時(shí)，相當(dāng)于抽取地下水，地下水水位下降?；诔毕盘栐诤畬又袀鞑ズ退p的過程，建立地下水模型求得解析解，提出用潮汐信號衰減率方法估計(jì)含水層的特性。潮汐信號衰減率r(x,t,ω)與海水振蕩的余切函數(shù)cotωt線性相關(guān)。

(2)線性關(guān)系中的斜率和截距都是壓力傳導(dǎo)系數(shù)的函數(shù)。提出了利用潮汐信號衰減率與余切函數(shù)的關(guān)系估計(jì)壓力傳導(dǎo)系數(shù)的方法(即用斜率和截距估計(jì)壓力傳導(dǎo)系數(shù))。數(shù)值仿真結(jié)果表明，該方法可以準(zhǔn)確地估算含水層的壓力傳導(dǎo)系數(shù)D。

(3)潮汐衰減率方法具有少打井，經(jīng)濟(jì)高效等優(yōu)點(diǎn)，但也存在一定的局限性。在一般情況下含水層長度l值的獲取需要一定的工作量，在實(shí)際應(yīng)用中需要根據(jù)情況權(quán)衡。