基于Richards方程切換的土壤水流及溶質(zhì)運移數(shù)值模擬

曾季才，查元源，楊金忠

（武漢大學(xué) 水資源與水電工程科學(xué)國家重點實驗室，湖北 武漢 430072）

1 研究背景

土壤飽和-非飽和帶是農(nóng)業(yè)、生態(tài)、水文地質(zhì)及陸面水文過程中的重要環(huán)節(jié)，研究飽和-非飽和土壤水分運動規(guī)律，對地下水資源管理、土壤水環(huán)境演變，及水文循環(huán)系統(tǒng)具有重要的科學(xué)意義。1930年代以來，基于質(zhì)量守恒及達西定律的Richards方程［1］被廣泛用于飽和-非飽和帶水分運動數(shù)值模擬［2］。原始的Richards方程以含水量θ為水量衡量變量，以壓力水頭h為勢能驅(qū)動變量，因而往往被稱為混合型Richards方程。為求解這種帶有兩個未知變量的二階偏微分方程，需要借助土壤水分特征曲線［3］來消除其中一個變量，因而得到分別含水量（θ）和壓力水頭（h）兩種形式的Richards方程。Richards方程同時具有退化橢圓-拋物型方程的雙重特征，因而難以獲取解析解［4］；與此同時，土壤水分特征曲線及邊界條件的高度非線性，使其數(shù)值解的穩(wěn)定性和質(zhì)量守恒問題難以完全解決［4-6］。

研究認為，θ型Richards方程嚴格遵循質(zhì)量守恒，并在干濕交替的大氣邊界條件下，相比h型Richards方程有更高的求解效率和數(shù)值穩(wěn)定性［7］，但這類方程難以準確刻畫非均質(zhì)及飽和土壤中的水流運動過程，限制了其發(fā)展；另一方面，由于水頭在飽和-非飽和非均質(zhì)土壤剖面中連續(xù)可導(dǎo)，h型Richards方程不存在上述困難，但這類方程往往難以保證質(zhì)量守恒［8-9］；直至Celia等［10］對質(zhì)量項進行修正后，傳統(tǒng)h型方程得以大量應(yīng)用于Picard迭代模型中，著名的例子有HYDRUS［11］和SWAP［12］軟件等。即便如此，h型方程的求解效率在處理干濕交替問題時仍然遠低于θ型方程［13］，且常常出現(xiàn)收斂困難等不穩(wěn)定現(xiàn)象，導(dǎo)致了水鹽運移模擬的低效率和低精度問題。為了結(jié)合兩種單一變量方程的優(yōu)勢，人們提出了傳統(tǒng)的主變量轉(zhuǎn)換方法［14］，它基于Newton-Raphson迭代，既能利用h型方程對飽和-非飽和非均質(zhì)問題的準確刻畫，又能保留θ型方程帶來的嚴格質(zhì)量守恒。但由于需要準備高階偏導(dǎo)的Jacobian矩陣，Newton-Raphson迭代算法在數(shù)學(xué)推導(dǎo)及代碼實現(xiàn)上更為困難；此外，在節(jié)點接近飽和狀態(tài)時，這類方法容易因主變量非平滑轉(zhuǎn)換出現(xiàn)錯誤解［15］，因而數(shù)值穩(wěn)定性有待提高。相比之下，Picard迭代算法具有代碼簡單、數(shù)值穩(wěn)定等特點［9，16］；同時，在該迭代算法框架下進行θ和h型方程（或主變量）的切換，能更為簡便地結(jié)合兩類Richards方程的數(shù)值優(yōu)勢，但目前少見報道。究其原因，一方面是廣義θ型Richards方程對非均質(zhì)及飽和土壤的處理［13］尚未廣泛得到應(yīng)用，另一方面是方程層面的切換難以避免在每個臨界面兩側(cè)新增兩個未知變量，從而無法簡單通過共同的方程矩陣進行求解。

