趙 潔,林 錦,吳劍鋒,吳吉春

(1.華北水利水電大學(xué)水資源學(xué)院，河南 鄭州 450046；2.表生地球化學(xué)教育部重點(diǎn)實(shí)驗(yàn)室/南京大學(xué)地球科學(xué)與工程學(xué)院水科學(xué)系，江蘇 南京 210023；3 南京水利科學(xué)研究院，江蘇 南京 210029)

海水入侵是濱海地區(qū)含水層的常見現(xiàn)象，嚴(yán)重影響了該地區(qū)地下水水質(zhì)[1-2]。了解濱海地區(qū)海水入侵的范圍、程度及未來發(fā)展趨勢[3-5]，可為區(qū)域地下水資源合理配置提供科學(xué)依據(jù)[6-8]。多項(xiàng)研究表明海平面升降、降水量變化等氣候因素對(duì)海水入侵程度均產(chǎn)生較大影響[9-14]。Xiao等[15]利用SEAWAT在Florida州沿海平原建立了一個(gè)三維變密度海水入侵模擬及預(yù)測模型，發(fā)現(xiàn)降水量變化和海平面升降對(duì)區(qū)域地下水位和水質(zhì)有重要影響。HUGMAN等[16]和Green等[17]分別在葡萄牙濱海地區(qū)和加拿大Atlantic濱海地區(qū)進(jìn)行了類似研究，后者研究表明降低降水補(bǔ)給量對(duì)咸淡水過渡帶附近的淺層至中層含水層影響較大。Carneiro等[18]和Unsal等[19]分別基于研究區(qū)水文地質(zhì)資料及政府間氣候變化專門委員會(huì)(IPCC)提供的降水量數(shù)據(jù)等資料，建立了不同地區(qū)的變密度地下水?dāng)?shù)值模擬及預(yù)測模型，研究結(jié)果均表明濱海含水層的海水入侵程度及地下水資源量受氣候變化因素影響較大，若減少補(bǔ)給且基于IPCC數(shù)據(jù)升高海平面，濱海含水層中海水入侵程度將會(huì)更加嚴(yán)重。此外，政府間氣候變化專門委員會(huì)IPCC第四次評(píng)價(jià)報(bào)告重申了未來氣候變化將會(huì)持續(xù)。目前，這些研究多集中在國外，國內(nèi)關(guān)于建立三維的實(shí)際海水入侵區(qū)地下水?dāng)?shù)值模擬模型及預(yù)測模型方面的研究較少。

自二十世紀(jì)八九十年代起，大連濱海地區(qū)海水入侵程度日趨嚴(yán)重(圖1)，本文基于校正好的三維變密度地下水?dāng)?shù)值模擬模型[20-21]，針對(duì)未來不同降水量情景設(shè)計(jì)了多種預(yù)測方案，包括降水頻率預(yù)測方案和CMIP5氣候模式下的未來降水量預(yù)測方案，對(duì)研究區(qū)海水入侵未來趨勢進(jìn)行了預(yù)測。

圖1 大連市歷年海水入侵變化趨勢圖

1 研究區(qū)概況

大連市地處遼東半島南端，屬于三面環(huán)海的丘陵地帶。周水子研究區(qū)是大連市海水入侵程度最嚴(yán)重地段之一，總面積約30 km2，研究區(qū)具體位置如圖2所示。研究區(qū)位于具有海洋性特點(diǎn)的大陸性季風(fēng)氣候區(qū)，雨季為6—9月。研究區(qū)第四系發(fā)育較差且局部缺失，之下為上元古界震旦系，巖性為泥粉晶灰?guī)r與白云巖，局部發(fā)育有輝綠巖條帶。在復(fù)雜的褶皺斷裂與負(fù)地形部位巖溶較發(fā)育，研究區(qū)主要的含水層由碳酸鹽巖裂隙巖溶含水層和基巖裂隙巖溶含水層組成。補(bǔ)給項(xiàng)包括降水入滲和充水?dāng)鄬拥膫?cè)向補(bǔ)給，排泄項(xiàng)包括人工開采、垂向蒸發(fā)和側(cè)向徑流排泄量。

圖2 研究區(qū)水文地質(zhì)圖

圖3 概念模型圖

