王?？频?/p>

摘 要：根據(jù)德州市區(qū)水文地質(zhì)條件，將含水層概化為非均質(zhì)、有越流向的二維承壓水流模型，根據(jù)水量均衡原理建立三角網(wǎng)格剖分節(jié)點的差分方程，采用逐次超松弛迭代方法解系數(shù)矩陣；根據(jù)實測資料對模型進行識別與調(diào)試，采用識別后的有關(guān)數(shù)據(jù)，進行地下水量均衡計算和資源量計算，發(fā)現(xiàn)補給量小于排泄量，年超采量為393萬m3/年，對未來地下水動態(tài)進行了預(yù)測，不同條件下地下水位降深達5.6～61.3m。研究發(fā)現(xiàn)隨著超采時間的延長及超采量的加大，地下水位將持續(xù)下降，地下水漏斗面積不斷擴大，導(dǎo)致地面沉降、水質(zhì)惡化等環(huán)境地質(zhì)問題進一步惡化。為防止地下水超采引起的環(huán)境地質(zhì)問題進一步惡化，提出了科學(xué)的防治措施及建議。

關(guān)鍵詞：地下水運動；數(shù)值模擬；地下水漏斗；地面沉降；環(huán)境地質(zhì)問題

1 數(shù)學(xué)模型

近年來，德州市地下水開采量大于補給量，地下水位處于持續(xù)下降狀態(tài)。根據(jù)深層水水資源概念，不僅將補給資源全部開采了出來，還利用了一部分儲存資源。本文在充分研究德州市水文地質(zhì)條件的基礎(chǔ)上進行了模型概化與計算。

1.1 數(shù)學(xué)模型的建立

根據(jù)研究區(qū)水文地質(zhì)條件，將主要開采含水層分為兩個含水巖組，上部稱為第一含水層（埋深300～500米以上），下部稱為第二含水層（埋深500～800米），中間為一越流層（見圖1）。第一、第二含水巖組均可概化為非均質(zhì)、有越流項的二維承壓水流模型，數(shù)學(xué)模型如下：

式中：H1、H2為第一、第二含水層水位標(biāo)高（m）；H10、H20為第一、第二含水層初始水位標(biāo)高（m）；T1、T2為第一、第二含水層導(dǎo)水系數(shù)（m2/d）；μ1*、μ2*為第一、第二含水層貯水系數(shù)（無量綱）；K 為弱透水層滲透系數(shù)（m/d）；M 為弱透水層滲透系數(shù)（m/d）；B= K/ M為越流系數(shù)（無量綱）；Γ2為滲流區(qū)域上的第二類邊界；q1、q2為第一、第二含水層為第二類邊界上的單寬流量（m3/d·m）；n1、n2為第一、第二含水層為第二類邊界某點的外法線方向；Q1、Q2為第一、第二含水層匯源項（m/d）；x、y為平面坐標(biāo)（m）；t為時間（d）；D為滲流區(qū)域（m2）。

1.2 數(shù)值解法

數(shù)學(xué)模型第一式的物理意義為：在區(qū)域面積和時刻均趨于無限小時含水層中地下水的均衡方程?；痉匠套髠?cè)的第一、二項表示側(cè)向流入和流出含水層流量的差值（單位時間單位面積含水層而言）；第三項為越流項；第四項為源匯項，表示單位面積含水層上垂向補給或排泄流量；右側(cè)表示單位時間、單位面積含水層貯水量的變化率。若面積和時段不取無限小而采用有限值，則可以用近似方法逼近方程各項的極限值，根據(jù)水量均衡原理可建立三角網(wǎng)格剖分節(jié)點的差分方程。節(jié)點i的差分方程如下：

第一含水巖組：

第二含水巖組：

式中：Cij為側(cè)向水量的水位系數(shù)；H1、H2為第一、第二含水層水位標(biāo)高；B 為越流系數(shù)；Q1（i）、Q2（i）為第一、第二含水層垂向水量；Fi為貯存量的水位系數(shù)；j為節(jié)點i周圍的節(jié)點號（i，j=1， …n，節(jié)點個數(shù)）；t為時間步長；k為計算時刻（k =1， …m，時段數(shù)）。

若用X表示未知數(shù)的列向量： ；用A表示方程組的系數(shù)矩陣，b表示常數(shù)項列陣，則上述數(shù)值方程組可寫為矩陣形式：AX=b。

