陽 玲,杜金月,王同科,趙志學,郝永紅
(1.天津師范大學數學科學學院,天津 300387;2.天津師范大學天津市水資源與水環境重點實驗室,天津 300387)
世界上有40%的人口居住在離海岸線100 km以內的海濱城市,給近海海域帶來了大量污染,破壞了生態環境[1-3]。源于陸地的海洋污染有44%是來自于地下水海底排泄和地表徑流[1]。雖然地下水海底排泄的流量是地表河流入海徑流量的40%[2],但是它的污染物濃度卻是地表徑流的10倍[3]。海底地下水排放被視為海洋污染最主要的方式之一[4-5]。預測海底地下水排放和輸送到海洋中的污染物對海洋生態和環境保護有著重要意義。
水文地質參數識別又稱水文地質逆問題求解,是建立地下水數值模型的關鍵步驟[6]。含水層的壓力傳導系數作為描述地下水溶質運移的關鍵參數被廣泛的研究[7-10]。傳統識別壓力傳導系數最廣泛的技術是抽水試驗和振蕩試驗[11-13],但是這些方法有一個共同的缺點:需要耗費大量的人力、財力和物力。為了更好地描述含水層的特性,有必要開發一種經濟高效的方法來估計含水層的壓力傳導系數。
海洋和河流的潮汐可以看作自然的振蕩抽水試驗,可通過分析潮汐和地下水響應關系估計壓力傳導系數。Ferris估計了潮汐河含水層的壓力傳導系數,他使用了時間滯后和潮汐與地下水位相關的潮汐效率因素的方法[14]。Erskine估計了非承壓和高度滲透沿海含水層的壓力傳導系數,該含水層位于英國的一個核電站附近,他用了時間滯后和潮汐效率因素方法[15]。Zhou利用觀測井觀測地下水水位振幅和潮汐確定了T/S(給水度與導水系數之比),寫出了T/S的表達式[16]。在貝爾港沿海巖溶含水層的背景下,類似的壓力傳導系數的估計表明,水動力環境包括基質、裂隙和管道流動系統。
然而,前人運用地下水水位對潮汐的響應識別含水層的參數多采用數值擬合的方法。在前人的研究基礎上,本文試圖獲得潮汐產生的壓力波在含水層中傳播的解析解,通過數值擬合估計含水層的壓力傳導系數。其能夠更加清楚地了解潮汐影響地下水的機理。
假設海濱區有一個水平、均質、等厚的含水層,它位于兩弱透水層之間(圖1)。因為地下水的壓力水頭高于含水層和海水的高度,所以地下水向海中排泄。
基于圖1的概念模型,建立笛卡爾直角坐標系,x軸的原點位于內陸,距海岸的水平距離為l m,且地下水水位是定水位H1。x軸是水平的,沿著含水層向海岸方向延伸,垂直方向為海岸線方向。在實際條件中左邊界為定水位邊界的可能性不大,更多的是流量邊界。假設有一條與海岸平行的運河,保證左邊界為定水位邊界。該模型需假設有一條與海岸平行的運河,且該運河與海不相通。

圖1 濱海區含水層概念模型
承壓含水層中地下水運動的一維偏微分方程為:
(1)
式中:H——水頭/m;
x——水平坐標/m;
t——時間/h;
T——含水層的滲透系數/(m·h-1);
S——含水層儲水系數;
l——含水層的長度/m。

(2)
邊界條件為:
H|x=0=H1,H|x=l=H2+Asinωt,t>0,
(3)
式中:H1——內陸含水層定水頭邊界/m;
H2——海水平均水位/m;
A——海水振蕩振幅/m;
ω——角速度/h-1。
將式(2)~(3)轉化到復數域上[17-18],即將有關物理量從實平面擴充到復平面[19]。變量x和t是相互獨立的。因此可以用分離變量法求解:
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)
因為式(7)關于時間是周期的,所以可以省略初始條件。把式 (4)~ (6)擴張到復數域,假設u2(x,t)=lm(U(x,t)),代入式(6)~(7),可以得到:
(8)
U(0,t)=0,U(l,t)=Aeiωt,t>0
(9)


(10)


(11)
最終求解結果:






(12)
式(12)可以簡化為:
(13)
其中:



(14)



(15)

(16)
(17)

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

表1 含水層參數

圖2 地下水對潮汐的響應
為更加清晰地了解不同位置地下水響應情況,圖3表達了x=30 m和x=80 m以及海水水位三種不同情況下地下水水位隨時間的變化。

圖3 x=30 m和x=80 m以及海水水位三種不同情況下地下水位隨時間的變化
為檢驗潮汐衰減率r(x,t,ω)與余切函數cotωt之間的線性關系。當x=80 m時地下水水位,r(x,t,ω)與cotωt呈線性相關(圖4)。

圖4 x=80 m時r(x,t,ω)與cotωt的線性關系
為了驗證潮汐信號衰減率方法的正確性,模擬野外地下水水位的觀測過程,對x=80 m處的地下水數據進行取樣,取樣時間步長為15 min,共獲得48個采樣點(圖5)。
由圖6可見,r(80,t,ω)與cotωt呈線性關系。由于采樣步長15 min,相對于樣本總體時間較小,圖6中部分點重合。

圖5 在x=80 m處,時間間隔為0.25 h,地下水對潮汐響應的合成數據

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

圖7 用α和β估計初值
(1)在濱海含水層中,地下水水位與海洋潮汐的交互作用類似于自然的振蕩抽水試驗。當漲潮時,相當于給地下水注水,地下水水位上升。當退潮時,相當于抽取地下水,地下水水位下降?;诔毕盘栐诤畬又袀鞑ズ退p的過程,建立地下水模型求得解析解,提出用潮汐信號衰減率方法估計含水層的特性。潮汐信號衰減率r(x,t,ω)與海水振蕩的余切函數cotωt線性相關。
(2)線性關系中的斜率和截距都是壓力傳導系數的函數。提出了利用潮汐信號衰減率與余切函數的關系估計壓力傳導系數的方法(即用斜率和截距估計壓力傳導系數)。數值仿真結果表明,該方法可以準確地估算含水層的壓力傳導系數D。
(3)潮汐衰減率方法具有少打井,經濟高效等優點,但也存在一定的局限性。在一般情況下含水層長度l值的獲取需要一定的工作量,在實際應用中需要根據情況權衡。