蔡松柏 邵明安 齊麗彬

摘 要：研究地下水中溶質的遷移規(guī)律可以為地下水環(huán)境質量預報、控制和管理服務。首次采用數(shù)值方法探討與位置坐標相關的溶質運移參數(shù)模型，采用即時時間域顯式差分法格式雙精度Fortran語言編程計算了變溶質運移參數(shù)模型的一維溶質運移問題，并進行了算法可靠性驗證;然后對本文首次提出的3種與位置坐標相關的變溶質運移參數(shù)模型進行了參數(shù)分析，獲得了這3種模型的初步認識。還將時間域顯式差分方程計算結果與解析解和實驗結果進行了驗證對比分析，說明算法可靠，顯式法列式簡單，更適于復雜模型計算。本文提出的算法和處理技巧可以應用到許多與地質相關的溶質運移問題的研究，通過建立不同模式的溶質運移參數(shù)模型，從而實現(xiàn)復雜條件下溶質運移過程的預測和模擬。

關鍵詞：一維溶質遷移;與深度相關;變溶質運移參數(shù);建模與模擬;即時顯式差分法

中圖分類號：O357.3

文獻標志碼：A

由于工農業(yè)生產、雨水徑流沖刷和人為因素等作用，造成水流中溶質瞬間注入或連續(xù)注入從而引起河流或地下水的水質改變。研究溶質瞬時注入或連續(xù)注入情況下不同時刻溶質濃度在河流或地下水流中的分布狀況、及不同位置時溶質濃度隨時間的變化情況，可通過數(shù)學模型來模擬或預報水質在時間或空間的變化，達到為水環(huán)境質量預報、控制和水資源管理服務的目的[1-22]。

