1. 引言
在我國許多地區(qū)，由于水資源嚴(yán)重短缺，超采地下水，導(dǎo)致大范圍的地下水位急劇下降，從而引起區(qū)域性的地面沉降、海水入侵以及濕地退化等生態(tài)環(huán)境問題，嚴(yán)重影響了區(qū)域經(jīng)濟(jì)社會可持續(xù)發(fā)展。如何合理地利用地下水資源，對區(qū)域和流域地下水、地面水進(jìn)行優(yōu)化配置，通過建立健全的流域水循環(huán)體系，保護(hù)地下水環(huán)境，恢復(fù)流域生態(tài)，實(shí)現(xiàn)水資源的可持續(xù)利用以及人與自然的協(xié)調(diào)，是當(dāng)今水資源研究領(lǐng)域十分重要的課題。
要合理地利用區(qū)域地下水資源，首要問題是要對區(qū)域地下水資源進(jìn)行符合實(shí)際的模擬，摸清區(qū)域地下水的運(yùn)動規(guī)律。然而，由于區(qū)域地下水模擬涉及大尺度地下水運(yùn)動，而傳統(tǒng)的地下水模擬手段和方法在處理大尺度地下水運(yùn)動模擬中所涉及的有限元計算中的時空尺度的選擇、水位和含水層信息不足等問題上方面存在困難，引起由于信息不足帶來的較大計算誤差。因此，針對尺度和信息不足問題，建立大區(qū)域地下水運(yùn)動模擬的理論和方法十分必要。
本文在融合地質(zhì)統(tǒng)計、逆問題理論和地下水運(yùn)動理論的基礎(chǔ)上，提出了建立大區(qū)域地下水運(yùn)動模擬的理論和方法。首先，提出了大區(qū)域地下水有限元計算的時空尺度理論公式；其次，針對信息不足，提出了大區(qū)域地下水位推定的方法；第三，根據(jù)推定的地下水位，運(yùn)用逆解析理論法對大區(qū)域地下水含水層透水系數(shù)進(jìn)行逆推定；最后，給出了運(yùn)用實(shí)例。
2. 大區(qū)域地下水有限元計算的時空尺度
在大區(qū)域地下水?dāng)?shù)值計算中，首先需要選擇研究對象所需劃分的計算網(wǎng)格大小或確定計算時間步長大小。尤其在模擬河道網(wǎng)基礎(chǔ)上的分布式水文模型和大區(qū)域地下水模型進(jìn)行地表水與地下水偶合計算時，同時兼顧河道網(wǎng)的數(shù)值地形網(wǎng)格大小和地下水網(wǎng)格大小，既滿足水文計算的精度同時使地下水計算不至于失穩(wěn)就顯得至關(guān)重要。然而，由于大區(qū)域地下水計算中，由于范圍較大，人們往往希望選擇較大的網(wǎng)格進(jìn)行計算，以減少計算工作量。這往往導(dǎo)致計算上的不穩(wěn)定。
對于地下水有限元計算的時空尺度（時空步長）的選擇，1977年，Newman等對二維地下水運(yùn)動方程的有限元解法中的不穩(wěn)定問題進(jìn)行了分析[1]，推測不合理的時間步長導(dǎo)致了混合型差分的出現(xiàn)，由此導(dǎo)致了解的不穩(wěn)定問題。針對Neman等的推測，Zhang（1992）[2] 和Wood（1996）[3]提出了二維地下水有限元計算的時間步長的條件。2002年，張祥偉等[4]運(yùn)用最大最小原理對大尺度二維和準(zhǔn)三維地下水有限元計算的不穩(wěn)定問題進(jìn)行了理論分析，提出了二維和準(zhǔn)三維地下水有限元計算的時空步長的理論公式。即對于準(zhǔn)三維地下水有限元計算的時空步長分別為：
空間步長條件：
3. 大區(qū)域地下水位的推定
在進(jìn)行大區(qū)域地下水計算中，需要根據(jù)實(shí)測的地下水位推定初始流場，以便檢驗(yàn)地下水?dāng)?shù)值計算精度、進(jìn)行非恒定地下水計算以及識別含水層參數(shù)。
對于初始流場的推定的方法，通常有Universal Kriging方法[5-6]，UK法的計算行列比較大，計算比較復(fù)雜。Newman等（1984）[7] 和Sun（1999）[8] 運(yùn)用Residual Kriging（RK）法進(jìn)行區(qū)域地下水位的推定，也就是將實(shí)測地下水趨勢面去除得到正態(tài)殘差，將正態(tài)殘差運(yùn)用Ordinary Kriging法進(jìn)行面上殘差的推定，再加上實(shí)測地下水面的趨勢得到三角形網(wǎng)格上各點(diǎn)的推定地下水位值。
上述方法在大區(qū)域地下水計算中遇到信息不足的問題，當(dāng)實(shí)測地下水位信息不足時，推定的初始流場會帶來較大的誤差，本文根據(jù)地形水文學(xué)的原理[9]，即地下水與地形之間存在的相關(guān)性，提出運(yùn)用數(shù)值地形模型（DEM）中的地形標(biāo)高作為輔助信息，修正實(shí)測地下水位得到的地下水趨勢面，然后，運(yùn)用Ordinary Kriging方法對修正后趨勢面得到的正態(tài)殘差對三角形網(wǎng)格點(diǎn)上的殘差進(jìn)行推定，每各網(wǎng)格點(diǎn)上推定的殘差再加上修正的趨勢值得到面上地下水位的推定值。本文將此方法定名為ROKMT（Residual Ordinary Kriging with Modified Trend）方法[10-11]，即地下水位由以下兩部分構(gòu)成：
大區(qū)域地下水位推定的方法如圖1。
4. 區(qū)域地下水透水系數(shù)的識別
在地下水流動模擬中，一般根據(jù)含水層的水文地質(zhì)條件和土壤特性，選定透水系數(shù)的取值。在大區(qū)域地下水計算中，含水層物理特性的信息常常十分有限，即使知道含水層的土壤類型，但由于同類土壤的透水系數(shù)的變化范圍很大，最小值與最大值之間甚至有100倍以上的變化幅度。這為正確選擇地下水參數(shù)帶來困難，本文針對式（5）的二維地下水運(yùn)動，根據(jù)推定的地下水位值，運(yùn)用Guass-Newton法對大區(qū)域地下水含水層透水系數(shù)進(jìn)行識別。
二維地下水運(yùn)動的基本方程為：
其中，h為地下水位，T為透水量系數(shù)，q為地下水涵養(yǎng)量或揚(yáng)水量。
方程式（5）的透水量系數(shù)的識別，有直接法和間接法。直接法是根據(jù)已知的地下水位值，直接從式（5）中反求透水量系數(shù)T。直接法的問題是，由于地下水位中含有誤差，細(xì)小的誤差將導(dǎo)致水位的偏微分很大的誤差，計算穩(wěn)定性差，而且往往引起不適定問題，大大降低了參數(shù)的計算精度。本文運(yùn)用間接法識別水文地質(zhì)參數(shù)。即給定含水層透水量T的初始值進(jìn)行迭代計算，根據(jù)3中推定的地下水位值和反復(fù)計算得到的計算水位值的殘差平方和最小作為目標(biāo)函數(shù)，計算最優(yōu)的透水量值T。目標(biāo)函數(shù)為：
其中，hobs為ROKMT法的推定地下水位值。Tk+1為第K次迭代的透水量系數(shù)。通過下式計算，
根據(jù)含水層的物理特性，透水量系數(shù)需滿足以下限制條件：
當(dāng)透水量系數(shù)識別完成后，運(yùn)用不規(guī)則三角形差分（TFDM）法進(jìn)行大區(qū)域地下水計算。綜合3和4的計算步驟，大區(qū)域地下水計算方法的如圖2。
5. 運(yùn)用實(shí)例-Sarobetsu濕地地下水模擬
本文將所提出的方法運(yùn)用到日本北海道Sarobetsu濕地地下水模擬中。
5.1 研究對象的概況
Sarobetsu濕地位于日本最北端（圖3），為日本最大的以水臺癬為主的高層濕地，面積635km2。近年來，由于農(nóng)田開墾等土地開發(fā)利用，地下水位降低，濕地出現(xiàn)退化，海水上朔，生態(tài)環(huán)境惡化。最主要的表現(xiàn)是：
5.1.1 濕地指示性植物減少。對地下水位敏感的濕地植物水臺癬面積逐漸減少，低竹類植物風(fēng)長，面積逐年擴(kuò)大，植物耗水量增加，植物葉面蒸散發(fā)增強(qiáng)。
5.1.2 湖泊面積減小。由于濕地干旱化導(dǎo)致1993年濕地內(nèi)三大湖泊兜湖、Paken湖和Peken湖的面積比1975年減少10%。三湖泊周邊的水臺癬急劇減少，生態(tài)環(huán)境惡化。
5.1.3 地面沉降。濕地中心區(qū)地面最大沉降達(dá)到50厘米。
5.1.4 海水上朔。由于地面沉降造成和Sarobetsu河的枯水期流量減少，海水沿Sarobetsu河口上朔流入Paken湖和Peken湖，破壞了湖周邊生態(tài)。
5.1.5 Sarobetsu河下游洪水泛濫。由于Sarobetsu濕地的地面沉降和海水的頂托，造成下游洪水泛濫。
因此，為恢復(fù)濕地和減少地面沉降，首要課題是研究該區(qū)域地下水位恢復(fù)的措施。前提是研究該區(qū)域地下水涵養(yǎng)機(jī)理、地下水的運(yùn)動規(guī)律。
5.2 地形和地質(zhì)
Sarobetsu濕地的地形為從北到南沿海岸線逐漸降低。Sarobetsu河和海岸之間為5-20米的沙丘地帶，濕地地面高程在2米至8米之間，周邊為起伏的臺地和丘陵（如圖4）。
該區(qū)域地質(zhì)上主要分為兩層，即更別層以上為泥巖層，為由南向北的向斜構(gòu)造，從東、北、西向日本海方向嵌入。更別層以上為30-40米的洪積層，主要以透水性強(qiáng)的泥炭和沙質(zhì)層為主。更別層以下為化石地下水，與上部地下水交換很少。本研究以更別層以上低地的非承壓地下水運(yùn)動為對象。
5.3 地下水位空間分布的推定
運(yùn)用本文提出的ROKMT法對Sarobetsu濕地地下水位進(jìn)行推定。
5.3.1 網(wǎng)格的劃分
基于250x250m的DEM地形標(biāo)高擴(kuò)展成500x500m的矩形網(wǎng)格，每個矩形網(wǎng)格分為兩格三角形。全對象區(qū)域分為1051個節(jié)點(diǎn)，1903個三角形網(wǎng)格（如圖5）。
5.3.2 地下水位的推定結(jié)果
根據(jù)1997年72個地下水監(jiān)測點(diǎn)同步地下水監(jiān)測值、5個Sarobetsu河的水位值和3個湖泊的水位值，運(yùn)用ROKMT法推定地下水位。另外的10個實(shí)測地下水位作為驗(yàn)證點(diǎn)（圖6）。
地下水位的推定結(jié)果如圖7。
10個驗(yàn)證點(diǎn)的平均誤差（ME），平均絕對誤差（MAE）和平方根二乘誤差（RMSE）分別為0.005m，0.437m和0.598m。從三種誤差看，推定的地下水位空間分布具有較好的精度，可以用來識別透水量系數(shù)。
5.3.3 透水系數(shù)的識別結(jié)果
�。�1）邊界條件的設(shè)定
海岸線一側(cè)和研究區(qū)內(nèi)部的河流、湖泊作為定水頭邊界。其他邊界位于丘陵的邊沿作為流量邊界。
�。�2）涵養(yǎng)率（入滲補(bǔ)給系數(shù)）和透水量系數(shù)的界限值的設(shè)定
運(yùn)用分布式Tank水文模型，對1997年降雨徑流的模擬，以及水平衡法水量平衡計算，降水對地下水的涵養(yǎng)率取0.23。根據(jù)研究區(qū)域的水文地質(zhì)條件，確定透水量系數(shù)的界限值為0.0043m2/d≤T≥6800m2/d。
采用Guass-Newton法和有限元方法，使用非融雪期1997年6月地下水位資料，進(jìn)行大區(qū)域地下水含水層透水系數(shù)進(jìn)行識別。識別的透水量除以含水層厚度得到含水層透水系數(shù)，結(jié)果如圖8。
5.3.4 98年-2000年Sarobetsu濕地地下水流動的非恒定模擬
為恢復(fù)濕地地下水位和生態(tài)，對1998年-2000年Sarobetsu濕地地下水流動進(jìn)行了非恒定模擬。
�。�1）涵養(yǎng)率的設(shè)定
由于濕地在每年的11月底次年的5月為降雪期，融雪對地下水有很大的補(bǔ)給作用，本文山本等人[12]的數(shù)字積雪-融雪模型進(jìn)行了模擬。模型中考慮地形標(biāo)高、氣溫、日照量，取2℃以下時的降水為降雪，同時根據(jù)風(fēng)速修正降雪量。1997年11月—2000年6月Sarobetsu濕地積雪和融雪模擬結(jié)果如圖9。
�。�2）積雪和融雪的計算
涵養(yǎng)率在非積雪期取0.23，在6月至9月期間，由于植物的蒸散發(fā)強(qiáng)烈，涵養(yǎng)率取0.21。在融雪期，由于地下水主要靠融雪補(bǔ)給，因此，1-3月融雪量的90%作為地下水涵養(yǎng)，4月-5月融雪量的60%作為地下水涵養(yǎng)。儲留系數(shù)S取0.2。
�。�3）模擬結(jié)果
運(yùn)用TFDM法對1998年1月—2000年12月Sarobetsu濕地地下水非恒定運(yùn)動進(jìn)行模擬，2000年的模擬結(jié)果如圖10。
（4）模擬結(jié)果的驗(yàn)證
對10個驗(yàn)證點(diǎn)的地下水位的非恒定計算結(jié)果進(jìn)行了評價，10個點(diǎn)的ME為-0.063m，0.691m和0.731m。其中，第7點(diǎn)和第10點(diǎn)的誤差較大，分別為1.67m和2.07m，原因是7點(diǎn)和第10點(diǎn)為丘陵山區(qū)，巖石較多，計算時涵養(yǎng)率設(shè)定過高。
通過地下水運(yùn)動的模擬，可以對濕地過去10-20年的地下水空間分布進(jìn)行再現(xiàn)，結(jié)合不同年份濕地植物遙感信息和地下水空間為濕地生態(tài)恢復(fù)找到最佳地下水調(diào)控面和濕地植物退化的臨界點(diǎn)，為濕地恢復(fù)尋求根據(jù)。其次，通過地下水運(yùn)動空間分布的模擬，可以為恢復(fù)濕地的水工程的方式、規(guī)模提供依據(jù)。
6. 結(jié)論
本文在融合地質(zhì)統(tǒng)計、逆問題理論和地下水運(yùn)動理論的基礎(chǔ)上，提出了建立大區(qū)域地下水運(yùn)動模擬的理論和方法，并運(yùn)用到日本Sarobetsu濕地地下水模擬中，為恢復(fù)濕地提供了科學(xué)依據(jù)。此外，本文提出的方法也可以運(yùn)用到其他大區(qū)域地下水的模擬和濕地保護(hù)中。
