淺水分潮對海岸潮致地下水波動特征的影響

樓悅安，宋志堯，孔 俊，程 晨，徐福敏

(1.中國人民解放軍92250部隊，浙江舟山 316041;2.南京師范大學虛擬地理環境教育部重點實驗室和大規模復雜系統數值模擬江蘇省重點實驗室，江蘇南京 210023;3.河海大學港口海岸與近海工程學院，江蘇南京 210098)

海岸帶是大陸和海洋的過渡地帶，受到大氣圈、水圈、巖石圈和生物圈等的共同作用，豐富的自然資源和良好的環境條件使之成為人口最為稠密、經濟最為發達的地區，在自然與人為災害因素的作用下，海岸帶已逐漸成為各種環境災害頻繁發生、疊合發展的地帶，如洪澇、風暴潮、海岸侵蝕、海水入侵、污染等，是典型的生態脆弱帶［1］。

我國擁有漫長的海岸線，隨著沿海地區社會經濟的不斷發展，使得因全球變暖而引發的海岸環境災害有進一步加劇的趨勢。已有的研究成果表明，海岸侵蝕、海水入侵及海岸污染等環境災害幾乎都與地下水運動有關，海岸潮致地下水波動的研究已日益為國內外學者所關注［2-8］，并形成了海灘地下水動力學的研究方向［9］。

Jacob［10］基于垂直灘面、順直岸線、單個承壓含水層中的一維流動，忽略非線性效應、毛細效應、垂向流等影響，海潮邊界條件為具有一定振幅的余弦函數的假設，與1950年首次提出了描述潮波在承壓潛水層傳播的最簡解析解。由于Jacob解基于上述的嚴格假設，因此灘面的坡度、岸線的形態特征及組成波的非線性作用等因素都可能造成Jacob解的不適用［11］。Dagan［12］針對垂直灘面，基于非線性地下水控制方程，采用攝動法得到了小振幅波情況下的地下水波動一階攝動解。Parlange［13］在此基礎上進一步推出了二階攝動解，并指出方程的非線性項在傳播過程中產生了高頻波動。Nielson［14］考慮了斜坡的影響，推導出了相應的攝動解，并與現場實測數據進行了比較分析。其后Li［15］引入等效邊界條件，采用攝動法得到了二階攝動解。但由于其采用的攝動參數與岸灘坡度有關，而攝動參數又要求小于1，所以該解不適用于緩坡情況。Teo等［16］綜合考慮方程的非線性和斜坡的影響，運用雙攝動參數，獲得了更高階的攝動解。理論上，攝動參數避免了坡度的影響，但由于長期項的存在，解的精度實際也只能夠保留至二階。對此，Kong et al［17］運用新的雙攝動參數，提出基于待定系數的新解析方法，獲得了斜坡海岸情況下非線性方程的高階攝動解，避免了長期項的出現，使攝動解在理論上可推至更高階。

上述攝動解的推導往往基于單個分潮的假設。Li［15］曾通過兩個半日潮研究過大小潮對海岸潛水層地下水波動的影響，指出大小潮引起的地下水波動比由單個半日潮引起的地下水波動向內陸傳播更遠的結論。實際上，海岸潮汐是十分復雜的，往往由多個分潮共同作用，其中淺水分潮作用顯著，所以考慮分潮及其淺水分潮綜合影響下的攝動解具有實際意義，但迄今尚無學者研究淺水潮波非線性作用產生的倍潮波即淺水分潮的影響。

為此，主要研究半日潮產生的淺水分潮對海岸潛水層地下水波動的影響，應用在處理地下水問題時常用的攝動法得到了包含半日潮M2和淺水分潮M4綜合影響下的地下水波動二階攝動解，并與現場實測數據和數值解進行了比較，分析了考慮淺水分潮后海岸潮致地下水波動特征的變化。

1 控制方程

針對一維海岸地下水運動，考慮半無限潛水層，底板水平不透水，不考慮毛細管壓力對地下水波動的影響，介質單一均質，忽略垂向流動，水位波動幅度相比于潛水層厚度較小(見圖1)，得到一維線性Boussinesq方程:

式中:h*(x*，t*)為相對于平均海平面的水位高度(m)，x*為水平方向坐標(m)，t*為時間(s)，z*為垂向坐標(m)，D為平均含水層厚度(m)，K為介質的滲透系數(m/s)，ne為有效孔隙度。

