胡 丹，伍靖偉，黃介生

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

0 引 言

土壤水和地下水運動是水文循環(huán)的重要環(huán)節(jié)，能夠準(zhǔn)確模擬土壤水、地下水運動對水資源高效利用具有重要意義[1,2]。地下水與土壤水之間有著直接的水力聯(lián)系和水通量交換，為了更好地模擬區(qū)域尺度土壤水、地下水的運動，近年來發(fā)展了一些飽和-非飽和耦合模型，這些模型的主要區(qū)別在于對土壤水運動的模擬，其主要方式有3種：第一種是水均衡模型，如MODFLOW中的REC和ET模塊[3]，以土壤水質(zhì)量平衡為原理，獨立考慮蒸發(fā)、根系吸水等水文過程，由于水均衡模型不涉及水頭，無法與地下水的計算直接聯(lián)系，過分簡化了土壤水運動，故會帶來較大誤差，但求解相對簡單；第二種是基于動力波方程，如MODFLOW-UZF1[4]，模擬土壤水運動時只考慮重力作用，而忽略毛管力作用；第三種是基于理查茲方程，如三維非飽和-三維飽和的耦合模型MODFLOW-VSF[5]、一維非飽和-三維飽和的耦合模型MODFLOW-HYDRUS[6]，理查茲方程是由達(dá)西定律和質(zhì)量守恒定律推導(dǎo)而得，對土壤水運動的描述更加準(zhǔn)確，但由于具有高度非線性，計算成本也更大。飽和-非飽和耦合模型能否做到模擬精度和計算成本的平衡直接關(guān)系到模型的運用推廣。Twarakavi等[7]對MODFLOW-UZF1、MODFLOW-VSF和MODFLOW-HYDRUS三種耦合模型進行了比較，發(fā)現(xiàn)進行區(qū)域模擬時，MODFLOW-HYDRUS和MODFLOW-VSF的模擬精度比MODFLOW-UZF1更高，在計算要求方面，MODFLOW-HYDRUS和MODFLOW-UZF1的計算量比MODFLOW-VSF更小，尤其對于淺地下水、蒸發(fā)強烈的地區(qū)，MODFLOW-HYDRUS進行土壤水-地下水區(qū)域模擬時，相比MODFLOW-UZF1和MODFLOW-VSF，在模擬精度和計算量方面能做到更好的平衡。

飽和-非飽和耦合模型進行區(qū)域水分運動模擬時，往往需要輸入較多的參數(shù)，且水文地質(zhì)參數(shù)和土壤水力參數(shù)空間變異性較大，導(dǎo)致模擬結(jié)果存在較大不確定性。模型參數(shù)的獲得目前主要有實測法和參數(shù)率定法，其中實驗法測得的數(shù)據(jù)由于尺度效應(yīng)的存在，很難直接運用到模型中。參數(shù)率定法包括試錯法、參數(shù)自動優(yōu)化方法以及試錯法與參數(shù)自動優(yōu)化法聯(lián)合優(yōu)化。傳統(tǒng)的試錯法依賴使用者對模型的了解及經(jīng)驗，主觀性太強，花費時間較長，且難以達(dá)到參數(shù)最優(yōu)化。參數(shù)自動優(yōu)化方法能夠較快得到優(yōu)化解，但是對初值的依賴較大，不同的初值可能導(dǎo)致差別較大的率定結(jié)果[8]。試錯法與參數(shù)自動優(yōu)化法聯(lián)合優(yōu)化既能發(fā)揮模型使用者的知識和經(jīng)驗，又能利用先進的優(yōu)化技術(shù)[9]。

PEST[10]是獨立于計算模型的參數(shù)優(yōu)化方法，已經(jīng)分別在地下水、土壤水運動模型中得到廣泛運用[11-13]，但是在耦合模型中運用效果和效率如何尚未有研究。

本文以內(nèi)蒙古河套灌區(qū)永聯(lián)試驗區(qū)地下水觀測井?dāng)?shù)據(jù)和土壤含水率實測數(shù)據(jù)為基礎(chǔ)，采用試錯法和PEST相結(jié)合的方法，對MODFLOW-HYDRUS模型的水文地質(zhì)參數(shù)和土壤水力參數(shù)進行優(yōu)化，以提升區(qū)域尺度上飽和-非飽和帶中水分運動的模擬效果。

1 材料與方法

1.1 研究區(qū)域及數(shù)據(jù)