針對以上兩個難題，本文提出一種通用的方程切換方法，它適用于Picard和Newton-Raphson等所有主流迭代求解算法，并能根據(jù)每個節(jié)點的水分狀態(tài)進行方程的自由選擇：當節(jié)點飽和度大于臨界閾值時，其控制方程選為h型Richards方程；否則切換為θ型。為實現(xiàn)本文方法，采用廣義θ型方程解決土壤非均質(zhì)問題，并采用h型方程解決土壤飽和問題；同時，采用隱式數(shù)值格式來進行交界面上不同控制方程的節(jié)點之間的水量平衡分析。本研究以Picard迭代算法為例，成功將該方法應(yīng)用到一維飽和-非飽和水流運動及溶質(zhì)運移數(shù)值模型方法中。通過室內(nèi)和數(shù)值試驗，驗證了本文模型相比傳統(tǒng)方法（HYDRUS-1D）的優(yōu)勢，并論述了該方法在區(qū)域飽和-非飽和土壤水分運動及溶質(zhì)運移模擬中所具有的應(yīng)用前景。

2 飽和-非飽和水流運動及溶質(zhì)運移模型

2.1 水流運動方程不考慮土壤骨架及水流的可壓縮性，一維飽和-非飽和土壤水流運動須滿足質(zhì)量守恒方程：?θ/?t=-?q/?z+s。式中，θ為土壤體積含水量，cm3· cm-3；t為時間，min；z軸以向下為正，cm；s為源匯項，min-1；q為達西通量，cm·min-1，可以用壓力水頭h或含水量θ兩種形式表達：q=-K·?h/?z+K=-D·?θ/?z+K，式中，K為土壤水力傳導(dǎo)度，cm ·min-1；D為土壤水力擴散度，cm2·min-1，D=K/C，其中C=?θ/?h為土壤容水度，cm-1。用Φ表示方程切換時采用的主變量（θ或h），則Richards方程的通用表達式為：

當Φ=h時，當Φ=θ時，。式（1）的初始條件可表示為Φt=0=Φ0，其一類（Dirichlet）邊界表示為同時，二類（Neumann）邊界表示為其中，Φ0為主變量在初始時刻的取值，Φ1為邊界上的指定主變量解；qw為邊界上的指定通量。

2.2 溶質(zhì)運移方程一維對流-彌散方程可以表示為：

式中，J=-θDz·?c/?z+qc為溶質(zhì)通量，g·min-1· cm-2；c為溶質(zhì)濃度，g·cm-3；R=1+ρκηcη-1/θ為遲滯系數(shù)，無量綱；Λ=μwθ+μsρκcη-1為一階反應(yīng)項，g · cm-3· min-1；Γ=γwθ+γsρ-scs為零階反應(yīng)項，min-1；Dz為縱向彌散系數(shù)，cm2· min-1，θDz=DL|q|+θDdτ。ρ為土壤容重，g · cm-3；μw、μs分別為一階反應(yīng)項的液相和固相反應(yīng)速率，min-1；γw、γs分別為零階反應(yīng)項的液相（g·cm-3·min-1）和固相（min-1）反應(yīng)速率；s為式（1）中的源匯項，min-1；cs為根系吸水帶走的溶質(zhì)濃度，mg·cm-3；DL為軸向機械彌散系數(shù)，cm；Dd為離子（或分子）擴散系數(shù)，cm2·min-1；τ為孔隙扭曲度，無量綱，τ=κ為用于線性擬合溶質(zhì)在溶解-結(jié)晶平衡相態(tài)的經(jīng)驗常數(shù)，cm3·g-1；η為Freundlich等溫吸附經(jīng)驗系數(shù)，當η=1時為線性吸附，否則為非線性。本文僅討論線性吸附的情況，因而式（2）為線性方程。初始條件表示為c（x，t）|t=0=c0（x），一類（Dirichlet）邊界表示為二類（Neumann）邊界表示為其中，c0為初始溶液濃度，c1為邊界Ss1上的指定溶質(zhì)濃度，Q為邊界Ss2上的指定溶質(zhì)通量。

3 方程切換方法

3.1 方程切換臨界閾值傳統(tǒng)主變量轉(zhuǎn)換方法往往采用若干個經(jīng)驗系數(shù)來判斷節(jié)點水分狀態(tài)是否滿足主變量轉(zhuǎn)換的要求［15］。本文提出用臨界有效飽和度Secrit作為閾值來切換每個節(jié)點的控制方程。土壤在Secrit時，具有最大容水度C，即當土壤水分變化Δθ時，其負壓變化Δh最小，由土壤水分條件強烈非線性變化而引入的數(shù)值不穩(wěn)定最小，此時滿足?C/?h=0。

