曾季才，查元源，楊金忠

（武漢大學(xué) 水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430072）

1 研究背景

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

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

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

2 飽和-非飽和水流運(yùn)動(dòng)及溶質(zhì)運(yùn)移模型

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

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

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

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

3 方程切換方法

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

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

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

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

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

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

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

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

反之，當(dāng)n+1節(jié)點(diǎn)為h型方程，n節(jié)點(diǎn)為θ型方程時(shí)，該通量表示為：

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

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

4 模型驗(yàn)證及分析

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

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

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

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

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

圖2 含水量剖面隨時(shí)間的變化

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

圖3 鎂離子穿透曲線

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

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

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

4.4 干濕交替的地表積鹽過程模擬非均質(zhì)土壤的水流運(yùn)動(dòng)過程具有強(qiáng)烈非線性，通常對(duì)數(shù)值工具的穩(wěn)定性及計(jì)算效率有更高的要求。傳統(tǒng)基于θ型Richards方程的方法往往因?yàn)闊o法準(zhǔn)確計(jì)算飽和土壤及非均質(zhì)界面處的節(jié)點(diǎn)間通量而少見應(yīng)用；同時(shí)，基于h型Richards方程的模型，例如HYDRUS和SWAP系列軟件，往往在干濕交替的大氣邊界上具有較低的計(jì)算效率甚至無法正確求解，此外容易造成較大的質(zhì)量誤差［13］。針對(duì)這個(gè)問題，本節(jié)參照Hills經(jīng)典問題［23］，進(jìn)行干濕交替的大氣邊界下，我國(guó)北方干旱-半干旱灌區(qū)非均質(zhì)土壤的積鹽過程模擬。將高度為1 m的土柱均勻分為5層（每層20 cm），自上而下，在第1、3、5層中填充均質(zhì)細(xì)砂土（見表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。算例對(duì)比相同網(wǎng)格密度（Δz=1 cm）下的HYDRUS-1D軟件及本文模型的計(jì)算效率與精度。模擬時(shí)長(zhǎng)365天。由于高密度網(wǎng)格時(shí)HYDRUS-1D軟件難以收斂，因而本算例以相同網(wǎng)格密度下的HYDRUS-1D軟件解作為參照解。

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

圖6 干濕交替的大氣邊界

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

圖7 鹽分濃度隨時(shí)間的變化

5 結(jié)論

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