陽 玲，張 銘，郝永紅

（1.天津師范大學 數(shù)學科學學院，天津 300387；2.天津師范大學 天津市水資源與水環(huán)境重點實驗室，天津 300387；3.日本地質(zhì)調(diào)查局，筑波茨城 日本 305-8567）

隨著中國沿海城市經(jīng)濟和社會的快速發(fā)展，陸源污染物無論從種類上還是數(shù)量上都呈上升趨勢，海洋污染物總量的85%以上來自陸源污染物.研究人員對大西洋海灣地下水中鐳同位素的示蹤研究表明，濱海地區(qū)地下水海底排泄量相當于入海河流水量的40%，但其攜帶的污染物濃度卻是地表水濃度的幾十倍甚至上百倍[1-2].地下水攜帶的大量陸源污染物（化肥、農(nóng)藥和重金屬等）通過海底排泄進入濱海水域，造成近岸海域水質(zhì)惡化，引發(fā)赤潮和綠潮等一系列海洋生態(tài)災害和環(huán)境突發(fā)事件.在地下水海底排泄過程中，含水層的水力特性起到關(guān)鍵性作用，因此，識別含水層水力參數(shù)是實現(xiàn)濱海及近海海域生態(tài)環(huán)保的關(guān)鍵.壓力傳導系數(shù)是地下水系統(tǒng)中溶質(zhì)運移和擴散的重要參數(shù)之一，通常由抽水試驗確定，即通過抽水和注水對地下水系統(tǒng)產(chǎn)生擾動（刺激），利用地下水位的響應計算壓力傳導系數(shù).此方法的不足之處在于抽水試驗需要耗費大量的人力、財力和物力[3].海水潮汐作為一種自然刺激，同樣會引起含水層地下水位的波動.因此，本研究將地下水水位對潮汐的響應看作是一種天然抽水試驗，以期利用地下水水位對潮汐的響應識別含水層壓力傳導系數(shù)，提出一種經(jīng)濟實用的識別濱海地區(qū)含水層壓力傳導系數(shù)的新方法.

1 模型建立和解析解

潮汐由萬有引力和地-月系統(tǒng)繞其共同質(zhì)心旋轉(zhuǎn)的綜合作用造成，通常由多個分潮構(gòu)成，因此又稱多頻潮汐.潮汐的多頻信號在含水層中傳播，只有最強的信號可以穿過含水層引起地下水水位的變化，并且能夠被地下水水位壓力感應器所捕獲.同時，由于含水層的過濾作用，較弱的潮汐信號會在經(jīng)過含水層時衰減消失.在以往地下水與海水交互作用的研究中，人們主要關(guān)注信號最強的單頻潮汐，相關(guān)研究如表1 所示.從數(shù)學角度考慮，多頻譜波動潮汐可以表達為

式（1）中：λ和φ分別為經(jīng)緯度；Ai（λ，φ）為振幅；Gi（λ，φ）為相位；χi為初始相位.

研究表明，式（1）的前2 項即可描述98%的潮汐特征[4-5].由此可見，雙頻潮汐可以更加準確地描述海水的周期波動.為了更加全面、準確地反映地下水水位對潮汐的響應，本研究在解析過程中選擇2 個分潮.

表1 各種含水層以及邊界條件Tab.1 Various aquifers and boundary conditions

本研究考慮的濱海地下水系統(tǒng)包含上、下2 個弱透水層及其中間夾層，即含水層.假設每一層均均勻且水平，系統(tǒng)可以概化為一維地下水流動概念模型，示意圖如圖1 所示，其中x 軸的起點位于垂直海灘與不透水層底板的交叉點，x 軸向陸地為正方向.以平均海水水位（mean sea level，MSL）為基準，雙頻潮汐振幅分別為A1和A2.

圖1 無限寬海岸含水層示意圖Fig.1 Schematic diagram of an infinite wide coastal aquifer

在這些假設下，承壓含水層中一維地下水流動的控制方程為

式（2）中：h 為地下水水位；x 為水平坐標；t 為時間；T為含水層的滲透系數(shù)；S 為含水層儲水系數(shù).式（2）嚴格應用于承壓含水層，但當?shù)叵滤袷幈群畬雍穸刃r，式（2）可用于非承壓含水層[5，14-16].

潮汐邊界條件為

式（3）中：h（0，t）為 x=0 處的地下水水位；Ai為潮汐振幅（i=1，2）；hMSL為平均海水水位；ωi為潮汐振蕩的角速度（i=1，2）；ci為潮汐振蕩的初始相位（i=1，2）.

無窮遠處，無流量邊界條件為