數(shù)值模擬水流溶質運移的主要方法有有限元法、邊界單元法、有限差分法和其他數(shù)值法等。由于進行水質模擬過程中涉及的參數(shù)在不同的時間和空間上變化較大，在水質模擬過程中，會造成數(shù)值彌散、數(shù)值波動和計算量巨大等問題。國內主要的研究有徐玉佩[2-3]用Galerkin法推導了二維地下水運動的有限元方程和溶質運移有限元方程，并編制了計算程序應用于明溝排水和田面保持淹灌水層情況下的沖洗種稻改良鹽堿地的分析，還與現(xiàn)有解析解和實驗資料進行了對照。徐玉佩還用水動力彌散理論進行了地下水中水鹽動態(tài)的研究，可直接用于鹽堿地的預防和改良。黃康樂[4]則提出了一種求解二維飽和-非飽和溶質運移問題的交替方向特征有限單元法（ADCG），在克服數(shù)值彌散、數(shù)值波動和提高計算速度方面具有顯著的優(yōu)越性。張效先等[5]提出了一種計算地下水中溶質運移的數(shù)值法，將多維問題分解成幾個一維問題， 減少了計算工作量， 并消除了數(shù)值波動現(xiàn)象，具有較高的精度。馮紹元[6]提出了特征有限元法的基本理論和實施過程，其計算實例表明能消除大部分的數(shù)值彌散而又不增加解的振蕩性。姚磊華[7]首先基于分步求解的思想， 用廣義迎風對偶單元均衡法求解對流定解問題， 對擴散定解問題采用一般的Galerkin有限元法求解， 不僅避免了數(shù)值彌散和過量問題， 而且避免了求節(jié)點速度這一步， 簡化了運算步驟。邱克儉等[8]結合拉格朗日觀點的運動坐標與歐拉觀點的固定坐標，分別求解對流-彌散方程的對流效應和彌散效應， 并通過與理論解和試驗結果做比較，驗證了模型的可靠性。任理等[9]采用混合拉普拉斯變換有限單元法求解對流占優(yōu)的地下水溶質運移問題，能有效地消除數(shù)值擴散和過量現(xiàn)象，特別適于大區(qū)域地下水污染的長期預報。李煥榮等[10-11]利用廣義差分法建立了一維非飽和水流問題的守恒形式的數(shù)值模型， 具有計算量小和穩(wěn)定性好的特點， 其提出的非粘性土壤水中溶質運移問題的守恒混合元格式，精度高且數(shù)值穩(wěn)定。段德宏等[14]以黃河蘭州段為例，采用一維穩(wěn)態(tài)河流水質模型，對岸邊19個污水排放口及共約32.9 km長的全河段污染狀況成功地進行了模擬。徐文彬等[15]根據(jù)數(shù)值隨時間步長變化的特點，提出了針對數(shù)值波動、彌散和過量現(xiàn)象的新改進方案，并用實例驗證了方法的可行性。何丕文[16]利用一維水質模型，采用有限差分法求解河流中污染物濃度，以淮河淮南段為例進行了驗證。馬東豪等[17]研究了土壤溶質遷移的兩流區(qū)與兩區(qū)模型，結果表明與兩區(qū)模型相比，兩流區(qū)模型可以更好地描述優(yōu)先流情況下的溶質穿透曲線。對于無優(yōu)先流情況， 連續(xù)流輸入情況下兩個模型的參數(shù)一致性較好;而在脈沖輸入時一致性較差。國外在溶質運移模型及其數(shù)值計算方法上的研究也很多，如JAISWAL等[18]運用Chebyshev配點法求解一維溶質運移問題的數(shù)值解。SILAVWE等[19]討論了某些數(shù)值方法用于一維溶質運移模型的參數(shù)估計，通過將數(shù)值解與試驗結果比較，反演溶質運移模型參數(shù)。JAISWAL等[20]采用算子矩陣法對多孔介質的非線性微分方程進行了數(shù)值分析，該法應用切比雪夫多項式并結合切比雪夫導數(shù)運算矩陣和譜配置法對解進行逼近。其優(yōu)點是可以將這些問題轉化為易于求解的代數(shù)方程組。CARR[21]利用拉普拉斯變換導出了層狀多孔介質溶質運移的新半解析解，該法在相鄰層的界面引入表示運移量的未知函數(shù)，將多層問題置于拉普拉斯域上分層解決，然后再將其數(shù)值反演回時域。ALLWRIGHT等[22] 對地下水溶質運移方程在對流明顯強于擴散時出現(xiàn)的數(shù)值波動提出了兩種新的數(shù)值近似計算方法，即只適用于對流項的迎風Crank-Nicolson格式和加權迎風-下風格式，并對這些新提出的格式進行了數(shù)值穩(wěn)定性分析和比較。發(fā)現(xiàn)如果Crank-Nicolson格式只用于對流項時很合宜。此外，顯式加權迎風-下風有限差分格式是對傳統(tǒng)顯式一階迎風格式的改進，而隱式加權一階迎風-下風有限差分格式在給定適當?shù)臋嘀匾蜃訒r無條件穩(wěn)定。

本文主要貢獻在于對一維溶質對流彌散基本方程，首先，建立了一維溶質運移3種變溶質運移參數(shù)模型，即線性變化模型/多項式變化模型、指數(shù)變化模型和冪變化模型;其次，采用顯式時域差分法成功地進行了溶質濃度時空分布的數(shù)值計算，有效克服了很多研究人員遭遇的數(shù)值波動、數(shù)值彌散和計算量大等難題，并與解析解和實驗結果做了對比，驗證了該數(shù)值方法的正確性;最后，對3種溶質運移參數(shù)模型進行了參數(shù)分析，獲得了3種模型的初步認識。

1 基本原理

研究溶質在水流中的遷移規(guī)律，試驗測定可采用電解質脈沖示蹤法等來進行。電解質脈沖在水流中的運移規(guī)律符合溶質對流彌散基本方程，根據(jù)夏衛(wèi)生等[12-13]的研究，其溶質運移參數(shù)隨水流位置而變化，在大規(guī)模的地下水流問題中，溶質運移參數(shù)亦會隨空間位置而變化。對其溶質運移規(guī)律的模擬可采用近似解析解和數(shù)值解。

