喬曉英，肖 平，王繼華，李俊亭，王 林

(1.長(zhǎng)安大學(xué)水利與環(huán)境學(xué)院，陜西 西安 710054；2.干旱區(qū)地下水文與生態(tài)效應(yīng)教育部重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710054；3.河南省地質(zhì)環(huán)境監(jiān)測(cè)院，河南 鄭州 450016)

國(guó)際、國(guó)內(nèi)已經(jīng)建立了多個(gè)包氣帶—地下水系統(tǒng)原位監(jiān)測(cè)與試驗(yàn)基地[1]，例如加拿大Borden 場(chǎng)地[2]，美國(guó)MADE 試驗(yàn)場(chǎng)[3]、Hanford 場(chǎng)地[4]、德國(guó)的Krauthausen場(chǎng)地[5]，國(guó)內(nèi)武漢大學(xué)建立的大型蒸滲儀系統(tǒng)[6]、長(zhǎng)安大學(xué)在渭河關(guān)中平原和鄂爾多斯風(fēng)沙灘地建成的多功能地表—地下水原位試驗(yàn)基地、中國(guó)科學(xué)院禹城綜合試驗(yàn)站、中國(guó)地質(zhì)科學(xué)院水文地質(zhì)與環(huán)境地質(zhì)研究所衡水原位試驗(yàn)場(chǎng)等[7]，實(shí)現(xiàn)了對(duì)地表—地下水系統(tǒng)狀態(tài)變量實(shí)時(shí)監(jiān)測(cè)，為開(kāi)展多尺度地下水文過(guò)程的機(jī)理研究提供了基礎(chǔ)[8]。

河南省鄭州地下水原位試驗(yàn)場(chǎng)改建后，成為國(guó)內(nèi)設(shè)備較全、觀測(cè)數(shù)據(jù)類(lèi)較多、數(shù)據(jù)采集自動(dòng)化程度較高的原位試驗(yàn)場(chǎng)之一。其最大特點(diǎn)是試驗(yàn)柱中的組合巖性，是依據(jù)河南省五個(gè)地貌單元典型剖面經(jīng)專(zhuān)家概化后的非均質(zhì)巖性仿真而構(gòu)成，經(jīng)過(guò)實(shí)時(shí)監(jiān)測(cè)，可獲得試驗(yàn)場(chǎng)尺度下充水(降雨入滲)與釋水(潛水蒸發(fā))水文過(guò)程、13種典型巖性水分特征曲線，繼而構(gòu)建降水入滲或潛水蒸發(fā)預(yù)報(bào)聯(lián)合模型[9]。但實(shí)際應(yīng)用中如何將所建場(chǎng)地尺度模型推廣應(yīng)用于復(fù)雜的野外狀況，探求天然狀態(tài)下潛水面、包氣帶介質(zhì)對(duì)大氣降水與蒸發(fā)(騰)的響應(yīng)機(jī)制，是水文地質(zhì)環(huán)境地質(zhì)工作者長(zhǎng)期面臨的難題之一，其核心問(wèn)題是包氣帶參數(shù)的空間變異問(wèn)題。本文試圖探討一種有較高精度的大氣降水對(duì)地下水補(bǔ)給與潛水蒸發(fā)影響的評(píng)估方法，對(duì)于正確評(píng)估黃河沖積平原地下水資源的可持續(xù)開(kāi)采與生態(tài)環(huán)境保護(hù)有重要的理論與實(shí)踐意義。

1 鄭州地下水均衡試驗(yàn)場(chǎng)數(shù)據(jù)監(jiān)測(cè)成果

