李選彧

(遼寧省遼陽(yáng)水文局，遼寧 遼陽(yáng) 111000)

0 引 言

研究表明[1-2]，海水入侵對(duì)濱海地區(qū)地下水水質(zhì)產(chǎn)生嚴(yán)重的影響，其含水層普遍存在海水入侵的現(xiàn)象。為更加科學(xué)合理的配置濱海地區(qū)地下水資源，有必要深入了解其海水入侵程度、范圍以及未來(lái)變化趨勢(shì)。目前，國(guó)內(nèi)外諸多學(xué)者探究了未來(lái)氣候變化與海水入侵程度的響應(yīng)關(guān)系，如海水入侵程度受降雨量變化、海平面上升等因素的影響，結(jié)果發(fā)現(xiàn)海水入侵過程受這些因素的影響較大；Green等基于當(dāng)?shù)厮牡刭|(zhì)資料和降水量數(shù)據(jù)，構(gòu)建了適用于加拿大Atlantic濱海地區(qū)的海水入侵模擬三維數(shù)學(xué)模型，并對(duì)海水入侵響應(yīng)氣候變化的程度進(jìn)行評(píng)價(jià)，結(jié)果表明咸淡水過渡帶附近的中層至淺層含水層受減少降水補(bǔ)給量的影響顯著；HUGMAN等利用溶質(zhì)與水流密度耦合的運(yùn)移模型，預(yù)測(cè)分析了2010-2099年葡萄牙濱海地區(qū)含水層海水入侵受水資源利用、氣候變化等因素的影響，研究表明地下水系統(tǒng)補(bǔ)給直接受到降雨量變化的影響，咸淡水界面的移動(dòng)與海水入侵程度有關(guān)，且水均衡能夠長(zhǎng)期的對(duì)海水入侵產(chǎn)生影響；Xiao等利用溶質(zhì)與水流三維變密度運(yùn)移模型，預(yù)測(cè)了Florida州沿海低洼沖擊平原區(qū)表層含水層響應(yīng)降雨量變化、海平面升降的關(guān)系，結(jié)果顯示氣候變化的主要形式為海平面升降和降雨量變化，并顯著影響著濱海地區(qū)地下水水質(zhì)與水位；Carneiro等基于IPCC提供的降雨量數(shù)據(jù)、海水入侵狀況和研究區(qū)水文地質(zhì)資料，利用變密度地下水?dāng)?shù)值模型預(yù)測(cè)了氣候變化對(duì)地下水資源量、海水入侵程度的影響作用；Unsal等對(duì)沿海含水層中海水入侵程度受氣候變化因素的影響利用三維數(shù)值模擬模型進(jìn)行評(píng)價(jià)，結(jié)果表明IPCC數(shù)據(jù)升高海平面及減少補(bǔ)給能夠進(jìn)一步加劇含水層中海水入侵現(xiàn)象。此外，未來(lái)氣候的持續(xù)變化被IPCC評(píng)價(jià)報(bào)告再次重申。由于枯水年降水量的減少使得沿海含水層補(bǔ)給量不斷下降，從而加重海水入侵程度和地下水咸化的現(xiàn)象。目前，國(guó)外關(guān)于三維數(shù)值模擬海水入侵的研究較多，而國(guó)內(nèi)涉及這些研究的減少。

近年來(lái)，地下水的大量開采使得丹東市濱海地帶海水入侵程度不斷加劇，有學(xué)者對(duì)研究區(qū)海水入侵程度、范圍利用SEAWAT構(gòu)建的三維變密度地下水?dāng)?shù)值模型進(jìn)行預(yù)測(cè)。鑒于此，文章對(duì)海水入侵程度受降水補(bǔ)給變化的影響運(yùn)用校正的數(shù)值模型研究，即探究未來(lái)氣候變化對(duì)海水入侵的影響。針對(duì)CMPI5氣候模式和多種降雨預(yù)測(cè)頻率，設(shè)計(jì)不同情景下未來(lái)降雨量預(yù)測(cè)方案，科學(xué)預(yù)測(cè)未來(lái)海水入侵變化趨勢(shì)，為區(qū)域水資源優(yōu)化配置和合理利用提供科學(xué)指導(dǎo)。丹東市歷年海水入侵變化特點(diǎn)同，見圖1。

