程博超，羅照陽，嚴佳偉，孔 俊，2

(1.河海大學港口海岸與近海工程學院，江蘇南京 210098;

2.河海大學水文水資源及水利工程科學國家重點實驗室，江蘇南京 210098)

目前，沿海地下水運動問題研究已成為國內(nèi)外的熱點課題之一［1-2］。傳統(tǒng)Boussinesq方程的理論研究［3-5］均假定地下水自由面是飽和含水層和不含水層間的突變界面，實際上這一突變界面并不存在，而是存在一過渡區(qū)(又稱毛細效應(yīng)區(qū))［6］。Millham等［7］曾經(jīng)指出，如果不考慮非飽和流的存在，會嚴重錯估地下水頭的時空特征。為了進一步分析2種不同區(qū)域的相互影響，Parlange等［8］引入Green等的假定，在方程中添加了高階項，提出考慮毛細效應(yīng)的Boussinesq方程。Barry等［9］研究表明，在一般情況下對于長周期波(如潮波)，緊鄰飽和區(qū)上層的毛細效應(yīng)作用不明顯，但對高頻波(波浪)來說，毛細效應(yīng)的誘發(fā)作用較明顯，在非飽和區(qū)作用下高頻波向潛水層傳播時，向陸一側(cè)衰減速度會有所減慢，位相偏移減小，傳播距離更遠。Li等［10-11］通過對高頻波動的數(shù)值模擬證實了這一現(xiàn)象的存在并將修正項引入中等水深Boussinesq方程，擴展了其在描述高頻波傳播方面的適用范圍。但上述研究考慮的毛細效應(yīng)區(qū)僅僅是位于非飽和區(qū)與飽和區(qū)交界處較窄的一段飽和負壓區(qū)，不能完全代表非飽和區(qū)的效用。為此，Kong等［12］針對上述不足，提出了一個在有限非飽和厚度區(qū)域內(nèi)、綜合考慮水平和垂向流的新型地下水波動控制方程。本文針對該方程通過攝動解析法剖析波動組成和波動的函數(shù)表達，進一步分析波動特征并對其影響因素展開分析。

1 控制方程

Kong等［12］基于Richards方程［13］，通過靜壓法和動壓校正法推求的地下水Boussinesq方程可表述為

其中

邊界條件為

式中:ne——土壤的有效孔隙率;h——地下水水頭高度;Ks——土壤的滲透系數(shù);α——與土壤非飽和特性成反比的參數(shù)，反映土壤的含水特性;Z0——土壤層頂面高度;D——平均海平面以下的潛水層厚度;A——地下水波動振幅;t——時間;ω——海潮波動頻率。

2 解析解

針對上述方程，Kong等［12］通過小參數(shù)攝動法推求出方程的解析解，得到1階主頻波和2階倍頻波(彌散波)的解析表達形式，考慮到1階主頻波和2階彌散波可基本反映潮波傳入潛水層后的波動特征，本文重點分析這2個頻率的波動特征和影響因素。1階主頻波表述為

2階彌散波的解析表達式為

其中

3 波動特征及影響因素分析

由于土壤的非飽和效應(yīng)，潮波自外海傳入潛水層后運動特征發(fā)生了明顯的改變，在傳播過程中，波的振幅出現(xiàn)了一定的衰減，同時產(chǎn)生了高頻波。本文基于上述1階主頻波和2階解析表達形式，分析各階組成波振幅沿程變化特征及其主要影響因素。

分別將α、Z0、Ks、ne作為研究變量，并對1階主頻波、2階彌散波波動方程(式(2)、式(3))中各常量參數(shù)賦值，得到5種取值方案，進而利用數(shù)值模擬方法繪出各取值情況下主頻波與彌散波的波幅曲線，從而分析各參數(shù)變量對波幅的影響機制。以下是方案的取值情況。方案1:Z0=5.5 m，Ks=4×10-5m/s，T=12 h，D=5 m，ne=0.3，α 為變量;方案2:Ks=4 ×10-5m/s，T=12 h，D=5 m，ne=0.3，α =0.2，Z0為變量;方案 3:Ks=4 ×10-5m/s，T=12 h，D=10 m，ne=0.1，α =0.2，Z0為變量;方案4:Z0=5.5 m，Ks=4 ×10-5m/s，T=12 h，D=5 m，α =0.2，ne為變量;方案5:Z0=5.5 m，T=12 h，D=5 m，ne=0.3，α =0.2，Ks為變量。