將研究區(qū)裂隙巖溶含水層概化為等效多孔介質(zhì)，建立了一個(gè)三維、不等厚、非均質(zhì)、各向異性的承壓—非承壓地下水流及溶質(zhì)運(yùn)移模擬模型[20-21]。根據(jù)含水層、隔水層的巖性及分布、地質(zhì)構(gòu)造、海水與含水層的水力聯(lián)系等信息確定了模型的邊界條件(圖3)，該海水入侵?jǐn)?shù)值模擬模型由一個(gè)水流模型和一個(gè)溶質(zhì)運(yùn)移模型組成[22]。對(duì)于水流模型來說，北部透水性差的輝綠巖條帶、東北部透水性差的平移斷層及南部水文地質(zhì)界線，可處理為隔水邊界；研究區(qū)西部邊界及中北部的充水?dāng)鄬?，定義為給定流量邊界；東部處理為海邊給定水頭邊界。對(duì)于溶質(zhì)運(yùn)移模型來說，東部海灣設(shè)為定濃度邊界；中部和西部的導(dǎo)水?dāng)鄬釉O(shè)為零濃度邊界；其余邊界均設(shè)置為零質(zhì)量通量邊界?；谝陨细拍钅Ｐ?，可建立相應(yīng)水流及溶質(zhì)運(yùn)移數(shù)學(xué)模型。Langevin等[23]基于質(zhì)量守恒定律及達(dá)西定律，推導(dǎo)得到的變密度地下水流運(yùn)動(dòng)控制方程式：

(1)

初始條件：

H(x,y,z,0)=H0(x,y,z) (x,y,z)∈Ω

(2)

第一類邊界條件：

H(x,y,z,t)|Γ1=H1(x,y,z,t) (x,y,z)∈Γ1

(3)

第二類邊界條件：

(4)

式中：x，y，z——與主滲透方向一致的坐標(biāo)軸；

Kf——等效淡水滲透系數(shù)；

Sf——等效淡水單位貯水系數(shù)；

c——溶質(zhì)濃度。

變密度地下水溶質(zhì)運(yùn)移模型的數(shù)學(xué)模型：

(5)

初始條件和邊界條件：

c(x,y,z,0)=c0(x,y,z)∈R

(6)

c(x,y,z,t)|Γ1=c1(x,y,z,t) (x,y,z)∈R

(7)

(8)

式中：R——計(jì)算域；

D——彌散系數(shù)張量；

c——溶質(zhì)濃度；

c*——井中鹽分的濃度；

c0——初始濃度；

c1——一類邊界濃度；

ux、uy、uz——x,y,z方向上的地下水實(shí)際平均流速。

由式(1)～(8)構(gòu)成了一個(gè)描述密度逐漸變化的海水入侵區(qū)地下水運(yùn)動(dòng)的完整數(shù)學(xué)模型。

采用變密度地下水?dāng)?shù)值模擬工具SEAWAT—2000進(jìn)行求解。平面上將研究區(qū)離散為81 行、166 列邊長為50 m的正方形網(wǎng)格，垂向上將模型劃分為5 層。模型中水文地質(zhì)參數(shù)的分區(qū)及初值來源于大連市地下水勘查相關(guān)資料[20-21]，模型校正后參數(shù)值可查閱文獻(xiàn)[20-21]。開采井采用Well程序包處理，中部、東部充水?dāng)鄬拥膫?cè)向補(bǔ)給采用Well程序包處理(設(shè)置為注水井)，降水入滲在模型中用Recharge程序包處理。模型初始條件由2007 年的水位觀測值及濃度觀測值獲得。時(shí)間模擬時(shí)長為2007年10月1日—2010年9月30日，共36 個(gè)應(yīng)力期，前兩年為校正期，后一年為驗(yàn)證期。采用試錯(cuò)法手動(dòng)調(diào)整對(duì)模型進(jìn)行校正，擬合結(jié)果良好。

2 氣候變化預(yù)測方案

2.1 預(yù)測方案構(gòu)建

基于校正模型，針對(duì)未來不同降水量情景設(shè)計(jì)了多種預(yù)測方案[20](表1)，包括降水頻率和CMIP5氣候模式下的未來降水量預(yù)測方案。預(yù)測時(shí)間為2010年9月30日—2040年9月30日，除降水條件發(fā)生變化外，其余所有條件與2010年保持一致。

降水頻率分析法將1964—2012年的年降水量由大到小進(jìn)行排序，并計(jì)算對(duì)應(yīng)年份的降水頻率，20 %為豐水年，50 %為平水年，80 %為枯水年。未來30 年的降水量由CMIP5(聯(lián)合模型對(duì)比項(xiàng)目第五階段)的7 種氣候模式下的預(yù)測降水量獲得。2008年9月，來自世界各地的20 個(gè)氣候模型小組召開會(huì)議，將地球系統(tǒng)項(xiàng)目的分析、集合與模擬的輸入相結(jié)合，設(shè)置了一套新的標(biāo)準(zhǔn)化氣候模型試驗(yàn)即聯(lián)合模型對(duì)比項(xiàng)目第五階段CMIP5，它提供氣候模型判斷、驗(yàn)證、比較、說明和數(shù)據(jù)存取的共享結(jié)構(gòu)[24-27]。本研究中采用了CMIP5中的7 個(gè)氣候模式下的未來降水?dāng)?shù)據(jù)，這7種氣候模式的簡稱分別為BCC(Beijing Climate Center)、BNU(Beijing Normal University Earth-System Models)、CNRM(Centre National de RecherchesMétéorologiques)、GISS(Goddard Institute for Space Studies)、MIROC(Model for Interdisciplinary Research on Climate)、MPI(Max Planck Institute)和 MRI(Meteorological Research Institute)。如BNU是以陸面模式CoLM為核心包含全碳循環(huán)過程的地球系統(tǒng)模式，將多種分量模式通過耦合器技術(shù)相耦合，并通過不同站點(diǎn)發(fā)布模式數(shù)據(jù)，實(shí)現(xiàn)數(shù)據(jù)共享[28]。

