霍思遠，靳孟貴

(1.中國地質(zhì)大學(xué)(武漢)生物地質(zhì)與環(huán)境地質(zhì)國家重點實驗室，湖北武漢 430074;2.中國地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430074)

降水或灌溉水通過包氣帶入滲到飽水帶的過程稱為降水或灌溉水的入滲補給。下滲補給地下水的部分灌溉水稱為灌溉回歸水［1］。降水灌溉入滲補給是華北平原淺層地下水的主要補給方式［2］。天然降水受氣候變化控制，具有一定的時空變異性。華北平原近50年氣候經(jīng)歷“冷濕-暖干”的變化過程，20世紀(jì)80年代以來由多雨期轉(zhuǎn)為少雨期［3］，降水量年際變化差異增大;同時，降水年內(nèi)分布極不均勻，全年降水量的75% ～85%分布在6—9月份;此外，次降水的變化呈現(xiàn)更為顯著的隨機性。降水在不同時間尺度上的變化特征造成了入滲補給規(guī)律的復(fù)雜性，而農(nóng)業(yè)灌溉又人為增加了入滲補給問題的不確定性，尤其是在干旱-半干旱農(nóng)業(yè)灌溉區(qū)，灌溉回歸水相對于降水入滲已經(jīng)成為主要的地下水補給來源，而傳統(tǒng)的補給量評價方法如地下水位波動法則無法適用灌溉區(qū)的地下水補給量評價［4］，因此灌溉回歸水的評價逐漸受到國內(nèi)外學(xué)者的關(guān)注。氚示蹤［5］、人工溴示蹤［6］、數(shù)值模擬［7］、地中滲透儀［8］等方法被用于灌溉回歸水的識別及評價，而在此基礎(chǔ)上，灌溉水的年內(nèi)分配對于入滲補給的影響需要進一步研究。

長期過量開采地下水已導(dǎo)致華北平原地下水位持續(xù)下降，并引發(fā)地面沉降、海水入侵等一系列地質(zhì)環(huán)境問題，詳細分析水量輸入項(降水、灌溉)的變化規(guī)律，對于合理評價入滲補給量和有效管理利用地下水資源具有重要意義。

本文以衡水地區(qū)為例，基于實測資料，利用HYDRUS［9］軟件建立一維變飽和水分運移模型，模擬分析了研究區(qū)近60年在天然降水條件下的垂向入滲補給規(guī)律以及次降水特征對于入滲補給的影響，同時評價了年周期內(nèi)各次灌溉的入滲補給量。

1 研究區(qū)概況

河北省衡水市位于華北平原中部的沖積湖積平原，屬溫帶大陸性季風(fēng)氣候，年均降水量約500 mm，降水多集中在7、8月份;年均水面蒸發(fā)量 976 mm(1985—2010年平均)［10］。研究區(qū)主要為冬小麥-夏玉米輪作，生育期(尤其是冬小麥)降水不能滿足作物生長需要，需開采地下水灌溉。本文選取河北省農(nóng)林科學(xué)院旱作農(nóng)業(yè)研究所的冬小麥-夏玉米實驗田為研究對象，實驗田位于河北省衡水市深州市護遲鎮(zhèn)(圖1)，不同試驗田塊進行不同的灌溉處理，其中2水田塊與當(dāng)?shù)貙嶋H灌溉條件較為接近。

圖1 實驗田示意圖及灌溉制度Fig.1 The experimental field and irrigation schedule

2 數(shù)值模擬

2.1 概念模型

根據(jù)實驗場地概況構(gòu)建一維變飽和水分運動模型，模擬天然降水及不同灌溉處理條件下水分垂向入滲補給淺層含水層的過程。如圖2所示，層狀非均質(zhì)包氣帶主要巖性為砂壤土，含有粉砂質(zhì)黏壤土夾層。上邊界為大氣邊界，考慮降水、灌溉、土壤蒸發(fā)、地表徑流和根系吸水作用。下邊界為自由排水邊界。

圖2 包氣帶巖性Fig.2 Lithology of the vadose zone

2.2 數(shù)學(xué)模型

式中:h——土壤負壓/cm;

c(h)——容水度/cm-1;

K(h)——土壤非飽和水力傳導(dǎo)度/(cm·d-1);

h0——初始土壤負壓值/cm;

