高翠萍，杜新強，楊悅鎖，賈思齊

(1.沈陽大學(xué)區(qū)域污染環(huán)境生態(tài)修復(fù)教育部重點實驗室，遼寧沈陽 110044；2.吉林大學(xué)地下水資源與環(huán)境教育部重點實驗室，吉林長春 130021)

?

基于參數(shù)靈敏度分析的地下水?dāng)?shù)值模擬研究
——以黑龍江樺南縣為例

高翠萍1，杜新強2，楊悅鎖1，賈思齊1

(1.沈陽大學(xué)區(qū)域污染環(huán)境生態(tài)修復(fù)教育部重點實驗室，遼寧沈陽 110044；2.吉林大學(xué)地下水資源與環(huán)境教育部重點實驗室，吉林長春 130021)

采用正交試驗法對地下水流數(shù)值模擬中的參數(shù)靈敏度進(jìn)行分析，以黑龍江省樺南縣為例，根據(jù)研究區(qū)實際的水文地質(zhì)條件選取了潛水的降水入滲補給系數(shù)、滲透系數(shù)、給水度、弱透水層的垂向滲透系數(shù)以及承壓水的滲透系數(shù)和彈性釋水系數(shù)6個參數(shù)進(jìn)行研究，每個參數(shù)選取5個水平值。結(jié)果表明：潛水含水層的滲透系數(shù)對模型計算水位的影響最大，因此在利用Visual MODFlow做水資源評價前查明研究區(qū)的水文地質(zhì)參數(shù)特別是潛水的滲透系數(shù)對于正確建立模型十分必要。

靈敏度分析；正交試驗法；數(shù)值模擬；地下水

地下水流數(shù)值模擬已經(jīng)成為地下水資源評價的重要方法，然而在建立數(shù)值模擬模型的過程中要對水文地質(zhì)條件進(jìn)行概化。概化使得模型參數(shù)具有不確定性，最終導(dǎo)致模型計算結(jié)果的不確定性。靈敏度分析能夠評價模型中各個參數(shù)的不確定性對模型運行結(jié)果的影響，因此，在做地下水流數(shù)值模擬的過程中，對模型參數(shù)進(jìn)行靈敏度分析非常必要[1]。參數(shù)靈敏度分析是基于數(shù)理統(tǒng)計原理分析的不確定性方法[2]，是地下水流數(shù)值模擬中不可或缺的一個步驟，目的是確定參數(shù)的不確定性對模型的影響程度[3]。在參數(shù)識別的過程中，若某個參數(shù)的微小變化能引起計算結(jié)果較大程度的改變，則表明該參數(shù)的靈敏度較高；反之，則較低。靈敏度分析方法種類很多，如多元回歸法、正交試驗法、區(qū)域靈敏度分析法、普適似然度不確定性估計法等[4]。筆者以黑龍江省樺南縣為例，利用Visual MODFlow軟件對地下水流進(jìn)行模擬，將正交試驗法[5-6]引入?yún)?shù)靈敏度分析中，并選取了潛水的降水入滲補給系數(shù)、滲透系數(shù)、給水度、弱透水層的垂向滲透系數(shù)以及承壓水的滲透系數(shù)和彈性釋水系數(shù)6個參數(shù)進(jìn)行靈敏度分析。

1　研究區(qū)概況

黑龍江省佳木斯市樺南縣城區(qū)(圖1)，東西長約12 km，南北寬約8 km，面積96 km2。多年平均氣溫3.1 ℃，多年平均降水量552 mm，多年平均水面蒸發(fā)量(E601)約1 340 mm。根據(jù)地下水的補給、徑流、排泄條件及水力聯(lián)系等特征將含水層劃分為潛水含水層和承壓含水層[7]。

2　模型的建立