以van Genuchten模型［3］為例，定義土壤參數(shù)p，它包括土壤殘余含水量θr（cm3cm-3），土壤飽和含水量θs（cm3·cm-3），表征土壤級配分布的參數(shù)α（cm-1）和無量綱參數(shù)n和m=1-1/n，以及飽和土壤水力傳導(dǎo)度ks（cm·min-1）。在節(jié)點飽和度為Secrit時，壓力水頭為hcrit|?C/?h=0=-m1-m/α，從而得到方程切換臨界閾值Secrit=（1+m）-m。其中，m和α對應(yīng)為節(jié)點上的土壤參數(shù)。在每一個時間步Δtj內(nèi)，當節(jié)點飽和度小于Secrit時，其控制方程為θ型方程；反之，則選擇h型方程。當不存在極大容水度時（例如Brooks-Corey模型［17］等），可以根據(jù)土壤類型進行經(jīng)驗取值，例如Secrit=0.4～0.9。

3.2 混合單元水量均衡分析對于相鄰節(jié)點n和n+1（如圖1），當（Sen-Secri）t（Sen+1-Secri）t＜0時，二者控制方程將不再一致，將該節(jié)點間的單元稱為“混合單元”，編號為n+1/2?；旌蠁卧目刂乒?jié)點控制方程分別為θ和h型。由于兩類Richards方程離散節(jié)點間水流通量存在內(nèi)在差異，本文采用和表示節(jié)點通量的兩種形式，將二者的加權(quán)平均作為混合單元內(nèi)的節(jié)點間通量

式中，Δzn+1/2為單元n+1/2的空間長度。式（4）是主變量轉(zhuǎn)換和方程切換算法中，單元內(nèi)節(jié)點間通量的通用表達式：當每個單元頂點為θ型方程時，ω=0；當均為h型方程時，ω=1；當兩個頂點控制方程不一致時，0＜ω＜1。本文采用算術(shù)平均方案，取ω=0.5。

圖1 節(jié)點n和n+1的方程切換原則及單元n+1/2內(nèi)節(jié)點間通量示意圖

特別地，基于Newton-Raphson迭代的傳統(tǒng)主變量轉(zhuǎn)換方法［14］可以認為是式（4）的一種特例，其滿足ω=1，且各水力參數(shù)及未知變量均用主變量的一階偏導(dǎo)來表示，例如其中，上標j和k為時間和迭代層數(shù)（下同），其數(shù)學(xué)推導(dǎo)不再贅述。為克服前述Newton-Raphson迭代在處理主變量轉(zhuǎn)換時存在的數(shù)值不穩(wěn)定，本文從Picard迭代出發(fā)，采用一種半隱格式來求解式（4）。便于簡潔，用來表示節(jié)點的導(dǎo)水度。

由于非均質(zhì)界面上節(jié)點含水量不連續(xù)且不可導(dǎo)，傳統(tǒng)θ型Richards方程及式（4）均不適用于處理節(jié)點n或n+1處于土壤分層界面的情況，因此需要在節(jié)點通量qn+1/2上引入一個非均質(zhì)修正項［13］：

當n節(jié)點為h型方程，n+1節(jié)點為θ型方程時，單元n+1/2內(nèi)的節(jié)點間通量表示為：

反之，當n+1節(jié)點為h型方程，n節(jié)點為θ型方程時，該通量表示為：

式（6）和式（7）的差異在于混合單元兩端未知變量h和θ的顯隱格式，可視為Picard迭代算法框架下的混合單元通量表達式。進而我們可以簡單地推導(dǎo)出，混合單元n+1/2的Richards方程的空間離散表達式為：為節(jié)點n的控制域長度。該式與無切換單元的空間離散格式（參見［15］）保持了一致。

而Richards方程的時間項離散則依據(jù)節(jié)點方程類型進行：當節(jié)點n為θ型方程時，當節(jié)點n為h型方程時，采用Celia格式［10］離散時間項：其中，Δtj+1=tj+1-tj為當前時間步長。此外，模型的時間步長自適應(yīng)和溶質(zhì)運移方程求解等數(shù)值算法均采納HYDRUS-1D［18］方案，代碼編寫基于HYDRUS 5.0軟件。

