基于Visual Modflow的秦皇島市海水入侵剖面

（河北省地礦局秦皇島礦產(chǎn)水文工程地質(zhì)大隊，河北 秦皇島066001）

海水入侵是現(xiàn)代社會所特有的環(huán)境問題，嚴重破壞了我國部分沿海地區(qū)的環(huán)境，阻礙了當?shù)毓まr(nóng)業(yè)和城市經(jīng)濟的發(fā)展。目前海水入侵問題已經(jīng)引起了國際社會的廣泛關(guān)注[1]。秦皇島是我國北方重要的能源輸出港口、著名的旅游城市。隨著社會經(jīng)濟建設(shè)地快速發(fā)展和城市化進程地加快，地下水開采強度不斷加大，打破了地下水的天然平衡，導(dǎo)致海水入侵，使地下水環(huán)境惡化，從而引起區(qū)域環(huán)境的破壞和生態(tài)系統(tǒng)的失衡，海水入侵地質(zhì)災(zāi)害頻發(fā)[2]。

Visual Modflow是由加拿大Waterloo水文地質(zhì)公司在原Modflow 軟件的基礎(chǔ)上應(yīng)用可視化技術(shù)開發(fā)研制的，其獨特的可視化菜單結(jié)構(gòu)允許操作者簡單快捷地定義模擬區(qū)范圍、選擇單位、方便地分配模型的屬性和邊界條件[3]。這個軟件包主要包括Modflow（水流評價）、Modpath（平面和剖面流線示蹤分析）、MT3D（溶質(zhì)運移評價）和Zone Budget（水量均衡計算）四大模塊。

1 海水入侵剖面選取

本研究選取秦皇島市海水入侵氯離子剖面—HQ6剖面作為示例剖面（見圖1）。

圖1 HQ6模擬區(qū)域示意圖

本區(qū)咸淡水界面存在一個較寬的過渡帶，濱海含水層中的海水入侵問題是可混液體間的水動力彌散問題[4]。鑒于本區(qū)海水入侵的特征和剖面觀測數(shù)據(jù)的限制，本研究建立了一個考慮密度差異的三維對流－彌散海水入侵模型[5]。

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

該模型的數(shù)學(xué)模型由地下水三維水流運移模型和污染質(zhì)運移模型耦合而成[6]。Modflow是一個三維有限差分地下水流動模型，方程如下：

?/?x（Kxx?h/?x）+?/?y（Kyy?h/?y）+?/?z（Kzz?h/?z）-W=Ss?h/?t

其中：Kxx，Kyy和Kzz為滲透系數(shù)在x，y和z方向上的分量。假定滲透系數(shù)的主軸方向與坐標軸方向一致，量綱為（LT-1）；h為水頭（L）；W為單位體積流量（T-1），用以代表流進匯或來自源的水量；Ss為空隙介質(zhì)的貯水率（L-1）；t為時間（T）。

邊界條件：

第一類邊界（水頭邊界）：

h—Γ1=h1（x，y，z）（x，y，z）∈Γ1

第二類邊界（流量邊界）：

K?h/?n—Γ2=q（x，y，z）（x，y，z）∈Γ2

第三類邊界（混合邊界）：

