未來氣候變化對丹東市海水入侵程度的影響預測

李選彧

(遼寧省遼陽水文局，遼寧 遼陽 111000)

0 引 言

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

近年來，地下水的大量開采使得丹東市濱海地帶海水入侵程度不斷加劇，有學者對研究區海水入侵程度、范圍利用SEAWAT構建的三維變密度地下水數值模型進行預測。鑒于此，文章對海水入侵程度受降水補給變化的影響運用校正的數值模型研究，即探究未來氣候變化對海水入侵的影響。針對CMPI5氣候模式和多種降雨預測頻率，設計不同情景下未來降雨量預測方案，科學預測未來海水入侵變化趨勢，為區域水資源優化配置和合理利用提供科學指導。丹東市歷年海水入侵變化特點同，見圖1。

圖1 丹東市歷年海水入侵變化特點

1 區域概況

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

研究區以局部缺失且發育較差的第四系為主，巖性多白云巖與泥粉晶灰巖，輝綠巖條帶局部發育較好。巖溶在負地形部位與復雜褶皺斷裂帶較發育，基巖裂隙巖溶含水層和碳酸鹽巖裂隙巖溶含水層組成研究區的主要含水層。其中，水斷層的側向補給和降水入滲為補給項，而側向徑流、垂向蒸發和人工開采為主要排泄項。

2 預測模型

2.1 建模流程

以概化的等效多孔介質替代研究區裂隙巖溶含水層，并構建一個非承壓-承壓、各向異性、非均質、不等厚、三維的溶質及地下水流運移模型。模型的邊界條件按照含水層與海水的水力聯系，隔水層和含水層的地質構造、分布及巖性等信息確定，將溶質運移模型與水流模型相耦合構成海水入侵數值模擬模型。設定水流模型的隔水邊界條件為研究區南部水文地質界線、東北部和北部透水性差的平移斷層、灰綠巖條帶，流量邊界條件為中北部的充水斷層及西部邊界，給定的水頭邊界為東部處理的海邊；設定溶質運移模型的濃度邊界為東部海灣，零濃度邊界為西部與中部的導水斷層，零質量通量邊界設置為其余各邊界，按照以上概念模型構建水流數學模型。以達西定律及質量守恒定律為基礎，Langevin等推導了地下水流運動變密度控制方程，其表達式為：

(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為等效淡水滲透系數和淡水單位貯水系數；θ、c為有效孔隙度和溶質濃度。推導的溶質運移數學方程，其表達式為：

(5)

其中，溶質運移數學模型的初始條件、邊界條件如下：

(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為彌散系數張量和計算域；c、c0、c1為溶質濃度、初始濃度和一類邊界濃度；c*為井中鹽分濃度；ux、uy、uz為地下水x、y、z方向上的實際平均流速。

2.2 校正擬合

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

2.3 預測方案構建

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

表1 海水入侵預測方案

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

長期試驗和近期試驗為CMIP5的兩種氣候實驗模式[5]，為改變長期試驗中地面覆蓋和大氣組成，系統模塊自帶反應不同氣候“強迫”的功能。BUN氣候模式在CMIP5試驗中有高排放RCP_8.5、中等排放RCP_4.5、低排放RCP_2.6三種溫室氣體排放場景，并設置2100年為預測終止時間。實驗表明，RCP_8.5、RCP_4.5、RCP_2.6三種排放情景下，BNU氣候模式的降水補給范圍為510.8-1006.5mm、608.2-1071.3mm、552.4-1226.7mm，平均值依次為772.5、821.6、761.8nn。

2.4 結果分析

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

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

表2 海水入侵程度對比

續表2 海水入侵程度對比

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

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

圖2 海水入侵復合區面積

3 結 論

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

2)人口的增長、經濟的發展以及對水資源需求的不斷增加，濱海地區海水入侵程度受地下水持續開發影響將更加嚴重，預測模型可為防止海水入侵的持續惡化和地下水資源優化管理提供決策依據。