P——降水量/(cm·d-1);

I——灌溉量/(cm·d-1);

E——蒸發(fā)量/(cm·d-1);

Roff——地表徑流/(cm·d-1);

t——時間/d;

B——模型下邊界;

S——根系吸水項/d-1，用式(2)表示:

式中:α(h，z)——水分脅迫函數(shù);

b(z)——根系分布函數(shù);

Tp——潛在蒸騰量/(cm·d-1)。

模型的上邊界通量計算:

式中:q——上邊界實際入滲量/(cm·d-1);

qp——上邊界潛在入滲量/(cm·d-1)，由表層巖性及含水率決定;

E——蒸發(fā)量/(cm·d-1);

EP——土壤潛在蒸發(fā)量/(cm·d-1);

h(0，t)——地表土壤負壓/cm;

hA——地表土壤允許的最小負壓值(hA=-10 000 cm)。

當(dāng)計算的地表土壤負壓值小于hA時，模型上邊界轉(zhuǎn)變?yōu)槎ㄋ^邊界計算(取h=hA);當(dāng)?shù)乇硗寥镭搲簼M足式(3)時，上邊界為變流量邊界，實際入滲量受潛在入滲量限制，當(dāng)降水強度大于潛在入滲能力時，地表徑流產(chǎn)生。

土壤水分特征曲線及非飽和水力傳導(dǎo)度采用van Genuchten 方程［11］表示:

式中:θr——土壤殘余含水率/%;

θs——土壤飽和含水率/%;

α、n、m——土壤水分特征曲線形狀參數(shù);

Ks——飽和水力傳導(dǎo)度/(cm·d-1);

Se——有效飽和度。

2.3 數(shù)值模型

2.3.1 離散化

采用等間隔剖分方式，以1 cm為間隔將包氣帶剖分為750個網(wǎng)格。天然降水模型模擬期從1957年1月1日—2013年12月31日，共計20 819d;從中選取典型年周期模擬分析人類灌溉活動的影響，典型年從2011年9月1日—2012年9月30日，共計396d。模型均采用變時間間隔剖分方式，初始時間間隔為1d，最小時間間隔為0.001d，最大時間間隔為2d。

2.3.2 模型參數(shù)獲取及識別

運用激光粒度儀將包氣帶不同深度所取環(huán)刀樣進行顆分定名，同時運用烘干法測定其干密度(表1)?；跍y試結(jié)果，運用Rosetta模型預(yù)測得到不同深度土樣的水力參數(shù)(表2)。

表1 顆粒分析及干密度Table 1 Particle-size analysis and bulk density

表2 Rosetta模型預(yù)測結(jié)果Table 2 Predicted soil hydraulic parameters by the Rosetta model

通過鋁盒烘干稱重法對試驗田塊的土壤重量含水率變化進行監(jiān)測，利用式(7)將其換算為土壤體積含水率。

式中:θv——體積含水率;

θω——重量含水率;

γ——干密度/(g·cm-3)。

根據(jù)模擬和實測土壤體積含水率對比(圖3)，對包氣帶表層砂壤土層及粉砂質(zhì)黏壤土層參數(shù)進行調(diào)整，得到識別后的土壤水力參數(shù)(表3)，用于模型計算。

圖3 模擬體積含水率與實測值對比Fig.3 Comparison between measured and simulated water contents

表3 識別后參數(shù)Table 3 Calibrated soil hydraulic parameters

2.3.3 邊界條件及初始條件

模型上邊界為大氣邊界，受降水、灌溉、蒸發(fā)、地表徑流及作物蒸騰作用影響，考慮作物根系吸水作用，最大根系深度為50 cm，根系分布密度從地表到最大根系深度線性遞減。根據(jù)實測氣象資料，采用Penman-Monteith公式計算參考作物騰發(fā)量。研究區(qū)采取小麥-玉米一年兩熟的種植制度，結(jié)合前人研究結(jié)果，得到作物系數(shù)［12］及葉面積指數(shù) LAI年內(nèi)變化規(guī)律［13］，利用式(9)計算作物潛在騰發(fā)量，利用式(10)計算土壤潛在蒸發(fā)量(Ep)及作物潛在蒸騰量(Tp)。

式中:ET0——參考作物騰發(fā)量/(mm·d-1);