3.1 α對波幅的影響

α是反映土壤非飽和持水特征的參數(shù)，其與土壤的非飽和程度成反比，其值越小，土壤的非飽和持水性越強。進一步結(jié)合已推求的解析解，對比不同α取值情況下各階組成波振幅的沿程變化趨勢，本文主要針對1階主頻波和2階彌散波進行討論。對于1階主頻波，沿程振幅呈指數(shù)型衰減。當α增大時，土壤的非飽和持水特性減弱，非飽和區(qū)持水能力降低，1階主頻波振幅的衰減速度加快，見圖1(a)(其中A1、A2分別表示1階主頻波與2階彌散波的波幅，下同)。對于2階彌散波，沿程的振幅表現(xiàn)為先遞增、后減小，反映了能量從低頻波向高頻波的遷移過程。當2階彌散波從1階主頻波獲得的能量大于其本身向3階高頻波的輸出時，表現(xiàn)為其振幅不斷增加;反之則表現(xiàn)為振幅不斷減小。從圖1(b)可以看出，當α取值越大時，2階彌散波所能達到的最大振幅越臨近海邊界，同時振幅也越小。這是由于非飽和區(qū)持水能力降低后增加了波動傳播的阻尼效應(yīng)，限制了2階彌散波的成長;當α減小時這種趨勢恰好相反，2階彌散波成長距離變長，最大波幅增大，表明土壤的非飽和效應(yīng)有利于讓能量更多地從低頻波向高頻波轉(zhuǎn)移，使高頻波的毛細效應(yīng)更加顯著。

圖1 α對波幅的影響(方案1)Fig.1 Impact ofα on amplitude in project 1

3.2 Z0對波幅的影響

Z0越大表明平均水面線高度以上的土壤層越厚。當保持外海潮汐振幅不變時，通過改變Z0，可分析不同非飽和區(qū)厚度影響下1階主頻波和2階彌散波的振幅沿程變化趨勢。從圖2(a)、圖3(a)可見，1階主頻波的振幅呈指數(shù)型衰減特征。當非飽和區(qū)厚度增大后，1階主頻波振幅沿程衰減速度變快。分析其原因在于隨著潛水面上層土壤層厚度的增大，非飽和區(qū)持水量也逐漸增大，上層非飽和區(qū)與下層飽和區(qū)之間垂向水量交換作用加強，由于土壤含水量與氣壓之間的關(guān)系，增厚的非飽和區(qū)需要獲取更多的水分得以維持，這就迫使1階主頻波在傳播過程中不僅需要在水平方向輸送水分，同時在垂向上還要不斷補給非飽和區(qū)，因而就促進了1階主頻波的能量衰減。此外從圖2(b)、圖3(b)可以看出，非飽和區(qū)厚度增大后，2階彌散波的成長受限，其所能達到的最大振幅也有所減小。這表明，在保持土壤持水特性不變的情況下，非飽和區(qū)厚度加大后，由于1階主頻波波能衰減速度加快，部分能量用于維持非飽和區(qū)內(nèi)水體的運動，限制了其本身能量向高頻波的遷移，導(dǎo)致高頻波的振幅呈現(xiàn)出減弱的趨勢。

圖2 Z0對波幅的影響(方案2)Fig.2 Impact of Z0 on amplitude in project 2

圖3 Z0對波幅的影響(方案3)Fig.3 Impact of Z0 on amplitude in project 3

3.3 n e對波幅的影響

從圖4可見，考慮非飽和區(qū)存在時，隨著ne增大，1階主頻波的振幅衰減速度增快，1階主頻波能量轉(zhuǎn)移到2階彌散波的速度加快，則2階彌散波達到最大的成長距離變短，所能達到的最大振幅減小。

