大小潮作用對潮灘沉積物層理影響的數值模擬研究

徐孟飄，東培華，馬駿，羅鋒，張長寬，范代讀，周曾,,5*

(1.河海大學 水文水資源與水利工程科學國家重點實驗室，江蘇 南京 210098；2.華設設計集團股份有限公司，江蘇 南京 210014；3.杭州市港航管理服務中心航道管理處，浙江 杭州 310016；4.同濟大學 海洋與地球科學學院，上海 200092；5.河海大學 江蘇省海岸海洋資源開發與環境安全重點試驗室，江蘇 南京 210098)

1 引言

潮灘是細顆粒沉積物在潮流為主導的水動力條件下形成的寬闊而平坦的淺灘[1]。潮灘上的水動力過程及沉積物輸移過程呈雙向反饋，集中表現為潮灘剖面形態的演變及潮灘上多組沉積物的分選現象。不同類型、屬性的沉積物顆粒在潮灘上呈現不同的空間分布特征，該現象通常被稱為沉積物分選。沉積物分選現象主要包括水平向分布和垂向分層兩個方面，常見的水平向分布現象例如潮灘沉積物粒徑的陸向細化，垂向分層現象一般通過潮灘底床沉積物層理來描述，沉積物在水平向和垂向的分選過程一般是緊密聯系的。

層理是描述沉積特征常用的術語[2]，潮汐沉積物中的潮汐層理展現出的周期韻律性可用于估算風暴潮事件、推算古地月軌道參數，有助于推算地史時期地?月系的演化歷史[3]，同時有助于分析潮灘沉積歷史及提高對海平面變化規律的認知等[3]，因此開展粉砂淤泥質潮灘上多組分泥沙的垂向分選機制的研究具有重要的實際意義。潮灘上潮汐的周期性作用（漲落潮周期、大小潮周期、季節性及更長時間尺度的周期）在沉積地貌上表現出的沖刷與淤積的交替促使了互層層理（砂質層和泥質層構成）的形成[4]。對于互層層理的形成機制及內在機理的認知仍未形成統一認識，許多學者從現場實際觀測入手分析其機制并取得一些進展。Deloffre等[5]在對法國、英國等河口潮間帶潮灘的研究中發現沉積物表現出大小潮周期性特征，王建等[6]對江蘇中部淤泥質潮灘進行實地觀測認為毫米級薄水平互層層理由半日潮產生，而厘米級厚砂泥互層層理則為半月天文潮（大小潮）的產物范代讀等[7]和Fan等[8]對長江口潮灘進行觀測，發現潮汐層偶的厚度與大小潮周期潮差旋回性有關，層偶厚度是潮差的函數，隨潮差由小潮?大潮?小潮而發生由薄?厚?薄的規律性變化。相較于現場實際觀測方面，數值模擬層面的研究相對較少，是有待進一步深入的方向。本文以潮流為主導作用的潮灘為例，建立了潮灘一維概化動力地貌模型，分析沉積物垂向分布規律，探究周期性潮汐條件及不同邊界含沙量對沉積韻律層結構（層理組成、層理厚度等）的影響。

2 研究方法

本文使用Delft3D模型系統中的水動力模塊、沉積物輸移模塊以及動力地貌演變模塊來進行一維水動力及地貌演變的數值模擬計算。Delft3D模型系統是一套可進行水流、波浪、水質、生態、泥沙輸移、床底地貌及各個過程之間相互作用數值模擬的強大開源軟件包，在內河、河口、海岸、沉積物、地貌演變水質等各方面被廣泛應用。

2.1 模型介紹

2.1.1 水動力模塊

假設潮灘沿岸均勻，可使用一維模式對潮灘剖面進行模擬，在靜壓假定和Boussinesq近似兩個假定的基礎上，對不可壓縮流體的Navier-Stokes方程進行求解。一維淺水方程包括連續性方程和動量守恒方程。

（1）基本方程

① 連續性方程：

② 動量守恒方程：

