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

(1.中國人民解放軍92250部隊(duì)，浙江舟山 316041;2.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實(shí)驗(yàn)室和大規(guī)模復(fù)雜系統(tǒng)數(shù)值模擬江蘇省重點(diǎn)實(shí)驗(yàn)室，江蘇南京 210023;3.河海大學(xué)港口海岸與近海工程學(xué)院，江蘇南京 210098)

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

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

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

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

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

1 控制方程

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

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

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

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

為考慮淺水分潮對地下水波動(dòng)的影響，假設(shè)外海潮汐過程為

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

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

2 二階攝動(dòng)解

以ε為攝動(dòng)參數(shù)，可設(shè)式(6)的攝動(dòng)解:

應(yīng)用Li［15］相同的方法，可解得二階攝動(dòng)解為

顯然，當(dāng)r=0時(shí)，上述二階解可退化為Li［15］的解。

3 解的驗(yàn)證

3.1 現(xiàn)場觀測

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

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

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

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

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

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

3.2 潮汐邊界條件

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

由于海岸實(shí)測潮位的觀測時(shí)間短，且有效時(shí)段更短，不宜應(yīng)用常用的調(diào)和分析方法分離各分潮，為此，假定海岸潮位為

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

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

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

圖4 實(shí)測水位與預(yù)測水位擬合ig.4 Tide level fitting result between observation and prediction

3.3 結(jié)果比較

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

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

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

4 分析討論

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

4.1 超高

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

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

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

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

4.2 振幅

圖7為兩種不同情況下二階攝動(dòng)解、數(shù)值解與1號測井的實(shí)測值過程線的比較圖。

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

4.3 相位

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

5 結(jié)語

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

［1］ 施雅風(fēng).我國海岸帶災(zāi)害的加劇發(fā)展及其防御方略［J］.自然災(zāi)害學(xué)報(bào)，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］.水文地質(zhì)工程地質(zhì)，1992，19(6):23-33.

［6］ 蔡祖煌，馬鳳山.海水入侵的基本理論及其在入侵發(fā)展預(yù)測中的應(yīng)用［J］.中國地質(zhì)災(zāi)害與防治學(xué)報(bào)，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］.地學(xué)前緣，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.