大區(qū)域地下水模擬的理論和方法

申請免費試用、咨詢電話：400-8352-114

　　1. 引言

　　在我國許多地區(qū)，由于水資源嚴(yán)重短缺，超采地下水，導(dǎo)致大范圍的地下水位急劇下降，從而引起區(qū)域性的地面沉降、海水入侵以及濕地退化等生態(tài)環(huán)境問題，嚴(yán)重影響了區(qū)域經(jīng)濟社會可持續(xù)發(fā)展。如何合理地利用地下水資源，對區(qū)域和流域地下水、地面水進行優(yōu)化配置，通過建立健全的流域水循環(huán)體系，保護地下水環(huán)境，恢復(fù)流域生態(tài)，實現(xiàn)水資源的可持續(xù)利用以及人與自然的協(xié)調(diào)，是當(dāng)今水資源研究領(lǐng)域十分重要的課題。

　　要合理地利用區(qū)域地下水資源，首要問題是要對區(qū)域地下水資源進行符合實際的模擬，摸清區(qū)域地下水的運動規(guī)律。然而，由于區(qū)域地下水模擬涉及大尺度地下水運動，而傳統(tǒng)的地下水模擬手段和方法在處理大尺度地下水運動模擬中所涉及的有限元計算中的時空尺度的選擇、水位和含水層信息不足等問題上方面存在困難，引起由于信息不足帶來的較大計算誤差。因此，針對尺度和信息不足問題，建立大區(qū)域地下水運動模擬的理論和方法十分必要。

　　本文在融合地質(zhì)統(tǒng)計、逆問題理論和地下水運動理論的基礎(chǔ)上，提出了建立大區(qū)域地下水運動模擬的理論和方法。首先，提出了大區(qū)域地下水有限元計算的時空尺度理論公式；其次，針對信息不足，提出了大區(qū)域地下水位推定的方法；第三，根據(jù)推定的地下水位，運用逆解析理論法對大區(qū)域地下水含水層透水系數(shù)進行逆推定；最后，給出了運用實例。

　　2. 大區(qū)域地下水有限元計算的時空尺度

　　在大區(qū)域地下水?dāng)?shù)值計算中，首先需要選擇研究對象所需劃分的計算網(wǎng)格大小或確定計算時間步長大小。尤其在模擬河道網(wǎng)基礎(chǔ)上的分布式水文模型和大區(qū)域地下水模型進行地表水與地下水偶合計算時，同時兼顧河道網(wǎng)的數(shù)值地形網(wǎng)格大小和地下水網(wǎng)格大小，既滿足水文計算的精度同時使地下水計算不至于失穩(wěn)就顯得至關(guān)重要。然而，由于大區(qū)域地下水計算中，由于范圍較大，人們往往希望選擇較大的網(wǎng)格進行計算，以減少計算工作量。這往往導(dǎo)致計算上的不穩(wěn)定。

　　對于地下水有限元計算的時空尺度（時空步長）的選擇，1977年，Newman等對二維地下水運動方程的有限元解法中的不穩(wěn)定問題進行了分析[1]，推測不合理的時間步長導(dǎo)致了混合型差分的出現(xiàn)，由此導(dǎo)致了解的不穩(wěn)定問題。針對Neman等的推測，Zhang（1992）[2] 和Wood（1996）[3]提出了二維地下水有限元計算的時間步長的條件。2002年，張祥偉等[4]運用最大最小原理對大尺度二維和準(zhǔn)三維地下水有限元計算的不穩(wěn)定問題進行了理論分析，提出了二維和準(zhǔn)三維地下水有限元計算的時空步長的理論公式。即對于準(zhǔn)三維地下水有限元計算的時空步長分別為：

　　空間步長條件：

　　3. 大區(qū)域地下水位的推定

　　在進行大區(qū)域地下水計算中，需要根據(jù)實測的地下水位推定初始流場，以便檢驗地下水?dāng)?shù)值計算精度、進行非恒定地下水計算以及識別含水層參數(shù)。

　　對于初始流場的推定的方法，通常有Universal Kriging方法[5-6]，UK法的計算行列比較大，計算比較復(fù)雜。Newman等（1984）[7] 和Sun（1999）[8] 運用Residual Kriging（RK）法進行區(qū)域地下水位的推定，也就是將實測地下水趨勢面去除得到正態(tài)殘差，將正態(tài)殘差運用Ordinary Kriging法進行面上殘差的推定，再加上實測地下水面的趨勢得到三角形網(wǎng)格上各點的推定地下水位值。