圖1 海岸潛水層潮致水位波動示意Fig.1 Sketch map of tide-induced groundwater fluctuations in coastal aquifer

在外海潮汐作用影響下，x*→∞時向陸一側水位梯度為0，即

為考慮淺水分潮對地下水波動的影響，假設外海潮汐過程為

其中，A1，A2分別為分潮M2，M4的振幅(m)，ω為分潮M2的角速度(rad/h)，δ為分潮M2與其淺水分潮M4之間的相位差(rad)。

引入新變量z*=x*-(t*)，將動邊界問題轉化為固定邊界問題，代入式(1)，并記r=A2/A1，得到的方程可寫為以下形式:

2 二階攝動解

以ε為攝動參數，可設式(6)的攝動解:

應用Li［15］相同的方法，可解得二階攝動解為

顯然，當r=0時，上述二階解可退化為Li［15］的解。

3 解的驗證

3.1 現場觀測

2010年7月1日至7月15日在蘇北弶港附近的潮灘區進行了一次現場觀測，潮灘剖面及測井布置如圖2所示。

圖2 測井布置Fig.2 Sketch map of arranged observation wells

弶港位于120°51'E，32°44'N，該海域潮汐變化顯著，外海屬正規半日潮，潮差最大超9 m，而近岸水深變淺，潮差減小，但淺水分潮作用顯著，主要是半日潮及其淺水分潮;同時海灘地質組成較單一，斜坡平緩，潛水層較淺，非常適合研究半日潮及其淺水分潮對地下水運動的影響。觀測儀器主要是測量水位的回聲探測器儀，測程0～12 m，精度0.1 cm，可自動記錄。根據閘外航道潮汐特征，測井長度為4 m，灘面最深處入土2.5 m。

各測井的地下水水位測量從7月1日上午08∶30開始，間隔為10 min;潮汐水位從7月6日15∶00開始觀測，間隔為1 h。

由圖3可以看出，各測井地下水位隨著潮汐水位的變化而變化。各測井地下水位的波動過程線具有非對稱性，水位上漲過程要明顯陡于下落過程，Knight［18］指出這是由于在漲潮過程中孔隙的透水性要明顯高于落潮造成的。各測井地下水位的波動振幅沿程減小，這是由于傳播過程中的能量耗散造成的;同時沿程各測井出現了相位滯后的現象，這是因為波動傳播過程阻尼作用導致的。

圖3 潮汐水位及各測井地下水位變化過程(基面為85高程)Fig.3 Time series of groundwater fluctuations in each well and tide level(85 datum)

3.2 潮汐邊界條件

為了充分研究淺水分潮對潮致地下水波動的影響，利用上述現場觀測數據中的某一天(2010年7月9日08:00至次日08:00)進行研究分析。

由于海岸實測潮位的觀測時間短，且有效時段更短，不宜應用常用的調和分析方法分離各分潮，為此，假定海岸潮位為

式中:h0為平均海平面高度(m);hM2，hM4分別為分潮M2及其淺水分潮M4的振幅(m);φ1，φ2分別為分潮M2，M4的初相位(rad)。

對表達式(7)，應用簡單實用的最小二乘法進行擬合，最終可得到h0=1.05 m，hM2=1.33 m，hM4=0.49 m，φ1=1.86 rad，φ2=3.30 rad。所以由分潮M2及其淺水分潮M4引起的閘外潮汐水位變化表達式為

圖4為預測水位與實測水位的擬合情況。由圖可以看出，兩者擬合良好，可以作為這里研究的海岸潮汐邊界條件。

圖4 實測水位與預測水位擬合ig.4 Tide level fitting result between observation and prediction

3.3 結果比較

選取在同一坡面上的1、2、4、5號測井的實測數據與攝動解及數值解進行比較分析。相關參數分別為:A1=1.33 m，A2=0.49 m，K=29.0 m/d，ne=0.5，tan(β)=0.20，D=5.0 m，δ= φ2-2φ1=-0.42 rad。

