地下水模擬研究管理論文

導語：地下水模擬研究管理論文一文來源于網友上傳，不代表本站觀點，若需要原創文章可咨詢客服老師，歡迎參考。

1.引言

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

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

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

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

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

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

空間步長條件：

3.大區域地下水位的推定

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

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

上述方法在大區域地下水計算中遇到信息不足的問題，當實測地下水位信息不足時，推定的初始流場會帶來較大的誤差，本文根據地形水文學的原理[9]，即地下水與地形之間存在的相關性，提出運用數值地形模型（DEM）中的地形標高作為輔助信息，修正實測地下水位得到的地下水趨勢面，然后，運用OrdinaryKriging方法對修正后趨勢面得到的正態殘差對三角形網格點上的殘差進行推定，每各網格點上推定的殘差再加上修正的趨勢值得到面上地下水位的推定值。本文將此方法定名為ROKMT（ResidualOrdinaryKrigingwithModifiedTrend）方法[10-11]，即地下水位由以下兩部分構成：

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

4.區域地下水透水系數的識別

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

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

其中，h為地下水位，T為透水量系數，q為地下水涵養量或揚水量。

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

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

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

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

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

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

5.1研究對象的概況

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

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

5.1.2湖泊面積減小。由于濕地干旱化導致1993年濕地內三大湖泊兜湖、Paken湖和Peken湖的面積比1975年減少10%。三湖泊周邊的水臺癬急劇減少，生態環境惡化。

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

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

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

因此，為恢復濕地和減少地面沉降，首要課題是研究該區域地下水位恢復的措施。前提是研究該區域地下水涵養機理、地下水的運動規律。

5.2地形和地質

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

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

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

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

5.3.1網格的劃分

基于250x250m的DEM地形標高擴展成500x500m的矩形網格，每個矩形網格分為兩格三角形。全對象區域分為1051個節點，1903個三角形網格（如圖5）。

5.3.2地下水位的推定結果

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

地下水位的推定結果如圖7。

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

5.3.3透水系數的識別結果

（1）邊界條件的設定

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

（2）涵養率（入滲補給系數）和透水量系數的界限值的設定

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

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

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

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

（1）涵養率的設定

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

（2）積雪和融雪的計算

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

（3）模擬結果

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

（4）模擬結果的驗證

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

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

6.結論

本文在融合地質統計、逆問題理論和地下水運動理論的基礎上，提出了建立大區域地下水運動模擬的理論和方法，并運用到日本Sarobetsu濕地地下水模擬中，為恢復濕地提供了科學依據。此外，本文提出的方法也可以運用到其他大區域地下水的模擬和濕地保護中。

參考文獻

[1]Neuman,S.P.andE.A.Jacobson（1984）：Analysisofnonintrinsicspatialvariabilitybyresidualkrigingwithapplicationtoregionalroundwaterlevels,Math.Geol,16(5),pp.499-521.

[2]Zhang,H.R.(1992):Theabnormalproblemfor2Dgroundwatersimulation(inChinese),J.HydroGeo.&Eng.Geo.No.1.

[3]Wood,W.L.(1996):Anoteonhowtoavoidspuriousoscillationinthefiniteelementsolutionoftheunsaturatedflowequation,J.Hydrol.176,pp.205-218.

[4]Zhang,X.W.,K.TakeuchiandH.Ishidaira(2001):StudyontheSpuriousOscillationandStabilityinQuasiThree-DimensionalGroundwaterSimulationUsingFEMandComparisonwithSFEMandTFDM.J.JapanSoci.Hydro.WaterResour.,14(7),pp.351-363。

[5]ASCEAmericanSocietyofCivilEngineeringTaskCommitteeongeostatisticstechniquesingeohydrolgy（1990）：Reviewofgeostatisticsingeohydrology1:Basicconcepts;2:Applications,ASCEJ.Hydraul.Eng.,116(5),pp.612-658.

[6]Aboufirassi,M.,andM.A.Marino（1983）：KrigingofwaterlevelsintheSoussaquifer,Morocco,Math.Geol.,15,pp.537-551.

[7]Neuman,S.P.andE.A.Jacobson（1984）：Analysisofnonintrinsicspatialvariabilitybyresidualkrigingwithapplicationtoregionalgroundwaterlevels,Math.Geol,16(5),pp.499-521.

[8]Sun,Ne-Zheng（1999）：Inverseproblemsingroundwatermodeling,KluwerAcademicPublishers,TheNetherlands,pp.161-210.

[9]ThomasD.（1994）：Hydroheomorphology-Anintroduction-,地形（日語）、第15巻A,pp.1-4．