平原區(qū)飽和-非飽和土壤水運動模型及數(shù)值算法研究

陳景波，王船海，杜世鵬，楊 海，張夢菲，朱立國

(1.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇南京210098；2.浙江省水利水電勘測設(shè)計院，浙江杭州310002；3.江蘇省水文水資源勘測局無錫分局，江蘇無錫214000)

?

平原區(qū)飽和-非飽和土壤水運動模型及數(shù)值算法研究

陳景波1，王船海1，杜世鵬2，楊 海1，張夢菲1，朱立國3

(1.河海大學水文水資源與水利工程科學國家重點實驗室，江蘇南京210098；2.浙江省水利水電勘測設(shè)計院，浙江杭州310002；3.江蘇省水文水資源勘測局無錫分局，江蘇無錫214000)

以土壤水運動理論為基礎(chǔ)，構(gòu)建了研究區(qū)三維飽和-非飽和土壤水運動模型，通過網(wǎng)格剖分模擬研究區(qū)各個位置土壤的體積含水率以及變換過程，并采用多個算例驗證模型的穩(wěn)定性和準確性。計算結(jié)果可知，土壤水分及地下水水位變化過程符合土壤水運動原理和日常經(jīng)驗，重現(xiàn)了特定土壤水運動的基本過程，為探究平原區(qū)水循環(huán)過程、水分轉(zhuǎn)化關(guān)系以及進行大尺度平原區(qū)水文預報提供理論基礎(chǔ)。

土壤水；水循環(huán)；飽和-非飽和流；數(shù)值模擬

0 引 言

平原區(qū)土壤水是水文循環(huán)中的關(guān)鍵紐帶，運移交換頻繁[1- 3]。現(xiàn)階段研究土壤水運動主要通過耦合現(xiàn)有的地表水和地下水模型來模擬計算，連接方式多以地下水模型進行三角形或矩形差分構(gòu)建基本單元[4]。在計算過程中，兩模型之間實際上是分離的，地表水和地下水模型分別獨立進行計算，模擬實際情況下的產(chǎn)流和地下水動態(tài)變化，模型之間只是通過下滲量的傳遞來進行簡單、松散的耦合[5]，如SWATMOD[6]模型、MODFLOW-DAFLOW[7]模型等。這就導致模型不能準確地模擬土壤水，特別是平原區(qū)飽和-非飽和帶土壤水的運動過程，對耦合模型要重點解決的諸如潛水蒸發(fā)、土壤水變化及產(chǎn)流量等問題不能很好地解決。程勤波等[8]在可變坡土槽進行了室內(nèi)人工降水入滲試驗，以此來分析MODFLOW-UZF1模型模擬降水入滲補給過程的精度。結(jié)果表明，模型不能準確地模擬入滲過程中土壤含水量的漸變過程，忽略了非飽和帶土壤水分運移中基質(zhì)勢的影響，在地下水淺埋地區(qū)模擬效果較差。

本文基于三維理查德方程推導并構(gòu)建了適用于平原區(qū)三維飽和-非飽和土壤水流動模型，并在太湖流域的金壇市建立野外實測試驗基地，以此來觀測分析和研究地勢平坦、地下水埋深淺的特殊條件下的土壤水運動和淺層地下水水位的變化，以及土壤水分對降水和蒸發(fā)的響應(yīng)。模型以土水勢為研究變量，首先對研究區(qū)域進行網(wǎng)格剖分，利用有限差分的格點法將研究對象離散化，進而模擬平原區(qū)淺層土壤含水率以及非飽和區(qū)與飽和區(qū)的變換過程。在模型求解過程中，運用“邊循環(huán)”的概念和矩陣標識法，進一步提高了模型的精確度和計算效率。

1 模型的求解

采用有限差分法進行求解。求解過程中，運用“邊循環(huán)”的概念和矩陣標識法達到快速準確求解模型的目的。采用有限差分方法將滲流區(qū)域劃分網(wǎng)格，并采用格點法，對方程中的左端項采用中心差分，右端時間項采用向前差分，方程離散如下