式（4）中：h（∞，t）為內(nèi)陸離海岸無窮遠處的地下水水位.式（4）表明，海水潮汐信號在含水層的傳播過程中具有衰減效應，當距離海岸線無窮遠時，潮汐對地下水水位的影響為0[7].

依據(jù)文獻[7]，將式（1）～式（3）擴張到復數(shù)域上求解，結(jié)合解的疊加原理，獲得地下水水位對潮汐響應的解析解

2 數(shù)值模擬

2.1 地下水對潮汐的響應

為了模擬地下水水位對海洋潮汐的響應，根據(jù)實際情況，在參數(shù)的合理取值范圍內(nèi)，對式（5）中參數(shù)進行賦值，具體參數(shù)值如表2 所示.

表2 含水層參數(shù)Tab.2 Parameter of aquifer

基于表2 中的參數(shù)，得到地下水水位變化對海洋潮汐的響應情況.圖2 描述了地下水水位h 隨x 和t變化而變化的過程.由圖2 可以看出，地下水水位h是關(guān)于時間t 和距離x 的二維函數(shù).x 越?。拷０叮?，地下水水位對潮汐的響應越大，x 越大（遠離海岸），地下水水位對潮汐的響應越小，即隨著距離的增加，地下水振幅逐漸衰減，至無窮遠處穩(wěn)定于平均海水水位.

為了更加清楚地了解地下水水位對潮汐的響應過程，分別給出x=0、100 和200 m 時，地下水水位對潮汐的響應情況，結(jié)果如圖3 所示.由圖3 可以看出，當x 較大（x=200 m），即遠離海岸線的地方，地下水水位較低，波動較小.相反，當x 較?。▁=100 m），即靠近海岸線的地方，地下水水位較高，波動較大，反映出潮汐信號在含水層中的傳播與位置x 和時間t 有關(guān).

圖2 地下水對潮汐的響應Fig.2 Response of groundwater to tides

圖3 當x=0、100 和200 m 時，地下水水位隨時間的變化Fig.3 Variation of groundwater level with time when x= 0，100 and 200 m

2.2 用最小二乘法識別壓力傳導系數(shù)

本研究的主要目的是運用解析解和數(shù)值模擬相結(jié)合的方法，驗證利用潮汐信號和地下水水位響應識別含水層的壓力傳導系數(shù)的可行性.首先，在x=200 m處選取t∈[96，121]時間段的地下水水位的101 個離散點，時間步長為0.25 h.然后，為了模擬實際地下水水位數(shù)據(jù)中由于設備和人為因素會造成的誤差，在離散點上分別疊加 N（0，0.12）和 N（0，0.052）的高斯白噪聲，生成地下水水位的仿真觀測數(shù)據(jù)，即獲得hm=h（x=200，t+ Δt，D=20 000）+N（0，0.12）和 hm=h（x=200，t+ Δt，D=20 000）+N（0，0.052）共 2 組仿真觀測數(shù)據(jù)，如圖4 所示.最后，基于仿真觀測數(shù)據(jù)，運用最小二乘法，估計壓力傳導系數(shù)，誤差函數(shù)為

式（6）中：hc為計算數(shù)據(jù)；hm為仿真觀測數(shù)據(jù).

圖4 仿真數(shù)據(jù)Fig.4 Synthetic data

為了識別壓力傳導系數(shù)D，使用Mathematica 軟件的內(nèi)部函數(shù)Findfit 計算誤差函數(shù)最小化時的滯后時間和壓力傳導系數(shù)D.Findfit 的默認值為線性最小二乘的Singular value decomposition 和非線性最小二乘的Levenberg-Marquardt[17].FindFit 在內(nèi)部建立了一個殘差函數(shù)和雅克比，進而用于高斯-牛頓法以尋找平方和的最小值或非線性最小二乘擬合[18].

在具體計算過程中，首先，給出初始值Δt=1.555 h，D=19 000 m2/h，第 1 種高斯白噪聲 N（0，0.12）條件下，識別結(jié)果為 Δt=-0.017 h，D=19 410 m2/h；第 2 種高斯白噪聲 N（0，0.052）條件下，識別結(jié)果為 Δt=-0.001 h，D=20 175 m2/h.2 種噪聲的滯后時間分別為0.017 h和-0.001 h，壓力傳導系數(shù)D 分別為19 410 m2/h 和20 176 m2/h，具體誤差如表3 所示.

表3 壓力傳導系數(shù)的誤差Tab.3 Errors in hydraulic diffusivity

在實際工程應用中，壓力傳導系數(shù)的誤差應該低于 5%，本研究加入高斯白噪聲 N（0，0.12）和 N（0，0.052）后，所得結(jié)果的相對誤差分別為2.949%和0.878%，均小于3%，符合工程誤差標準，充分證明運用本研究方法識別壓力傳導系數(shù)D 是可行的.