q（x，y，z）—Γ3= k′（h-h0）/(B′）（x，y，z）∈Γ3

式中Γ1、Γ2、Γ3分別表示一、二、三類邊界。

初始條件：

h|t=0=h0（x，y，z，t）Kx、Ky、Kz分別為沿X、Y、Z坐標軸方向的水力傳導(dǎo)率；h（x，y，z，t）為水頭；K為滲透系數(shù)。

污染質(zhì)運移模型（MT3DMS）的基本方程為：（?（nC））/?t=?/?x[n（Dxx?C/?x+Dxy?C/?y+Dxz?C/?z-CVx）]+?/?y[n（Dxy?C/?x+Dyy?C/?y+Dyz?C/?z-CVy）]+?/?z[n（Dxz?C/?x+Dyz?C/?y+Dzz?C/?z-CVz）]+GCe-WC，（x，y，z）∈Ω，t＞0

C（x，y，z，t）—t=0=C0（x，y，z），（x，y，z）∈Ω，t=0

C（x，y，z，t）—Г1=C1（x，y，z，t），（x，y，z）∈Г1，t＞0

D·C/（?n→）—Г2=f2（x，y，z，t），（x，y，z）∈Г2，t＞0

式中：n為有效孔隙度；C（x，y，z，t）為污染溶質(zhì)濃度；Ce為污染源濃度；D為彌散系數(shù)張量，用矩陣表示為：

矩陣中各個分量分別為：

V為地下水滲透速率，V=（Vx2+Vy2+Vz2）Vx、Vy、Vz為V在x、y、z軸上的投影Vx=-1/nKx?h/?x Vy=-1/nKy?h/?y Vz=-1/nKz?h/?z

3 海水入侵剖面模型建立

3.1 概念模型

實例區(qū)南部近海，北部為洋河，東西向剖面不做研究，可設(shè)定為人為邊界。自上而下可分為8個含水層。其中，亞粘土透水性較差，為弱透水層；而砂卵石層、砂礫石層、細砂層為透水層；下伏基巖為隔水底板。

根據(jù)Visual Modflow 用戶手冊的格式和操作步驟錄入數(shù)據(jù)，將模擬區(qū)平面剖分為10行、95列，垂向分為5層，共剖分8 550個六面體單元。本次建模選擇的長度單位為m，時間單位為d，滲透系數(shù)單位為m/s，抽水率單位為m/d，垂直補給單位為mm/a，質(zhì)量單位為kg，濃度單位為mg/L。在軟件中生成地形、淺層含水層底板、弱透水層底板和深層含水層底板數(shù)據(jù)數(shù)字高程圖形，導(dǎo)入滲透系數(shù)等參數(shù)。

3.2 模型邊界條件

3.2.1 地下水流場邊界

模擬區(qū)南部近海，北部為洋河，故可以設(shè)定為水頭邊界；下伏基巖為隔水邊界；剖面的東西向不考慮，可設(shè)定為人為邊界。

3.2.2 濃度邊界

剖面北向南向海逐漸接近，故將剖面北側(cè)（即模型左側(cè)）設(shè)定為淡水定濃度，模型南側(cè)（即模型右側(cè)）設(shè)定為近海水定濃度。

3.3 時間離散

本模型總模擬時段為10 a（2008年5月1日至2018年5月1日），考慮到實際情況，不同的剖面劃分不同的應(yīng)力期，起始運移步長設(shè)定為0.001 d，遞增系數(shù)為1.1，最大運移步長設(shè)定為200 d。

3.4 模型初始條件

3.4.1 初始流場

采用2008年5月常觀井和統(tǒng)測井實際監(jiān)測數(shù)據(jù)，結(jié)合給定的邊界條件作為已知水位值的結(jié)點，通過Kriging插值方法生成模型的初始水位場。

3.4.2 初始濃度場

由于氯離子是地下水系統(tǒng)中穩(wěn)定的離子，不易被含水層的孔隙介質(zhì)所吸附，而且其濃度與礦化度呈線性關(guān)系，因此將其作為評價研究區(qū)海水入侵標志。以剖面計算區(qū)內(nèi)2008年5月統(tǒng)測井和常觀井的實測濃度值和模型中第一類濃度邊界（北部邊界和海岸邊界）上的結(jié)點的濃度邊界值作為已知濃度值的結(jié)點，通過二次Akima插值法生成模型的初始濃度場。

3.5 模型的驗證

本模型選取2008年5月至2012年5月枯水期（5月份取樣）作為模型的檢驗期，選取濃度觀測井的實測值與模型計算值進行對比，以驗證模型的精度。觀測數(shù)據(jù)和模型運算數(shù)據(jù)比對見圖2。

從觀測井和模型運算的數(shù)據(jù)比對可以看出，模型的運算數(shù)據(jù)都比較平滑，這主要是由于模型建立對各種實際情況的概化引起的，但是就整體趨勢而言，模型和實際觀測數(shù)據(jù)比較切合。

圖2 Hq6-6 觀測數(shù)據(jù)和模型運算數(shù)據(jù)比對曲線

4 海水入侵監(jiān)測剖面預(yù)測

模型充分考慮了研究區(qū)的實際情況。模型的模擬再現(xiàn)了2008年5月至2012年5月海水入侵的全過程。由于該段時間內(nèi)的工況條件沒有發(fā)生改變，模型運算的結(jié)果也與實際情況比較擬合，該過程也完整地反映了Cl-濃度場的變化過程。模型運行初期的濃度變化較大，但是總體趨勢與實際觀測數(shù)據(jù)相符，截止到模型運行1 825 d，模型已經(jīng)完全平衡，為了盡可能準確地預(yù)測不同工況條件下的海水入侵變化趨勢，按照不同的工況條件給予預(yù)測（圖3）。

圖3 2 800 m3/d 工況條件下1 825 d 氯離子變化

以往的預(yù)測針對降水概率等進行不同工況的劃分，從模型的輸入源進行調(diào)整，從而改變了地下水的均衡條件。本次模擬采用改變抽水井抽水量的方式，從排泄項方面來改變地下水的平衡，從而模擬不同工況條件下的海水入侵狀況。

模型0～1 825 d的穩(wěn)定抽水最終確定為2 800 m3/d（圖3），分兩種工況條件以1 825 d作為預(yù)測的起點，分別以2 400 m3/d、3 100 m3/d，進行模擬預(yù)測2015年5月、2018年5月、2020年5月的入侵情況（見圖4、圖5）。

圖4 2 400 m3/d 工況條件下氯離子濃度變化

圖5 3 100 m3/d 工況條件下氯離子濃度變化

5 模擬結(jié)果分析

由模擬結(jié)果可以看出，海水入侵剖面在抽水量2 800 m3/d工況條件下，模型運行1 825 d基本區(qū)域穩(wěn)定。以1 825 d為預(yù)測起點，在2 800 m3/d工況下，海水明顯呈海退現(xiàn)象，而在3 100 m3/d工況下，海水入侵繼續(xù)發(fā)生。

因此秦皇島海水入侵災(zāi)害最主要是由不合理的人類活動（超采地下水，沿海工程）引起的，地下水補給量的減少則對海水入侵起到加速作用。因此，適當?shù)販p少開采量，防止地下水位降低，可以有效地緩解海水入侵，實現(xiàn)地下水資源的可持續(xù)開發(fā)利用以支持該區(qū)經(jīng)濟社會的可持續(xù)發(fā)展。