式中，x為空間水平坐標；u為垂線平均流速在x方向的分量（單位：m/s）；h為全水深，即水面到水底的距離（單位：m）；η為相對于某標準（例如：平均海平面，MSL）的水位值（單位：m）；v為平面紊動黏性系數（單位：m2/s）；g為重力加速度（單位：m/s2）；C為謝才系數（單位：m1/2/s）。

模型中設置某一時間步長中當計算單元的水深小于 0.1 m（臨界值）時，處于“干”狀態，該計算單元不參與計算。當該計算單元水深大于0.2 m（臨界值）時，處于“濕”狀態，該計算單元重新被激活參與模型計算。

（2）由潮流引起的底部切應力

計算公式為： τc= ρgu2/C2，式中，ρ為水的密度（單位：kg/m3）。

（3）波浪模塊

Delft3D默認的波浪模塊為SWAN（Simulating Waves Nearshore）波浪模型，由于長周期模擬需要的計算時長較長因此本文未使用SWAN模塊，而是采用了基于 Young 和 Verhagen[9]以及 van Rijn[10]提出的公式開發的簡單的風浪模型，該公式于2009年被證實適用于模擬淺灘上風浪產生的作用[11]，該公式將無量綱波能量（ε =g2E/UW4）及無量綱譜峰波頻率（?=與無量綱風區及無量綱水深聯系起來，對已知風場（例如：風區長度、風速及水深）的波高及波周期進行估算。

式中，A1=0.493δ0.75，B1=3.13×10?3χ0.57，A2=0.331δ1.01，B2=5.215×10?4χ0.73；UW為海拔10 m處的風速；F為風區長度；TP為峰波周期；E為波能，其公式為E=ρgHS2/8（HS為有效波高）。當風條件給定時，解以上方程可估算有效波高HS和峰波周期TP。由于波浪在傳播過程中波陡增大，波浪在淺灘上會破碎，衡量波浪破碎的典型破碎指標 γ為有效波高和水深的比值（γ =HS/h），在本文中，參考Roberts等[12]參數設置，設定 γ =0.5。在波浪破碎后，波高變為受水深限制，即HS=γh。當波浪破碎時，波浪引起的底部切應力達到最大值且隨著水深的向陸減小而減小[13]。根據線性波理論，底部最大軌道速度為

且由波浪引起的底床切應力為

式中，k為波數；fw為由Soulsby[14]定義的波浪摩擦因子；為非黏性沉積物的中值粒徑。綜合底床最大切應力是引起沉積物輸移的主要因素，由潮流和風浪引起的底床切應力的疊加通常大于兩者的線性總和，在本模型中參考文獻[14]，綜合底床最大切應力的公式為

2.1.2 沉積物輸運模塊

模塊中的沉積物分為黏性和非黏性兩種，在計算沉積物的輸運時，黏性部分與非黏性部分使用不同方法進行計算，即該模塊忽略了沉積物混合物的臨界再懸浮條件與單個種類沉積物的臨界條件的差別。黏性沙和非黏性沙的物理過程（例如侵蝕、淤積和輸移）均分為不同部分考慮。

黏性沉積物的輸運使用帶源匯項的平流方程來描述黏性沙的侵蝕與淤積過程

式中，Qmud,e與Qmud,d分別代表侵蝕量與淤積量，該沖淤量使用被廣泛運用的經典的Patheniades-Krone公式來計算，公式如下：

式中，Me為沖刷系數（單位：kg/（m2·s））；ωs為沉降速度（單位：m/s）；c是垂向平均濃度（單位：kg/m3）；τcr,e和τcr,d（ 單位：Pa）分別表示黏性沉積物的臨界起動切應力及臨界沉降切應力。該模型中使用Winterwerp[15]提出的“連續沉積”的概念，因此黏性沉積物的臨界沉降切應力 τcr,d取 較大值（默認為1 000 Pa），則上公式可近似等于Qmud,d=ωsc。

對于非黏性沉積物，使用Sousby-van Rijn公式來計算潮流與波浪共同作用下的輸移量。總沉積物輸移量分為推移質和懸移質兩類來計算

其中，

式中，Acal是 校準系數；CD是僅有潮流產生的阻力系數；Urms為 波質點軌道流速的均方根ucr為泥沙顆粒起動的臨界流速；s為沉積物的相對密度 ；D?無量綱粒徑大小