由圖5比較攝動解、數值解和實測值，可以看出:1)對于某一測井來說，波動振幅隨著攝動解階數的增加與實測值愈接近，由于高階攝動解越趨復雜，一般只能推導至二階解，如果推至三階甚至更高階的解，攝動解與實測值應該越來越接近。當然，在水位谷值處，攝動解與實測值之間還可能存在較大的偏差，這主要是由于坡面溢出點的形成造成的，此時作為潮致地下水運動控制方程的邊界條件已發生了根本的變化，漲潮時臨海一側地下水頭為潮汐水位，而落潮時由于潛水層內落潮較慢，溢出點與潮汐水位不再同步(如圖1)。同時數值解相比二階解具有更高的精度，與實測值也更為接近;2)1號測井由于最靠近外海，所以由攝動解反映的地下水位波動過程線與實測地下水波動的相位基本一致，而2、4、5號測井由攝動解反映的地下水波動相比實測地下水位過程線會整體向右略有偏移，這主要是由于攝動解是在十分概化的條件下獲得的，如斜坡坡度是規則的、潛水層是等厚度的、介質參數(滲透系數和孔隙度等)是均勻的等等，又沒有考慮溢出面、垂直流、毛細作用及組成波的相互作用等的影響，造成計算的波長要比實際波長偏大。這一現象與Nielson［14］的研究結果基本一致。而數值解較二階攝動解在相位偏移上有一定程度的改善。綜上所述，各階攝動解基本上能反映海岸潮致地下水波動的趨勢性特征，攝動解的階數越高，其所反映的波動特征越充分。

圖5 各測井各階攝動解與數值解及實測過程線比較(相對于平均海平面)Fig.5 Time series comparison of perturbation solutions，numerical results and observed data in each well(relative to averaged sea level)

需要強調的是，即使采用實測的邊界潮位進行數值模擬，計算值與實測值在潮致地下水波下落過程的后期仍存在一定的誤差，表明目前常用的Boussinesq方程，是基于潛水層等深和介質參數(滲透系數和孔隙度等)是均勻的等假定條件得到的，又沒有考慮溢出面、垂直流、毛細作用等的影響，所以只能反映海岸潮致地下水波動的基本特征，如超高、振幅和相位等。由于這里只是討論淺水分潮對這些基本特征的影響，而且這種影響也是相對的，因而基此進行分析研究目前仍是可行的。

4 分析討論

為了研究淺水分潮對海岸地下水波動基本特征的影響，將從超高、振幅、相位等潮致地下水波動特征進行比較分析。

4.1 超高

在潮汐作用下，海岸潛水層內的平均地下水位會產生相對于平均海平面的雍高。當沙灘傾角較小時，這種由于潮汐作用而產生的超高非常明顯。同時地下水超高也是海岸地下水波動非線性作用的一個基本特征。Philip［19］和Smiles［20］曾基于一維非線性Boussinesq方程計算了直墻海岸內的地下水超高;Parlange在其二階攝動解和實驗的基礎上對地下水位超高進行了分析，證實了潮致地下水水位超高的存在。根據之前獲得的二階攝動解，在周期平均后，忽略其中的周期項，得到地下水超高沿程二階解析解:

相應不考慮淺水分潮的解即Li的解反映的地下水超高沿程變化:

圖6 超高沿程分布比較Fig.6 Comparison of water table over-heights along inland direction

1、2號測井由于會被海水淹沒，所以其分析計算得到的超高值并不能真實反映其本質特征。選取離外海較遠的4、5號測井的實測超高與解析解及數值解進行比較(見圖6)。考慮淺水分潮后，解析解反映的地下水位超高有了明顯的提升，但是與實測地下水位超高還是存在一定的差值。一方面是由于忽略了方程自身的非線性作用;另一方面是由于在水位谷值處攝動解與實測值之間存在一定的偏差。而數值解反映的地下水位超高相比解析解要更接近實測值。由此可知，如果忽略淺水分潮的影響將會極大地低估地下水位超高值。

4.2 振幅

圖7為兩種不同情況下二階攝動解、數值解與1號測井的實測值過程線的比較圖。

圖7 不同情況下二階攝動解、數值解與實測過程線比較Fig.7 Time series comparison of second-order perturbation solutions，numerical results and observed data in No.1 well

4.3 相位

同樣由圖7可以看出，僅考慮半日潮M2的攝動解反映的地下水位波動曲線較之實測地下水位波動曲線存在向右偏移的趨勢，而考慮淺水分潮M4后的攝動解反映的地下水位波動曲線與實測地下水位波動曲線的相位基本一致。這表明，考慮淺水分潮后可以在一定程度上改善攝動解與實測值之間的相位偏離。

