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

申請(qǐng)免費(fèi)試用、咨詢電話：400-8352-114

　　1. 引言

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

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

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

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

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

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

　　空間步長(zhǎng)條件：

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

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

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

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

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

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

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

　　二維地下水運(yùn)動(dòng)的基本方程為：

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

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

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

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

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

　　5. 運(yùn)用實(shí)例-Sarobetsu濕地地下水模擬

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

　　5.1 研究對(duì)象的概況

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

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

　　5.1.2 湖泊面積減小。由于濕地干旱化導(dǎo)致1993年濕地內(nèi)三大湖泊兜湖、Paken湖和Peken湖的面積比1975年減少10%。三湖泊周邊的水臺(tái)癬急劇減少，生態(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)動(dòng)規(guī)律。

　　5.2 地形和地質(zhì)

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

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

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

　　運(yùn)用本文提出的ROKMT法對(duì)Sarobetsu濕地地下水位進(jìn)行推定。

　　5.3.1 網(wǎng)格的劃分

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

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

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

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

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

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

　　（1）邊界條件的設(shè)定

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

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

　　運(yùn)用分布式Tank水文模型，對(duì)1997年降雨徑流的模擬，以及水平衡法水量平衡計(jì)算，降水對(duì)地下水的涵養(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í)別。識(shí)別的透水量除以含水層厚度得到含水層透水系數(shù)，結(jié)果如圖8。

　　5.3.4 98年-2000年Sarobetsu濕地地下水流動(dòng)的非恒定模擬

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

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

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

　?。?）積雪和融雪的計(jì)算

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

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

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

　　（4）模擬結(jié)果的驗(yàn)證

　　對(duì)10個(gè)驗(yàn)證點(diǎn)的地下水位的非恒定計(jì)算結(jié)果進(jìn)行了評(píng)價(jià)，10個(gè)點(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ū)，巖石較多，計(jì)算時(shí)涵養(yǎng)率設(shè)定過(guò)高。

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

　　6. 結(jié)論

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