2.1水文地質(zhì)概念模型的建立研究區(qū)潛水的天然水力坡度為1.1×10-3～1.5×10-2，承壓水的天然水力坡度為1.3×10-3～7.4×10-3，流場較為平緩，滲流基本符合達(dá)西定律，潛水和承壓水之間通過越流進(jìn)行水量交換，且水流各要素隨時間變化，因此整個滲流區(qū)域地下水的運動特征為三維非穩(wěn)定流；水文地質(zhì)參數(shù)隨空間變化，概化為非均質(zhì)各向同性。綜上所述，該研究將水文地質(zhì)概念模型概化為非均質(zhì)各向同性的三維非穩(wěn)定流模型。潛水含水層北部邊界為河流，因此，將北部邊界概化為水頭邊界，其他邊界根據(jù)地下水流場概化為流量邊界；承壓含水層所有的邊界均根據(jù)地下水流場概化為流量邊界，弱透水層側(cè)向邊界概化為隔水邊界。模型頂部為潛水面邊界，主要接受大氣降水補給，同時以垂直蒸發(fā)的形式排泄；模型底部有相對隔水的泥巖，概化為不透水邊界。

2.2數(shù)學(xué)模型的建立根據(jù)水文地質(zhì)概念模型，建立如下地下水流數(shù)學(xué)模型：

圖1　研究區(qū)區(qū)位示意Fig.1　The location of study area

2.3數(shù)值模型的建立及求解該模型應(yīng)用了Visual MODFlow軟件，方法是在計算區(qū)域內(nèi)采用矩形剖分和線性插值，應(yīng)用伽遼金有限差分法進(jìn)行求解。網(wǎng)格長120 m，寬80 m，網(wǎng)格剖分情況見圖2。該水位統(tǒng)測時間是2009年11月11日至2010年10月31日，選取2009年11月11日至2010年5月6日為識別期，共分為36個時段，2010年5月6日至10月31日為驗證期，共分為36個時段。經(jīng)過模型的識別和驗證，最后確定水文地質(zhì)參數(shù)見表1，水位擬合情況見圖3、4。

圖2　模型網(wǎng)格平面Fig.2　Model grid plane

3　靈敏度分析

根據(jù)研究區(qū)的實際水文地質(zhì)條件，同時綜合數(shù)值模型，選取了潛水的降水入滲補給系數(shù)、滲透系數(shù)、給水度、弱透水層的垂向滲透系數(shù)、承壓水的滲透系數(shù)和彈性釋水系數(shù)進(jìn)行分析，選取的水平數(shù)為5，即第3個是已識別的參數(shù)，其余為已識別的參數(shù)上下浮動10%和20%(表2)。

表1　識別后水文地質(zhì)參數(shù)值

注：Ⅰ、Ⅱ、Ⅲ分別表示潛水含水層、承壓水含水層及弱透水層的3個分區(qū)。

Note：Ⅰ，Ⅱ，Ⅲ stands for phreatic aquifer, confined aquifer and weak permeable layer.

圖3　93號井識別期水位擬合情況Fig.3　The fitting of water level in the identification stage of No.93 well

圖4　93號井驗證期水位擬合情況Fig.4　The fitting of water level in the verification stage of No.93 well

根據(jù)表2中的影響因素和水平，選擇了正交試驗L25(56)(表3)。根據(jù)表3中的因素及水平值，在Visual MODFlow中輸入相應(yīng)參數(shù)的對應(yīng)值進(jìn)行試驗，共25次，將每次試驗結(jié)果的水位值輸出，將水位值代入式(2)，計算每個試驗的敏感度指數(shù)[3]。

(2)

式中，L為敏感度指數(shù)，m2；m為長觀孔個數(shù)(有10個長觀孔，m取10)；n為識別階段各長觀孔水位觀測次數(shù)(觀測次數(shù)為35，n取35)；h＇ij為敏感度試驗時計算的與各長觀孔觀測時刻對應(yīng)的水位，m；hij為識別階段與各長觀孔觀測時刻對應(yīng)的計算水位，m。計算得到的敏感度指數(shù)見表3。由表3可知，對地下含水系統(tǒng)計算水位的影響從大到小依次為B、F、E、D、C、A。

表2　各參數(shù)不同水平值

注：Ⅰ、Ⅱ、Ⅲ分別表示潛水含水層、承壓水含水層及弱透水層的3個分區(qū)。