3 計算結果與討論

根據(jù)上面的理論基礎，利用雙精度FORTRAN語言編寫程序，分別采用方程（19）所代表的數(shù)值方法和方程（18）所代表的解析方法求得溶質濃度分布，以檢驗數(shù)值方法的可靠性。并對不同溶質運移參數(shù)變化模式下瞬時注入溶質后水流中溶質濃度的時空分布進行了數(shù)值模擬和參數(shù)分析。

3.1 一維溶質運移參數(shù)線性變化時方程數(shù)值解與解析解比較

由夏衛(wèi)生等[12-13]的試驗結果得知，坡面薄層水流的溶質擴散系數(shù)及水流速度均隨位置而變化。本計算的變溶質運移參數(shù)模型采用D（x）=D0（1+bx/l），u（x）=u0（1+bux/l）。根據(jù)夏衛(wèi)生等[12-13]的試驗數(shù)據(jù)，取D0=5.1×10-2m2/s，b=-0.55，u0=1.02 m/s和bu=0.042。l=3.5 m，t0=0.6 s， 注入溶液濃度c1=0.63，水流起始溶質濃度和平衡濃度c0=0.20，Δt=0.04 sec和Δx=0.025 m，分別求得了水流方向上不同位置處溶質濃度的數(shù)值解和近似解析解。

圖1給出了數(shù)值解、近似解析解與文獻[12]實測結果的比較，從圖中可以看出兩種模擬結果同夏衛(wèi)生等[12-13]的結果一致，測量曲線和模擬曲線的溶質濃度最大值均呈逐漸減小的趨勢;在距離較

大處其模擬結果和實測結果吻合較好。近似解析解和數(shù)值解能較好地吻合，特別是距離越大，二者吻合得越好。這說明顯式差分法求解溶質運移方程的精度已經能夠達到要求。

表1給出了一維常溶質運移參數(shù)模型在脈沖注入后溶質濃度隨深度、時間變化的數(shù)值解式（20）與解析解式（18）的比較，其中，常溶質運移參數(shù)D0=5.1×10-2m2/s，u0=1.02 m/s，其余計算參數(shù)同圖1取值，從表中可以看出二者吻合非常好，說明了數(shù)值方法是可靠的。

由圖1和表1可知，常溶質運移參數(shù)一維溶質運移基本方程的解析解與變溶質運移參數(shù)一維溶質運移基本方程的數(shù)值解基本一致。因此，求解溶質運移基本方程可以采用本文數(shù)值方法。

3.2 一維溶質運移參數(shù)多項式變化時的數(shù)值結果

給定D0=0.03 m2/s，u0=1.00 m/s，計算水流位置l=3.5 m，脈沖注入溶質濃度時間t0=0.6 s，注入溶液濃度c1=1.00，水流起始溶質濃度和平衡濃度c0=0.20，計算步長Δt=0.04 sec和Δx=0.025 m，當式（8）和（9）中模型參數(shù)c=0.1，b=±0.6，±0.4，±0.2，對cu，bu和c，b同步一樣取完全相同值時，其數(shù)值解法求出的溶質濃度沿水流位置的分布規(guī)律如圖2所示，圖2中c和b同時代表cu，bu和c，b的取值。從圖2可以看出，c，cu的零值附近微小變化對溶質濃度分布影響較計算方法帶來的差別要小。c，cu值一定時，隨著b，bu從-0.6到0.6的增加，某時刻溶質濃度位置上的分布，其峰值和峰寬增加，波峰位置滯后。時間越長，這種濃度分布的波峰變化規(guī)律越明顯。在1 s時間以內，本文對b，bu 和c，cu 的不同取值對溶質濃度分布影響不大。

3.3 一維溶質運移參數(shù)冪函數(shù)變化時的數(shù)值結果