Rn——凈太陽輻射/(J·m-2·d-1);

G——土壤熱通量/(J·m-2·d-1);

ea——實際水汽壓/kPa;

es——飽和水汽壓/kPa;

γ——濕度計常數(shù)/(kPa·℃-1);

Δ——飽和水汽壓曲線斜率/(kPa·℃-1);

T——溫度/℃;

u2——高度2 m處的平均風(fēng)速/(m·s-1)。

式中:ETp——作物潛在騰發(fā)量/(mm·d-1);

Kc——作物系數(shù)。

式中:LAI——作物葉面積指數(shù);

Kgr——消光系數(shù)，取值 0.39［24］。

圖4 衡水地區(qū)年降水量及年潛在騰發(fā)量(1957—2013)Fig.4 Annual precipitation and potential evapotranspiration from 1957 to 2013 in the Hengshui area

下邊界為自由排水邊界。根據(jù)野外實測剖面含水率分布，通過循環(huán)迭代計算得到的剖面穩(wěn)定含水率分布作為模型的初始條件。

3 結(jié)果及討論

3.1 降水變化對入滲補給的影響

降水是研究區(qū)淺層地下水的天然補給來源，其大小及變化是決定入滲補給量的關(guān)鍵。包氣帶是調(diào)節(jié)降水入滲補給地下水過程的樞紐，一次降水過后，入滲水量首先要補充包氣帶前期水分虧缺并滿足作物根系吸水需求，剩余水量受重力勢和基質(zhì)勢的共同作用向下運動，入滲水流穿過根系區(qū)域并達到極限蒸發(fā)深度以下時，不再被蒸騰蒸發(fā)所消耗，此時的水量均能到達潛水面形成補給，即潛在補給量。分析模型不同埋深的通量變化，當(dāng)埋深大于3 m時沒有向上的水流通量產(chǎn)生，因此選取3 m處的水流通量表征潛在補給量的變化。

如圖6所示，降水的年內(nèi)變化存在季節(jié)性差異，研究區(qū)降水集中在6—8月，雨季豐沛的降水有利于補給的形成，并產(chǎn)生相應(yīng)的補給峰值;旱季降水較少，土壤較為干旱，次數(shù)不多的少量降水主要補充包氣帶的水分虧缺，難以形成有效的補給。

圖5 2011—2012年潛在蒸發(fā)Ep和潛在蒸騰TpFig.5 Potential evaporation Epand potential transpiration Tpfrom 2011 to 2012

圖6 潛在補給量與降水的響應(yīng)關(guān)系Fig.6 Relationship between precipitation and recharge rate

降水年際變化主要用年降水量差異表征。天然降水條件下的包氣帶水均衡:

式中:P——降水量/cm;

E——土壤蒸發(fā)量/cm;

T——蒸騰量/cm;

Roff——地表徑流量/cm;

D——補給量/cm;

ΔW——土壤水總量變化量/cm。

可以看出，年降水量的變化會導(dǎo)致水均衡發(fā)生改變?；谀Ｐ陀嬎愕慕Y(jié)果，得出模擬期內(nèi)(1957—2013年)逐年的水均衡情況，并與當(dāng)年的降水量進行相關(guān)分析，得到年降水量與補給量、根系吸水量、蒸發(fā)量、地表徑流量的相關(guān)系數(shù)分別為:0.71、0.74、0.59、0.44。其中補給量與降水量的顯著正相關(guān)性表明年降水量的變化是導(dǎo)致補給量改變的主要因素。如圖4所示，研究區(qū)近60年的年降水量介于262.8～934.5 mm，平均值為514 mm，變差系數(shù)0.28。將系列年降水資料按降序排列，并逐個累加次數(shù)，計算得到了研究區(qū)年降水量的累積頻率曲線(圖7)。參照累積頻率曲線對研究區(qū)進行降水年型劃分，代表年份的年降水量及對應(yīng)的模型計算得出的年補給量如表4所示，在極端豐水條件下，研究區(qū)補給量達到151.5 mm/a，在相對枯水的年份補給量僅有30 mm/a，補給量隨著降水量的減少而減少。研究區(qū)由天然降水形成的多年平均補給量為66.6 mm/a，平均入滲補給系數(shù)為0.13。

圖7 年降水量累計頻率曲線Fig.7 Cumulative frequency curve of annual precipitation