用水量均衡原理推導(dǎo)出的不規(guī)則網(wǎng)格有限差分方程，其系數(shù)矩陣A是一個具有對角線優(yōu)勢的高度稀疏的對稱正定矩陣，對于大型稀疏矩陣（A的階數(shù)n較大，但零元素較多），利用迭代法求解是較合適的。由逐次超松弛迭代方法SOR（Successive Over Relaxation Method）收斂的判別定理知：如果A為對稱正定矩陣且0<ω<2，則解AX=b的SOR方法收斂。因此，本文選擇SOR方法，它是高斯—塞德爾迭代方法的一種加速方法，是解大型稀疏矩陣方程組的有效方法之一，它具有計算公式簡單，程序設(shè)計容易，占用計算機內(nèi)存少等優(yōu)點。下面對其迭代公式作一簡述。

現(xiàn)狀開采量條件下，隨開采時間的增加，地下水水位降深愈來愈大，t=10年時，第一含水層水位降深s為15.06～33.66m，第二含水層水位降深s為22.38～31.32m；第二含水層開采量增加20%條件下，相同時間與節(jié)點的降深大于現(xiàn)狀開采量下的降深，t=10年時，第一含水層水位降深s為17.75～36.62m，第二含水層水位降深s為23.12～36.24m。

3結(jié)論

①根據(jù)德州市區(qū)地下水運動數(shù)值模擬，發(fā)現(xiàn)現(xiàn)狀地下水年超采量為393萬m3/年，連年超采已引起區(qū)域地下水位降落漏斗，并已引起地面沉降、水質(zhì)惡化等一系列的環(huán)境地質(zhì)問題，如1991年至2006（15年）德州市累計地面沉降總量為710.5mm，年均沉降量為47.4mm/a，沉降中心累計沉降量為992mm，沉降范圍已與周邊的河北省衡水、滄州地面沉降連為一片。

②為防止地下水超采引起的環(huán)境地質(zhì)條件進一步惡化，建議總體規(guī)劃、合理利用水資源，施行最嚴格的水資源管理制度，建立地下水環(huán)境監(jiān)測網(wǎng)絡(luò)，建立有效的節(jié)水措施及地下水保護措施，實現(xiàn)水資源的可持續(xù)利用。

參考文獻

[1] 陳崇希，唐仲華.地下水流動問題數(shù)值方法[M].中國地質(zhì)大學(xué)出版社.2002.

[2] 孫訥正.地下水流的數(shù)學(xué)模型和數(shù)值方法[M].地質(zhì)出版社.1981.

[3] 朱學(xué)愚，錢孝星，劉新仁.地下水資源評價[M].南京大學(xué)出版社.1987.

[4] 李慶揚，王能超，易大義.數(shù)值分析[M].華中科技大學(xué)出版社.2005.

[5] 薛禹群.地下水動力學(xué)[M].地質(zhì)出版社.1986.

[6] 張蔚榛.地下水非穩(wěn)定流計算和地下水資源評價[M].科學(xué)出版社.1983.

[7] 德州市水文地質(zhì)調(diào)查研究報告[R].德州市水文地質(zhì)勘察院.2006.

[8] 陸金甫，關(guān)治.偏微分方程數(shù)值解法[M].清華大學(xué)出版社.1987.

[9] 王大純，張權(quán)，史毅虹等.水文地質(zhì)學(xué)基礎(chǔ)[M].地質(zhì)出版社.1994.

[10] 蔡文曉.德州市深層地下水開采與地面沉降關(guān)系研究[D].吉林大學(xué).2009.

[11] 張纓，周家權(quán)，任緒偉.德州市地下水漏斗引發(fā)生態(tài)問題及應(yīng)對措施[C].

[12] 趙全升，馮娟，安樂生.德州市深層地下水水質(zhì)演化研究[J].地理科學(xué)，2009，29（5）.

[13] 馮娟.德州市深層地下水動態(tài)變化與模擬研究[D].青島大學(xué).2008.

[14] 盧文喜.地下水運動數(shù)值模擬過程中邊界條件問題探討[J].水利學(xué)報，2003，（3）.

[15] 王元行，王吉良，王華敏.華北地區(qū)滄縣出現(xiàn)的地下水環(huán)境問題與對策[J].水資源保護，2003，（4）.