未来气候变化对丹东市海水入侵程度的影响预测

李选彧

(辽宁省辽阳水文局，辽宁 辽阳 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)人口的增长、经济的发展以及对水资源需求的不断增加，滨海地区海水入侵程度受地下水持续开发影响将更加严重，预测模型可为防止海水入侵的持续恶化和地下水资源优化管理提供决策依据。