4 模型驗證及分析

為體現(xiàn)本文方法的數(shù)值優(yōu)勢，以經(jīng)典室內(nèi)試驗來驗證本文模型的效率與精度；以土壤水鹽運移數(shù)值模擬的兩個傳統(tǒng)方法難以解決的普遍難題為例，證明了本文方法的優(yōu)勢。所有算例以HYDRUS-1D軟件的高密度時間和空間網(wǎng)格解Φref為參照解（Φ表示含水量θ、壓力水頭h或溶質(zhì)濃度c），得到模型解Φ的均方根誤差：

以Em表示不同模型的相對質(zhì)量誤差：

以矩陣求解次數(shù)作為通用指標來對比本文方法相比HYDRUS-1D商業(yè)軟件的計算成本。所有土壤水流運動及溶質(zhì)運移參數(shù)見表1。

4.1 入滲-排水試驗為測試本文模型求解水流過程的有效性及模型穩(wěn)定性，采用經(jīng)典的Abeele［19］土柱入滲-排水試驗進行數(shù)值驗證。試驗采用的一維土柱直徑為300 cm，高度為600 cm，填充均質(zhì)的Bandelier火山凝灰?guī)r（見表1（土壤#1））。試驗包括兩步：（1）以定水頭（htop=0）上邊界進行極度干燥條件的土壤強入滲，初始條件為h（z，0）=－729.7 cm，模擬時長5 d；（2）對灌滿后的土柱進行為期100 d的自由排水試驗觀測，上邊界密封，設(shè)為零通量，初始條件為h（z，0）=0，下邊界為（?h/?z）bot=0。以網(wǎng)格密度Δz=5 cm、Δtmax=0.5 d的HYDRUS-1D和本文模型解進行對比。其中，HYDRUS-1D軟件求解方程 次，本文模型求解方程 次，計算成本節(jié)約19%。分析其壓力水頭和含水量剖面（見圖2a），認為模型水流計算的求解精度較高且數(shù)值穩(wěn)定。試驗排水階段水流條件較為溫和，采用本文模型及參照模型解對比觀測值，得到圖2（b）所示土壤含水量剖剖面隨時間的變化。本文模型求解非穩(wěn)定水流過程的精度較高，對比高密度時空網(wǎng)格參照解θref，其RMSE（θ，t）＜0.002 cm3cm-3。

表1 土壤水流運動及溶質(zhì)運移參數(shù)

圖2 含水量剖面隨時間的變化

4.2 鎂離子（Mg2+）穿透試驗為驗證本文模型溶質(zhì)運移模塊的計算精度，以經(jīng)典的（a）Selim非線性吸附試驗［20］和（b）Lai&Jurinak線性吸附試驗［21］為參照，進行穩(wěn)定流狀態(tài)下的溶質(zhì)運移模擬。相關(guān)參數(shù)見表1（土壤#2和#3）。以HYDRUS-1D軟件的高密度網(wǎng)格（Δz=0.025 cm、Δt=0.1 h）作為參照解。試驗（a）采用0.271 cm/h穩(wěn)定通量的飽和CaCl2溶液（5 mmol/L）的10.75 cm土柱，在模擬初始的358.05小時內(nèi)，注入14.26倍孔隙體積量的MgCl2溶液脈沖（5 mol/L），隨后恢復(fù)注入原始濃度的CaCl2溶液。采用非線性吸附函數(shù)模擬鎂離子的交換過程，整個試驗持續(xù)600小時。試驗（b）則采用2.64 cm/h穩(wěn)定通量的CaCl2溶液（125 mmol/L）的25 cm土柱，在初始的3.9小時內(nèi)，同時注入兩份各62.5 mmol/L的MgCl2和CaCl2溶液，隨后恢復(fù)注入原始濃度的CaCl2溶液。采用線性吸附函數(shù)模擬鎂離子的交換過程，整個試驗持續(xù)30 h。兩個鎂離子穿透試驗的數(shù)值模型上邊界為定溶質(zhì)通量邊界（二類邊界）：下邊界為零梯度自由邊界：在土柱底部進行鎂離子濃度檢測，得到如圖3所示，本文模型的溶質(zhì)運移模塊計算結(jié)果準確，RMSE（c，t）＜0.005 mmol/L。