研究區(qū)域為內(nèi)蒙古河套灌區(qū)五原縣永聯(lián)村境內(nèi)的永聯(lián)試驗區(qū)(見圖1)，東經(jīng)108°00′～108°01′，北緯41°04′～41°05′，南北長約13 km，東西寬約5 km，總面積約2 875 hm2，其中灌溉面積約1 453 hm2，土地利用類型主要是耕地、裸地(包括荒地、渠道、排水溝等)和村莊。試驗區(qū)地勢自南向北稍有傾斜，中間較高，東西兩側(cè)略低，整體較為平坦，地面高程1 028.6～1 025.6 m。試驗區(qū)邊界條件較為清楚，北邊是六排干排水邊界，南邊是皂火渠，東西兩邊分別是永什分干溝和乃永分干溝(見圖1)。

試驗區(qū)內(nèi)布置有11口地下水觀測井(見圖1)，每5 d觀測1次地下水位，每個地下水觀測井附近布置有剖面土壤含水率監(jiān)測點，每10 d進行1次取樣測定，取樣深度為0～100 cm，分為0～10，10～30，30～50，50～70，70～100 cm 5個層次。地下水位觀測采用測鐘法，土壤含水率采用烘干法測定。

圖1 研究區(qū)示意圖Fig.1 Sketch map of study area

1.2 飽和-非飽和耦合模型

非飽和帶水流運動考慮為一維垂向運動，用一維Richards方程來描述：

(1)

式中：θ為土壤的體積含水量，cm3/cm3；t為時間，d；h為壓力水頭，m；z為空間坐標(biāo)，向上為正，m；S為匯項，m3/(m3·d)；K(h)為非飽和水力傳導(dǎo)度，m/d。

飽和帶水流運動為三維運動，其控制偏微分方程為：

(2)

式中：h為水頭，m；x，y，z為空間坐標(biāo)，m；Kxx，Kyy，Kzz為滲透系數(shù)在x，y，z方向上的分量，這里假定滲透系數(shù)主軸方向與坐標(biāo)軸方向一致，m/d；W為匯項，m3/(m3·d)；t為時間，d；SS為給水度，無量綱。

非飽和帶與飽和帶的耦合過程參考Seo等[5]的介紹。

本文采用2007年5月1日至11月30日(共214 d)期間實測的地下水位和土壤含水率數(shù)據(jù)對模型參數(shù)進行率定。

(1)網(wǎng)格劃分與土壤剖面設(shè)置。將研究區(qū)域劃分為25×50×2網(wǎng)格，其中有效網(wǎng)格為742個，水平方向上網(wǎng)格大小為160 m×264 m，垂直方向上，第一、二層底部高程分別為1 020.29 m和973.75 m。根據(jù)地面高程、地下水位、土地利用類型等情況，采用k均值聚類分析方法，在研究區(qū)內(nèi)設(shè)置25個土壤剖面(見圖2)，其中1～7號土壤剖面為耕地，8～14號為排水溝，15～16號為人民支渠，17號為皂火渠，18～21號為荒地，22～25號為村莊。土壤剖面的高度為地面到1 020.29 m高程處。

圖2 網(wǎng)格與土壤剖面Fig.2 Grids and soil profiles

(2)初始輸入數(shù)據(jù)。以2007年5月1日實測的地下水位數(shù)據(jù)插值得到整個區(qū)域的地下水位作為初始地下水位[見圖3(b)]；土壤剖面的初始壓力水頭分布按照Seo等[5]的方法計算得到。

圖3 地面高程和初始地下水位Fig. 3 Ground surface elevation and initial groundwater table elevation

(3)邊界條件。上邊界條件為變化的大氣邊界和灌溉補給。從五原縣氣象站獲得該研究區(qū)域氣象資料，包括降雨量、氣溫、濕度、風(fēng)速、日照時數(shù)等，采用FAO-56[14]介紹的方法計算每日的潛在蒸發(fā)量和蒸騰量，將人民支渠口部的凈引水量作為區(qū)域內(nèi)的灌溉水量，而灌溉滲漏量概化為人民支渠上的線性補給。耕地(1～7號土壤剖面)的上邊界條件包括降雨、蒸發(fā)、蒸騰和灌溉[見圖4(a)]，荒地(18～21號土壤剖面)、排水溝(8～14號土壤剖面)和皂火渠(17號土壤剖面)的上邊界條件包括降雨和蒸發(fā)[見圖4(b)]，人民支渠(15～16號土壤剖面)的上邊界條件包括降雨、蒸發(fā)和灌溉滲漏量[見圖4(c)]，村莊(22～25號土壤剖面)的上邊界條件只有降雨[見圖4(d)]。第二含水層下是約1 m的黏土層，故將下邊界條件視為零通量邊界。將六排干、永什分干和乃永分干三條排水溝設(shè)置為排水邊界，皂火渠灌溉渠道設(shè)置為河流邊界。

