基于Modflow的虎林灌區(qū)地下水開采三維數(shù)值模擬

韓　雪，　馬大男

（黑龍江科技大學(xué) 建筑工程學(xué)院，哈爾濱 150022）

基于Modflow的虎林灌區(qū)地下水開采三維數(shù)值模擬

韓雪，馬大男

（黑龍江科技大學(xué) 建筑工程學(xué)院，哈爾濱 150022）

針對地下水開采引起的地面沉降問題，以雞西市虎林灌區(qū)為例，采用國際通用的地下水三維滲流與地面沉降耦合計算軟件——Processing Modflow，建立虎林灌區(qū)地下水-地面沉降耦合模型。以2014年年均開采量為依據(jù)，對2015年的地下水位及沉降進(jìn)行數(shù)值模擬計算，并通過鉆孔及室內(nèi)土工實驗分析其沉降機理。結(jié)果表明：2015年水位呈現(xiàn)北部及東南部水位高、西部較低，沉降量由中心區(qū)向四周逐漸增大的趨勢；沉降嚴(yán)重區(qū)域土樣的天然含水率、密度隨深度增加呈現(xiàn)減小的趨勢，相對密度隨深度變化不明顯；土樣的塑性指數(shù)與其相對密度呈現(xiàn)正相關(guān)，變形量與其壓力呈正相關(guān)，土體的壓縮系數(shù)和滲透系數(shù)與壓力呈負(fù)相關(guān)。

地下水開采；室內(nèi)土工實驗；沉降機理；Processing Modflow；耦合模型；相關(guān)性

虎林灌區(qū)是三江平原的糧食主產(chǎn)區(qū)，為保證黑龍江省糧食安全做出了顯著貢獻(xiàn)。長期以來，由于灌區(qū)地表水灌溉設(shè)施不完善、水利工程老化情況較為嚴(yán)重，灌區(qū)出現(xiàn)地下水超采現(xiàn)象，導(dǎo)致地面沉降，對區(qū)域農(nóng)業(yè)生產(chǎn)和地區(qū)經(jīng)濟(jì)的可持續(xù)發(fā)展造成了不利影響。因此，建立地下水-地面沉降水土耦合模型，研究地下水超采引起的地面沉降問題及其沉降機理，為有關(guān)部門的科學(xué)決策提供依據(jù)，顯得十分必要［1-4］。

為此，筆者以虎林灌區(qū)為研究對象，應(yīng)用軟件Processing Modflow建立地下水-地面沉降耦合模型，并對沉降嚴(yán)重區(qū)域利用現(xiàn)場鉆孔及室內(nèi)土工實驗進(jìn)行沉降機理分析。

1　水流模型與土力學(xué)模型

1.1水流模型

對于非均質(zhì)、空間三維非穩(wěn)定流系統(tǒng)，可用地下水流的三維運動偏微分方程來描述［5］：

式中：kxx、kyy、kzz——含水層坐標(biāo)軸方向的滲透系數(shù)，md-1；

h——水頭值，L；

W——匯源項，d-1；

Ss——儲水率，m-1；

t——時間，d；

Ω——立體計算域。

1.2土力學(xué)模型

假定地層總應(yīng)力不發(fā)生變化，且土層以豎向變形為主，水平方向的應(yīng)變變化微小，可忽略不計。則土體有效應(yīng)力變化值，可由式（2）求得［5］:

式中：ρw——水的密度，g·cm-3；

g——重力系數(shù)，N·kg-1；

Δh——水頭，L。

土體有效應(yīng)力的增加或減小，使得黏性土層在垂直方向上的壓縮量或回彈量呈線性增長。土層彈性變形量可通過式（3）計算:

Δb——壓縮量，m；

Δt——土層壓縮經(jīng)歷時間，s；

dA——土層厚度，m；

Ss'k——彈性骨架釋水系數(shù)，壓縮為正，回彈為負(fù)。

將式（3）中的qi代入到方程（1）中，給定初始條件和邊界條件就可以得到地下水-地面沉降耦合模型，其方程為

初始條件：

邊界條件：

式中：n——外法線單位向量；

Г0——水流區(qū)域上的第二類邊界；

q——邊界上的流量。

2　數(shù)值模擬

2.1模擬區(qū)邊界特性及地層巖性

灌區(qū)的邊界范圍：北以阿布沁河為界，南至穆棱河，西邊與慶豐農(nóng)場接壤，東至烏蘇里江沿岸的珍寶島濕地自然保護(hù)區(qū)。模擬區(qū)地層巖性見圖1。