原鄭州地下水均衡試驗(yàn)場(chǎng)始建于1981年，1983年開(kāi)始監(jiān)測(cè)。試驗(yàn)柱是地下水均衡試驗(yàn)場(chǎng)的核心構(gòu)筑物[10]。改建后的試驗(yàn)柱介質(zhì)來(lái)自黃河北沖積平原、黃河南沖積平原、淮河沖積平原、南陽(yáng)盆地、豫西黃土丘陵等五個(gè)地貌單元的包氣帶(約7 m深度)巖性及其迭置關(guān)系，共計(jì)25個(gè)試驗(yàn)柱，仍采用原有地下水位埋深，分別為1，2，3，5，7 m五種深度。每個(gè)試驗(yàn)柱都設(shè)計(jì)有獨(dú)立供水系統(tǒng)，以保證“水位控制埋深線”處的水均衡。試驗(yàn)柱亦呈環(huán)形布置，放置在以排氣通道為中心5個(gè)地下監(jiān)測(cè)室中，共安裝設(shè)置了140個(gè)土壤負(fù)壓傳感器[11-12]，140個(gè)土壤含水率與溫度傳感器。針對(duì)鄭州地下水均衡試驗(yàn)場(chǎng)周邊高樓林立，失去了應(yīng)用氣象部門(mén)信息的可能。因而安裝了3 m高主桿的HOBO野外氣象站。采用U30-NRC采集器，采集間隔最小可自定義為1 s，另外還有2套自行建立蒸發(fā)(含降水)對(duì)比試驗(yàn)觀測(cè)。一套是E—601，另一套是自行設(shè)計(jì)的雙圈水面蒸發(fā)筒[13-14]。

2 均衡試驗(yàn)場(chǎng)與野外現(xiàn)場(chǎng)的包氣帶建模方法對(duì)比

首先是研究目的各不相同。前者是結(jié)合擬解決的實(shí)際問(wèn)題，比如針對(duì)包氣帶的多層結(jié)構(gòu)，探討包氣帶在某些特定條件下的水分運(yùn)移規(guī)律；后者則是解決區(qū)域潛水面蒸發(fā)與入滲量等實(shí)際問(wèn)題。因而前者是研究基礎(chǔ)，后者是研究目標(biāo)。二者最大差別是包氣帶參數(shù)的選取和土質(zhì)結(jié)構(gòu)的不同。在試驗(yàn)場(chǎng)尺度下，試驗(yàn)土柱巖性結(jié)構(gòu)是仿真了野外實(shí)際條件，單就某一層來(lái)講是均一的。根據(jù)試驗(yàn)柱中觀測(cè)數(shù)據(jù)(含水率、負(fù)壓、溫度)確定的參數(shù)可能與野外實(shí)際(按顆分結(jié)果分類(lèi)定名的同一巖性)參數(shù)值有較大的差異。這種差異主要來(lái)源于包氣帶土質(zhì)沉積過(guò)程(沉積環(huán)境)的隨機(jī)性，諸如沉積過(guò)程中水流速度、水中泥沙結(jié)構(gòu)及影響沉積過(guò)程的其他自然環(huán)境等因素。因而實(shí)際上對(duì)于同一種巖性參數(shù)來(lái)說(shuō)，在時(shí)空分布上也會(huì)有變異。二者相同之處在于構(gòu)建的數(shù)學(xué)模型、邊界條件基本一致，初始條件的確定需結(jié)合野外條件選擇有所不同。國(guó)內(nèi)曾在河南省商丘大吳莊進(jìn)行過(guò)較大地塊(長(zhǎng)70 m、寬40 m)參數(shù)變異性研究，研究目標(biāo)是土壤，取樣深度最大也只有40 cm[15]。其研究結(jié)論對(duì)于解決區(qū)域性(千米級(jí))包氣帶參數(shù)變異性問(wèn)題具有一定的參考價(jià)值。本次研究側(cè)重于模型的空間尺度、原狀土樣測(cè)試、包氣帶參數(shù)個(gè)數(shù)等方面做進(jìn)一步改進(jìn)。

圖1 改建后的鄭州地下水均衡試驗(yàn)場(chǎng)俯視圖[6]Fig.1 Top view of the Zhengzhou groundwater balance experiment site after reconstruction