圖1 丹東市歷年海水入侵變化特點(diǎn)

1 區(qū)域概況

丹東市地處東北亞中心地帶，屬于一座沿海、沿江、沿邊城市，其核心區(qū)位置處E124°23′、N40°07′，總面積1.52萬(wàn)km2。境內(nèi)水系發(fā)育良好，江河密布，其中流域面積超過4983km2的有渾江、大洋河、鴨綠江和叆河，其它大小河流944條，水源涵養(yǎng)能力好，林草覆蓋率61.6%。該區(qū)域?qū)贉嘏瘞啙駶?rùn)季風(fēng)氣候，由于地貌形態(tài)差異其氣候環(huán)境變化明顯，北部氣溫6-7℃，南部8-9℃，降水量為648.4-1761.7mm，暴雨集中期多7月—8月中旬，夏季降水約占年的2/3，地形以丘陵山地、平原谷底為主。借鑒水文資料，丹東地區(qū)的平均徑流深81mm，地下水補(bǔ)給量16.82億m3，平均水資源量84.88億m3/a，區(qū)域產(chǎn)水模數(shù)57.48萬(wàn)m3/km2，屬于遼寧地區(qū)水資源最豐富且降雨量最多的地區(qū)[3]。

研究區(qū)以局部缺失且發(fā)育較差的第四系為主，巖性多白云巖與泥粉晶灰?guī)r，輝綠巖條帶局部發(fā)育較好。巖溶在負(fù)地形部位與復(fù)雜褶皺斷裂帶較發(fā)育，基巖裂隙巖溶含水層和碳酸鹽巖裂隙巖溶含水層組成研究區(qū)的主要含水層。其中，水?dāng)鄬拥膫?cè)向補(bǔ)給和降水入滲為補(bǔ)給項(xiàng)，而側(cè)向徑流、垂向蒸發(fā)和人工開采為主要排泄項(xiàng)。

2 預(yù)測(cè)模型

2.1 建模流程

以概化的等效多孔介質(zhì)替代研究區(qū)裂隙巖溶含水層，并構(gòu)建一個(gè)非承壓-承壓、各向異性、非均質(zhì)、不等厚、三維的溶質(zhì)及地下水流運(yùn)移模型。模型的邊界條件按照含水層與海水的水力聯(lián)系，隔水層和含水層的地質(zhì)構(gòu)造、分布及巖性等信息確定，將溶質(zhì)運(yùn)移模型與水流模型相耦合構(gòu)成海水入侵?jǐn)?shù)值模擬模型。設(shè)定水流模型的隔水邊界條件為研究區(qū)南部水文地質(zhì)界線、東北部和北部透水性差的平移斷層、灰綠巖條帶，流量邊界條件為中北部的充水?dāng)鄬蛹拔鞑窟吔纾o定的水頭邊界為東部處理的海邊；設(shè)定溶質(zhì)運(yùn)移模型的濃度邊界為東部海灣，零濃度邊界為西部與中部的導(dǎo)水?dāng)鄬?，零質(zhì)量通量邊界設(shè)置為其余各邊界，按照以上概念模型構(gòu)建水流數(shù)學(xué)模型。以達(dá)西定律及質(zhì)量守恒定律為基礎(chǔ)，Langevin等推導(dǎo)了地下水流運(yùn)動(dòng)變密度控制方程，其表達(dá)式為：

(1)

其中，變密度地下水流的初始條件、第一類和第二類邊界條件為：

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

(2)

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

(3)

(4)

式中：Kf、Sf為等效淡水滲透系數(shù)和淡水單位貯水系數(shù)；θ、c為有效孔隙度和溶質(zhì)濃度。推導(dǎo)的溶質(zhì)運(yùn)移數(shù)學(xué)方程，其表達(dá)式為：

(5)