2.1.3 地貌演變模塊

床面變形方程為

式中，D是沉積物的沉積通量；E是沉積物的侵蝕通量；ρs為沉積物的密度。地貌演變的計算是基于黏性沙與非黏性沙這兩種沉積物在每個水動力計算步長上 的侵蝕或淤積量的更新。

2.1.4 底床分層模塊

模型中將初始底床分為若干薄層，分為3大類：沖淤層、交換層以及基準層。沖淤層的厚度 δa是一個用戶根據模擬需要選取的參數，在模型計算時沖淤層中的沉積物發生侵蝕或者淤積與水體中的沉積物進行交換，但厚度始終保持不變來保證模型的穩定性，交換層在計算過程中厚度可變化，但存在一個最大厚度 δu。

模型開始計算之前，不同種類沉積物（黏性沙及非黏性沙）由初始設定的比例進行均勻混合且被分為若干層，當模型開始計算時，最頂層的沖淤層會發生侵蝕或淤積，且沖淤層下的交換層數量及厚度也會隨之發生變化。底床分層模型及模型中的侵蝕或淤積的計算過程見圖1。當床面發生侵蝕時，沖淤層會發生沖刷且部分沉積物會被帶走，當侵蝕量大于 δa厚的沖淤層時，其下的交換層會向沖淤層進行沉積物補充使侵蝕過程繼續執行，在下一個計算步長開始之前交換層會對被侵蝕過的沖淤層執行補充過程，使沖淤層厚度再次回到定值 δa，此時交換層厚度減小，以此反映床面的侵蝕過程。當床面發生淤積時，沖淤層會變厚，在下一個計算步長開始之前會執行補充過程，沖淤層厚度將保持定值 δa將多余的沉積物向下傳遞并入交換層，若交換層達到所設定的厚度值 δu，則底床中 將會產生一層新的交換層，以此反映床面的淤積過程。

圖1 底床分層模型及其計算過程Fig.1 Bed stratigraphy modeling and schematic diagram of calculation process

2.2 模型建立

本文基于我國沿海典型粉砂淤泥質潮灘概化模型對潮灘沉積物垂向分層機制開展研究，研究使用的底床分層模型由于模型設置垂向層數較多，長期沉積過程計算時間很長，且本研究主要關注的是垂向的沉積物分布特征，因此選擇一維模型更為合適，在不影響研究需求的前提下大大提高了模型計算效率。

粉砂淤泥質潮灘坡度較緩，約為1‰～3‰[16]，本文基于Zhou等[17]的研究進行模型設定，采用一維模型，初始床面為高程從–14 m 至–1 m（平均海平面以下）的緩坡，且垂向分層層厚1 mm。模型使用50 m的等距矩形網格，共計280個網格，時間步長Δt=0.25 min，滿足模型計算穩定性及精度要求。圖2為初始地形剖面示意圖，假設初始底床由50%的黏土、25%粉砂和25%的細砂組成，且這3種沉積物沿深度方向均勻混合。

圖2 初始床面組成Fig.2 Initial bed composition

模型的東側為開邊界其余均為固邊界，開邊界水動力條件主要使用M2及S2兩個分潮，M2為太陰主要半日分潮，S2為太陽主要半日分潮。為研究垂向分層結構對不同潮流條件的響應，設置了6組工況進行計算，工況設置如表1所示。

表1 水動力條件工況設置Table 1 Hydrodynamic condition setting