從均衡試驗(yàn)場(chǎng)監(jiān)測(cè)成果到實(shí)際應(yīng)用，不可回避的核心技術(shù)是區(qū)域條件下包氣帶參數(shù)的變異性問(wèn)題。國(guó)際學(xué)術(shù)界自20世紀(jì)70年代提出包氣帶物理參數(shù)的空間變異性以來(lái)，多數(shù)學(xué)者認(rèn)為，包氣帶參數(shù)變異性是指包氣帶介質(zhì)的物理參數(shù)變異是空間的函數(shù)[15]，而同一巖性在不同深度的物理參數(shù)變異則是時(shí)間的函數(shù)。以往土壤特性空間變異性研究基于觀測(cè)或取樣資料，分析土壤各特性參數(shù)的空間變化特征、參數(shù)間的空間關(guān)系，以確定合理的取樣點(diǎn)數(shù)目，從而對(duì)未測(cè)點(diǎn)的參數(shù)進(jìn)行最優(yōu)估值等。進(jìn)一步結(jié)合標(biāo)定理論的應(yīng)用來(lái)分析和預(yù)測(cè)狀態(tài)變量的空間分布[16]。例如在歐美國(guó)家，基本上是采用野外采樣，通過(guò)室內(nèi)試驗(yàn)測(cè)試包氣帶物理參數(shù)，建立仿真大氣降水與蒸發(fā)(騰)數(shù)學(xué)模型，預(yù)測(cè)大氣降水對(duì)地下水的補(bǔ)給與蒸發(fā)(騰)。這種方法理論嚴(yán)謹(jǐn)，但還是沒(méi)有突破需要現(xiàn)場(chǎng)參數(shù)測(cè)試的局限，因而解決實(shí)際問(wèn)題尚有一定差距。而且研究成果多集中在表層土壤水分的空間結(jié)構(gòu)分析，對(duì)包氣帶剖面土壤水分的空間結(jié)構(gòu)特性研究不足，缺乏系統(tǒng)分析實(shí)測(cè)土壤水分空間結(jié)構(gòu)變化規(guī)律的實(shí)際應(yīng)用[9]。

3 包氣帶參數(shù)變異性研究思路與方法

針對(duì)包氣帶參數(shù)變異性問(wèn)題，筆者認(rèn)為有兩種研究思路可借鑒：(1)用概率模型表述；(2)數(shù)學(xué)地質(zhì)模型表述。基于此，本文試圖通過(guò)原狀土采樣、分析測(cè)試、數(shù)學(xué)建模等手段，探討均衡試驗(yàn)場(chǎng)監(jiān)測(cè)數(shù)據(jù)推廣到實(shí)際應(yīng)用的一種方法(圖2)。其中，構(gòu)建包氣帶野外參數(shù)模型是關(guān)鍵步驟，包括按照一定精度和置信水平，確定取樣數(shù)目；根據(jù)半方差和自相關(guān)圖分析土壤特性的空間結(jié)構(gòu)(方向性和相關(guān)距離)；應(yīng)用Kriging法進(jìn)行內(nèi)插計(jì)算。

圖2 野外監(jiān)測(cè)與均衡試驗(yàn)場(chǎng)監(jiān)測(cè)集成流程圖Fig.2 Flow chart of integration of field monitoring and balanced experiment site monitoring

包氣帶參數(shù)變異性研究方法主要有傳統(tǒng)的統(tǒng)計(jì)方法、時(shí)間序列分析方法、地統(tǒng)計(jì)學(xué)方法、隨機(jī)模擬方法、分?jǐn)?shù)維方法以及GIS的空間變異分析方法[18]。地統(tǒng)計(jì)學(xué)方法由于它注重變量因子的空間過(guò)程，考慮其空間分布特征和空間自相關(guān)而得到廣泛應(yīng)用[19]。它是通過(guò)變異函數(shù)可以確定和比較變量因子的空間變異程度及空間變異尺度, 以提供地理學(xué)、生態(tài)學(xué)和土壤學(xué)對(duì)自然現(xiàn)象及過(guò)程的空間變異特征解釋。目前將地統(tǒng)計(jì)學(xué)方法和GIS 結(jié)合起來(lái),一方面利用GIS 的空間分析功能,利用計(jì)算機(jī)的先進(jìn)技術(shù)方便地實(shí)現(xiàn)地統(tǒng)計(jì)學(xué)的計(jì)算內(nèi)插和制圖要求;另一方面能夠很好地描述因子的空間結(jié)構(gòu)特征及其時(shí)間變化規(guī)律，是分析土壤水分空間特征及其變異規(guī)律最為有效的方法之一[17,20-22]。