　　上述方法在大區(qū)域地下水計算中遇到信息不足的問題，當(dāng)實測地下水位信息不足時，推定的初始流場會帶來較大的誤差，本文根據(jù)地形水文學(xué)的原理[9]，即地下水與地形之間存在的相關(guān)性，提出運用數(shù)值地形模型（DEM）中的地形標(biāo)高作為輔助信息，修正實測地下水位得到的地下水趨勢面，然后，運用Ordinary Kriging方法對修正后趨勢面得到的正態(tài)殘差對三角形網(wǎng)格點上的殘差進行推定，每各網(wǎng)格點上推定的殘差再加上修正的趨勢值得到面上地下水位的推定值。本文將此方法定名為ROKMT（Residual Ordinary Kriging with Modified Trend）方法[10-11]，即地下水位由以下兩部分構(gòu)成：

　　大區(qū)域地下水位推定的方法如圖1。

　　4. 區(qū)域地下水透水系數(shù)的識別

　　在地下水流動模擬中，一般根據(jù)含水層的水文地質(zhì)條件和土壤特性，選定透水系數(shù)的取值。在大區(qū)域地下水計算中，含水層物理特性的信息常常十分有限，即使知道含水層的土壤類型，但由于同類土壤的透水系數(shù)的變化范圍很大，最小值與最大值之間甚至有100倍以上的變化幅度。這為正確選擇地下水參數(shù)帶來困難，本文針對式（5）的二維地下水運動，根據(jù)推定的地下水位值，運用Guass-Newton法對大區(qū)域地下水含水層透水系數(shù)進行識別。

　　二維地下水運動的基本方程為：

　　其中，h為地下水位，T為透水量系數(shù)，q為地下水涵養(yǎng)量或揚水量。

　　方程式（5）的透水量系數(shù)的識別，有直接法和間接法。直接法是根據(jù)已知的地下水位值，直接從式（5）中反求透水量系數(shù)T。直接法的問題是，由于地下水位中含有誤差，細小的誤差將導(dǎo)致水位的偏微分很大的誤差，計算穩(wěn)定性差，而且往往引起不適定問題，大大降低了參數(shù)的計算精度。本文運用間接法識別水文地質(zhì)參數(shù)。即給定含水層透水量T的初始值進行迭代計算，根據(jù)3中推定的地下水位值和反復(fù)計算得到的計算水位值的殘差平方和最小作為目標(biāo)函數(shù)，計算最優(yōu)的透水量值T。目標(biāo)函數(shù)為：

　　其中，hobs為ROKMT法的推定地下水位值。Tk+1為第K次迭代的透水量系數(shù)。通過下式計算，

　　根據(jù)含水層的物理特性，透水量系數(shù)需滿足以下限制條件：

　　當(dāng)透水量系數(shù)識別完成后，運用不規(guī)則三角形差分（TFDM）法進行大區(qū)域地下水計算。綜合3和4的計算步驟，大區(qū)域地下水計算方法的如圖2。

　　5. 運用實例-Sarobetsu濕地地下水模擬

　　本文將所提出的方法運用到日本北海道Sarobetsu濕地地下水模擬中。

　　5.1 研究對象的概況

　　Sarobetsu濕地位于日本最北端（圖3），為日本最大的以水臺癬為主的高層濕地，面積635km2。近年來，由于農(nóng)田開墾等土地開發(fā)利用，地下水位降低，濕地出現(xiàn)退化，海水上朔，生態(tài)環(huán)境惡化。最主要的表現(xiàn)是：

　　5.1.1 濕地指示性植物減少。對地下水位敏感的濕地植物水臺癬面積逐漸減少，低竹類植物風(fēng)長，面積逐年擴大，植物耗水量增加，植物葉面蒸散發(fā)增強。

　　5.1.2 湖泊面積減小。由于濕地干旱化導(dǎo)致1993年濕地內(nèi)三大湖泊兜湖、Paken湖和Peken湖的面積比1975年減少10%。三湖泊周邊的水臺癬急劇減少，生態(tài)環(huán)境惡化。

　　5.1.3 地面沉降。濕地中心區(qū)地面最大沉降達到50厘米。

　　5.1.4 海水上朔。由于地面沉降造成和Sarobetsu河的枯水期流量減少，海水沿Sarobetsu河口上朔流入Paken湖和Peken湖，破壞了湖周邊生態(tài)。

　　5.1.5 Sarobetsu河下游洪水泛濫。由于Sarobetsu濕地的地面沉降和海水的頂托，造成下游洪水泛濫。

　　因此，為恢復(fù)濕地和減少地面沉降，首要課題是研究該區(qū)域地下水位恢復(fù)的措施。前提是研究該區(qū)域地下水涵養(yǎng)機理、地下水的運動規(guī)律。

　　5.2 地形和地質(zhì)

　　Sarobetsu濕地的地形為從北到南沿海岸線逐漸降低。Sarobetsu河和海岸之間為5-20米的沙丘地帶，濕地地面高程在2米至8米之間，周邊為起伏的臺地和丘陵（如圖4）。

　　該區(qū)域地質(zhì)上主要分為兩層，即更別層以上為泥巖層，為由南向北的向斜構(gòu)造，從東、北、西向日本海方向嵌入。更別層以上為30-40米的洪積層，主要以透水性強的泥炭和沙質(zhì)層為主。更別層以下為化石地下水，與上部地下水交換很少。本研究以更別層以上低地的非承壓地下水運動為對象。