圖1　模擬區(qū)地層巖性Fig.1 Simulation area’s formation lithology

灌區(qū)地層共為五層，第一層土層為0～5 m灰黃色的亞黏土，土層厚度為5.0 m；第二層土層為5.0～7.6 m灰綠黃色的淤泥質(zhì)黏土，厚度為2.6 m，為不含水土層；第三層土層為7.6～17.4 m灰黃色的粉質(zhì)黏土，粉質(zhì)黏土夾亞黏土，土層厚度為9.8 m，土層含水飽和；第四層土層為17.4～40.2 m灰白及灰色的粉質(zhì)黏土，土層厚度為22.8 m，土層含水飽和；第五層土層為40.2～51.5 m灰綠灰黑色的淤泥質(zhì)黏土，土層致密，不含水。

2.2模型邊界概化與網(wǎng)格劃分

采用Processing Modflow進(jìn)行孔隙介質(zhì)中地下水運動的數(shù)值模擬［6-7］，首先建立離散化的三維模型（網(wǎng)格剖分、層數(shù)等）［7］。根據(jù)水文地質(zhì)特征和模擬計算需求，將地層概化為三層，第一層為淺部弱潛水層，位于承壓含水層頂部，0～7.6 m，為亞黏土和淤泥質(zhì)黏土層，厚度為7.6 m，不含水，其滲透性差，概化為隔水邊界；第二層為承壓含水層，7.6～40.2 m，為粉質(zhì)黏土層，厚度為32.6 m，含水飽和；第三層為不透水層，位于承壓含水層底部，40.2～51.5 m，為淤泥質(zhì)黏土層，厚度為11.3 m，致密不含水，其滲透性差，概化為隔水邊界。

由于灌區(qū)北部、南部是阿布沁河和穆棱河，西鄰慶豐農(nóng)場，東至烏蘇里江，因此，北部、南部均為虎林灌區(qū)的分水嶺，可以概化為定水頭邊界；西部為慶豐農(nóng)場平原區(qū)，淺層地下水不與外界進(jìn)行水流交換，可概化為隔水邊界；東部以烏蘇里江為界，常年地下水位變化不大，可視為定水頭邊界；中部七虎林河河流穿越，系統(tǒng)分析表明，七虎林河為模擬區(qū)的優(yōu)勢河流，可視為河流邊界?；⒘止鄥^(qū)概化圖見圖2。

圖2　虎林灌區(qū)邊界概化Fig.2 Hulin irrigation district’s generalized boundary

灌區(qū)共14口井，自上而下依次編號為1～14號。在模擬范圍為479 km2平面內(nèi)，將模擬區(qū)剖分為不同尺寸大小的矩形網(wǎng)格，一般區(qū)域的單元網(wǎng)格劃分成500 m×500 m，井群重點區(qū)域網(wǎng)格加密，單元網(wǎng)格劃分為85 m×85 m。垂直方向上剖分三層，加密后網(wǎng)格共分154行和156列，共剖分24 024個單元，模擬區(qū)網(wǎng)絡(luò)剖分見圖3。

2.3模型識別與驗證

2.3.1模型基本參數(shù)設(shè)定

三維數(shù)值模型采用分區(qū)賦值方法，其主要參數(shù)包括水平向滲透系數(shù)、垂直滲透系數(shù)、貯水系數(shù)、有效孔隙度、降水入滲系數(shù)及灌溉入滲系數(shù)［8］。

圖3　模擬區(qū)網(wǎng)格剖分Fig.3 Simulation area mesh

水平方向上，模擬區(qū)潛水含水層滲透系數(shù)和給水度可以分為兩個區(qū)，七虎林河以北為Ⅰ區(qū)，七虎林河以南為Ⅱ區(qū)（見圖4）。各分區(qū)初值參數(shù)如表1所示。

表1　虎林灌區(qū)含水層參數(shù)初值Table 1　Initial parameters of aquifer parameters in Hulin irrigation district

圖4　虎林灌區(qū)含水層參數(shù)分區(qū)Fig.4　Parameter zoning map of aquifer in Hulin irrigation district

垂直方向上，淺部弱潛水層為亞黏土，滲透系數(shù)（kxx、kyy）為6 m/d，kzz為0.6 m/d，給水度0.06，彈性釋水系數(shù)為0.000 27；中層承壓水層為粉質(zhì)黏土，滲透系數(shù)（kxx、kyy）為20.76 m/d，kzz為2.1 m/d，給水度為0.21，彈性釋水系數(shù)為0.000 40；底部不透水層為淤泥質(zhì)黏土，滲透系數(shù)（kxx、kyy）為5 m/d，kzz為0.5 m/d，給水度0.01，彈性釋水系數(shù)為0.000 21。