2.3 潮汐振幅在含水層中的衰減

潮汐信號在含水層中傳播時，隨著傳播距離的增加，振幅逐漸衰減，且壓力傳導系數(shù)越大，衰減速度越小，潮汐信號傳播得越遠.因此，在實際操作中，需要依據(jù)含水層的透水特性確定地下水水位觀測井的位置.如果觀測井的位置超出潮汐信號所能到達的距離范圍，則觀測不到地下水水位對潮汐的響應.為此，有必要研究潮汐信號在不同壓力傳導系數(shù)的含水層中的衰減過程及其傳播距離，為野外觀測井定位提供依據(jù).

為了比較精確地確定壓力傳導系數(shù)，地下水水位對潮汐的響應不宜過小.基于實際情況，本研究設定觀測井位置處的地下水水位波動不得小于潮汐振幅的0.1 倍.即針對不同壓力傳導系數(shù)的含水層，計算地下水水位衰減至潮汐振幅的0.1 倍時的x 值.根據(jù)表2中的參數(shù)帶入解析解式（5）得到潮汐最大振幅為Am=3.198 m，在含水層壓力傳導系數(shù)為4 000～200 000 m2/h的范圍內(nèi)選取具有代表性的12 個數(shù)值，得到其衰減曲線如圖5 所示.由圖5 可知，壓力傳導系數(shù)D 越大，衰減距離x 越大，不同傳導系數(shù)情況下，潮汐信號衰減至0.1 Am時的衰減距離如表4 所示.

圖5 潮汐振幅信號在不同壓力傳導系數(shù)含水層中傳播時的衰減曲線Fig.5 Attenuation decaying of tide amplitude in aquifers with different hydraulic diffusivity

表4 在不同的壓力傳導系數(shù)D 的條件下地下水水位振幅衰減到0.1 Am 的距離Tab.4 Distance of groundwater amplitude decaying to 0.1 Am under different hydraulic diffusivity D

根據(jù)實際情況，壓力傳導系數(shù)D 的最大值約為2×105m2/h.由表 4 可以看出，當 D=2×105m2/h 時，振幅衰減到0.1 Am時x 值為2 590.4 m，不同的壓力傳導系數(shù)對應的x 值不同.在實際應用中，為了得到更準確的數(shù)據(jù)，應依據(jù)表4 中提供的位置鉆取觀測井.

3 討論

本研究利用數(shù)值模擬的方法基于地下水水位對潮汐響應的解析解得到相關(guān)結(jié)果，而實際問題更為復雜，本研究尚存在諸多未考慮因素，如海灘是有斜坡的，含水層底部與水平方向也往往存在一定傾斜角，但本研究模型中假設海灘和含水層均為水平，且一般非承壓含水層可能通過滲漏影響地下蓄水層的潮汐漲落，即實際問題中可能還含有滲漏項.這些因素都會影響地下水對潮汐的響應，也會影響在含水層中振幅的衰減.假設海灘斜坡為β，含水層底部與水平方向夾角為 θ，則當 β=90°，θ=0°時，就退化成為本研究中模型.后續(xù)工作可以考慮滲漏和斜坡混合問題，建立新的模型，用該方法來識別壓力傳導系數(shù)D，且可以討論不同的D 與振幅衰減的關(guān)系，從而更好地確定打井位置.

4 結(jié)論

本研究建立了濱海地區(qū)地下水水位對多頻潮汐的響應模型，將模型擴展到復數(shù)域，并結(jié)合解的疊加原理求得解析解，模擬了地下水水位對潮汐的響應過程，并在此基礎上，通過數(shù)值實驗采用最小二乘法識別了含水層的壓力傳導系數(shù)，得到以下結(jié)論：

（1）由海水潮汐引起的地下水水位波動可以作為天然的抽水試驗確定含水層水力傳導系數(shù).

（2）通過數(shù)值仿真，采用最小二乘擬合函數(shù)，借助Mathematica 軟件識別出壓力傳導系數(shù)D 和滯后時間，相對誤差小于3%，滿足工程地質(zhì)要求，也說明該方法可以準確識別壓力傳導系數(shù)D.

（3）對壓力傳導系數(shù)D 和距離x 進行敏感性分析，把0.1 倍最大振幅作為一個標準，一般打井位置應小于振幅衰減至0.1 倍時對應的x.

（4）本研究具有一定的局限性，如沒有考慮實際問題中普遍存在滲漏和斜坡，但對研究有滲漏和斜坡的問題具有參考價值.