圖4 不同土壤剖面的上邊界條件Fig.4 Upper boundary conditions of different soil profiles

(4)模型參數(shù)。研究區(qū)域兩個含水層分別是沙壤土和砂土，水文地質(zhì)參數(shù)的初始值根據(jù)五原縣永聯(lián)村鉆孔水文地質(zhì)報告取值(見表1)，van Genuchten[15]模型中土壤水力參數(shù)的初始值根據(jù)Carsel and Parrish[16]取值(見表1)。然后采用試錯法對模型各參數(shù)進行調(diào)整(見表1)，最后進行參數(shù)自動優(yōu)化。

表1 土壤水力參數(shù)和水文地質(zhì)參數(shù)取值Tab.1 Values of soil hydraulic parameters and hydrogeological parameters

1.3 PEST參數(shù)優(yōu)化算法

PEST算法采用高斯牛頓法與梯度下降法結(jié)合快速收斂性，又具備梯度下降法的全局搜索性，它能夠使目標(biāo)函數(shù)快速有效地收斂，進行參數(shù)優(yōu)化時模型運行的次數(shù)比其他算法更少[17]。

PEST算法的核心是對目標(biāo)函數(shù)求解最小值，目標(biāo)函數(shù)是實驗觀測值與計算值的帶權(quán)重差異函數(shù)，根據(jù)參數(shù)優(yōu)化的需求和實驗觀測值的種類，可將觀測值分組并賦予相應(yīng)的權(quán)重，以達(dá)到模擬值與實測值的最佳匹配。本研究包含地下水位、土壤含水率2組目標(biāo)函數(shù)，其公式如下：

(3)

式中：Φ為目標(biāo)函數(shù)；Oi和Si分別表示第i個實驗觀測值和第i個模擬值；p表示實驗觀測值的個數(shù)；ω表示權(quán)重系數(shù)；下標(biāo)wt和sm分別對應(yīng)地下水位和土壤含水率。

1.4 評價方法

采用決定系數(shù)(R2)、納什系數(shù)(NSE)、均方根誤差(RMSE)、標(biāo)準(zhǔn)化均方根誤差(NRMSE)4種指標(biāo)來評價模型的模擬效果：

(7)

決定系數(shù)(R2)用來衡量模擬值和實測值的相關(guān)密切程度，R2越接近1表示相關(guān)性越好；納什系數(shù)(NSE)可用來評價模型的模擬質(zhì)量，NSE越接近1說明模擬質(zhì)量越好；均方根誤差(RMSE)用來衡量模擬值與實測值之間的偏差，RMSE越接近0，表示模擬值與觀測值的差異越小，模擬效果越好；標(biāo)準(zhǔn)化均方根誤差(NRMSE)也是用來衡量模擬誤差的大小，NRMSE值越小，說明模擬效果越好。

2 結(jié)果與分析

參數(shù)自動優(yōu)化共選取19個參數(shù)作為優(yōu)化對象，包括15個土壤水力參數(shù)和4個水文地質(zhì)參數(shù)，參數(shù)優(yōu)化初值為試錯法調(diào)整后的結(jié)果，各個參數(shù)的取值范圍和優(yōu)化結(jié)果見表1。

將試錯法和PEST參數(shù)優(yōu)化后的地下水位模擬值和觀測井實測值進行比較(見圖5)，同時將模擬開始后第46、86、128、168天土壤剖面含水率的模擬值和實測值進行比較(見圖6)。將所有的地下水位觀測值與對應(yīng)時刻的模擬值繪制在一張圖中(見圖7)，將所有的土壤含水率實測值與對應(yīng)時刻的模擬值繪制在一張圖中(見圖7)，并計算評價模型效果的4個指標(biāo)(見表2)。