Note:Ⅰ,Ⅱ,Ⅲ stands for phreatic aquifer, confiend aquifer and weak permeable layer.

表3　L25(56)正交試驗結(jié)果

接下表

續(xù)表3

4　結(jié)論

參數(shù)靈敏度分析選取的6個介質(zhì)參數(shù)的敏感度指數(shù)大小代表其影響數(shù)值模擬結(jié)果的強弱，研究結(jié)果表明，在用Visual MODFlow建立的研究區(qū)地下水流數(shù)值模型中，分別將以上6個介質(zhì)參數(shù)做相同幅度的改變，重新運行計算出的觀測孔處的水位動態(tài)，與改變參數(shù)前的相比，變化幅度不同，其中改變潛水含水層的滲透系數(shù)變化幅度最大，因此在利用Visual MODFlow做水資源評價前應(yīng)查明研究區(qū)的水文地質(zhì)參數(shù)特別是潛水的滲透系數(shù)，這對于正確建立模型十分必要。

[1] 束龍倉，劉佩貴，劉波，等.傍河水源地數(shù)學(xué)模型的參數(shù)靈敏度分析：以遼寧省北票市某傍河水源地為例[J].工程勘察，2006(8)：29-31.

[2] 薛禹群，謝春紅.地下水?dāng)?shù)值模擬[M].北京：科學(xué)出版社，2007.

[3] 翟遠(yuǎn)征.鄂爾多斯白堊系盆地烏蘭陶勒蓋水源地地下水流數(shù)值模擬中介質(zhì)參數(shù)的敏感度分析[D].長春：吉林大學(xué)，2008.

[4] 張麗華，蘇小四，孟祥菲，等.地下水流數(shù)值模擬參數(shù)全局靈敏度分析[J].中國農(nóng)村水利水電，2014(8)：92-97.

[5] 苗勝軍，李長洪，文俊，等.基于正交試驗設(shè)計的滑帶土參數(shù)敏感性分析[J].中國礦業(yè)2007，16(9)：76-79.

[6] 郝靜，賈仰文，張永祥，等.應(yīng)用正交試驗法分析地下水流模型參數(shù)靈敏度[J].人民黃河，2015(9)：66-68.

[7] 高翠萍.黑龍江省樺南縣地下水資源評價與合理開發(fā)利用研究[D].長春：吉林大學(xué)，2012.

Study on Groundwater Numerical Simulation Based on Parameter Sensitivity Analysis—A Case Study of Huanan County in Heilongjiang Province

GAO Cui-ping1, DU Xin-qiang2, YANG Yue-suo1et al

(1. Key Laboratory of Regional Environment & Eco-remediation of Ministry of Education, Shenyang University, Shenyang, Liaoning 110044; 2. Key Laboratory of Groundwater Resources and Environment of Ministry of Education, Jilin University, Changchun, Jilin 130021)

The orthogonal experimental method was used to do the parameter sensitivity analysis of groundwater numerical simulation in Huanan County of Heilongjiang Province. 6 parameters including feed coefficient of precipitation infiltration, permeability coefficient, water supply degree, vertical permeability coefficient of weak permeable layer, permeability coefficient of confined water, elastic storativity, were considered to do sensitivity analysis according to the hydrogeological conditions, and 5 level values of each parameter were selected. The results showed that the permeability coefficient of unconfined aquifer has the greatest impact to the water level calculated by the numerical simulation model. It is necessary to find out the hydrogeological parameters especially the permeability coefficient of unconfined aquifer when doing water resource evaluation by use of Visual MODFlow.

Sensitivity analysis; Orthogonal experiment method; Numerical simulation; Groundwater

遼寧省創(chuàng)新團隊項目(LT2015017)；沈陽市科技創(chuàng)新專項資金(F15-113-9-00)。

高翠萍(1986- )，女，內(nèi)蒙古通遼人，助理工程師，從事地下水?dāng)?shù)值模擬及地下水污染修復(fù)研究。

2016-06-08

S 273.4

A

0517-6611(2016)20-053-03