陳景波，王船海，杜世鵬，楊 海，張夢(mèng)菲，朱立國(guó)

(1.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇南京210098；2.浙江省水利水電勘測(cè)設(shè)計(jì)院，浙江杭州310002；3.江蘇省水文水資源勘測(cè)局無(wú)錫分局，江蘇無(wú)錫214000)

?

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

陳景波1，王船海1，杜世鵬2，楊 海1，張夢(mèng)菲1，朱立國(guó)3

(1.河海大學(xué)水文水資源與水利工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，江蘇南京210098；2.浙江省水利水電勘測(cè)設(shè)計(jì)院，浙江杭州310002；3.江蘇省水文水資源勘測(cè)局無(wú)錫分局，江蘇無(wú)錫214000)

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

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

0 引 言

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

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

1 模型的求解

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

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

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

2 模型驗(yàn)證

2.1 垂向一維降雨入滲模擬

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

圖1 不同降雨歷時(shí)土壤水分剖面

從圖1可知，降雨開始后，在積水點(diǎn)之前，降雨全部下滲到土柱中，濕潤(rùn)峰面下移；當(dāng)?shù)乇砗蔬_(dá)到相對(duì)飽和，濕潤(rùn)峰面全部進(jìn)入土體中，即t=3 h時(shí)，土體開始積水入滲，此后濕潤(rùn)峰繼續(xù)下移；t=13 h時(shí)，入滲率開始出現(xiàn)拐點(diǎn)，入滲率明顯下降，此后當(dāng)濕潤(rùn)峰全部通過(guò)底部出溢面時(shí)，入滲率達(dá)到飽和點(diǎ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è)土柱的初始?jí)毫λ^為-5 m。z方向的飽和導(dǎo)水率為0.003 44 cm/s，x、y方向的飽和導(dǎo)水率是z方向的1/4。40 h降雨計(jì)算結(jié)果見圖2。

圖2 土壤含水率等值線

從圖2可知，降雨落到網(wǎng)格后，向下和四周擴(kuò)散。當(dāng)x、y方向上的導(dǎo)水率值相等時(shí)，由中心網(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。對(duì)比可知，降雨落到網(wǎng)格上后，10 cm處受降雨影響的面積大約為1 m2，20 cm處大約為0.49 m2，由上到下，受影響面積逐漸減小。對(duì)比模型初始和計(jì)算完成后的水量與降雨入滲及產(chǎn)流量，模型計(jì)算符合水量平衡要求。

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

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

圖3 土壤含水率等值線

2.3 蒸發(fā)模擬

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

圖4 土壤含水率等值線

2.4 淺層地下水水位模擬

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

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

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

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

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

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

2.4 實(shí)際應(yīng)用

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

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

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

3 結(jié) 語(yǔ)

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

圖8 土壤含水率計(jì)算和實(shí)測(cè)對(duì)比

過(guò)程模型，構(gòu)建了三維飽和-非飽和土壤水運(yùn)動(dòng)模型，運(yùn)用有限差分法推導(dǎo)出三維飽和-非飽和流方程的差分方程。運(yùn)用“邊循環(huán)”和矩陣標(biāo)識(shí)法，縮短模型計(jì)算時(shí)間，求解速度大大提高。實(shí)際應(yīng)用表明，該模型符合土壤水運(yùn)動(dòng)原理和日常經(jīng)驗(yàn)，重現(xiàn)了特定土壤水運(yùn)動(dòng)的基本過(guò)程，為探究平原區(qū)水循環(huán)過(guò)程、水分轉(zhuǎn)化關(guān)系以及進(jìn)行大尺度平原區(qū)水文預(yù)報(bào)提供理論基礎(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]劉新仁， 楊海艦. 土壤水動(dòng)力學(xué)在平原水文模擬中的應(yīng)用[J]. 河海大學(xué)學(xué)報(bào)， 1989， 17(4)： 9- 15．

[3]張蔚榛. 地下水與土壤水動(dòng)力學(xué)[M]. 北京： 中國(guó)水利水電出版社， 1996．

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

[5]雷志棟， 楊詩(shī)秀， 謝森傳. 土壤水動(dòng)力學(xué)[M]. 北京： 清華大學(xué)出版社， 1988．

[6]張雪剛， 毛媛媛， 董家瑞， 等. SWAT模型與MODFLOW模型的耦合計(jì)算及應(yīng)用[J]. 水資源保護(hù)， 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]程勤波， 陳喜， 趙玲玲， 等. 飽和與非飽和帶土壤水動(dòng)力耦合模擬及入滲試驗(yàn)[J]. 河海大學(xué)學(xué)報(bào)， 2008， 37(3)： 284- 289．

(責(zé)任編輯 楊 健)

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

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

陳景波(1990—)，男，河北邯鄲人，碩士研究生，研究方向?yàn)榱饔蛩哪M及預(yù)報(bào)．

S152.7(253)

A

0559- 9342(2016)09- 0013- 04