其中，溶質(zhì)運(yùn)移數(shù)學(xué)模型的初始條件、邊界條件如下：

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

(6)

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

(7)

(8)

式中：D、R為彌散系數(shù)張量和計(jì)算域；c、c0、c1為溶質(zhì)濃度、初始濃度和一類邊界濃度；c*為井中鹽分濃度；ux、uy、uz為地下水x、y、z方向上的實(shí)際平均流速。

2.2 校正擬合

采用公式(1)-(8)可以完整的描述地下水運(yùn)動(dòng)模型，該數(shù)學(xué)模型能夠反映海水入侵區(qū)的密度變化特征。然后利用SEAWAT-2000軟件完成變密度地下水的數(shù)值模擬，將研究區(qū)離散成平面上邊長(zhǎng)為50m的正方形網(wǎng)格，該網(wǎng)格行、列數(shù)為158和80，將模型沿垂向分為5層。丹東市地下水勘查相關(guān)資料提供模型中水文地質(zhì)參數(shù)初始值及其分區(qū)，校正模型時(shí)可適當(dāng)調(diào)整分區(qū)數(shù)值及其大小，在水平范圍內(nèi)各分區(qū)水文地質(zhì)存在較大差異，且各水文地質(zhì)參數(shù)值隨著垂向距離的增加而減少，查閱文獻(xiàn)確定具體參數(shù)值。通過設(shè)置注水井(Well程序包實(shí)現(xiàn))處理東部、中部充水?dāng)鄬拥膫?cè)向補(bǔ)給，同時(shí)完成開采井的處理，利用Recharge程序包對(duì)模型的降水入滲過程處理，以2015年的濃度及水位觀測(cè)值獲取模型初始條件。設(shè)置36個(gè)應(yīng)力期(2015.10.1-2018.09.30)作為模擬時(shí)長(zhǎng)，模型校正期為前2a，驗(yàn)證期為后1a，模型的校正利用手動(dòng)調(diào)整試錯(cuò)法，并取得了良好的擬合效果[4]。

2.3 預(yù)測(cè)方案構(gòu)建

考慮CMPI5氣候模式和多種降雨頻率，利用校正的模型設(shè)計(jì)不同情景下未來(lái)降雨量預(yù)測(cè)方案，在此基礎(chǔ)上分析未來(lái)海水入侵的變化趨勢(shì)。模擬時(shí)間2015.10.1-2018.09.30，降水條件為模擬預(yù)測(cè)的唯一變量，即源匯項(xiàng)等其余所有條件均保持不變。海水入侵預(yù)測(cè)方案，見表1。

表1 海水入侵預(yù)測(cè)方案

按照從大到小的順序排列1968-2015年的年降水量，并將每年對(duì)應(yīng)序的列號(hào)利用降水頻率分析法加1生成一組新的數(shù)列，然后求解該數(shù)值的倒數(shù)即可獲取此年份的降水頻率，設(shè)置80%、50%、20%為枯水、平水和豐水年。根據(jù)7種氣候模式(聯(lián)合模型對(duì)比項(xiàng)目第5階段CMIP5)下的預(yù)測(cè)降水量確定未來(lái)30a降水量，CMIP5提供了數(shù)據(jù)存儲(chǔ)、模型比較、驗(yàn)證分析、判斷說(shuō)明的共享結(jié)構(gòu)。未來(lái)降雨數(shù)據(jù)來(lái)源于CMIP5中MRI、MPI、MIROC、GISS、CNRM、BNU、BCC7個(gè)氣候模式的預(yù)測(cè)降水量。例如，覆蓋全碳循環(huán)過程且以陸面模式Co LM為核心的BNU地球系統(tǒng)模式，通過耦合器技術(shù)將海冰、海洋、大氣、陸地分量模式相耦合，并向CMIP5、PCMDI、WDCC站點(diǎn)運(yùn)用地球系統(tǒng)網(wǎng)絡(luò)(ESG)站點(diǎn)發(fā)布模式數(shù)據(jù)，從而實(shí)現(xiàn)數(shù)據(jù)信息的全球共享[4]。