表4 不同降水年型對應(yīng)降水量及補給量Table 4 Rainfall and recharge in different typical years

3.2 次降水特征對入滲補給的影響

降水量的年際和季節(jié)變化可以從宏觀角度判斷地下水的入滲補給量，但是單次降水的特征，如次降水量、降水歷時、降水強度等因素對于實際入滲補給過程的影響更為顯著。因此，將典型年模擬的第一天進行以小時為時間間隔的加密處理，并賦予不同雨量、不同降水歷時(降水強度)的單次降水(表5)，以分析次降水特征對于入滲補給規(guī)律的影響。

模擬結(jié)果表明，次降水量及降水強度共同決定了入滲補給和地表徑流的變化規(guī)律。如圖8所示，降水入滲補給系數(shù)隨降水強度的增大而減小，當(dāng)降水強度小于50 mm/h時，單次降水可以形成明顯的補給，在同一雨強下入滲補給系數(shù)隨降水量增加而增大，雨量很小的降水(小于10 mm)在補充前期包氣帶水分虧缺、同時被蒸發(fā)蒸騰消耗之后，無法形成補給;而當(dāng)降水強度大于50 mm/h時，即使降水量很大，形成的有效補給量也是有限的，補給系數(shù)則隨著降水量的增加逐漸趨近于零。地表徑流的產(chǎn)生是導(dǎo)致入滲補給系數(shù)減小的主要原因，表層土壤滲透性與降水強度的相對大小決定了地表徑流的形成。如表3所示，試驗場地表層土壤的飽和滲透系數(shù)為25 mm/h，當(dāng)降水強度小于土壤的滲透能力時，沒有地表徑流產(chǎn)生，地表入滲強度受降水控制;而當(dāng)降水強度大于25 mm/h時，地表徑流系數(shù)與降水強度呈現(xiàn)正相關(guān)性，地表入滲強度受土壤滲透性限制。

表5 次降水特征及模擬結(jié)果Table 5 Characteristics of single rainfall and corresponding simulation results

圖8 入滲補給系數(shù)、地表徑流系數(shù)與次降水強度的關(guān)系Fig.8 Relationship between single rainfall intensity and recharge coefficient and runoff coefficient

3.3 灌溉水對入滲補給的影響

作物種植期間的灌溉活動人為增加了表層的水量輸入，改變天然條件下的垂向入滲補給過程，如圖9所示，設(shè)計灌溉量最大的4水田塊年周期內(nèi)的各次灌溉均會形成相應(yīng)的補給峰值，旱季少量的降水難以形成有效的補給，而水量較大的灌溉則較易形成補給，小麥底墑水以及前兩次春灌均使?jié)撛谘a給通量增加，并在灌溉20～25d后達到峰值，此后緩慢減小直至補給結(jié)束。雨季降水使土壤保持較高的含水率，此時對玉米進行灌溉會形成明顯的補給。

圖9 灌溉對入滲補給的影響Fig.9 Relationship between irrigation and recharge rate

補給過程曲線表明各次灌溉對補給的影響各不相同。分析不同灌溉處理模型補給總量之間的差值，可以得到各次灌溉所形成的入滲補給量，進而計算單次灌溉的入滲補給系數(shù)。

各次灌溉會對淺層地下水形成不同程度的補給，這主要受氣候條件、土壤前期含水率以及作物生長周期等因素的控制(表6)。研究區(qū)7、8月份雨水豐沛，玉米生長期內(nèi)包氣帶水分虧缺較少，大部分灌溉水穿過根區(qū)入滲補給地下水。9月底小麥播種前澆底墑水時土壤含水率仍然較高，底墑水灌溉在補足包氣帶水分虧缺后會對淺層地下水形成補給。小麥生長階段，較少的降水造成土壤前期含水率較低，一部分灌溉水被用來補充包氣帶水分虧缺，另一方面，小麥生長需要大量水分(圖10)，每次春灌均有明顯的根系吸水響應(yīng)，伴隨春灌次數(shù)增加，根系吸水量明顯增多(圖11)，因此春灌入滲補給系數(shù)較小。較大的入滲補給系數(shù)表明灌溉量偏大，在實際條件下，應(yīng)基于當(dāng)年的降水情況適時決定灌溉次數(shù)及灌水量的大小。