3.1 包氣帶參數(shù)的概率模型

用經(jīng)典統(tǒng)計(jì)學(xué)的基本理論，對(duì)不同位置、深度土壤含水量、干容重、顆分(黏粒和粉粒組成)、飽和滲透系數(shù)、孔隙率以及含水率、負(fù)壓等參數(shù)統(tǒng)計(jì)其特征值。例如均值、中值、最大值、最小值、極差、標(biāo)準(zhǔn)差、離散系數(shù)。當(dāng)離散系數(shù)Cv≤0.1為弱變異性, 0.1

現(xiàn)以干容重(λ；g/cm3)為例分析其統(tǒng)計(jì)分布特征。同一巖性在不同位置(平面與剖面)所取樣品干容重最大值、最小值不同。并且這種變化是隨機(jī)的，可依據(jù)統(tǒng)計(jì)方法獲得均方差(σ)與數(shù)學(xué)期望(α)。將干容重視為隨機(jī)變量δ，若

p(α-σ<δ<α+σ)=0.68

(1)

p(α-2σ<δ<α+2σ)=0.956

(2)

p(α-3σ<δ<α+3σ)=0.997

(3)

就認(rèn)為干容重γ服從正態(tài)分布。

同理，對(duì)于包氣帶介質(zhì)顆粒分析的不均勻系數(shù)(Cu)或曲率系數(shù)(Cc)、飽和滲透系數(shù)(K)及給水度(μ)的變異性都可以如此分析。

但是含水率與負(fù)壓之間為函數(shù)關(guān)系，不能用上述分析方法，表述二者關(guān)系稱(chēng)為土壤水分特征曲線。目前統(tǒng)計(jì)模型有Van Geunchten模型、Brooks-Corey模型、Gardner-Russo模型、Campbell模型、Williams模型、Mckee和Bumb模型、Frdlund和Xing模型、Broadbridge-White模型、Burdine模型、Mualem模型等十余種[23]，以Van Geunchten模型為例，土壤含水率與負(fù)壓經(jīng)驗(yàn)方程：

θ=θγ+(θS-θγ)[(1+α|h|)n]-m

m=1-1/n

(4)

其中，α、n和m為統(tǒng)計(jì)經(jīng)驗(yàn)參數(shù)；θs、θr、αw和nw為擬合吸濕過(guò)程參數(shù)；θs、θr、αd和nd為擬合脫濕過(guò)程參數(shù)；θs為飽和含水率；θr為殘留含水率。

鄭州地下水均衡試驗(yàn)場(chǎng)的巖性概化為：細(xì)砂、粉細(xì)砂、粉質(zhì)黏土、粉土、黃土狀粉土、黃土狀粉質(zhì)黏土夾粉質(zhì)黏土、黏土7大類(lèi)13種。經(jīng)過(guò)一年以上含水率、負(fù)壓監(jiān)測(cè)就可以獲得13種巖性建立脫濕與吸濕兩種狀態(tài)下的水分特征曲線模型。

滲透系數(shù)與負(fù)壓關(guān)系的Van Geunchten經(jīng)驗(yàn)方程為：

(5)

式中：K(h)——滲透系數(shù)/ (m·d-1)；

Ks——飽和滲透系數(shù)/ (m·d-1)，非飽和區(qū)則為壓力水頭的函數(shù)。

K(S)=KsS0.5[1-(1-S1/m)m]2

(6)

由式(4)～(6)可以看出，控制含水率或滲透系數(shù)經(jīng)驗(yàn)方程都分別有1～2個(gè)獨(dú)立待定常數(shù)。分析水分特征曲線的變異性本質(zhì)是分析待定常數(shù)的變異性。

3.2 區(qū)域尺度包氣帶參數(shù)變異地學(xué)統(tǒng)計(jì)模型

3.2.1合理選擇場(chǎng)地取樣數(shù)量