由于粉砂淤泥質潮灘主要由黏土、粉砂及細砂構成，因此模型同時考慮黏性沙與非黏性沙，設定3種沉積物：黏土、粉砂和細砂。粉砂介于黏性和非黏性之間，處理時有時被當做黏性砂，但也有一些學者發現由于粉砂的成分與砂相似，透水性強、黏性差，一些國內外同行在研究中也有將之處理為非黏性砂。本研究中的粉砂粒徑設置為30 μm，在粉砂界定的3.9～62.5 μm范圍內，顆粒相對較大，將其考慮為非黏性沙。其中，黏土設置為黏性沙且以懸移質的方式輸移，而粉砂與細砂均設置為非黏性沙，既可以懸移質的方式輸移也以推移質的方式輸運。初始底床由2 m的黏土、1 m的粉砂和 1 m的細砂組成，3種沉積物均勻混合。經過調試，黏土沉速取為0.5 mm/s，臨界起動切應力為 0.2 N/m2，沖刷率系數取為 5×10?4kg/(m2·s)，為了使黏性沙始終處于沉降狀態，因此黏性沙的沉降臨界切應力取較大值設定為 1 000 N/m2（默認值），黏性沉積物干密度為500 kg/m3。屬于非黏性沙的細砂的中值粒徑D50取為90 μm，粉砂的中值粒徑D50取為30 μm，由于Delft3D軟件內默認的非黏性沙計算公式不能設置粒徑小于 63 μm的顆粒，因此本文調用Van Rijn[10]公式進行非黏性沙的計算，非黏性沉積物干密度為1 600 kg/m3。在東側開邊界給予泥沙含沙量邊界條件：黏土為 0.006 kg/m3，粉砂為 0.004 kg/m3，細砂為0.001 kg/m3。由于本研究是概化模型，側重于揭示一些機制，若邊界含沙量設置過大，易形成邊界的劇烈淤積造成模型無法長期運行，因此選取的工況邊界含沙量相對較小，但不會影響所得結論。

3 結果與分析

3.1 垂向分層結構對潮流的響應

為研究潮流對潮灘底床沉積物垂向結構的影響，模型采用了單M2分潮作用及M2、S2共同作用兩種工況，其模擬的沉積物分選、分層特征結果如圖3所示，計算結果是10 a接近平衡態。由于邊界含沙量相對較低，潮間帶發育相對較慢，且表現出不完全，但這對本文側重說明的泥沙垂向層理現象影響較小。

圖3 單M2分潮作用下（a）以及大小潮作用下（b）潮灘沉積物分選、分層特征，代表點A、B、C將用作后續分析Fig.3 Results of sediment sorting and layering under only M2tide (a) and spring/neap tidal cycles (b),representative points A,B and C are used for subsequent analyses

在僅M2分潮作用下，選取潮間帶上的3點A、B、C點來進行對比分析（A點距離岸線0.5 km，B點距離岸線 1 km，C 點距離岸線 1.5 km），圖 4為 3個代表點僅在M2分潮作用下垂向沉積物占比圖。

圖4 僅 M2分潮作用下 A、B、C 3 點處垂向沉積物占比Fig.4 Percentage of vertical sediments at points A,B and C under the action of M2 tidal constituent

由圖4可見，潮間帶上A、B、C 3點，由于底床高程的不斷抬高，同一位置的水動力條件發生改變，各類沉積物垂向含沙量也發生改變。A、B、C 3點細砂垂向含沙量占比幾乎為0。A點位于高程最大處，潮間帶上部黏性沉積物占主導，因此垂向沉積物含沙量中黏土占比均大于粉砂占比。B點位于黏土與粉砂分布的過渡區域，垂向上在高程較小處粉砂占比較多，但隨著沉積物淤積底床抬升，動力條件減弱，在高程較小處沉積物的垂向含沙量發生變化（黏土占比大于粉砂占比），而C點處于粉砂主導區域，因此粉砂占比整體大于黏土占比。盡管3個代表點在垂向占比上有所不同，但在垂向上含沙量曲線較為平滑沒有明顯波動，且均未見到明顯的層理現象，可能原因包括：（1）模型中底床分層層厚設置較大（1 mm）限制了單個潮周期內潮汐層偶現象的形成；（2）模型中將粉砂當作一個組分，未考慮粗、細粉砂的不同行為過程；（3）模型未考慮絮凝沉積作用和細粒沉積物在灘上沉積后暴露空氣發生脫水固結作用對再侵蝕的影響等。

圖5 大小潮作用下 A、B、C 3 點處垂向沉積物占比Fig.5 Percentage of vertical sediments at points A,B and C under the action of the superposition of M2and S2constituents