式中，K為x，y，z方向上的導水率；H為總水頭；Sf為土壤空隙中含水百分比；us為給水度；C為比水容量；h為壓力水頭；W為源匯項的水量。

方程聯(lián)系著格點(i,j,k)在t時刻與t+Δt時刻的水頭，同時還聯(lián)系著與這個格點相鄰的6個格點的水頭。對每個未知水頭的格點都列出上述方程，當格點靠近邊界時用給定水頭或給定流量的邊界條件，這樣可以組成一個線性代數(shù)方程組。其中，方程左端的水頭全部采用t+Δt時刻的值即采用全隱式解法；將上述方程化簡為m×x=y方程的格式，并與其周圍相鄰網(wǎng)格的水頭關(guān)系方程式進行聯(lián)立，構(gòu)成如MX=Y形式的矩陣以方便求解。系數(shù)矩陣M有以下特點：每個網(wǎng)格最多會產(chǎn)生7個非零元素，對于寬廣計算區(qū)域，網(wǎng)格很多，方程階數(shù)很高，對這樣的高稀疏矩陣，計算的效率就顯得尤為重要。本文采用矩陣標識法求解方程，使網(wǎng)格編碼具有任意性，避免非零元素的存儲和計算，又不用考慮收斂問題，較好的提高了模型的精確度和計算效率。

2 模型驗證

2.1 垂向一維降雨入滲模擬

1 m×1 m×2.3 m的垂向土柱，垂向上剖分成230層，dx=dy=1 m、dz=0.01 m。土柱的初始壓力水頭為-0.8 m，上邊界給4 mm/h的降雨，下邊界、四周邊界沒有水量交換。模型計算結(jié)果見圖1。

圖1 不同降雨歷時土壤水分剖面

從圖1可知，降雨開始后，在積水點之前，降雨全部下滲到土柱中，濕潤峰面下移；當?shù)乇砗蔬_到相對飽和，濕潤峰面全部進入土體中，即t=3 h時，土體開始積水入滲，此后濕潤峰繼續(xù)下移；t=13 h時，入滲率開始出現(xiàn)拐點，入滲率明顯下降，此后當濕潤峰全部通過底部出溢面時，入滲率達到飽和點，土柱開始飽和入滲。

2.2 降雨模擬

2.2.1 單網(wǎng)格降雨入滲

2.5 m×2.5 m×1 m的土柱剖分成25×25×100的網(wǎng)格，即dx=dy=0.1 m、dz=0.01 m。選上邊界中任意網(wǎng)格，給予連續(xù)10 mm/h 的18 h降雨，下邊界、四周邊界都沒有水量交換，各層的土壤水力參數(shù)設(shè)為相同，設(shè)土柱的初始壓力水頭為-5 m。z方向的飽和導水率為0.003 44 cm/s，x、y方向的飽和導水率是z方向的1/4。40 h降雨計算結(jié)果見圖2。

圖2 土壤含水率等值線

從圖2可知，降雨落到網(wǎng)格后，向下和四周擴散。當x、y方向上的導水率值相等時，由中心網(wǎng)格流入周邊網(wǎng)格的流速、水量都是一致的。圖2a中，有降雨網(wǎng)格的土壤含水率由0.257上升到0.295，降雨影響到的外圍網(wǎng)格的土壤含水率上升到0.260；圖2b中，有降雨網(wǎng)格的土壤含水率由0.257上升到0.268，降雨影響到的外圍網(wǎng)格的土壤含水率由0.257上升到0.258。對比可知，降雨落到網(wǎng)格上后，10 cm處受降雨影響的面積大約為1 m2，20 cm處大約為0.49 m2，由上到下，受影響面積逐漸減小。對比模型初始和計算完成后的水量與降雨入滲及產(chǎn)流量，模型計算符合水量平衡要求。

2.2.2 多網(wǎng)格降雨入滲

