波浪非線性對溶質進入海底沉積物運移過程影響研究?

劉小麗，賀建檁

(1.山東省海洋環境地質工程重點實驗室，山東 青島 266100；2.中國海洋大學 環境科學與工程學院，山東 青島 266100)

大量污染物的排海以及海上溢油事件的頻發，使得海洋環境污染問題受到越來越多的關注。懸浮于水體中的污染物質，在水動力的作用下會進入到海底的沉積物中，從而對海洋生態環境產生進一步影響。波浪作用下污染物質(溶質)進入沉積物遷移過程的研究，對于海洋環境的評估具有重要意義。

作為近岸海域中壓力梯度的來源，波浪會引起海底沉積物中孔隙水的滲流，其對溶質輸運的影響不容忽視[1]。對于平底海床，由于規則波的周期循環特征，一個波浪周期內上覆海水與沉積物床之間的凈交換水量基本為0，因此，當海水與沉積物界面上的溶質濃度均勻分布時，規則波作用下對流作用對溶質在沉積物中的運移過程的貢獻可以忽略[2-4]。Harrison等[2]對波浪作用下溶質進入平海床運移特征的理論分析表明，波浪通過促進溶質的機械彌散作用而加速溶質的遷移，波浪作用下溶質在沉積物中的有效擴散率可比分子擴散速率大4個數量級。Clark等[5]通過波浪水槽試驗研究了行進波對溶質進入礫石床運移特征的影響，結果顯示，行進波極大程度地提高了溶質進入底床的運移速率和運移深度，而進一步的理論分析表明，行進波主要通過增強機械彌散作用促進溶質運移。波浪作用下上覆水體與沉積物床之間的溶質交換過程可以通過溶質的對流-彌散方程進行描述。Qian等[6]基于二維溶質對流-彌散方程，對線性行進波和駐波作用下溶質進入湖床沉積物的溶質擴散過程進行了數值模擬，計算結果表明，波浪作用可以使湖床沉積物對氧的吸收率提高2個數量級。Geng和Boufadel[7]利用溶質對流-彌散方程，通過數值模擬研究了波浪和潮汐作用下淺灘含水層溶質運移特征，結果表明，海洋水動力改變了溶質運移的路徑，波浪增加了灘面的交換通量，并拓寬了溶質交換區域。劉小麗和劉明珠[8]基于二維Biot固結理論和對流-彌散方程，考慮了波致海底沉積物的變形效應，對行進波作用下海水中的溶質進入沉積物中的運移過程進行了數值模擬，計算結果表明，沉積物的變形效應會促進溶質的運移速率，增大溶質的運移深度。

以上研究表明，波浪對溶質進入沉積物的運移過程具有明顯的促進作用，利用溶質的對流-彌散方程可以較好地模擬水動力作用下溶質的運移過程，但目前絕大多數的研究都是基于線性波浪作用的分析。當波浪由深水區傳至近岸淺水區域時，隨著水深的減小，波浪的非線性特征不斷增強，當近岸區的相對水深(即水深與波長之比)小于1/8～1/10時，線性波浪理論甚至Stokes非線性波浪理論已經不再適用，而需要采用橢圓余弦波理論對近岸淺水區的非線性波浪特征進行描述[9-10]。

因而，對于近岸淺水區，有必要對非線性波浪作用下溶質進入沉積物中的運移特征進行分析，以揭示波浪非線性特征對溶質運移過程的影響。本文基于Biot多孔介質理論和溶質運移的二維對流-彌散方程，建立了波浪作用下考慮沉積物變形效應的溶質遷移數值計算模型，通過數值模擬，對線性波，非線性二階Stokes波以及一階橢圓余弦波作用下的溶質遷移過程進行對比分析，以揭示波浪非線性對溶質進入沉積物遷移特征的影響。

1 數值計算模型

如圖1所示，為波浪作用下海床溶質運移計算的幾何模型示意圖，坐標軸原點o位于海床面的左端點，x軸位于海床面，y軸以向上為正，波長為L，波高為H，水深為d，海床厚度h。海水-沉積物交界面處的溶質濃度為c0=1 mg·L-1；沉積物中的初始溶質濃度為cs0= 0。

1.1 海床表面波壓力

(1)線性波

海床表面波壓力方程為

(1)

ω2=gktanh(λd)。

(2)

(2)二階Stokes波

二階Stokes非線性波在海床表面波壓力方程式為[12]：

Pb(x)=P0(C11cos(λx-ωt)+C02+C20cos2(λx-ωt)+C22cos2(λx-ωt))。

(3)

式中：

(3)一階橢圓余弦波

一階橢圓余弦波在海床表面的波壓力解析式為[13]：

(4)

1.2 滲流場控制方程

基于Biot彈性多孔介質固結理論[14]，考慮沉積物的變形與孔隙水滲流的耦合效應，控制方程如下：

(5)

(6)

(7)

(8)

β為孔隙流體壓縮系數，表達式為[15]:

(9)

式中：Kw為孔隙水體積壓縮模量，數值為2×109Pa；Sr為飽和度；pw0表示絕對孔隙水壓力。