長(zhǎng)期試驗(yàn)和近期試驗(yàn)為CMIP5的兩種氣候?qū)嶒?yàn)?zāi)Ｊ絒5]，為改變長(zhǎng)期試驗(yàn)中地面覆蓋和大氣組成，系統(tǒng)模塊自帶反應(yīng)不同氣候“強(qiáng)迫”的功能。BUN氣候模式在CMIP5試驗(yàn)中有高排放RCP_8.5、中等排放RCP_4.5、低排放RCP_2.6三種溫室氣體排放場(chǎng)景，并設(shè)置2100年為預(yù)測(cè)終止時(shí)間。實(shí)驗(yàn)表明，RCP_8.5、RCP_4.5、RCP_2.6三種排放情景下，BNU氣候模式的降水補(bǔ)給范圍為510.8-1006.5mm、608.2-1071.3mm、552.4-1226.7mm，平均值依次為772.5、821.6、761.8nn。

2.4 結(jié)果分析

為更加直觀的對(duì)比分析，定義總的海水入侵復(fù)合區(qū)面積S*為水平方向上所有模型層Cl-濃度超過250mg/l的單元格面積之和。海水入侵程度對(duì)比，見表2。在水面上對(duì)含水層250mg/L Cl-等值線沿所有垂向上的投影即為海水入侵區(qū)投影面積線，通過疊加生成最大區(qū)域內(nèi)海水的入侵邊界線。

從表2可知，枯水年S1-Y的海水入侵程度稍嚴(yán)重于豐水年及平水年的S2-Y、S3-Y，并且越接近底部含水層該變化趨勢(shì)就越明顯。此外，考慮潮汐作用下的海水入侵復(fù)合區(qū)面積與現(xiàn)狀年(2015.10.1)海水入侵程度相比，枯水年S1-Y、平水年S2-Y、豐水年S3-Y模擬方案增加了20.00%、18.05%、13.17%。結(jié)合表1中，豐水、平水和枯水年的年降雨量為775.5mm、614.0mm、445.1mm，相對(duì)于其它年份降水量最小的為枯水年。進(jìn)一步分析引起該模擬結(jié)果的原因發(fā)現(xiàn)，地下水位在較少降水補(bǔ)給的枯水年相對(duì)較低，由此加劇了海水位與淡水位間的不平衡程度，更易發(fā)生海水入侵的現(xiàn)象，而豐水年的變化趨勢(shì)相反。所以，區(qū)域降雨量與海水入侵程度之間存在負(fù)相關(guān)性，即海水入侵程度隨著降雨量的減少而加劇。地下水在海水入侵最嚴(yán)重的1990s受到嚴(yán)重破壞，并對(duì)當(dāng)?shù)禺a(chǎn)生生活極大的不利影響，為有效控制地下水的開采當(dāng)?shù)卣扇×嗽S多工程措施，在一定程度上減弱了海水入侵程度，但未來(lái)仍不可避免的出現(xiàn)海水入侵現(xiàn)象[6]。

表2 海水入侵程度對(duì)比

續(xù)表2 海水入侵程度對(duì)比

通過對(duì)比3種排放情景下BNU氣候模式的海水入侵程度，可以模擬平原地區(qū)含水層系統(tǒng)受不同排放情境下預(yù)降雨的影響。海水入侵程度在S4-1-Y、S4-2-Y、S4-3-Y平原地區(qū)含水層系統(tǒng)模擬預(yù)測(cè)方案中幾乎相同，幾乎可以忽略S4-1-Y預(yù)測(cè)方案的海水入侵程度。從表1可知，3種排放情景下BNU氣候模式的未來(lái)降雨量為776.2mm、821.0mm、728.3mm，其中未來(lái)降水量最小的為S4-1-Y預(yù)測(cè)方案，各方案相差較小。進(jìn)一步分析引起該預(yù)測(cè)結(jié)果的原因發(fā)現(xiàn)，未來(lái)降雨量相差不大使得地下水位和海水入侵程度相差較小，未來(lái)降雨量相對(duì)最小的S4-1-Y預(yù)測(cè)方案呈現(xiàn)出微小的海水入侵嚴(yán)重趨勢(shì)，由于3種方案的降水量相差較少使得海水入侵程度預(yù)測(cè)結(jié)果也差別不大。所以，可認(rèn)為未來(lái)降雨量與海水入侵程度存在近似的負(fù)相關(guān)性。此外，BNU氣候模式與其它6種模式存在基本相似的預(yù)測(cè)結(jié)果，即海水入侵受不同排放場(chǎng)景的影響可忽略不計(jì)，所以海水入侵預(yù)測(cè)結(jié)果可隨機(jī)選取BNU氣候模式來(lái)展示。