表6 灌溉入滲補給系數(shù)Table 6 Irrigation recharge coefficient

圖10 根系吸水量變化Fig.10 Variation in root water uptake rate

3.4 研究區(qū)入滲補給評價結(jié)果對比

林丹［14］利用溴示蹤法對不同實驗田的入滲補給量進行評價，得出入滲補給量隨灌溉次數(shù)增加而增加的結(jié)論。其中3水田塊補給量出現(xiàn)異常小值，可能是由于優(yōu)先流入滲或示蹤劑投放、溴離子測定等實驗過程中的誤差造成的。與示蹤試驗相比，數(shù)值模型得到了相同的結(jié)果，其中1水和2水田塊模擬得到的補給量與示蹤試驗結(jié)果基本一致，0水田塊模擬值偏大，4水田塊則偏小。

利用數(shù)值模擬以及示蹤試驗評價入滲補給量均存在局限性。示蹤法受野外實際條件影響較大，評價結(jié)果具有一定的不確定性;而模型模擬法則是將野外條件簡化處理，需要盡可能詳細的資料來減小與實際條件的偏差。兩種方法得到的相似結(jié)果在一定程度上相互印證了評價結(jié)果的可靠性。

圖11 不同灌溉處理年根系吸水量Fig.11 Annual root water uptake quantity under different irrigation fields

圖12 數(shù)值模擬與Br示蹤結(jié)果對比Fig.12 Comparison of the annual recharge rates between numerical modeling and bromine tracing

近年來，國內(nèi)外學(xué)者針對華北平原入滲補給評價問題開展了大量理論和實例研究，筆者總結(jié)前人基于不同評價方法得出的結(jié)果，與本文模型計算結(jié)果進行對比。如表7所示，基于不同方法計算得到的華北平原降水灌溉入滲補給系數(shù)介于0.08～0.24之間，其中有灌溉情況下的入滲補給系數(shù)大于無灌溉條件，本文模擬計算得到的降水入滲補給系數(shù)為0.13、降水和灌溉綜合入滲補給系數(shù)為0.20(2水)，與前人研究結(jié)果接近，具有較高的可靠性。

表7 評價結(jié)果對比Table 7 Comparison of recharge coefficients

4 結(jié)論

(1)研究區(qū)近60年年降水量差異明顯(262.8～934.5 mm，變差系數(shù)0.28)，年入滲補給量與降水量顯著正相關(guān)，枯水年份降水入滲補給量為30 mm/a，豐水年入滲補給量為120～150 mm/a;多年平均降水入滲補給量為66.6 mm/a，平均降水入滲補給系數(shù)為0.13。

(2)次降水特征及表層包氣帶滲透性是決定實際入滲補給過程的主要因素，較小的次降水在補充包氣帶水分虧缺后無法形成補給，當(dāng)次降水強度大于飽和滲透系數(shù)(25 mm/h)時，入滲補給主要受表層土壤滲透性限制，入滲補給系數(shù)隨降水強度增大而減小。

(3)農(nóng)業(yè)灌溉增加了地表輸入水量，導(dǎo)致入滲補給量隨灌溉量增加而增加;小麥底墑水及玉米灌溉對應(yīng)入滲補給系數(shù)較大，表明實驗田灌溉量偏大，實際條件下應(yīng)基于當(dāng)年降水情況適時確定合理的灌水定額及灌溉次數(shù)。

［1］ 張人權(quán)，梁杏，靳孟貴，等.水文地質(zhì)學(xué)基礎(chǔ)［M］.6版.北京:地質(zhì)出版社，2011.［ZHANG R Q，LIANG X，JIN M G，et al.Fundamental of hydrogeology［M］. 6th ed. Beijing:Geology Publishing House，2011.(in Chinese)］

［2］ 李金柱.降水入滲補給系數(shù)綜合分析［J］.水文地質(zhì)工程地質(zhì)，2009，36(2):29-33.［LI J Z.An analysis of the coefficient of replenishment from infiltration of precipitation［J］.Hydrogeology ＆ Engineering Geology，2009，36(2):29-33.(in Chinese)］