溶質運移參數(shù)采用式（10）和（11）的冪函數(shù)變化模型，給定模型參數(shù)b=2且c=±0.9，±0.6，±0.3下，cu，bu和c，b取值完全相同，其余計算參數(shù)D0，u0，l，t0， c1， c0，Δt和Δx同以上多項式模式中參數(shù)的取值，一維溶質運移濃度分布數(shù)值解如圖3所示。圖3表明，參數(shù)中指數(shù)b/c，bu/cu較大時，c，cu值的大小對溶質濃度分布的影響較小，時間越短，這種影響越小。

3.4 一維溶質運移參數(shù)指數(shù)變化時的數(shù)值結果

圖4所示為溶質運移參數(shù)模型采用式（12）和（13）模型的數(shù)值計算結果，計算中取模型參數(shù)b=0.125且c=±0.9，±0.6，±0.3下，cu，bu和c，b同步取完全相同的值，其余計算參數(shù)D0，u0，l，t0，c1，c0，Δt和Δx同以上多項式模式中參數(shù)的取值。從圖4看出模型參數(shù)的變化對濃度分布形態(tài)的影響較小，對于本例，當b，bu一定時，c，cu的增大會導致濃度分布峰值位置滯后，濃度峰值增加。c，cu一定時，b，bu在一定范圍內的變化對溶質濃度分布的影響很小。

4 結論

本文采用顯式時間差分格式計算變溶質運移參數(shù)模型一維水流溶質運移過程，計算結果通過解析解進行驗證，并和試驗結果進行對比分析，說明算法可靠，顯式法列式簡單，更適于復雜模型計算[1];模型越復雜，本法的計算效率越高。本文首次進行了變溶質運移參數(shù)模型溶質運移問題的數(shù)值計算，并建立不同的溶質運移參數(shù)模型，對該模型進行了參數(shù)分析。提出的系列算法和處理技巧可以應用到溶質注入水流后，變溶質運移參數(shù)模型的溶質遷移規(guī)律研究，通過建立不同模式下溶質運移參數(shù)模型，從而實現(xiàn)復雜條件下溶質運移過程的預測和模擬。

參考文獻：

[1] 邵明安， 王全九， 黃明彬.土壤物理學[M]. 北京： 高等教育出版社， 2006.

[2] 徐玉佩. 地下水中溶質運移數(shù)值計算（有限元法）及其在沖洗種稻改良鹽堿地中的應用[J]， 水利學報， 1983，14（4）：1-15.

[3] 徐玉佩. 地下水中溶質運移數(shù)值計算的沿流線分段追蹤法[J]， 水利學報， 1986， 17（12）：38-47.

[4] 黃康樂. 求解二維飽和-非飽和溶質運移問題的交替方向特征有限單元法[J].水利學報， 1988，19（7）：1-7.

[5] 張效先， 張志剛， 李法虎. 計算地下水中溶質運移的時間分裂Eulerian-Lagrangian有限差分法[J].水利學報，1991，22（8）：35-45.

[6] 馮紹元. 解地下水中溶質運移問題的特征有限元法[J].武漢水利電力學院學報， 1991，24（4）：410-417.

[7] 姚磊華. 三維溶質運移問題的分步廣義迎風解法[J].長春地質學院學報，1997，27（4）：424-428.

[8] 邱克儉， 戚隆溪， 秦寧. 飽和非飽和土壤中溶質運移的數(shù)值模擬[J].力學學報， 1993，25（1）：40-47.

[9] 任理，李春友，李韻珠. 粘性土壤溶質運移新模型的應用[J].水科學進展，1997，8（4）：321-328.

[10]李煥榮， 羅振東， 謝正輝， 等. 非飽和土壤水流問題的廣義差分法及其數(shù)值模擬[J].計算數(shù)學，2006，28（3）：321-336.

[11]李煥榮， 羅振東. 非粘性土壤中溶質運移問題的守恒混合有限元法及其數(shù)值模擬[J].計算數(shù)學， 2010，32（2）： 183-194.