4.3 砂土淋鹽過程模擬本節(jié)參考Miller經(jīng)典算例［22］，模擬干旱灌區(qū)厚包氣帶淋鹽過程，設(shè)置均質(zhì)土柱高度為1000 cm，土壤參數(shù)見表1（#2）。初始條件為靜水壓狀態(tài)，潛水面位于柱子底部（即h（z，0）=z-1000 cm）。在地表0～10 cm埋深處，存在均勻分布的高濃度鹽分20 g/L，其它埋深鹽分背景值為0.5 g/L。模擬過程中，上邊界以定水頭htop=0 cm進行淋鹽，灌溉水質(zhì)0.5 g/L；下邊界為定水頭hbot=10 cm，地下水鹽分濃度為2 g/L，全部淋鹽過程持續(xù)260 min。主要溶質(zhì)運移參數(shù)見表1中土壤#4。算例在t=80 min、160 min、260 min處獲取觀測剖面，對比高密度網(wǎng)格（Δz=0.25 cm）的HY?

圖3 鎂離子穿透曲線

圖4 壓力水頭、含水量和溶質(zhì)濃度剖面隨時間的變化

DRUS-1D模型解，得到如圖4的（a）壓力水頭、（b）含水量和（c）鹽分濃度剖面隨時間的變化情況。

采用網(wǎng)格分辨率為Δz=0.5 cm、1.0 cm、1.25 cm和2.0 cm的本文模型及相應(yīng)網(wǎng)格密度的HY?DRUS-1D解對比高密度網(wǎng)格參照解，得到圖5所示鹽分濃度、土壤含水量剖面解的RMSE隨時間的變化?？傮w上，本文模型能保證不同網(wǎng)格分辨率上更高的求解精度。對于網(wǎng)格密度為Δz=0.5 cm、1.0 cm、1.25 cm和2.0 cm的情況，本文模型相比HYDRUS-1D軟件，分別節(jié)省了16.8%、21.2%、12.8%和7.1%的計算成本，壓力水頭解的求解誤差降低了8%～79%，含水量求解誤差降低了3.5%～60%，溶質(zhì)求解誤差降低了9.6%～43%。因而在保證更高計算精度的同時，能實現(xiàn)更高效率的求解。以Δz=1.0 cm為例，HYDRUS-1D軟件則需要求解方程10499次，總體質(zhì)量誤差為0.234%；而本文模型求解方程次數(shù)為8270，質(zhì)量誤差為7.25×10-5%，節(jié)省了21%以上的計算成本，降低了99.97%的質(zhì)量誤差。從算法穩(wěn)健性來看，本文模型能適應(yīng)更大的時間步長，及更高密度的網(wǎng)格。自適應(yīng)時間步長方案下，本文模型最大時間步長是HYDRUS-1D軟件的1.20～1.44倍，最小時間步長則為后者102～104倍，從而明顯節(jié)省計算工作量；而在更高密度網(wǎng)格下，本文模型能更大幅度地提高計算精度。

4.4 干濕交替的地表積鹽過程模擬非均質(zhì)土壤的水流運動過程具有強烈非線性，通常對數(shù)值工具的穩(wěn)定性及計算效率有更高的要求。傳統(tǒng)基于θ型Richards方程的方法往往因為無法準確計算飽和土壤及非均質(zhì)界面處的節(jié)點間通量而少見應(yīng)用；同時，基于h型Richards方程的模型，例如HYDRUS和SWAP系列軟件，往往在干濕交替的大氣邊界上具有較低的計算效率甚至無法正確求解，此外容易造成較大的質(zhì)量誤差［13］。針對這個問題，本節(jié)參照Hills經(jīng)典問題［23］，進行干濕交替的大氣邊界下，我國北方干旱-半干旱灌區(qū)非均質(zhì)土壤的積鹽過程模擬。將高度為1 m的土柱均勻分為5層（每層20 cm），自上而下，在第1、3、5層中填充均質(zhì)細砂土（見表1中土壤參數(shù)#5），在2、4層中填充均質(zhì)黏壤土（見表1中土壤參數(shù)#6）。在初始狀態(tài)下，土柱水頭均勻分布，h（z，0）=-100 cm，初始鹽分背景濃度為c（z，0）=1.5 g/L。下邊界指定水頭和鹽分濃度：hbot=0、cbot=1.5 g/L；考慮復(fù)雜的干濕交替大氣邊界（見圖6）。上邊界溶質(zhì)輸入量近似設(shè)為0.5 g/L。算例對比相同網(wǎng)格密度（Δz=1 cm）下的HYDRUS-1D軟件及本文模型的計算效率與精度。模擬時長365天。由于高密度網(wǎng)格時HYDRUS-1D軟件難以收斂，因而本算例以相同網(wǎng)格密度下的HYDRUS-1D軟件解作為參照解。