［3］ 譚方穎，王建林，宋迎波，等.華北平原近45年氣候變化特征分析［J］.氣象，2010，36(5):40-45.［TAN F Y，WANG J L，SONG Y B，et al.Characteristics of climate change in the North China Plain for recent 45 years［J］.Meteorological Monthly，2010，36(5):40-45.(in Chinese)］

［4］ Scanlon B R，Reedy R C，Gates J B.Effects of irrigated agroecosystems ［J］. Water Resource Research，2010，46:1-14.

［5］ Qin D J，Qian Y P，Han L F，et al.Assessing impact of irrigation water on groundwater recharge and quality in arid environment using CFCs，tritium and stable isotopes，in the Zhangye Basin，Northwest China［J］.Journal of Hydrology，2011:194-208.

［6］ 汪丙國.地下水補給評價方法研究［D］.武漢:中國地質(zhì)大學(xué)(武漢)，2008.［WANG B G.Research on estimating methods of groundwater recharge:a case in North China Plain［D］.Wuhan:China University of Geosciences，2008.(in Chinese)］

［7］ Lu X H，Jin M G，Van Genuchten M Th，et al.Groundwater recharge at five representative sites in the Hebei Plain，China［J］.Ground Water，2011，49(2):286-294.

［8］ 王洪超，聶慶林，聶秋月，等.魯北平原區(qū)井灌回歸系數(shù)分析［J］.地下水，2008，30(3):17-23.［WANG H C，NIE Q L，NIE Q Y，et al.Analyzing the coefficient of irrigation return flow in Lubei Plain［J］.Groundwater，2008，30(3):17-23.(in Chinese).］

［9］ SimunekJ， van Genuchten M Th， Sejna M.Development and applications of the HYDRUS and STANMOD software package sand related codes［J］.Vadose Zone Journal，2008，7(2):587-600.

［10］ 張彥增，秦建文，喬光建.河北省平原區(qū)水面蒸發(fā)量變化趨勢及影響因素［J］.南水北調(diào)及水利科技，2011，9(4):63-65.［ZHANG Y Z，QIN J W，QIAO G J. Changing trends ofwatersurface evaporation and its influencing factors in the plain area ofHebeiprovince［J］. South-to-North Water Diversion and Water Science＆ Technology，2011，9(4):63-65.(in Chinese)］

［11］ 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:892-898.

［12］ 韓淑敏，程一松，胡春勝.太行山山前平原作物系數(shù)與降水年型關(guān)系探討［J］.干旱地區(qū)農(nóng)業(yè)研究，2005，23(5):152-158.［HAN S M，CHENG Y S，HU C S.Relationship between crop coefficient and precipitation pattern in the piedmont of Mt.Taihang［J］.Agricultural Research in the Arid Area，2005，23(5):152-158.(in Chinese)］

［13］ 付雪麗.冬小麥——夏玉米產(chǎn)量性能動態(tài)特征及其主要栽培措施效應(yīng)［D］.北京:中國農(nóng)業(yè)科學(xué)院作物科學(xué)研究所，2009.［FU X L.Quantitative analysis on yield performance dynamic process of winter wheat and summer maize and high yielding approaches［D］.Beijing:Institute of Crop Sciences of CAAS，2009.(in Chinese)］

［14］ 林丹.包氣帶變化及其對地下水補給的影響［D］.武漢:中國地質(zhì)大學(xué)(武漢)，2014.［LIN D.The changes of vadose zone and its impacts on groundwater recharge［D］. Wuhan: China University of Geosciences，2014.(in Chinese)］

［15］ 譚秀翠，楊金忠，宋雪航，等.華北平原地下水補給量計算分析［J］.水科學(xué)進展，2013，24(1):76-84.［TAN X C，YANG J Z，SONG X H，et al.Estimating of groundwater recharge in North China Plain［J］.Advances in Water Science，2013，24(1):77-84.(in Chinese)］

［16］ Huo S Y，Jin M G，Liang X，et al.Changes of vertical groundwater recharge with increase in thickness of vadose zone［J］.Journal of Earth Science.2014，25(6):1043-1050.

［17］ 朱長坤.河北平原三維水文地質(zhì)建模與降水入滲補給研究［D］.武漢:中國地質(zhì)大學(xué)(武漢)，2014.［ZHU C K.Three-dimensional hydrogeological modeling and precipitation recharge research in Hebei Plain［D］.Wuhan:China University of Geosciences，2014.(in Chinese)］