2.3.2源匯項處理

灌區(qū)地下水補給來源，主要接受大氣降水入滲補給、農(nóng)田灌溉回滲補給及阿布沁河、穆棱河、七虎林河和烏蘇里江滲漏補給。工農(nóng)業(yè)及生活用水開采、淺層含水層的蒸發(fā)，越流排泄及植物蒸騰排泄是本區(qū)地下水的主要排泄方式。其中，入滲和蒸發(fā)利用recharge程序包和蒸發(fā)包處理，人工開采在模型中利用程序包well處理。

2.3.3模型驗證

模擬的時間范圍為2014年1月1日運行到2015年1月1日，共365個時段，時間歩長為365 d。將初始水位和其他要素輸入模型后，進(jìn)行水位輸出，對比模擬流場和初始流場，反復(fù)調(diào)試參數(shù)，進(jìn)行模型識別和校正。初始流場與實際末流場擬合圖如圖5所示。總體來看，模型的總體擬合程度較好。

圖5　初始流場與實際末流場擬合Fig.5 Fitting of initial and actual flow field

2.4結(jié)果分析

基于Processing Modflow建立虎林灌區(qū)2015年全年的地下水流場與地面沉降耦合三維模型，如圖6～8所示。

圖6　含水層水位線等值線云圖（單位：m）Fig.6 Contour lines of water level line（unit:m）

從圖6、7可以看出，含水層在靠近烏蘇里江，穆棱河和阿布沁河的水位變化比較明顯，水位變化量為20 cm，中心區(qū)水位變化較初始水位有明顯回升趨勢。一方面，研究區(qū)北部、東部以及南部受烏蘇里江和阿布沁河的入滲補給明顯；另一方面含水層受到降水等面狀補給，中部水位出現(xiàn)回升現(xiàn)象。平原區(qū)地下水位變化由北向南逐漸變淺，地下水位變化明顯。

圖7　365 d含水層水位線三維立體圖（單位：m）Fig.7 365 d three-dimensional of water level（unit:m）

圖8　365 d含水層沉降量等值線平面云圖（單位：m）Fig.8 365 d equivalent line of water bearing stratum（unit:m）

從圖8中可以看出，由于含水層補給、排泄方式不同，該層各區(qū)域的沉降量也不同，模擬區(qū)沉降量呈現(xiàn)由中心區(qū)向周邊江河逐步增大的總體趨勢，最大沉降值2.5 cm。

3　模擬區(qū)沉降機理分析

3.1綜合鉆孔的選取

研究模擬區(qū)淺層地下水開采，分別取不同深度含水層的原狀土進(jìn)行研究。利用DC-1水文地質(zhì)雙層抽水鉆孔在沉降嚴(yán)重區(qū)鉆孔取樣。該孔位于黑龍江省農(nóng)墾總局牡丹江分局八五四農(nóng)場的大王家村西北，揭穿整個模擬區(qū)的淺層含水層，具有典型性和代表性。文中通過研究該孔土樣的物理性指標(biāo)分析虎林灌區(qū)地面沉降區(qū)的沉降機理。

3.2黏性土物理特性

此次DC-1孔內(nèi)取樣深度范圍0～52 m，其中，5.5～5.8 m，37.5～37.8 m，49.1～49.4 m深度的土樣編號分別為原001、原005、原009，實驗項目均為土的基本物理參數(shù)測試、排水固結(jié)及滲透實驗。為充分了解土體的物理性質(zhì)，測定抽水鉆孔原狀土樣密度、天然含水率、相對密度、液塑限四項指標(biāo)，并總結(jié)實驗所得的土體物理參數(shù)和壓縮系數(shù)之間存在的規(guī)律。

3.2.1基本物理參數(shù)測試結(jié)果

土樣基本物理參數(shù)與壓縮系數(shù)隨深度變化曲線如圖9所示。從圖9可以看出，測試的原狀土土樣的天然含水率在24.86%～25.71%；密度值是在1.73～1.95 g/cm3；相對密度在2.71～2.79。天然含水率值隨深度增加呈現(xiàn)減小的趨勢，同時呈現(xiàn)上下波動的變化特點，埋深在27.5 m以上的土樣天然含水率隨深度的增加而逐漸減小，埋深27.5～41.5 m之間土樣的天然含水率有逐漸增大的趨勢：密度隨深度的增加也呈現(xiàn)減小的趨勢，其波動并不明顯，在埋深11.3 m處，密度突然變大，其值為1.95 g/cm3；其他深度值趨于平穩(wěn)。相對密度隨埋深增加的變化不明顯，總體呈現(xiàn)上下波動的趨勢，其值總體上在2.71～2.79之間變化。當(dāng)然，隨深度變化，土樣的土體結(jié)構(gòu)和所處應(yīng)力環(huán)境對土體基本物理參數(shù)產(chǎn)生重要的影響。