運用上述土柱條件，給上層網(wǎng)格中相互對稱的4個角落網(wǎng)格連續(xù)10 mm/h的降雨。計算結(jié)果見圖3。從圖3可知，降雨分布在平面對稱的4個網(wǎng)格中。在網(wǎng)格土壤特性一致的前提下，土壤水的運動也保持一致，土壤含水率的值相同，且4個角上的等值線保持對稱，降雨所在網(wǎng)格下20 cm處的含水率為0.445，降雨影響到的最大區(qū)域為0.9 m，該處的土壤含水率為0.426。

圖3 土壤含水率等值線

2.3 蒸發(fā)模擬

土柱信息同上。在最上層的中心網(wǎng)格給1 mm/h的蒸發(fā)，其他網(wǎng)格沒有降雨和蒸發(fā)。計算結(jié)果見圖4。從圖4可知，當中間網(wǎng)格由于蒸發(fā)導致土壤含水率降低時，四周網(wǎng)格的水勢平衡被破壞，土壤水在土水勢的作用下，由土水勢的高處向低處運動。因此，中心網(wǎng)格(有蒸發(fā))的土壤含水率減小幅度最大，由0.258減小到0.253，由中心網(wǎng)格擴散，四周的網(wǎng)格土壤含水率減小幅度逐步減小。

圖4 土壤含水率等值線

2.4 淺層地下水水位模擬

平原區(qū)地下水埋深較淺，一次降雨對潛水位影響較大。在此，模擬降雨條件下的淺層地下水水位選取3 m×3 m×2 m的一個大土柱，土柱上部高程一致，將土柱剖分成30×30×200的網(wǎng)格，即dx=dy=0.1 m、dz=0.01 m，相當于分成了90個2 m的小土柱，在整個大土柱的左上角的0.5 m×0.5 m的區(qū)域內(nèi)給148 mm/h的降雨。地下埋深3 m處的不同歷時x方向上淺層地下水水位隨時間變化見圖5。從圖5可知，土柱含水率在時間上成上升趨勢，但在x方向上從左向右有減小的趨勢。

圖5 x方向上淺層地下水水位隨時間變化

40 h土柱地下水水位渲染見圖6。從圖6可知，降雨發(fā)生前，大土柱的地下水水位保持穩(wěn)定，且各網(wǎng)格的地下水水位持平。降雨發(fā)生在某一片區(qū)域后，經(jīng)土壤水垂直下滲和水平運移，降雨所在網(wǎng)格的小土柱的地下水水位率先上升，由降雨中心向外，各小土柱的地下水水位依次增大。而離有降雨網(wǎng)格最遠的對角網(wǎng)格的地下水水位保持不變，從地下水水位最高處到最低處，有規(guī)律的下降，下降曲線應(yīng)是光滑的。由于模型計算誤差和網(wǎng)格劃分不夠密，導致圖形存在鋸齒。

圖6 40 h土柱地下水水位渲染

不同地下水水位土壤水分剖面見圖7。從圖7可知，地下水水位上升過程中，土壤水分剖面曲線先是迅速由扁平變成尖瘦；在地下水水位由0.695 m上升到0.855 m的過程中，土壤水分剖面變動很小。

圖7 不同地下水水位土壤水分剖面

2.4 實際應(yīng)用

試驗區(qū)選址在常州市金壇朱林鎮(zhèn)紅旗圩村。該地區(qū)地形平坦、坡度小、河網(wǎng)交織、湖塘遍布。土壤多為第四系沖、洪積松散沉積和湖相沉積的粉粘土，顆粒細微，滲透性差。該地區(qū)地下水年平均埋深1 m左右，地下水水位對降雨輸入的響應(yīng)迅速。該區(qū)域能較好地反映太湖流域平原河網(wǎng)區(qū)的基本水文特征。

降雨資料由Campbell TE525MM系列翻斗式雨量筒測定，精度為0.1 mm；土壤含水率采用TDR傳感器觀測。應(yīng)用上述模型進行模擬，對比模擬計算含水率和實測含水率。10、20 cm處土壤含水率計算和實測對比見圖8。