圖5為M2和S2分潮疊加的大小潮作用下潮間帶上離岸 0.5 km、1 km、1.5 km 的 A、B、C 3點處的垂向沉積物占比分布，由于潮間帶上的沉積物主要由粉砂與黏土組成，因此圖中只對比了粉砂與黏土的垂向占比分布情況。由圖可見，相較于圖4，垂向上黏土與粉砂的含量分布線呈鋸齒狀表現出明顯的波動且具有一定的規律性，這種周期性的規律主要是由大小潮的周期旋回造成的。A點處高程較小處的鋸齒狀較高程大處更明顯，是由于隨著沉積物的不斷淤積地貌不斷變化使得水深在減小，潮動力相對的越來越弱，相對較粗顆粒的沉積物（粉砂）更難以被搬運至近岸處。從水平向沉積物的分布規律來看，B點位于黏土區域與粉砂區域的過渡區域，在此區域中黏土與粉砂對于潮汐能量的變化響應較好沉積物占比變化劇烈，因此鋸齒狀明顯。而C點位于粉砂與細砂的過渡帶，該點垂向上高程較小處細砂含量相對較大，因此由黏土與細砂形成的層理結構表現較差。

為研究沉積物垂向分布與潮汐周期的對應關系，選取離岸0.6 km、高程約為–1.4 m處的P1點以及離岸0.2 km、高程約為–1.1 m處的P2點進行分析（即模型的一個網格50 m），圖6為黏土與粉砂的垂向分布情況，圖中的沉積物于29 d內形成，包含兩個大小潮周期，該段時間內的潮位過程、P2點累計層數及含沙量變化如圖7所示。

圖6a代表黏土的占比分布，顏色越紅代表黏土占比越多，越藍代表占比越少，圖6b代表粉砂的占比分布，同理越紅表示粉砂越多，越藍代表粉砂越少，由圖可見，P1點與P2點均包含完整的潮汐韻律層，韻律層由位于上部的砂質層及位于下部的泥質層構成，且韻律層的個數與大小潮周期數量相一致。P2點垂向層理結構中黏土占比均大于P1點，這是由于P2點較P1點高程更大，水動力條件相對弱，粗顆粒沉積物更難輸移至較高處。由圖7可見，該潮位過程呈現“大潮–小潮–大潮–小潮”的變化趨勢，與潮位過程相對應，垂向分層結構中沉積物顆粒均呈現“較粗–較細–較粗–較細”的變化規律。同時潮汐大小潮周期性的變化規律也體現于層理厚度的旋回性變化之中，由于大潮期間潮汐能量大，沉積物在較強的動力條件下更易起動并發生輸運，潮差大時層理較厚而潮差小時能量衰減，因此層理較薄。

圖6 P1 點（a，b）和 P2 點（c，d）兩個大小潮周期內垂向沉積物占比Fig.6 Percentage of vertical sediments in two spring-neap tidal cycles at poiots P1 and P2

圖7 潮位過程（a），P2 點累計層數（b）和含沙量變化（c）Fig.7 Tidal level (a)， cumulative layers (b) and changes of percentage of vertical sediments (c)

3.2 邊界含沙量對垂向分層結構的影響

現場觀測表明，水體中懸浮沉積物的濃度也會影響大小潮周期層偶結構的組成。因此，為研究水體含沙量對于沉積層偶結構特征的影響，在相同水動力條件的基礎上設置3種邊界含沙量工況，如表2所示。3種工況中，工況三的含沙量量級與通常認為的相對接近，工況一及工況二的設置實際上是對層理結構中“粗–細”形成機制的對比試驗及研究，工況二及工況三的設置實際上是對層理結構中“厚–薄”形成機制的對比試驗及研究。此外，由于本研究是概化模型側重于揭示一些機制，在邊界含沙量設置過大時，容易形成邊界的劇烈淤積造成模型無法長期運行，因此選取的工況邊界含沙量相對較小，但不會影響所得結論。

表2 3種工況的邊界含沙量設置Table 2 Three settings of sediment concentration