圖9　土樣基本物理參數(shù)與壓縮系數(shù)隨深度變化曲線Fig.9 Soil sample’s basic physical parameters and compression coefficient change with depth

3.2.2液塑限實驗結(jié)果

由實驗結(jié)果可知，實驗土樣的液限（ωL）在29%左右，塑限（ωp）在20%左右，其值均隨著深度的增加而逐漸增大，從而說明土體中黏粒的含量在逐步增加。從圖9中可以看出，土樣的塑性指數(shù)與其相對密度呈現(xiàn)正相關(guān)性。相對密度的大小與其礦物組成成分有關(guān)，隨著土樣黏粒含量的增加而增大；土樣的塑性指數(shù)和相對密度又與土體的壓縮系數(shù)成一定的負(fù)相關(guān)性，也就是說，塑性指數(shù)越大，相對密度越大，土體越不容易壓縮。

3.3土樣壓縮變形規(guī)律

由實驗可知，土體的初始孔隙比（e0）隨著埋深的增加而呈現(xiàn)逐漸減小的趨勢。壓縮系數(shù)在深度范圍11.3～21.6 m，27.0～36.0 m兩處有峰值，在22.0～26.0 m，36.0～46.0 m處有兩個谷值，其他深處的值變化不大，這與土樣所處的地質(zhì)環(huán)境和土體結(jié)構(gòu)密切相關(guān)。

3.3.1變形量隨壓力的變化特征

不同深度土樣的變形量隨壓力變化曲線如圖10所示。由圖10可以看出，土樣的變形量隨壓力的增大而增大，在壓力100～400 kPa荷載作用下，土樣的變形量變化較為劇烈，當(dāng)壓力達(dá)到1 000 kPa時，土樣的變形量趨于穩(wěn)定，尤其對于原009土樣，這種變化趨勢更為明顯。綜合分析，相同結(jié)構(gòu)的土體，其壓縮變形量隨土體埋深的增加而呈現(xiàn)逐漸減小的趨勢。

圖10　不同深度土樣變形量隨壓力變化曲線Fig.10 Varying depth of soil sample deformation with pressure

3.3.2土體壓縮系數(shù)隨壓力變化特征

不同深度土樣的壓縮系數(shù)隨壓力變化曲線如圖11所示。由圖11可以看出，土樣的壓縮系數(shù)隨著壓力的增加而呈現(xiàn)逐漸減小的趨勢。原001和原005土樣在壓力范圍100～250 kPa之間，波動較大，出現(xiàn)谷峰值。三種土樣在前期壓應(yīng)力較小的情況下，壓縮系數(shù)變化劇烈，其增大值和減小值變化幅度非常大。埋深在49.1～49.4 m的原009土樣為低壓縮性土，埋深在37.5～37.8 m的原005土樣為中壓縮性土，埋深在5.5～5.8 m的原001土樣為中壓縮性土。

3.3.3土體滲透系數(shù)隨埋深變化特征

不同深度土樣的滲透系數(shù)變化曲線如圖12所示。由圖12可知，0～50 m土樣的滲透系數(shù)數(shù)量級為101～10-2，且隨著土樣埋深的增大，滲透系數(shù)逐漸減小。

圖11　不同深度土樣的壓縮系數(shù)隨壓力變化曲線Fig.11 Compression coefficient of different depth soil samples with pressure change curve

圖12　不同深度土樣的滲透系變化曲線Fig.12 Change curve of soil samples with different depth

綜合上述分析，在研究區(qū)內(nèi)伴隨著地下水的開采，地下水位降低，土體孔隙水壓力減小，而有效應(yīng)力增大，土體的變形量與其壓力呈正相關(guān)性，土體的壓縮系數(shù)和滲透系數(shù)與壓力呈負(fù)相關(guān)性。這些結(jié)果顯示，由于土樣中黏性土和淤泥質(zhì)黏土占的比例比較大，所以容易形成永久性地面沉降。該研究結(jié)果對揭示研究區(qū)地面沉降機理有重要意義。

4　結(jié)　論