圖5 對比砂土淋鹽數(shù)值求解剖面溶質(zhì)濃度、含水量的RMSE

圖6 干濕交替的大氣邊界

在干旱地區(qū)，大量的潛在騰發(fā)將造成土壤積鹽（見圖7（a）），尤其是地下水位較淺的情況。本文以虛擬算例為例，模擬了這種積鹽過程。圖7（b）為位于z=10 cm（砂土層）的觀測點上，鹽分的變化過程。當出現(xiàn)短時間的降雨或灌溉入滲時，土壤鹽分得到迅速淋洗，但由于強烈、持續(xù)的騰發(fā)作用，這種積鹽過程再次發(fā)生。通常認為，干濕交替的大氣邊界模擬容易造成基于h型Richards方程的數(shù)值模型（如HYDRUS、SWAP等）出現(xiàn)較大幅度的數(shù)值振蕩。尤其當土壤較干燥時，少量的水分輸入往往導(dǎo)致壓力水頭在多個數(shù)量級上發(fā)生變化，因而模型需要更多的迭代來完成一次求解。但對于θ型方程而言，干濕交替使得含水量θ在一個時間步內(nèi)的變化幅度不超過θs～θr。本文模型以10 s完成計算（方程求解9293次），而HYDRUS-1D則耗時606 s（方程求解699698次），本文模型節(jié)省98%以上的計算成本。值得一提的是，為獲取更高精度的數(shù)值解，本文方法能高效求解更高密度的數(shù)值網(wǎng)格，例如，Δz=0.1 cm時，計算耗時107 s，而HYDRUS-1D的計算則因不收斂而中斷。對比相對質(zhì)量守恒誤差，本文模型（1.4×10-3%）相比HYDRUS-1D軟件（7.5%）而言，質(zhì)量誤差減小99.98%。

圖7 鹽分濃度隨時間的變化

5 結(jié)論

針對傳統(tǒng)土壤水分運動數(shù)值模型存在數(shù)值振蕩、質(zhì)量誤差大等潛在問題，本文提出一種通用的控制方程切換技術(shù)，并成功應(yīng)用于一維非均質(zhì)飽和-非飽和水流運動及溶質(zhì)運移數(shù)值模擬中。通過一系列室內(nèi)及數(shù)值試驗，認為本文模型能獲得有效的數(shù)值解。對比通用軟件HYDRUS-1D，本文模型經(jīng)論證實現(xiàn)了數(shù)值精度的普遍提升。采用經(jīng)典數(shù)值算例，討論了傳統(tǒng)數(shù)值方法難以解決的兩個普遍問題：（1）在干燥砂土淋鹽問題中，本文模型能在更大的時間步長下克服傳統(tǒng)水流模型難收斂的問題，并證明了基于方程切換技術(shù)的數(shù)值模型在計算精度上有顯著提升。其中，計算成本節(jié)省7.1%～21.2%，壓力水頭誤差減小8%～79%，含水率誤差減小3.5%～60%，溶質(zhì)誤差減小9.6%～43%，質(zhì)量誤差減小99.98%。（2）在大氣邊界干濕交替條件的非均質(zhì)土壤積鹽問題中，本文模型能大幅降低計算成本，并能在更高密度的網(wǎng)格下快速獲取高精度解（此時HYDRUS-1D軟件難以收斂）。本文提出的通用方程切換技術(shù)能應(yīng)用于二維、三維水流運動及溶質(zhì)運移模型中；對于大區(qū)域水流運動與溶質(zhì)運移模擬，本文模型能實現(xiàn)高效率數(shù)值求解。