未来气候变化对大连周水子地区海水入侵程度的影响预测

赵 洁,林 锦,吴剑锋,吴吉春

(1.华北水利水电大学水资源学院，河南 郑州 450046；2.表生地球化学教育部重点实验室/南京大学地球科学与工程学院水科学系，江苏 南京 210023；3 南京水利科学研究院，江苏 南京 210029)

海水入侵是滨海地区含水层的常见现象，严重影响了该地区地下水水质[1-2]。了解滨海地区海水入侵的范围、程度及未来发展趋势[3-5]，可为区域地下水资源合理配置提供科学依据[6-8]。多项研究表明海平面升降、降水量变化等气候因素对海水入侵程度均产生较大影响[9-14]。Xiao等[15]利用SEAWAT在Florida州沿海平原建立了一个三维变密度海水入侵模拟及预测模型，发现降水量变化和海平面升降对区域地下水位和水质有重要影响。HUGMAN等[16]和Green等[17]分别在葡萄牙滨海地区和加拿大Atlantic滨海地区进行了类似研究，后者研究表明降低降水补给量对咸淡水过渡带附近的浅层至中层含水层影响较大。Carneiro等[18]和Unsal等[19]分别基于研究区水文地质资料及政府间气候变化专门委员会(IPCC)提供的降水量数据等资料，建立了不同地区的变密度地下水数值模拟及预测模型，研究结果均表明滨海含水层的海水入侵程度及地下水资源量受气候变化因素影响较大，若减少补给且基于IPCC数据升高海平面，滨海含水层中海水入侵程度将会更加严重。此外，政府间气候变化专门委员会IPCC第四次评价报告重申了未来气候变化将会持续。目前，这些研究多集中在国外，国内关于建立三维的实际海水入侵区地下水数值模拟模型及预测模型方面的研究较少。

自二十世纪八九十年代起，大连滨海地区海水入侵程度日趋严重(图1)，本文基于校正好的三维变密度地下水数值模拟模型[20-21]，针对未来不同降水量情景设计了多种预测方案，包括降水频率预测方案和CMIP5气候模式下的未来降水量预测方案，对研究区海水入侵未来趋势进行了预测。

图1 大连市历年海水入侵变化趋势图

1 研究区概况

大连市地处辽东半岛南端，属于三面环海的丘陵地带。周水子研究区是大连市海水入侵程度最严重地段之一，总面积约30 km2，研究区具体位置如图2所示。研究区位于具有海洋性特点的大陆性季风气候区，雨季为6—9月。研究区第四系发育较差且局部缺失，之下为上元古界震旦系，岩性为泥粉晶灰岩与白云岩，局部发育有辉绿岩条带。在复杂的褶皱断裂与负地形部位岩溶较发育，研究区主要的含水层由碳酸盐岩裂隙岩溶含水层和基岩裂隙岩溶含水层组成。补给项包括降水入渗和充水断层的侧向补给，排泄项包括人工开采、垂向蒸发和侧向径流排泄量。

图2 研究区水文地质图

图3 概念模型图

将研究区裂隙岩溶含水层概化为等效多孔介质，建立了一个三维、不等厚、非均质、各向异性的承压—非承压地下水流及溶质运移模拟模型[20-21]。根据含水层、隔水层的岩性及分布、地质构造、海水与含水层的水力联系等信息确定了模型的边界条件(图3)，该海水入侵数值模拟模型由一个水流模型和一个溶质运移模型组成[22]。对于水流模型来说，北部透水性差的辉绿岩条带、东北部透水性差的平移断层及南部水文地质界线，可处理为隔水边界；研究区西部边界及中北部的充水断层，定义为给定流量边界；东部处理为海边给定水头边界。对于溶质运移模型来说，东部海湾设为定浓度边界；中部和西部的导水断层设为零浓度边界；其余边界均设置为零质量通量边界。基于以上概念模型，可建立相应水流及溶质运移数学模型。Langevin等[23]基于质量守恒定律及达西定律，推导得到的变密度地下水流运动控制方程式：