由達西定律，孔隙水在水平和豎向的實際滲流速度u、v分別為：

(10)

1.3 溶質運移方程

依據質量守恒定律，多孔介質中溶質二維對流-彌散方程為[16]：

(11)

(12)

(13)

(14)

1.4 邊界條件與初始值

結合圖1所示的幾何模型示意圖，數值模型的邊界條件如下。

在海床表面y=0位置，土體的剪切力和有效正應力均為零；海床面的孔隙水壓力為波浪在海床表面的波壓力Pb(x)，根據需要可分別為線性波、二階Stokes波和一階橢圓余弦波的波壓力，從而可以對不同類型波浪場作用下的海床滲流場進行計算；溶質在海床面上的濃度c(x,0,t)=c0=1 mg·L-1。

在海床底部y=-h處，為不透水邊界，土體位移為0，即：

在海床兩側x=0和x=L處，為周期性邊界，即：

p|x=0=p|x=L,wx|x=0=wx|x=L，

wy|x=0=wy|x=Lc|x=0=c|x=L。

此外，由于砂土海床的滲透系數較小，沉積物中的滲流符合達西定律，其對波浪場的影響較小[18]，故不考慮多孔介質海床對波浪場的影響。

2 數值模型驗證

2.1 非線性波浪計算結果驗證

通過與波浪試驗結果的對比，驗證本文非線性波計算結果的正確性。首先對非線性Stokes波的結果進行驗證，如圖2(a)中所示，為解析解計算的海床面波壓力與Qi等[19]波浪水槽試驗測試的海床表面波壓力對比圖，相應的波浪參數為H=0.095 m，T=1.2 s，d=0.5 m，由圖可知，Stokes非線性波的計算結果與試驗結果吻合較好，驗證了本文Stokes波計算結果的正確性。其次，對一階橢圓余弦波的計算結果進行驗證，如圖2(b)中所示，為一階橢圓余弦波理論計算得到的波浪自由面高程與Chang等[20]試驗記錄的橢圓余弦波自由面數據對比圖，其中波浪周期T=2.0 s，相對波高H/d=0.15，水深d=24 cm，波長L=2.97 m；圖中試驗數據為距離造波板4.8 m位置處測量的橢圓余弦波的自由面高程，可以看出本文解析解計算結果和試驗結果一致性較好。綜上可知，本文基于非線性波浪理論計算得到的結果正確，能夠反映波浪的非線性特征。

2.2 滲流計算結果驗證

Hsu和Jeng[15]基于Biot固結理論，考慮了海床土體的變形效應，推導了波浪作用下海床瞬態響應的解析解。將本文數值計算結果與該解析解進行對比。基本參數為：波高H=5 m，波長L=121.19 m，水深d=20 m，周期T=10 s，孔隙率n=0.333，泊松比μ=0.35，剪切模量G=5×106Pa，海床厚度h=25 m。對比結果如圖3所示，圖3(a)為不同滲透系數下孔壓響應對比；圖3(b)為不同剪切模量下孔壓響應對比。從圖3中對比結果可以看出，數值計算結果與解析解一致性較好，表明當前模型能夠較準確地模擬波致海床的滲流場。

2.3 溶質運移計算結果驗證

Packman等[21]基于礫石底床，通過水槽試驗探究了不同水力條件下上覆水體中的溶質進入礫石床的情況，得出了上覆水體中溶質濃度隨時間的變化曲線。將本文數值模型計算的溶質運移結果與試驗結果進行對比。分別選取試驗工況7和工況9的結果進行比較，計算參數為：孔隙率n=0.38，平均顆粒粒徑dg=0.6 cm，滲透系數k=0.15 m·s-1，沉積物厚度為0.19 m，工況7中水深d=20.4 cm，波長L=32 cm，波高H=2 cm；工況9中水深d=20.4 cm，波長L=32 cm，波高H=3.76 cm。對比結果如圖4所示，可以看出，數值模型的計算結果和試驗測試結果較為一致，表明了本模型用于溶質運移計算結果的可靠性。

3 結果分析與討論

隨著波浪向岸傳播，水深逐漸變淺，波浪非線性特征逐漸增強。當波浪Ursell數(其計算表達式為HL2/d3)大于40時[22]，橢圓余弦波理論能更有效地表達淺水非線性波浪的特征。為了更好地分析波浪非線性對海底沉積物中溶質運移特征的影響，同時對線性波浪、二階Stokes波和一階橢圓余弦波作用下溶質進入沉積物的過程進行數值模擬，具體計算參數見表1。波浪的Ursell數為43，屬于典型的淺水非線性波浪。

3.1 波浪非線性對滲流場特征的影響

溶質運移特征主要取決于滲流場特征，因此，首先對不同類型波浪作用下沉積物中的滲流場特征進行分析。