包氣帶參數(shù)變異性是一個(gè)隨機(jī)問(wèn)題，首要考慮樣本容量的選擇。如果取樣點(diǎn)數(shù)目過(guò)少，所得結(jié)果缺乏代表性，甚至是錯(cuò)誤的結(jié)論，但觀測(cè)點(diǎn)數(shù)目過(guò)多，則需要消耗較多的人力、物力。用有限觀測(cè)值去估計(jì)該參數(shù)的均值〈或期望)時(shí)，為保證足夠的精度，取樣或觀測(cè)點(diǎn)的數(shù)目應(yīng)合理確定[15]。若認(rèn)為所采樣品獨(dú)立，而且取樣容量足夠，中心極限定理成立，可由置信區(qū)間滿足取樣數(shù)目合理性。

(7)

式中：Pl——置信水平，一般取值為95%；

Δ——估值精度。

由中心極限定理隨機(jī)變量為標(biāo)準(zhǔn)正態(tài)分布：

(8)

結(jié)合式(7)、(8)，可知取樣點(diǎn)數(shù)目N為：

N=3.84(σ/Δ)2

(9)

若取Δ=κμ，則可表述為：

N=3.84(CV/k)2

(10)

當(dāng)置信水平為95%，k=10%時(shí)，則取樣數(shù)目對(duì)應(yīng)于弱變異性，N<4 ，中等變異性，N=4～400，強(qiáng)變異性，N>400。

實(shí)際應(yīng)用中，用樣本方差代替總體方差，由統(tǒng)計(jì)學(xué)原理：

(11)

其中，λα·f為t分布的特征值，可以查表得到，這樣滿足要求的取樣數(shù)目即為：

N=λ2α·f(S/Δ)2

(12)

根據(jù)式(12)可知，取樣數(shù)目和所取置信水平及精度要求有關(guān)。

以河南省鄭州地下水均衡試驗(yàn)場(chǎng)的推廣應(yīng)用為例，擬選擇黃河南岸沖洪積平原某一試驗(yàn)場(chǎng)地為長(zhǎng)9 km，寬6 km，布設(shè)網(wǎng)度為500 m，共計(jì)24 km2，117個(gè)采樣點(diǎn)，見(jiàn)圖3。試驗(yàn)場(chǎng)地在地下水位埋深7米內(nèi)有7種巖性組合，預(yù)計(jì)一種巖性的采樣容量大概是100，一個(gè)采樣點(diǎn)不同深度采樣5～7個(gè)，共計(jì)采樣約600個(gè)。這個(gè)樣本容量是區(qū)別于以往研究的特點(diǎn)之一，對(duì)于描述區(qū)域尺度包氣帶參數(shù)變異性的模型基本滿足。

圖3 原位采樣點(diǎn)的布設(shè)Fig.3 Layout of in situ sampling points

此種原位采樣點(diǎn)布置方案的特點(diǎn)是可組成500 m×500 m的網(wǎng)度，從而為獲得包氣帶參數(shù)變異性模型與取樣網(wǎng)度尺寸的關(guān)系提供借鑒。包氣帶樣品為原狀樣，區(qū)別與以往擾動(dòng)樣，直徑95 mm，高100 mm，擬選用取樣效率較高的直推履帶式取樣鉆(Geoprobe，6620DT)取樣。并采用自行研制測(cè)量原狀樣的滲透系數(shù)與孔隙率及水分特征曲線測(cè)試儀器進(jìn)行測(cè)試。

3.2.2包氣帶參數(shù)空間分布相關(guān)性的半方差分析

(13)

式中：N——所取測(cè)點(diǎn)的“對(duì)”數(shù)；

(14)

表1 常用半方差函數(shù)擬合方程

3.2.3最優(yōu)內(nèi)插的Kriging法

在半方差分析的基礎(chǔ)上，可以對(duì)未知點(diǎn)的參數(shù)值進(jìn)行最優(yōu)內(nèi)插估計(jì)，即所謂的Kriging法。Kriging最優(yōu)內(nèi)插法利用了區(qū)域化變量的原始數(shù)據(jù)和半方差函數(shù)的結(jié)構(gòu)性，對(duì)未知采樣點(diǎn)的區(qū)域化變量的取值進(jìn)行線性無(wú)偏最優(yōu)估計(jì)，最終生成研究對(duì)象的空間分布。相比于一般線性?xún)?nèi)插方法，由于其方差較小，因而估值精度較高[16]。