（1）基于Processing Modflow軟件建立三維地下水-地面沉降耦合模型，模擬結(jié)果顯示：2015年水位呈現(xiàn)北部及東南部水位高、西部較低、沉降量由中心區(qū)向四周逐漸增大的趨勢。

（2）沉降嚴(yán)重區(qū)的現(xiàn)場鉆孔與室內(nèi)實驗結(jié)果表明，沉降重災(zāi)區(qū)綜合鉆孔土樣的天然含水率、密度隨深度增加呈現(xiàn)減小的趨勢 ，相對密度隨深度變化不明顯。

（3）土樣的壓縮變形規(guī)律為：土樣的塑性指數(shù)與其相對密度呈現(xiàn)正相關(guān)，變形量與其壓力呈正相關(guān)，土體的壓縮系數(shù)和滲透系數(shù)與壓力呈負(fù)相關(guān)。這一規(guī)律對揭示研究區(qū)地面沉降機理有重要意義。

［1］吳昌友，李中才，王福林.三江平原撓河流域地下水流數(shù)學(xué)模型的建立及仿真［J］.水動力學(xué)研究與進(jìn)展，2011，26（3）:384-389.

［2］LIN LIN，YANG JIN ZHONG，ZHANG BIN.A simplified numerical model of 3-d groundwater and solute transport at large scale area［J］.Journal of Hydrodynamics，2010，22（3）:319-328.

［3］羅美芳，陳遠(yuǎn)法，杜麗坷，等.浙江溫州市永強平原地面沉降分析與趨勢預(yù)測［J］.中國地質(zhì)災(zāi)害與防治學(xué)報，2012，23（3）:51-56.

［4］孫香泰，艾曉燕，張力春.三江平原地下水生態(tài)水位初步研究［J］.黑龍江水利技，2012，40（8）:166-168.

［5］付延玲，郭正法.Processing Modflow在地下水滲流與地面沉降研究中的應(yīng)用［J］.勘察科學(xué)技術(shù)，2006，24（4）:19-23.

［6］陶月贊，姚 梅.地下水滲流力學(xué)的發(fā)展進(jìn)程與動向［J］.吉林大學(xué)學(xué)報:地球科學(xué)版，2007，37（2）:221-230.

［7］CHIANG WENHSING，WOLFGANG KINZELBACH.3D-Groundwatermodeling with PMWIN-A Simulation S-ystem f or Modeling Ground water Flow and Pollution［M］.［S.l.］：Springer，2001.

［8］路瑞利，方樹星，王紅雨.基于Modflow的某水源區(qū)地下水開采三維數(shù)值模擬［J］.武漢大學(xué)學(xué)報：工學(xué)版，2011，44（5）:618-623.

［9］孫袖正.地下水流的數(shù)學(xué)模型和數(shù)值方法［M］.北京:地質(zhì)出版社，1981.

（編輯荀海鑫）

MODFLOW -based 3D numerical simulation of groundwater exploitation in Hulin irrigation area

HAN Xue，MA Da’nan
（School of Civil Engineering，Heilongjiang University of Science＆Technology，Harbin 150022，China）

This paper is focused on addressing ground settlement resulting from groundwater exploitation.The research exemplified by a case study of irrigation areas in Hulin in Jixi city involves using Processing Modflow-an internationally accepted coupling calculation software for 3D groundwater seepage and surface subsidence and thereby developing the coupling-model between groundwater and ground subsidence typical of Hulin irrigation area.The study describes a simulation calculation of the groundwater level and subsidence in 2015，according to the annual mean production in 2014 and an analysis of the mechanism behind the settlement by drilling and laboratory soil test.Results show that，in 2015 there is a higher water level in the north and southeast，but a lower water level in the west，with the settlement tending to have a gradual increase from the central area to surrounding areas；the natural moisture content and density tend to decrease due to an increasing depth in the soil samples taken from the area with severe subsidence，but there is no obvious change in relative density；a positive correlation exists between the plastic index and relative density of the soil sample，and between the deformation amount and the pressure，but a negative correlation occurs between the compression coefficient and the osmotic coefficient of the soil on the one hand and the pressure on the other hand.

groundwater mining；laboratory soiltestin；mechanism of settlement；Processing Modflow；coupling-model；relevance

10.3969/j.issn.2095-7262.2015.06.016

P641.8

2095-7262（2015）06-0654-06

A

2015－10－10

黑龍江省自然科學(xué)基金項目（E201462）

韓雪（1969－），男，吉林省榆樹人，教授，博士，研究方向:邊坡穩(wěn)定技術(shù)及其理論，E-mail:898531033@qq.com。