表1 海水入侵程度預(yù)測方案

注：S1-Y中S1代表未來降水方案1，Y表示考慮了潮汐作用；S4-1-Y中S4代表特定氣候模式下的降水方案4，1代表排放情景，Y表示該預(yù)測模型考慮了潮汐作用；除歷史降水外，其余方案均為為不同氣候模式下的未來降水?dāng)?shù)據(jù)。

CMIP5包含不同氣候模式試驗(yàn)[29-30]，在試驗(yàn)中模型提供了不同的氣候“強(qiáng)迫”反應(yīng)改變大氣組成和地面覆蓋。如在CMIP5試驗(yàn)中BNU氣候模式有3種溫室氣體濃度排放情景(低排放情景RCP_2.6，中等排放情景RCP_4.5，高排放情景RCP_8.5)，預(yù)測結(jié)果截止到2100 年，試驗(yàn)結(jié)果顯示BNU氣候模式3 種排放情景下的未來年降水量范圍分別為550.2～1 228.8 mm，607.1～1 069.1 mm，和511.6～1 007.6 mm，年平均降水量為762.7，820.3，773.8mm。

2.2 預(yù)測結(jié)果與分析

為方便對(duì)比，將所有模型層在水平方向上Cl-濃度大于250mg/L的單元格面積總和定義為總的海水入侵復(fù)合區(qū)面積(即表2 中S*)。圖4中顯示的海水入侵區(qū)投影面積線為所有垂向上的含水層250 mg/L Cl-等值線在平面上進(jìn)行投影，疊加而形成的海水入侵范圍最大區(qū)域的邊界線。

圖4(a)、圖5及表2為降水頻率分析得到的預(yù)測降水量對(duì)未來海水入侵程度影響的結(jié)果。從圖4(a)可看出枯水年(S1-Y)較平水年及豐水年(S2-Y 和S3-Y)海水入侵程度稍微嚴(yán)重。從圖5(垂向1-5層海水入侵程度圖)還可看出，越接近含水層底部上述結(jié)果越明顯。從表2 可看出，與現(xiàn)狀年(2010年10月1日)的海水入侵程度相比，考慮潮汐作用下30 年后總的海水入侵復(fù)合區(qū)面積在枯水年(S1-Y)、平水年(S2-Y)及豐水年(S3-Y)預(yù)測方案下分別增加了20.66 %、18.83 %、16.91 %，又三種預(yù)測方案下的年平均降水量分別為444.6、614.9、774.3 mm(表1)，枯水年年平均降水量最少，故推測出現(xiàn)該模擬結(jié)果的可能原因是，枯水年較少的降水補(bǔ)給導(dǎo)致了相對(duì)較低的地下水位，淡水位與海水位之間的不平衡加劇，海水入侵更容易發(fā)生，豐水年則完全相反，故未來海水入侵程度與未來降水量呈負(fù)相關(guān)關(guān)系。二十世紀(jì)八九十年代大連海水入侵程度最嚴(yán)重，政府高度重視并采取了“引碧入連”工程等多種措施來控制地下水開采，故海水入侵程度雖已相對(duì)減弱，但未來海水入侵趨勢不可避免。