[12]夏衛(wèi)生， 雷廷武， 趙軍，等. 薄層水流速度測量系統(tǒng)的研究[J]. 水科學進展， 2003， 14（6）： 781-784.

[13]夏衛(wèi)生， 雷廷武， 張晴雯，等. 坡面薄層水流中電解質脈沖遷移模型[J].水利學報，2003，34 （11），90-95.

[14]段德宏， 王根霞， 王萍. 一維穩(wěn)態(tài)河流水質的計算機模擬[J]. 水資源與水工程學報， 2005， 16（3）： 54-56.

[15]徐文彬， 徐維生， 方濤. 溶質運移數(shù)值模擬彌散問題再探[J]. 災害與防治工程， 2007， 1（6）： 56-61.

[16]何丕文. 有限差分法在河流水質預測中的應用[J]. 長江大學學報（自然科學版），2006，3（1）： 38-40.

[17]馬東豪，王全九. 土壤溶質遷移的兩區(qū)模型與兩流區(qū)模型對比分析[J]. 水利學報，2004，35（6）：92-97.

[18]JAISWAL S，CHOPRA M，SENG-HUAT O，et al. Numerical solution of one-dimensional finite solute transport system with first type source boundary condition[J].International Journal of Applied and Computational Mathematics，2016，3（4）：1-11.

[19]SILAVWE D D，BRINK I C，WALLIS S G. Assessment of some numerical methods for estimating the parameters of the one-dimensional advection-dispersion model[J]. Acta Geophysica，2019，67（3）：999-1016.

[20]JAISWAL S，CHOPRA M. Numerical solution of non-linear partial differential equation for porous media using operational matrices[J]. Mathematics and Computers in Simulation，2019，160：138-154.

[21]CARR E J. New semi-analytical solutions for advection-dispersion equations in multilayer porous media[J].Transport in Porous Media，2020，135（1）：39-58.

[22]ALLWRIGHT A， ATANGANA A. Augmented upwind numerical schemes for the groundwater transport advection-dispersion equation with local operators[J].International Journal for Numerical Methods in Fluids， 2018，87（9）：437-462.

（責任編輯：曾 晶）

A Numerical Study on Depth-Relevant Solute

Advection-dispersion Model

CAI Songbai1， SHAO Mingan*2， QI Libin3

（1.College of Resources and Environmental Science， Hunan Normal University， Changsha 410081， China;2.Institute of Soil and Water Conservation， Chinese Academy of Science， Yangling 712100， China;3.Department of Biochemical Engineering， Sanmenxia Polytechnic， Sanmenxia 472000， China）

Abstract：

Research on modelling of solute transport behaviour in groundwater could serve for environmental quality prediction， controlling and management of groundwater. This work has for the first time studied depth-relevant solute transport parameter model by adopting numerical method. By coding a double precision FORTRAN program， this work has calculated the solution of one-dimensional solute transport model of changing solute transport parameter by instant explicit time domain difference method. Then parameter analysis has been done for the 3 changing solute transport models of groundwater and some conceptual knowledge of the model has been obtained.The results of explicit difference equation in time domain are verified and compared with the analytical solution and experimental results. The results show that the algorithm is reliable， the explicit formula is simple， and it is more suitable for complex model calculation. Finally， the 3 changing solute transport models and its numerical analysis methods presented can be applied for many fields such as 2-D and/or 3-D solute transport problem for modelling and prediction of complicated geological condition of solute transport in groundwater.

Key words：

one dimension solute transport; depth-relevant; changing solute transport coefficient; modelling and simulation; instant explicit difference method

收稿日期：2020-11-11

基金項目：國家自然科學基金資助項目（40371060）

作者簡介：蔡松柏（1962—），男，副教授，博士，博士后，研究方向：土壤物理學，E-mail：songbai_cai@sina.com.

通訊作者：邵明安，E-mail：mashao@ms.iswc.ac.cn.