采用2008年5月1日至11月31日期間的地下水位和土壤含水率實測值對模型進行驗證，將地下水位觀測值與對應(yīng)時刻的模擬值繪制在一張圖中(見圖8)，將土壤含水率實測值與對應(yīng)時刻的模擬值繪制在一張圖中(見圖8)，并計算評價模型效果的4個指標(biāo)(見表2)。

圖5 地下水位模擬值與實測值對比Fig.5 Computed vs. observed groundwater table

圖6 不同時刻土壤剖面含水率模擬值與實測值對比Fig.6 Computed vs. observed soil moisture of different time

模擬期觀測變量試錯法R2NSERMSENRMSEPESTR2NSERMSENRMSE率定期地下水位0．7650．7310．5210．1040．8080．7860．4640．093土壤含水率0．4160．4080．0640．1370．5290．5280．0570．122驗證期地下水位0．8040．6030．5800．1290．8240．6410．5520．123土壤含水率0．3250．0320．0660．1740．4380．3950．0530．140

圖7 率定期地下水位與土壤含水率散點圖Fig. 7 Scatter plots of groundwater table and soil moisture for calibration period

圖8 驗證期地下水位與土壤含水率散點圖Fig. 8 Scatter plots of groundwater table and soil moisture for validation period

從圖5可以看出，地下水位模擬值與觀測值的變化趨勢基本一致；從圖6可以看出，模擬開始后不同時刻土壤剖面含水率模擬值與實測值較為接近；從圖7、圖8可以發(fā)現(xiàn)，地下水位和土壤含水率散點都較好地集中在1:1直線附近。

表2的模型效果評價指標(biāo)顯示，PEST算法進行參數(shù)優(yōu)化，提高了模型的模擬效果。

在率定期，地下水位的決定系數(shù)由0.765提高到0.808，提高了6%；納什系數(shù)由0.731提高到0.786，提高了8%；均方根誤差由0.521減小到0.464，減小了11%；標(biāo)準(zhǔn)化均方根誤差由0.104減小到0.093，減小了11%。

在率定期，土壤含水率的決定系數(shù)由0.416提高到0.529，提高了27%；納什系數(shù)由0.408提高到0.528，提高了29%；均方根誤差由0.064減小到0.057，減小了11%；標(biāo)準(zhǔn)化均方根誤差由0.137減小到0.122，減小了11%。

在驗證期，地下水位的決定系數(shù)由0.804提高到0.824，提高了2%；納什系數(shù)由0.603提高到0.641，提高了6%；均方根誤差由0.580減小到0.552，減小了5%；標(biāo)準(zhǔn)化均方根誤差由0.129減小到0.123，減小了5%。

在驗證期，土壤含水率的決定系數(shù)由0.325提高到0.438，提高了35%；納什系數(shù)由0.032提高到0.395，提高了11.3倍；均方根誤差由0.066減小到0.053，減小了20%；標(biāo)準(zhǔn)化均方根誤差由0.174減小到0.140，減小了20%。

同時從表2可以看出，耦合模型對地下水位的模擬效果要優(yōu)于土壤含水率。從試錯法的結(jié)果看，地下水位在率定期和驗證期決定系數(shù)分別為0.765和0.804，納什系數(shù)分別為0.731和0.603；而土壤含水率在率定期和驗證期決定系數(shù)分別為0.416和0.325，納什系數(shù)分別為0.408和0.032。通過PEST參數(shù)優(yōu)化后，地下水位在率定期和驗證期決定系數(shù)分別為0.808和0.824，納什系數(shù)分別為0.786和0.641；而土壤含水率在率定期和驗證期決定系數(shù)分別為0.529和0.438，納什系數(shù)分別為0.528和0.395。土壤含水率的各項評價指標(biāo)都低于地下水位的相應(yīng)指標(biāo)。對土壤水運動模擬效果較差，可能的原因是本次模擬將整個土壤剖面視為均質(zhì)土壤，與實際土壤剖面為非均質(zhì)土壤的情況不符。

3 結(jié) 語

(1)本文采用試錯法和PEST優(yōu)化算法相結(jié)合的方法對MODFLOW-HYDRUS耦合模型的水文地質(zhì)參數(shù)和土壤水力參數(shù)進行優(yōu)化，提高了模擬精度。與試錯法的結(jié)果相比，PEST參數(shù)優(yōu)化后決定系數(shù)提高了2%～35%，納什系數(shù)提高了6～11.3倍，均方根誤差減小了5%～20%，標(biāo)準(zhǔn)化均方根誤差減小了5%～20%，其中驗證期土壤含水率的各項指標(biāo)提高或減小的幅度均為最大，決定系數(shù)、納什系數(shù)、均方根誤差和標(biāo)準(zhǔn)化均方根誤差變化幅度分別為35%、11.3倍、20%和20%。