從圖8可以看出，模型較好地模擬了平原區(qū)飽和-非飽和帶土壤含水率的變化情況，特別是在表層10 cm處，計算擬合數(shù)據(jù)與實際觀測數(shù)據(jù)的相關(guān)性R為0.64，個別場次的降雨模擬誤差低于10%，能夠反映土壤水的整體變化趨勢。相對于實測值，模擬結(jié)果整體較低，需對模型進行改進。

3 結(jié) 語

本文以土壤水動力學理論為基礎(chǔ)，輔以考慮土壤蒸發(fā)、降雨土壤入滲等水文循環(huán)

圖8 土壤含水率計算和實測對比

過程模型，構(gòu)建了三維飽和-非飽和土壤水運動模型，運用有限差分法推導出三維飽和-非飽和流方程的差分方程。運用“邊循環(huán)”和矩陣標識法，縮短模型計算時間，求解速度大大提高。實際應(yīng)用表明，該模型符合土壤水運動原理和日常經(jīng)驗，重現(xiàn)了特定土壤水運動的基本過程，為探究平原區(qū)水循環(huán)過程、水分轉(zhuǎn)化關(guān)系以及進行大尺度平原區(qū)水文預報提供理論基礎(chǔ)。

[1]DEMETRIOU C， PUNTHAKEY J F. Evaluating sustainable groundwater management options using the MIKE SHE integrated hydrogeological modelling package[J]. Environmental Modelling & Software， 1998， 14(2- 3)： 129- 140．

[2]劉新仁， 楊海艦. 土壤水動力學在平原水文模擬中的應(yīng)用[J]. 河海大學學報， 1989， 17(4)： 9- 15．

[3]張蔚榛. 地下水與土壤水動力學[M]. 北京： 中國水利水電出版社， 1996．

[4]陳喜， 程勤波， 張志才. 飽和-非飽和水流數(shù)值模擬[M]. 北京： 中國水利水電出版社， 2011．

[5]雷志棟， 楊詩秀， 謝森傳. 土壤水動力學[M]. 北京： 清華大學出版社， 1988．

[6]張雪剛， 毛媛媛， 董家瑞， 等. SWAT模型與MODFLOW模型的耦合計算及應(yīng)用[J]. 水資源保護， 2010， 26(3)： 49- 52．

[7]XU X， HUANG G， ZHAN H， et al. Integration of SWAP and MODFLOW- 2000 for modeling groundwater dynamics in shallow water table areas[J]. Journal of Hydrology， 2012， 412(1)： 170- 181．

[8]程勤波， 陳喜， 趙玲玲， 等. 飽和與非飽和帶土壤水動力耦合模擬及入滲試驗[J]. 河海大學學報， 2008， 37(3)： 284- 289．

(責任編輯 楊 健)

Study on Flow Model and Numerical Simulation of Unsaturated and Saturated Soil Water in Plain Area

CHEN Jingbo1, WANG Chuanhai1, DU Shipeng2, YANG Hai1, ZHANG Mengfei1, ZHU Liguo3

(1. State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, Jiangsu, China; 2. Water Resources and Hydropower Survey and Design Institute of Zhejiang Province, Hangzhou 310002, Zhejiang, China; 3. Wuxi Branch, Hydrology and Water Resources Survey Bureau of Jiangsu Province, Wuxi 214000, Jiangsu, China)

Based on the theory of soil water movement, a 3D model is constructed to simulate the movement of unsaturated and saturated water in plain area. The finite difference method is used to simulate soil moisture content and transformation process in each location of study area. Some case examples are used to verify the stability and accuracy of model. The calculation results show that the process of soil water moisture and groundwater level change is reasonable in theory, the simulation results of actual soil water situation is relatively reasonable, which can reflect the hydrological cycle in plain area. The results provide theoretical basis for exploring water cycle process, water transformation relationship and large-scale hydrological forecasting in plain area．

soil water; water cycle; unsaturated and saturated flow; numerical simulation

2016- 04- 06

江蘇省水利科技項目(2015007)；國家水體污染控制與治理科技重大專項(2014ZX07101- 011)

陳景波(1990—)，男，河北邯鄲人，碩士研究生，研究方向為流域水文模擬及預報．

S152.7(253)

A

0559- 9342(2016)09- 0013- 04