圖4 參數(shù)n e對波幅的影響(方案4)Fig.4 Impact of n e on amplitude in project 4

3.4 Ks對波幅的影響

Ks越大，土壤的滲透性越強。從圖5可見，當考慮土壤非飽和效應(yīng)時，隨著Ks增大，1階主頻波的振幅沿程衰減速率減慢。對于2階彌散波，隨著Ks增大，能量從1階主頻波轉(zhuǎn)移到2階彌散波的速度減慢，2階彌散波所能達到的最大振幅也變大。

圖5 參數(shù)Ks對波幅的影響(方案5)Fig.5 Impact of Ks on amplitude in project 5

4 結(jié) 論

a.潛水層內(nèi)1階主頻波的振幅沿程呈指數(shù)形式衰減;2階彌散波的振幅沿程變化特征為先增大、后減小，體現(xiàn)出能量從低頻波向高頻波遷移的過程。

b.各因素對比表明，非飽和區(qū)土壤持水能力越強，土壤滲透系數(shù)越大，越有利于高頻波的產(chǎn)生和傳播。

c.考慮到土壤層的非飽和持水效應(yīng)，2階彌散波所能達到的最大振幅比不考慮土壤非飽和效應(yīng)時大，說明考慮土壤非飽和效應(yīng)時2階彌散波所獲得的能量較多。因此土壤的非飽和效應(yīng)對于誘發(fā)非飽和區(qū)含水層內(nèi)的高頻波起到促進作用，而傳統(tǒng)的地下水波動方程忽略了這一物理現(xiàn)象，導(dǎo)致對高階倍頻波的形成認識存在不足。

d.Kong等［12］提出的新型地下水Boussinesq方程中，關(guān)于土壤含水量和土壤滲透性的參數(shù)指標過于單一，為進一步深入剖析不同土壤環(huán)境下的地下水波動特征，還需要對上述控制方程做進一步完善，以期更深入地認識潛水層內(nèi)的水波運動理論。

［1］陳吉余，陳沈良.中國河口海岸面臨的挑戰(zhàn)［J］.海洋地質(zhì)動態(tài)，2002，18(1):1-5.(CHEN Jiyu，CHEN Shenliang.Estuarine and coastal challenges in China［J］.Marine Geology Letters，2002，18(1):1-5.(in Chinese))

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

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

［4］LI Ling，BARRY D A，PATTIARATCHI C B.Numerical modeling of tide-induced beach water table fluctuations［J］.Coastal Engineering，1997，30:105-123.

［5］PARLANGE JY，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:251-263.

［6］TURNER I.Water table outcropping on macro-tidal beaches:a simulation model［J］.Marine Geology，1993，115:227-238.

［7］MILLHAM N P，HOWESB L.A comparison of methods to determineKin a shallow coastal aquifer［J］.Ground Water，1995，33(1):49-57.

［8］PARLANGE J Y，BRUTSAERT W.A capillary correction for free surface flow of groundwater［J］.Water Resource Research，1987，23:805-808.

［9］BARRY D A，BARRY SJ，PARLANGE J Y.Capillarity correction to periodic solutions of the shallow flow approximation［M］.Washington，D.C.:American Geophysical Union，1996:496-510.

［10］LI Ling，BARRY D A，PATTIARATCHI C B.Numerical modeling of tide-induced beach water table fluctuations［J］.Coastal Engineering，1997，30:105-123.

［11］LILing，BARRY D A，CUNNINGHAM C，et al.A two-dimensional analytical solution of ground-water responses to tidal loading in an estuary and ocean［J］.Advance of Water Resource，2000，23:825-833.

［12］ KONG Jun，SHEN Chengji，XIN Pei，et al.Capillary effect on water table fluctuations in unconfined aquifers［J］.Water Resource Research，2013，49:1-6

［13］GARDNER W R.Some steady-state solutions of the unsaturated moisture flow equations with application to evaporation from a water table［J］.Soil Science，1958，85(4):228-232.