目前在ArcGIS的Geostatistical Analyst模塊采用kringing插值。需要注意的是標(biāo)值前將特異性值剔除，才能避免插值結(jié)果偏離實(shí)際值，可以利用該模塊的數(shù)據(jù)檢查工具(ESDA)來(lái)完成。

4 區(qū)域潛水面入滲或蒸發(fā)量估算

獲得野外包氣帶參數(shù)變異特征后，結(jié)合均衡試驗(yàn)場(chǎng)監(jiān)測(cè)資料，構(gòu)建降水入滲與蒸發(fā)條件下的數(shù)學(xué)模型：

(15)

式中：θ——含水量/( cm3· cm-3)；

T——時(shí)間/d；

h——壓力水頭/cm；

K(θ)——非飽和滲透系數(shù)/(cm·h-1)；

R——植物根系吸水源匯項(xiàng) /(cm3·cm-3·h-1)；

Z——垂直坐標(biāo)軸，將坐標(biāo)原點(diǎn)選在地面，向下設(shè)為正。

試驗(yàn)柱數(shù)學(xué)模型下邊界取定水頭邊界(壓力水頭為零)，用實(shí)測(cè)下邊界流量進(jìn)行校驗(yàn)。

試驗(yàn)柱上邊界有三種狀態(tài)：有壓入滲、無(wú)壓入滲和蒸發(fā)。有壓和無(wú)壓入滲可分別取變水頭邊界和流量邊界；蒸發(fā)邊界較復(fù)雜，下邊界取定水頭值，各類(lèi)巖性水分特征曲線取實(shí)際測(cè)定值，以試驗(yàn)柱監(jiān)測(cè)數(shù)據(jù)為基礎(chǔ)，用數(shù)值模型方法反求上界蒸發(fā)強(qiáng)度EvtTop，結(jié)合表層(如5 cm) 實(shí)測(cè)平均含水率值θTop(或飽和度STOP)，統(tǒng)計(jì)表層經(jīng)驗(yàn)蒸發(fā)規(guī)律函數(shù)EvtTop(θTop)，作為上部非線性通量邊界條件：

EvtTop(θTop)=Evt0f((θTop-θγ)/(θs-θγ))=

Evt0f(STOP)

STOP=(θTop-θγ)/(θs-θγ)

(16)

式中：EvtTop——試驗(yàn)柱表層(如5 cm)日均蒸發(fā)強(qiáng)度/(cm·h-1)；

θTop——實(shí)驗(yàn)柱表層日均含水率/(cm3·cm-3)；

STOP——試驗(yàn)柱表層日均飽和度；

Evt0——?dú)庀笥^測(cè)日均蒸發(fā)強(qiáng)度/(cm·h-1)。

試驗(yàn)柱初始條件根據(jù)所測(cè)的土壤含水量或負(fù)壓確定。

另外，當(dāng)試驗(yàn)柱數(shù)學(xué)模型建立之后，可在試驗(yàn)柱上邊界處置不同溶質(zhì)濃度，探討在降水與蒸發(fā)條件下的溶質(zhì)運(yùn)移特征，并建立包氣帶水分、熱量、溶質(zhì)運(yùn)移耦合模型，為土壤面狀污染預(yù)測(cè)提供信息。

黃河南沖洪積區(qū)野外試驗(yàn)場(chǎng)地推廣應(yīng)用實(shí)踐，可為黃河北沖洪積區(qū)、淮河沖積平原、南陽(yáng)盆地、豫西黃土丘陵區(qū)潛水面蒸發(fā)與入滲補(bǔ)給的現(xiàn)場(chǎng)模擬提供示范作用，對(duì)于探討黃河沖積平原降雨入滲補(bǔ)給或潛水蒸發(fā)特征，評(píng)價(jià)區(qū)域地下水資源的開(kāi)發(fā)潛力提供一定的科學(xué)依據(jù)。