(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——与主渗透方向一致的坐标轴；

Kf——等效淡水渗透系数；

Sf——等效淡水单位贮水系数；

c——溶质浓度。

变密度地下水溶质运移模型的数学模型：

(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——计算域；

D——弥散系数张量；

c——溶质浓度；

c*——井中盐分的浓度；

c0——初始浓度；

c1——一类边界浓度；

ux、uy、uz——x,y,z方向上的地下水实际平均流速。

由式(1)～(8)构成了一个描述密度逐渐变化的海水入侵区地下水运动的完整数学模型。

采用变密度地下水数值模拟工具SEAWAT—2000进行求解。平面上将研究区离散为81 行、166 列边长为50 m的正方形网格，垂向上将模型划分为5 层。模型中水文地质参数的分区及初值来源于大连市地下水勘查相关资料[20-21]，模型校正后参数值可查阅文献[20-21]。开采井采用Well程序包处理，中部、东部充水断层的侧向补给采用Well程序包处理(设置为注水井)，降水入渗在模型中用Recharge程序包处理。模型初始条件由2007 年的水位观测值及浓度观测值获得。时间模拟时长为2007年10月1日—2010年9月30日，共36 个应力期，前两年为校正期，后一年为验证期。采用试错法手动调整对模型进行校正，拟合结果良好。

2 气候变化预测方案

2.1 预测方案构建

基于校正模型，针对未来不同降水量情景设计了多种预测方案[20](表1)，包括降水频率和CMIP5气候模式下的未来降水量预测方案。预测时间为2010年9月30日—2040年9月30日，除降水条件发生变化外，其余所有条件与2010年保持一致。

降水频率分析法将1964—2012年的年降水量由大到小进行排序，并计算对应年份的降水频率，20 %为丰水年，50 %为平水年，80 %为枯水年。未来30 年的降水量由CMIP5(联合模型对比项目第五阶段)的7 种气候模式下的预测降水量获得。2008年9月，来自世界各地的20 个气候模型小组召开会议，将地球系统项目的分析、集合与模拟的输入相结合，设置了一套新的标准化气候模型试验即联合模型对比项目第五阶段CMIP5，它提供气候模型判断、验证、比较、说明和数据存取的共享结构[24-27]。本研究中采用了CMIP5中的7 个气候模式下的未来降水数据，这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为核心包含全碳循环过程的地球系统模式，将多种分量模式通过耦合器技术相耦合，并通过不同站点发布模式数据，实现数据共享[28]。

表1 海水入侵程度预测方案

注：S1-Y中S1代表未来降水方案1，Y表示考虑了潮汐作用；S4-1-Y中S4代表特定气候模式下的降水方案4，1代表排放情景，Y表示该预测模型考虑了潮汐作用；除历史降水外，其余方案均为为不同气候模式下的未来降水数据。

CMIP5包含不同气候模式试验[29-30]，在试验中模型提供了不同的气候“强迫”反应改变大气组成和地面覆盖。如在CMIP5试验中BNU气候模式有3种温室气体浓度排放情景(低排放情景RCP_2.6，中等排放情景RCP_4.5，高排放情景RCP_8.5)，预测结果截止到2100 年，试验结果显示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 预测结果与分析

为方便对比，将所有模型层在水平方向上Cl-浓度大于250mg/L的单元格面积总和定义为总的海水入侵复合区面积(即表2 中S*)。图4中显示的海水入侵区投影面积线为所有垂向上的含水层250 mg/L Cl-等值线在平面上进行投影，叠加而形成的海水入侵范围最大区域的边界线。

图4(a)、图5及表2为降水频率分析得到的预测降水量对未来海水入侵程度影响的结果。从图4(a)可看出枯水年(S1-Y)较平水年及丰水年(S2-Y 和S3-Y)海水入侵程度稍微严重。从图5(垂向1-5层海水入侵程度图)还可看出，越接近含水层底部上述结果越明显。从表2 可看出，与现状年(2010年10月1日)的海水入侵程度相比，考虑潮汐作用下30 年后总的海水入侵复合区面积在枯水年(S1-Y)、平水年(S2-Y)及丰水年(S3-Y)预测方案下分别增加了20.66 %、18.83 %、16.91 %，又三种预测方案下的年平均降水量分别为444.6、614.9、774.3 mm(表1)，枯水年年平均降水量最少，故推测出现该模拟结果的可能原因是，枯水年较少的降水补给导致了相对较低的地下水位，淡水位与海水位之间的不平衡加剧，海水入侵更容易发生，丰水年则完全相反，故未来海水入侵程度与未来降水量呈负相关关系。二十世纪八九十年代大连海水入侵程度最严重，政府高度重视并采取了“引碧入连”工程等多种措施来控制地下水开采，故海水入侵程度虽已相对减弱，但未来海水入侵趋势不可避免。

图4(b)为BNU气候模式3种温室气体浓度排放情景下的未来海水入侵程度对比图，可看出不同温室气体排放情景下的预测降水量对滨海地区含水层系统海水入侵程度的影响。从图4(b)可看出，模拟预测方案(S4-1-Y、S4-2-Y及S4-3-Y)中海水入侵程度几乎相同，从图上放大区域可以看出，红色线条即方案S4-1-Y下海水入侵程度有微小减弱(程度太小，几乎可忽略)，又由表1可得，BNU气候模式3种温室气体浓度排放情景下的未来年平均降水量分别为726.7、820.3、773.8 mm，相差不太大，方案S4-1-Y下未来年平均降水量相对最小。推测出现该预测结果可能原因是，相差不太大的未来降水量导致未来地下水位也差别不大，故未来海水入侵程度也几乎无差别，方案S4-1-Y下未来年平均降水量相对最小，故该方案下未来海水入侵程度较其它方案呈现微小的严重趋势，但三种方案下预测结果基本差别不大。因此，可近似认为未来海水入侵程度与未来降水量呈负相关关系。值得注意的是，其它6种气候模式下的海水入侵程度预测结果与此预测结果基本相似，即不同的温室气体浓度排放情景对海水入侵的结果几乎没有影响，故随机选择BNU气候模式下的预测结果作为展示。

表2 不同降水方案下海水入侵程度对比

图4 不同预测方案下海水入侵复合区面积对比

图5 方案S1-Y, S2-Y和S3-Y下1～5层海水入侵程度对比

图6 7种气候模式下未来总的海水入侵复合区面积对比

图6为7种气候模式下30 年后的总的海水入侵复合区面积对比柱状图。图4(c)为CNRM、MIROC、MPI这3 种气候模式在温室气体浓度排放情景RCP_4.5下，未来30年海水入侵程度的对比。从图7、图4(c)及表2可看出，MPI气候模式(S7-2-Y)下的海水入侵程度最严重，总的海水入侵复合区面积为2.41 km2，CNRM气候模式(S5-2-Y)下的海水入侵程度最不严重，总的海水入侵复合区面积为2.38 km2，MIROC气候模式(S6-2-Y)下的海水入侵程度适中，总的海水入侵复合区面积为2.39 km2。由表1可得，CNRM、MIROC、MPI三种气候模式温室气体浓度为RCP_4.5的排放情景下，未来年平均降水量分别为804.4，687.0，645.7 mm。推测可能原因是，3 种预测方案中CNRM气候模式下未来年平均降水量相对最大，因此该模式下预测模型中地下水位相对最高，最不容易发生海水入侵，模型预测结果也印证了这个观点，故可认为未来海水入侵程度与未来降水量呈负相关关系，降水量越大，海水入侵程度越不严重，MPI气候模式下预测结果也与该观点吻合。

随着经济发展和人口增长，对地下水资源量的需求会源源不断地增加，未来海水入侵程度较现状会更加严重，故本研究可为合理配置地下水可开采量来控制海水入侵提供了理论基础。

3 结论

(1)设计了两种未来降水量预测方案对未来海水入侵趋势进行预测：未来海水入侵程度与未来降水量呈近似负相关关系，即降水量越少，海水入侵程度越严重。不同降水补给量对模型最后运行结果影响不同，枯水年下的未来海水入侵程度较丰水年、平水年更加严重。7 种气候模式中，MPI气候模式(S7-2-Y)下的未来海水入侵程度最严重，CNRM气候模式(S5-2-Y)下的未来海水入侵程度最不严重。每种气候模式不同温室气体浓度排放情景下的预测降水量对未来海水入侵程度几乎无影响。与CMIP5不同气候模式下的预测降水量相结合的海水入侵预测模型，对未来海水入侵趋势的预测较降水频率预测方案精度更高。两种预测方案下，未来海水入侵程度将会进一步加重。

(2)人口增长和经济发展将导致用水量逐步增加，滨海地区海水入侵程度也将会逐步严重，模拟及预测模型可为未来的地下水资源管理工作提供理论基础，防止海水入侵程度继续恶化。