圖4(b)為BNU氣候模式3種溫室氣體濃度排放情景下的未來海水入侵程度對(duì)比圖，可看出不同溫室氣體排放情景下的預(yù)測降水量對(duì)濱海地區(qū)含水層系統(tǒng)海水入侵程度的影響。從圖4(b)可看出，模擬預(yù)測方案(S4-1-Y、S4-2-Y及S4-3-Y)中海水入侵程度幾乎相同，從圖上放大區(qū)域可以看出，紅色線條即方案S4-1-Y下海水入侵程度有微小減弱(程度太小，幾乎可忽略)，又由表1可得，BNU氣候模式3種溫室氣體濃度排放情景下的未來年平均降水量分別為726.7、820.3、773.8 mm，相差不太大，方案S4-1-Y下未來年平均降水量相對(duì)最小。推測出現(xiàn)該預(yù)測結(jié)果可能原因是，相差不太大的未來降水量導(dǎo)致未來地下水位也差別不大，故未來海水入侵程度也幾乎無差別，方案S4-1-Y下未來年平均降水量相對(duì)最小，故該方案下未來海水入侵程度較其它方案呈現(xiàn)微小的嚴(yán)重趨勢，但三種方案下預(yù)測結(jié)果基本差別不大。因此，可近似認(rèn)為未來海水入侵程度與未來降水量呈負(fù)相關(guān)關(guān)系。值得注意的是，其它6種氣候模式下的海水入侵程度預(yù)測結(jié)果與此預(yù)測結(jié)果基本相似，即不同的溫室氣體濃度排放情景對(duì)海水入侵的結(jié)果幾乎沒有影響，故隨機(jī)選擇BNU氣候模式下的預(yù)測結(jié)果作為展示。

表2 不同降水方案下海水入侵程度對(duì)比

圖4 不同預(yù)測方案下海水入侵復(fù)合區(qū)面積對(duì)比

圖5 方案S1-Y, S2-Y和S3-Y下1～5層海水入侵程度對(duì)比

圖6 7種氣候模式下未來總的海水入侵復(fù)合區(qū)面積對(duì)比

圖6為7種氣候模式下30 年后的總的海水入侵復(fù)合區(qū)面積對(duì)比柱狀圖。圖4(c)為CNRM、MIROC、MPI這3 種氣候模式在溫室氣體濃度排放情景RCP_4.5下，未來30年海水入侵程度的對(duì)比。從圖7、圖4(c)及表2可看出，MPI氣候模式(S7-2-Y)下的海水入侵程度最嚴(yán)重，總的海水入侵復(fù)合區(qū)面積為2.41 km2，CNRM氣候模式(S5-2-Y)下的海水入侵程度最不嚴(yán)重，總的海水入侵復(fù)合區(qū)面積為2.38 km2，MIROC氣候模式(S6-2-Y)下的海水入侵程度適中，總的海水入侵復(fù)合區(qū)面積為2.39 km2。由表1可得，CNRM、MIROC、MPI三種氣候模式溫室氣體濃度為RCP_4.5的排放情景下，未來年平均降水量分別為804.4，687.0，645.7 mm。推測可能原因是，3 種預(yù)測方案中CNRM氣候模式下未來年平均降水量相對(duì)最大，因此該模式下預(yù)測模型中地下水位相對(duì)最高，最不容易發(fā)生海水入侵，模型預(yù)測結(jié)果也印證了這個(gè)觀點(diǎn)，故可認(rèn)為未來海水入侵程度與未來降水量呈負(fù)相關(guān)關(guān)系，降水量越大，海水入侵程度越不嚴(yán)重，MPI氣候模式下預(yù)測結(jié)果也與該觀點(diǎn)吻合。

隨著經(jīng)濟(jì)發(fā)展和人口增長，對(duì)地下水資源量的需求會(huì)源源不斷地增加，未來海水入侵程度較現(xiàn)狀會(huì)更加嚴(yán)重，故本研究可為合理配置地下水可開采量來控制海水入侵提供了理論基礎(chǔ)。

3 結(jié)論

(1)設(shè)計(jì)了兩種未來降水量預(yù)測方案對(duì)未來海水入侵趨勢進(jìn)行預(yù)測：未來海水入侵程度與未來降水量呈近似負(fù)相關(guān)關(guān)系，即降水量越少，海水入侵程度越嚴(yán)重。不同降水補(bǔ)給量對(duì)模型最后運(yùn)行結(jié)果影響不同，枯水年下的未來海水入侵程度較豐水年、平水年更加嚴(yán)重。7 種氣候模式中，MPI氣候模式(S7-2-Y)下的未來海水入侵程度最嚴(yán)重，CNRM氣候模式(S5-2-Y)下的未來海水入侵程度最不嚴(yán)重。每種氣候模式不同溫室氣體濃度排放情景下的預(yù)測降水量對(duì)未來海水入侵程度幾乎無影響。與CMIP5不同氣候模式下的預(yù)測降水量相結(jié)合的海水入侵預(yù)測模型，對(duì)未來海水入侵趨勢的預(yù)測較降水頻率預(yù)測方案精度更高。兩種預(yù)測方案下，未來海水入侵程度將會(huì)進(jìn)一步加重。

(2)人口增長和經(jīng)濟(jì)發(fā)展將導(dǎo)致用水量逐步增加，濱海地區(qū)海水入侵程度也將會(huì)逐步嚴(yán)重，模擬及預(yù)測模型可為未來的地下水資源管理工作提供理論基礎(chǔ)，防止海水入侵程度繼續(xù)惡化。