鄭州地下水均衡試驗場的改建工程
——試驗場監(jiān)測資料的推廣應(yīng)用核心問題

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

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

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

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

1 鄭州地下水均衡試驗場數(shù)據(jù)監(jiān)測成果

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

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

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

(3)

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

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

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

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

m=1-1/n

(4)

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

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

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

(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)驗方程都分別有1～2個獨立待定常數(shù)。分析水分特征曲線的變異性本質(zhì)是分析待定常數(shù)的變異性。

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

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

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

(7)

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

Δ——估值精度。

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

(8)

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

N=3.84(σ/Δ)2

(9)

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

N=3.84(CV/k)2

(10)

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

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

(11)

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

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

(12)

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

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

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

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

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

(13)

式中：N——所取測點的“對”數(shù)；

(14)

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

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

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

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

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

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

(15)

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

T——時間/d；

h——壓力水頭/cm；

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

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

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

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

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

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

Evt0f(STOP)

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

(16)

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

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

STOP——試驗柱表層日均飽和度；

Evt0——氣象觀測日均蒸發(fā)強度/(cm·h-1)。

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

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

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