圖8為不同工況下潮間帶區域（0～4 km）黏土與粉砂占比分布對比，黏土占比圖中顏色越紅代表黏土占比越多，顏色越藍代表黏土占比越少，粉砂占比圖中顏色含義同理。在工況一情況下，黏土邊界含沙量略大于粉砂含沙量，潮汐韻律層結構于離岸0～2 km處較為明顯，主要是由黏土與粉砂組成，而2～4 km處主要是由粉砂與細砂組成的層理。在粉砂邊界含沙量大于黏土邊界含沙量的工況二的情況中，可見在0.5～4 km范圍內粉砂占比明顯增大，在離岸0～2 km處雖然仍是黏土與粉砂組成的層理結構，但不論是大潮期形成的砂質層還是小潮期形成的泥質層層理中的沉積物顆粒一定程度上均發生粗化，而在離岸2～4 km范圍內由粉砂與細砂構成的韻律層中的沉積物顆粒在一定程度上均發生了細化。工況三中由于黏土與細砂的含沙量是工況二中的10倍，由圖8明顯可見潮間帶區域內整體淤積且最高處達到1 m左右，且離岸2～4 km范圍內黏土與粉砂占比均有所增大，其中粉砂占比增大更明顯，該范圍內變為黏土、粉砂與細砂3種沉積物構成的層理結構，韻律層內的沉積物顆粒發生了明顯的細化。并且，當邊界含沙量濃度增大時，整個潮間帶區域內大小潮對應的層理更加完整。

圖8 不同工況下潮間帶區域黏土與粉砂占比分布Fig.8 Distribution of clay and silt fractions in intertidal flat under different sediment concentration conditions

圖9為代表點（離岸1 km處）不同工況下兩個大小潮周期內黏土與粉砂的垂向占比分布，由圖可見，當粉砂邊界含沙量大于黏土邊界含沙量時，潮汐韻律層中的沉積物均發生粗化。而當黏土與粉砂邊界含沙量整體增大時（工況三）該代表點處兩個潮周期內產生的潮汐韻律層明顯增厚，但整體沉積物顆粒粒徑變小，這主要是由于粉砂的邊界含沙量雖然較大，但水體的挾沙能力有限，很大一部分粉砂較難隨潮流輸移至高灘處而是停留在高程較小處，而較大部分黏土卻能輸移至高灘處并沉積，因此粉砂占比反而減小。

圖9 代表點（離岸1 km處）不同工況下兩個大小潮周期內沉積物垂向分布Fig.9 Vertical distribution of sediments in two tidalcycles under different conditions at representative points

4 結論

潮灘垂向沉積韻律層的形成是多種動力因子共同作用的結果，其中大小潮周期性變化是韻律層形成的重要原因之一，本文著重從數值模擬角度探討了大小潮周期的影響，分析了潮灘垂向沉積物的分選與分層結構形成的機理。具體結論如下：

（1）大小潮作用下，潮間帶區域可見明顯的分層結構（主要由黏土與粉砂形成），黏土與粉砂含沙量曲線沿垂向表現出明顯的波動呈鋸齒狀。

（2）結果表明一個層理結構由泥質層與砂質層構成，其中泥質層在小潮期間形成，而砂質層于大潮期間形成，同時層理的厚度也與大小潮周期相關，大潮時層理較厚而小潮時層理較薄。

（3）除了周期性的潮汐條件這個主要因素，懸浮泥沙濃度也是影響大小潮周期性層理結構的重要因子。邊界含沙量中粉砂占比增大會使得大小潮周期性層理結構中泥質層與砂質層均粗化且砂質層厚度增大，當邊界含沙量整體顯著增大時，潮灘上的垂向潮汐韻律層會更加完整且厚度明顯增大。

在本文的研究中，在計算沉積物的輸運時，黏性與非黏性沙使用不同方法進行計算，忽略了泥沙混合物由于黏性沙的存在呈現出的一定的黏性，可能會使得泥沙輸運通量略有增大。同時，在潮灘底床垂向沉積物觀測中，除了大小潮表現出的韻律性層理結構，漲落潮期間形成的潮汐層偶（基本沉積單元）以及波浪、風暴（風暴來臨的水位、風暴歷時）等動力因子作用下的層理重塑或破壞也常被觀測到的現象，是未來數值模擬研究的重點，進而可以更加深入地掌握潮汐層理的形成與演化機制。