(2)MODFLOW-HYDRUS耦合模型模擬區(qū)域尺度上地下水、土壤水運動時，在地下水運動模擬上的表現(xiàn)要優(yōu)于土壤水運動的模擬。從試錯法與PEST參數(shù)優(yōu)化后結(jié)果來看，地下水位的決定系數(shù)在0.765～0.824之間，納什系數(shù)在0.603～0.786之間；而土壤含水率的決定系數(shù)在0.325～0.529之間，納什系數(shù)在0.032～0.528之間，明顯低于地下水位相對應(yīng)的指標(biāo)。導(dǎo)致土壤水運動模擬效果不佳的原因可能是將整個土壤剖面視為均質(zhì)土壤，與實際土壤剖面非均質(zhì)土壤的情況不符，在今后運用模型時可以考慮按照實際的非均質(zhì)土壤剖面來設(shè)置。

□

[1] 楊欣坤, 王 宇, 趙蘭坡, 等. 土壤水動力學(xué)參數(shù)及其影響因素研究進展[J]. 中國農(nóng)學(xué)通報, 2014,30(3):38-43.

[2] 薛禹群. 中國地下水?dāng)?shù)值模擬的現(xiàn)狀與展望[J]. 高校地質(zhì)學(xué)報, 2010,16(1):1-6.

[3] Harbaugh A W, E R Banta, M C Hill, et al. MODFLOW-2000, the U.S. Geological Survey modular ground-water model user guide to modularization concepts and the ground-water fl ow process.USGS, Denver, CO,2000.

[4] Niswonger R G, D E Prudic, R S Regan. Documentation of the Unsaturated-Zone Flow (UZF1) package for modeling unsaturated flow between the land surface and the water table with MODFLOW-2005. Techniques and Methods 6-A19. Available at http:∥pubs.usgs.gov/tm/2006/tm6a19/ (verified 14 Apr. 2008). USGS, Denver, CO,2006.

[5] Thoms R B, R L Johnson, R W Healy. User’s guide to the Variably Saturated Flow (VSF) process for MODFLOW. Techniques and Methods 6-A18. Available at http://pubs.usgs.gov/tm/2006/ tm6a18/ (verified 15 Apr. 2008). USGS, Reston, VA, 2006.

[8] 舒曉娟, 陳洋波, 黃鋒華, 等. PEST在WetSpa分布式水文模型參數(shù)率定中的應(yīng)用[J]. 水文, 2009,29(5):45-49.

[9] 章四龍, 劉九夫. 通用模型參數(shù)率定技術(shù)研究[J]. 水文, 2005,25(1):9-12.

[10] Doherty J. PEST Model-Independent Parameter Estimation User Manual[M]. 5th Edition. Brisbane, Australia: Watermark Numerical Computing, 2004.

[11] Zyvoloski G, E Kwicklis, A A. Eddebbarh, et al. The site-scale saturated zone flow model for Yucca Mountain: calibration of different conceptual models and their impact on flow paths[J]. Journal of Contaminant Hydrology, 2003,62-63(2):731-750.

[12] 董艷輝, 李國敏, 郭永海, 等. 應(yīng)用并行PEST 算法優(yōu)化地下水模型參數(shù)[J]. 工程地質(zhì)學(xué)報, 2010,18(1):140-145.

[13] Fang Q X, Green T R, Ma L, et al. Optimizing soil hydraulic parameters in RZWQM2 under fallow conditions[J]. Soil Science Society of America Journal, 2010,74(6):1 897-1 913.

[14] Allen R G, Pereira L S, Raes D, et al. Crop Evapotranspiration: Guidelines for Computing Crop Water Requirements[D]. Rome, Italy: United Nations Food and Agriculture Organization, 1998.

[15] van Genuchten, M.Th. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980,44(44):892-898.

[16] Carsel R F, Parrish R S. Developing joint probability distributions of soil water retention characteristics[J]. Water Resources Research, 1988,24(5):755-769.

[17] Bahremand A, De Smedt F. Distributed hydrological modeling and sensitivity analysis in Torysa Watershed, Slovakia[J]. Water Resources Management, 2008,22(3):393-408.