6 結(jié)論

（1）應(yīng)用Visural Modflow軟件建立海水入侵剖面模型可行，該軟件可有效地將海水入侵可視化，以直觀反映海水入侵問題，所建模型和預(yù)警系統(tǒng)可信度高，為進一步的研究提供了可靠依據(jù)。

（2）本區(qū)的地下水位下降是由于大量開采地下水所造成的，地下水位的快速下降引起海水入侵，使地下水中Cl-含量急劇增加。從模型的預(yù)測結(jié)果看，適當?shù)販p少開采量，防止地下水位降低，可以有效地緩解海水入侵，實現(xiàn)地下水資源的可持續(xù)開發(fā)利用以支持該區(qū)經(jīng)濟社會的可持續(xù)發(fā)展。保護好秦皇島地下水資源，防治海水入侵，必須從本區(qū)的具體地質(zhì)，水文地質(zhì)條件，地下水開采現(xiàn)狀及水資源狀況等實際情況出發(fā)，確定有效的保護方法，才能防止海水入侵的發(fā)展。

[1]高學(xué)平，楊建華，涂向陽，等.帷幕防治海水入侵的數(shù)值模擬研究[J].海洋環(huán)境科學(xué)，2006，25（4）：55-58.

[2]孫娟，楊燕雄.秦皇島市海水入侵特征[J].中國環(huán)境管理干部學(xué)院學(xué)報，2007，17（2）：51-53，68.

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

[4]艾康洪，陳崇希.漫尾島咸淡水界面運移剖面二維流水質(zhì)模型模擬研究[J].勘探科學(xué)技術(shù)，1994（6）：3-9.

[5]成建梅，陳崇希，吉孟瑞，等.山東煙臺夾河中、下游地區(qū)海水入侵三維水質(zhì)數(shù)值模擬研究[J].地質(zhì)前緣，2001，8（1）：179-184.

[6]李俊亭.地下水流數(shù)值模擬[M].北京：地質(zhì)出版社，1989.