如圖5所示，為3種波浪理論下海床表面的波壓力分布圖，從中可以看出，線性波在波峰和波谷處的波壓力絕對值相等，相對于靜水位面對稱；非線性波的波峰和波谷相對于靜水位面呈現明顯的非對稱性，波峰處波壓力增加，波谷位置的波壓力降低。一階橢圓余弦波的非線性特征相對于二階Stokes波更加明顯，一階橢圓余弦波和二階Stokes波在波峰處的波壓力分別為線性波波壓力的165%和120%，而在波谷位置二者的波壓力基本一致，均為線性波波壓力的80%。

表1 基本計算參數Table 1 Parameters for the base case

如圖6所示為線性波作用下一個周期內不同時刻的孔隙水流速矢量圖和等值面圖。由于海床面不同位置存在水頭差(壓力差)，因此在海床內部會發生滲流，如圖6中所示，波峰處發生向下的滲流，波谷處發生向上的滲流；對于同一位置，其滲流流速的數值和方向隨著波壓力的變化呈現出周期性變化。以海床表面任意點處的滲流為例，如圖7中所示，豎向和水平向的滲流流速均隨時間呈現出周期性變化，受波浪非線性的影響，一階橢圓余弦波作用下波峰處的豎向滲流流速與波節點附近的水平向滲流流速均顯著大于線性波和二階Stokes波；總體而言，線性波和二階Stokes波的流速數值相差較小。對3種波浪作用下沉積物和上覆海水之間的水量交換進行分析可知，無論是哪種波浪，一個周期內流入和流出沉積物的水量基本相同，即沉積物和上覆海水之間水的凈交換量為0，因此，正如Harrison等[2]所述，在周期性循環波浪荷載作用下，對流作用對溶質質量運移的貢獻基本可以忽略，波浪主要通過機械彌散作用促進溶質在沉積物中的運移。

3.2 波浪非線性對溶質運移特征的影響

一階橢圓余弦波作用下周期平均縱向水動力彌散系數的增大，促進了溶質在沉積物中的運移。如圖11中所示，在同一深度，一階橢圓余弦波作用下的溶質濃度大于線性波作用下的溶質濃度，表明一定條件下的淺水非線性波可以促進溶質在沉積物中的運移。進入沉積物中的溶質質量隨波浪作用時間的變化如圖12中所示，從中可以更為清楚地看出，一階橢圓余弦波作用下進入沉積物中的溶質質量明顯大于線性波，3 000個波浪周期作用后，質量增加量可達線性波的7.2%。另外，從圖11和12均可以看出，二階Stokes波作用下溶質濃度的分布以及進入沉積物中的溶質質量均與線性波作用下的結果一致，表明二階Stokes波的非線性效應對溶質運移過程的影響很小，基本可以忽略。上述模擬結果的分析表明，一階橢圓余弦波非線性特征較強，能顯著促進溶質在海底沉積物中的運移；相比之下，二階Stokes波的非線性特征相對較弱，其對溶質運移過程的影響相對小。

3.3 波浪參數影響分析

3.3.1 波高的影響 根據前述內容可知，可以通過分析不同波浪參數對孔隙水滲流流速和溶質縱向水動力彌散系數的影響規律，來反映其對溶質運移的影響規律。Ursell數的大小表征了波浪非線性特征的強弱[13]，因此，基于Ursell數的變化來分析波浪非線性特征強弱的變化。

對一階橢圓余弦波和線性波作用下波高對孔隙水豎向滲流流速和水動力彌散系數的影響進行對比，結合波浪的非線性特征，以分析不同波浪理論下波高對溶質運移過程的影響。波高H的取值分別為2、2.5和3 m，其余參數見表1中所示。

4 結論

以二階Stokes波和一階橢圓余弦波為例，通過與線性波作用下沉積物中孔隙水滲流場和溶質濃度場計算結果的對比，分析了波浪非線性對非反應性溶質從上覆海水進入沉積物中運移特征的影響，主要得到以下結論。

(1) 在周期性循環波浪荷載作用下，沉積物中的滲流場呈現出周期性循環特征，對于平底海床，波浪主要通過機械彌散作用促進溶質沿沉積物深度方向的遷移，一個波浪周期內的平均縱向水動力彌散系數可以用來表征溶質運移速率。

(2) 二階Stokes波的非線性特征對溶質運移的影響基本可以忽略；一階橢圓余弦波的非線性特征顯著，能在一定程度上促進溶質在沉積物中的運移，與線性波相比，其周期平均縱向水動力彌散系數的增加量可達28%。

(3) 線性波和一階橢圓余弦波作用下，波浪參數對溶質運移速率的影響不盡相同。隨著波高增大，兩種波浪作用下溶質運移速率均有所增加；隨著水深的增大，線性波作用下溶質運移速率降低，而一階橢圓余弦波作用下溶質運移速率增加；隨著周期的增大，線性波作用下溶質運移速率增加，而一階橢圓余弦波作用下溶質運移速率降低。

(4) 一階橢圓余弦波的非線性特征隨著水深的減小、波高和周期的增大而增強；與線性波相比，一階橢圓余弦波的非線性特征對溶質運移速率的促進程度隨著水深的減小、波高和周期的增大而降低。