研究區(qū)混合海水入侵面積在7種氣候模式下的對(duì)比情況，海水入侵復(fù)合區(qū)面積，見圖2。根據(jù)RCP_4.5溫室氣體濃度排放場(chǎng)景下3種氣候模式MPI、MIROC、CNRM的海水入侵程度，海水入侵最嚴(yán)重的為S7-2-Y方案(MPI氣候模式)，最不嚴(yán)重的為S5-2-Y方案(CNRM氣候模式)，海水入侵程度適中的為S6-2-Y(MIROC氣候模式)，所對(duì)應(yīng)的總的海水入侵復(fù)合區(qū)面積為2.416km2、2.351km2、2.396km2。由表1可知，RCP_4.5溫室氣體濃度排放場(chǎng)景下3種氣候模式MPI、MIROC、CNRM的未來(lái)降雨量為642.8mm、682.1mm、806.7mm。進(jìn)一步分析引起該預(yù)測(cè)結(jié)果的原因，年降雨量最大的為CNRM氣候模式，所以3種預(yù)測(cè)方案中該模式具有相對(duì)最高的地下水位，其發(fā)生海水入侵的難度也最大，這個(gè)觀點(diǎn)與預(yù)測(cè)結(jié)果相符，由此判定為未來(lái)降雨量與海水入侵程度存在負(fù)相關(guān)性，即海水入侵程度隨降雨量的增大而減弱，該觀點(diǎn)與CNRM氣候模式預(yù)測(cè)保持較高的一致性。隨著人口的增長(zhǎng)、經(jīng)濟(jì)的發(fā)展以及對(duì)水資源需求的不斷增加，未來(lái)海水入侵程度受地下水持續(xù)開發(fā)影響將更加嚴(yán)重，本研究可為海水入侵控制和地下水資源優(yōu)化配置提供科學(xué)指導(dǎo)[29-36]。

圖2 海水入侵復(fù)合區(qū)面積

3 結(jié) 論

1)對(duì)未來(lái)海水入侵趨勢(shì)利用兩種降雨量預(yù)測(cè)方案進(jìn)行模擬：未來(lái)降雨量與海水入侵程度之間近似存在負(fù)相關(guān)性，即海水入侵程度隨降雨量的減少而增強(qiáng)。模型運(yùn)行結(jié)果受不同降水補(bǔ)給量的影響顯著，未來(lái)平水年、豐水年的海水入侵程度低于枯水年；海水入侵程度在S5-2-Y(CNRM氣候模式)中最不嚴(yán)重，而在S7-2-Y(MPI氣候模式)中最嚴(yán)重，未來(lái)海水入侵程度受不同排放情景下每種氣候模式的影響可不略不計(jì)。海水入侵預(yù)測(cè)模型與CMPI5預(yù)測(cè)降水量相結(jié)合，能夠更加精準(zhǔn)的預(yù)測(cè)海水入侵變化趨勢(shì)，基于兩種方案的海水入侵預(yù)測(cè)結(jié)果均顯示呈加重趨勢(shì)。

2)人口的增長(zhǎng)、經(jīng)濟(jì)的發(fā)展以及對(duì)水資源需求的不斷增加，濱海地區(qū)海水入侵程度受地下水持續(xù)開發(fā)影響將更加嚴(yán)重，預(yù)測(cè)模型可為防止海水入侵的持續(xù)惡化和地下水資源優(yōu)化管理提供決策依據(jù)。