　　5.3 地下水位空間分布的推定

　　運用本文提出的ROKMT法對Sarobetsu濕地地下水位進行推定。

　　5.3.1 網(wǎng)格的劃分

　　基于250x250m的DEM地形標(biāo)高擴展成500x500m的矩形網(wǎng)格，每個矩形網(wǎng)格分為兩格三角形。全對象區(qū)域分為1051個節(jié)點，1903個三角形網(wǎng)格（如圖5）。

　　5.3.2 地下水位的推定結(jié)果

　　根據(jù)1997年72個地下水監(jiān)測點同步地下水監(jiān)測值、5個Sarobetsu河的水位值和3個湖泊的水位值，運用ROKMT法推定地下水位。另外的10個實測地下水位作為驗證點（圖6）。

　　地下水位的推定結(jié)果如圖7。

　　10個驗證點的平均誤差（ME），平均絕對誤差（MAE）和平方根二乘誤差（RMSE）分別為0.005m，0.437m和0.598m。從三種誤差看，推定的地下水位空間分布具有較好的精度，可以用來識別透水量系數(shù)。

　　5.3.3 透水系數(shù)的識別結(jié)果

　?。?）邊界條件的設(shè)定

　　海岸線一側(cè)和研究區(qū)內(nèi)部的河流、湖泊作為定水頭邊界。其他邊界位于丘陵的邊沿作為流量邊界。

　?。?）涵養(yǎng)率（入滲補給系數(shù)）和透水量系數(shù)的界限值的設(shè)定

　　運用分布式Tank水文模型，對1997年降雨徑流的模擬，以及水平衡法水量平衡計算，降水對地下水的涵養(yǎng)率取0.23。根據(jù)研究區(qū)域的水文地質(zhì)條件，確定透水量系數(shù)的界限值為0.0043m2/d≤T≥6800m2/d。

　　采用Guass-Newton法和有限元方法，使用非融雪期1997年6月地下水位資料，進行大區(qū)域地下水含水層透水系數(shù)進行識別。識別的透水量除以含水層厚度得到含水層透水系數(shù)，結(jié)果如圖8。

　　5.3.4 98年-2000年Sarobetsu濕地地下水流動的非恒定模擬

　　為恢復(fù)濕地地下水位和生態(tài)，對1998年-2000年Sarobetsu濕地地下水流動進行了非恒定模擬。

　?。?）涵養(yǎng)率的設(shè)定

　　由于濕地在每年的11月底次年的5月為降雪期，融雪對地下水有很大的補給作用，本文山本等人[12]的數(shù)字積雪-融雪模型進行了模擬。模型中考慮地形標(biāo)高、氣溫、日照量，取2℃以下時的降水為降雪，同時根據(jù)風(fēng)速修正降雪量。1997年11月—2000年6月Sarobetsu濕地積雪和融雪模擬結(jié)果如圖9。

　?。?）積雪和融雪的計算

　　涵養(yǎng)率在非積雪期取0.23，在6月至9月期間，由于植物的蒸散發(fā)強烈，涵養(yǎng)率取0.21。在融雪期，由于地下水主要靠融雪補給，因此，1-3月融雪量的90%作為地下水涵養(yǎng)，4月-5月融雪量的60%作為地下水涵養(yǎng)。儲留系數(shù)S取0.2。

　?。?）模擬結(jié)果

　　運用TFDM法對1998年1月—2000年12月Sarobetsu濕地地下水非恒定運動進行模擬，2000年的模擬結(jié)果如圖10。

　?。?）模擬結(jié)果的驗證

　　對10個驗證點的地下水位的非恒定計算結(jié)果進行了評價，10個點的ME為-0.063m，0.691m和0.731m。其中，第7點和第10點的誤差較大，分別為1.67m和2.07m，原因是7點和第10點為丘陵山區(qū)，巖石較多，計算時涵養(yǎng)率設(shè)定過高。

　　通過地下水運動的模擬，可以對濕地過去10-20年的地下水空間分布進行再現(xiàn)，結(jié)合不同年份濕地植物遙感信息和地下水空間為濕地生態(tài)恢復(fù)找到最佳地下水調(diào)控面和濕地植物退化的臨界點，為濕地恢復(fù)尋求根據(jù)。其次，通過地下水運動空間分布的模擬，可以為恢復(fù)濕地的水工程的方式、規(guī)模提供依據(jù)。

　　6. 結(jié)論

　　本文在融合地質(zhì)統(tǒng)計、逆問題理論和地下水運動理論的基礎(chǔ)上，提出了建立大區(qū)域地下水運動模擬的理論和方法，并運用到日本Sarobetsu濕地地下水模擬中，為恢復(fù)濕地提供了科學(xué)依據(jù)。此外，本文提出的方法也可以運用到其他大區(qū)域地下水的模擬和濕地保護中。