5 結語

針對淺水分潮M4對海岸潮致地下水波動特征的影響進行了研究。通過對考慮淺水分潮M4后的二階攝動解與實測值的比較可以看出，兩者之間的峰值較為吻合，但因目前研究所應用的控制方程即Boussinesq方程忽略了溢出面、垂直流、毛細作用等因素，包括數值解在內計算的谷值誤差仍較大;同時對某一測井來說，攝動解階數越高，越能充分反映地下水波動特征，相對而言數值解能較好地反映實測的潮致地下水波動特征。概括地說，淺水分潮對潮致地下水波動的振幅及其超高值有較顯著的影響，由于斜坡的非線性作用，這種影響明顯大于海岸邊界淺水分潮與半日潮的振幅比;同時可較好地改善相位上出現的偏移，表明在研究海岸潮致地下水波動特征時淺水分潮作用是不能忽略的。

［1］ 施雅風.我國海岸帶災害的加劇發展及其防御方略［J］.自然災害學報，1994，3(2):3-14.

［2］ Chappell J，Eliot I G，Bradshaw M P.Experimental control of beach face dynamics by water table pumping［J］.Eng.Geol.，1979，14:29-41.

［3］ Fragoso G，Sepencer T.Physiographic control on the development of spartina marshes［J］.Science 322，2008，1064.

［4］ Pinder G F，Cooper H H.A numerical technique for calculating the transient position of the saltwater front［J］.Water Resources Research，1970，6(4):875-882.

［5］ 薛禹群，謝春梅，吳吉春.海水入侵研究［J］.水文地質工程地質，1992，19(6):23-33.

［6］ 蔡祖煌，馬鳳山.海水入侵的基本理論及其在入侵發展預測中的應用［J］.中國地質災害與防治學報，1996，7(3):1-9.

［7］ Moore W S.Large groundwater inputs to coastal waters revealed by 226Ra enrichment［J］.Nature，1996，380(18):612-614.

［8］ 賈國東，黃國倫.海底地下水排放:重要的海岸帶陸海相互作用過程［J］.地學前緣，2005，12(特刊):29-35.

［9］ Horn D P.Beach groundwater dynamics［J］.Geomorphology，2002，48:121-146.

［10］ Jacob C E.Flow of Groundwater in Engineering Hydraulics［M］.New York:［s.n.］，1950:321-386.

［11］ Trefry M G，Johnston C D.Pumping test analysis for a tidally forced aquifer［J］.Ground Water，1998，36(3):427-433.

［12］ Dagan D.Second-order theory of shallow free-surface flow in porous media［J］.The Quarterly Journal of Mechanics and Applied Mathematics，1967，20(4):517-526.

［13］ Palange J Y，Stagnitti F，Starr J L，et al.Free-surface flow in porous media and periodic solution of the shallow-flow approximation［J］.Journal of Hydrology，1984，70(1-4):251-263.

［14］ Nielsen P.Tidal dynamics of the water table in beaches［J］.Water Resources Research，1990，26(9):2127-2134.

［15］ Li L，Barry D A，Stagnitti F，et al.Beach water table fluctuations due to spring-neap tides:moving boundary effects［J］.Adv Water Resour，2000，23(8):817-824.

［16］ Teo H T，Jeng D S，Seymour B R，et al.A new analytical solution for water table fluctuations in coastal aquifers with sloping beaches［J］.Adv Water Resour，2003，26(12):1239-1247.

［17］ Kong J，Song Z Y，Xin P，et al.A new analytical solution for tide-induced groundwater fluctuations in an unconfined aquifer with a sloping beach［J］.China Ocean Eng，2011，25(3):479-494.

［18］ Knight J H.Steady period flow through a rectangular dam［J］.Water Resources Research，1981，4(17):1222-1224.

［19］ Philip J R.Periodic nonlinear diffusion:An integral relation and its physical consequences［J］.Aust J Phys，1973，26:513-519.

［20］Smiles D E，Stokes A N.Periodic solutions of a nonlinear diffusion equation used in groundwater flow theory:Examination using a Hele-Shaw model［J］.Journal of Hydrology，1976，31(1-2):27-35.