基于离散裂隙基质模型的水平双裂隙含水层海水入侵机制

谢一凡 曾祎芃 杨杰 叶逾 吴吉春 鲁春辉

摘要：探究水平双裂隙含水层中的海水入侵机制，对理解裂隙含水层中的海水入侵机制具有重要意义。本文构建了大型室内砂槽模型，并通过离散裂隙基质数值模型分析了不同分布形式的水平双裂隙对盐水楔、混合带等海水入侵特征的影响。结果表明：① 水平双裂隙能显著增强盐水楔及混合带的发展，裂隙间的咸淡水混合尤为显著;2条主干裂隙会限制彼此的作用，位于近海侧、含水层中部的裂隙作用更加显著。② 位于近海侧的裂隙高度越高，则盐水楔发展程度和裂隙间混合作用越显著，而近陆侧裂隙高度则仅对盐水楔尾部有明显影响。③ 位于含水层中部、底部的裂隙，越靠近内陆则盐趾长度更大、盐水体积更小，而远离盐水楔时则会使混合带和盐趾长度急剧变小。

关键词：海水入侵;地下水动力;水平双裂隙;室内试验;数值模拟

中图分类号：P641

文献标志码：A

文章编号：1001-6791（2023）03-0454-11

收稿日期：2022-11-24;

网络出版日期：2023-04-14

网络出版地址：https：∥kns.cnki.net/kcms/detail/32.1309.P.20230413.1832.003.html

基金项目：国家重点研发计划资助项目（2021YFC3200500）;国家自然科学基金资助项目（42277190）

作者简介：谢一凡（1987—），男，江苏南京人，副教授，博士，主要从事海岸带水动力方面研究。

E-mail：yfxie@hhu.edu.cn

通信作者：魯春辉，E-mail：clu@hhu.edu.cn

海水入侵是海岸带一种常见的环境地质灾害，严重制约着沿海经济和社会的可持续发展［1-5］。裂隙含水层是全球分布最广泛的含水层类型之一，也是中国较常见的滨海含水层类型［6］。相比其他类型的含水层，裂隙含水层更易受海水入侵的影响［7-8］，多个国家的原位观测数据和场地案例研究均证实了这一点［5，7，9-10］。目前，裂隙含水层是海水入侵研究领域的前沿方向与研究热点［8］，但尚在起步阶段，主要集中在静态边界下裂隙含水层海水入侵的数值模拟研究，缺乏基于室内试验的海水入侵物理机制研究。因此，通过构建室内物理模型深入探索裂隙含水层中的海水入侵过程与规律，不仅能够丰富、完善海水入侵机制，也对环境地质灾害防治等实际需求具有重要实践意义。

虽然裂隙含水层中的海水入侵机制尚不清晰，但其他类型非均质含水层的研究已表明非均质性可影响地下水流动路径，进而影响海水入侵［1，11］。Lu等［12］研究表明，分层非均质含水层对咸淡水混合带位置有决定性作用，且流线的折射和分散作用会令高渗透层上方的低渗透层中的混合带显著变宽。随机非均质含水层则会导致混合体积的增加［13］。同时，国内外关于单裂隙、裂隙网络含水层中水流和溶质运移的研究，揭示了裂隙产生的优先流动路径对于溶质运移的作用，以及裂隙的物理与几何特征（尺寸、分布、数量等）对于动力学系统的重要性［14-18］，但大部分研究尚未涉及海岸带复杂水盐运移过程。

现有裂隙含水层中海水入侵研究已初步给出研究所用的基本假设与关键因子。现场勘测研究发现裂隙的位置及其与抽水井、海洋的水力连接决定了海水入侵程度［5，7］。裂隙含水层的现有数值模型主要有隐式模型和显式模型2类［19］，分别具有高效率和高精度的特点，2类模型的典型方法分别为等效多孔介质模型（EPM）、双重孔隙度模型（DPM）和离散裂隙基质模型（DFM）、离散裂隙网络模型（DFN）［19］。其中，DFM仅保留了部分裂隙，其余裂隙则通过提升尺度技术与基质一同被等效为多孔介质，并同时考虑裂隙与等效多孔介质中的流动［19］。同时，本领域的数值模拟研究已发现裂隙的密度、方向等特征对海水入侵有宏观影响，且显式模型比隐式模型更适合研究仅含少量主干裂隙的含水层中的海水入侵机制［20-21］。基于达西定律和DFM，Etsias等［19］于2021年建立了裂隙含水层海水入侵物理模型，并对裂隙特征进行了数值敏感性分析，发现裂隙的位置和长度会显著影响盐水楔的分布，且对抽水引发的盐水上升过程也有重要影响［22］。然而，Etsias等［19］仅考虑了含单条裂隙的含水层，而实际滨海裂隙含水层常发育多条平行的主干裂隙［23］。

本文建立水平双裂隙非承压含水层海水入侵的大型砂槽物理模型，应用SUTRA（saturated-unsaturated transport）构建DFM，以厘清不同位置组合的2条水平主干裂隙对海水入侵的宏观影响规律和裂隙间的水盐运移模式。通过构建4种不同类型的裂隙分布组合（含5组试验工况，21组数值模拟工况），详细研究不同水平双裂隙分布对11个典型海水入侵特征值的影响。本研究有助于理解裂隙含水层中海水入侵的基本机制，为实际滨海裂隙含水层的管理、利用和保护提供重要科学依据。

1 模型设置

1.1 室内物理试验设置

试验在长4 m、高0.8 m、宽0.02 m的大型砂槽中进行，砂槽系统如图1（a）所示，砂槽左右两侧设置了盐水和淡水水槽，盐水水槽内设的栅栏式溢流槽以及与淡水水槽相接的蠕动泵可分别精确控制海、陆边界所需的水位（0.54 m）与流量（0.62 mL/s）。参考Etsias等［19］的试验设置，主干裂隙由长1 m、高0.01 m、宽0.02 m且表面分布0.8 mm孔径细孔的中空不锈钢冲孔管模拟（图1（b））;等效多孔介质由粒径为1 mm的石英砂模拟。经过测定，试验砂的饱和渗透率为1.33×10-9 m2，孔隙度为0.45，与相关研究的参数值相近［19，22］。试验砂的纵向弥散度（αL）为0.005 m，横向弥散度（αT）设为αL的1/10［12］，即0.000 5 m。裂隙的孔隙度设为1，其渗透率则通过数值敏感性分析确定［19］，为多孔介质渗透率的300倍，满足Matthi等［24］工作中的裂隙-基质渗透率比范围（102～104）。

配置密度为1 025 kg/m3的盐水模拟海水并加入诱惑红作为染色剂［19］。试验期间每30 min使用Atago PAL-SALT盐度计监测盐度，并缓慢加入浓缩盐水以确保试验期间海水边界的盐度和水位恒定。当30 min内的盐水盐度和盐趾位置变化不超过1%时，认为试验已经达到稳态。砂槽装填高度为0.7 m，并在海水边界设置长1.2 m，坡度为1∶4的海滩斜坡。本文共进行5组试验工况，工况1为均质情况;以砂槽左下角为原点，工况2—5中2条裂隙的底部中点坐标分别为（1.0 m，0.3 m）和（2.0 m，0.1 m）、（1.0 m，0.2 m）和（2.0 m，0.2 m）、（1.7 m，0.3 m）和（2.0 m，0.1 m）、（1.0 m，0.3 m）和（1.5 m，0.1 m）。

1.2 数值模型设置

采用由美国地质调查局开发的SUTRA进行二维数值模拟。SUTRA是基于有限差分法和有限元法模拟变密度、变饱和度水流运动及溶质运移的程序，已广泛应用于海水入侵模拟研究［19，22，25］，所用控制方程和van Genuchten经验公式参见文献［25］。本文假设流体密度（ρ）随质量浓度（c）线性变化，变化率ρ/c=714.3 kg/m3［25］，并控制海、陆边界的水位和流量以保证水流为层流［19］。模拟的21个工况中，最大流速为0.005 7 m/s，雷诺数为75.49，远小于Etsias等［19］工作中层流的雷诺数值（256.3），验证了达西定律在本文工作中的适用性。

数值模型示意如图1（c），模型不考虑降水和蒸散发，所有设置均与室内试验一致（表1）。研究区被剖分为500×120个网格，共60 621个节点，其中裂隙被显式剖分为二维网格。离散网格的最大Péclet数Pe=1.2，满足Pe≈ΔL/αL≤4的数值稳定性判据，其中ΔL为沿地下水流动方向单元两侧的传输距离［25］。边界ABC上，对流入含水层的节点设为海水质量浓度（cs=35 kg/m3），而对流出含水层的节点设为零浓度梯度［25］。模拟时间步长为1 s、总时间为1 d，所有模拟工况均可达到稳态。

2 结果与分析

2.1 室内试验与数值模拟结果

图2展示了5组试验工况达到稳态的结果和对应的数值模拟结果，显示裂隙改变了其附近流体流动方向、加速了其内的流体流动，并在其上、下方形成低流速区。工况1为均质情况。工况2中盐水楔占领了裂隙下方的低流速区，且裂隙中的高速水流限制了盐水楔的高度，故盐水楔呈现阶梯状；同时，咸淡水在低速区流速接近，近海侧裂隙的上、下方均形成较大的混合区。工况3中2条裂隙形成1条长裂隙，促进了盐趾长度。对比工况2、工况3的结果，可以发现裂隙高度是盐水楔高度的重要影响因素。工况4中含水层中部裂隙对盐水楔形态的影响较小，但其下方低流速区使盐水突破了含水层底部裂隙的限制，在裂隙间形成溶质浓度较小的混合区。工况5中含水层底部裂隙完全处于盐水楔内部，故盐水楔特征与Etsias等［19］工作中的单裂隙情况类似。同时，图2也显示数值结果与试验结果吻合较好，仅有的细微差别可能来源于填砂均匀程度、裂隙位置等试验中难以准确控制的部分。

2.2 敏感性分析

根據均质情况的海水入侵水平和垂直范围，构建了4组数值模拟敏感性分析，共21个工况（表2），以研究水平双裂隙含水层中裂隙的高度、水平位置对海水入侵的影响。所分析的海水入侵特征包括盐趾趾长（海水50%浓度等值线与含水层底部的交点位置［19］）、盐水楔质心位置、盐水体积、混合体积（以海水10%和90%浓度等值线为界限［20］）、混合带质心位置、混合带垂直厚度和水平宽度。

2.2.1 近海侧裂隙的高度变化对海水入侵特征的影响

图3展示了水平双裂隙含水层中近海侧裂隙位于不同高度时的海水入侵，相关数值总结于表3。

如图3（a）所示，近海侧裂隙对淡水的汇聚作用随其高度上升变强，裂隙汇聚了大量淡水，减少了淡水对盐水楔的挤压，促进了盐水楔的扩张。同时，当近海侧裂隙升高时，近陆侧裂隙对淡水的汇聚随着其与淡水接触面积的减小而递减，对底层盐水的抑制作用逐步减弱。相较均质情况，A1—A5的盐趾长度分别增加了29.11%、30.61%、32.57%、36.01%和40.52%，且增量逐渐递增（图3（b）、表3）;盐水楔质心位置逐渐向右上方迁移。其中，A1—A3中质心垂直位置略低于均质情况，展现了裂隙对盐水楔高度的限制作用。此外，5组工况的盐水体积均大于均质情况。其中，A2、A3中盐水体积较小，此时因近海侧裂隙汇聚淡水而减少的盐水体积（裂隙上方）大于因其高度上升而增加的盐水体积（裂隙下方）。随后，盐水体积随近海侧裂隙高度上升递增。

结合图3和表3可以看到，A1—A5中混合体积以及混合带的厚度、宽度均大于均质情况，裂隙间产生了大面积的混合带，其混合程度随近海侧裂隙高度上升而增强。A1—A5中混合体积的变化趋势为先增大后减小;峰值于A3发生（约为均质情况的5.14倍），此时近海侧裂隙位于淡水排泄出口和近陆侧裂隙中间，在其上、下方均有较大面积的低流速区，为混合带的产生和发展提供了条件。随着近海侧裂隙的位置更高，主要混合带逐渐从其上部转移到下部，表现为混合带质心水平位置的右移。A1—A3中最大混合带垂直厚度出现在近海侧裂隙上方，随着上方混合带的变小而递减（A2—A3）;A4、A5中则转移至裂隙之间，并因高浓度盐水体积扩张而递减。A1、A2中最大混合带水平宽度出现在近陆侧裂隙下方，A3—A5中则出现在近海侧裂隙下方。混合带平均厚度变幅不大，但其平均宽度随近海侧裂隙高度上升递增。

2.2.2 近陆侧裂隙的高度变化对海水入侵特征的影响

图4展示了水平双裂隙含水层中近陆侧裂隙位于不同高度时的海水入侵，相关数值总结于表4。

近陆侧裂隙的主要影响范围是盐水楔尾部（图4（a）），故其对盐水楔形态及溶质混合程度的影响小于近海侧裂隙。海水入侵程度随着近陆侧裂隙高度上升呈先增大后减小的趋势，具体体现在盐水楔尾部形状、盐趾长度和盐水体积的变化上。B2中2条裂隙首尾相连，此时盐趾长度、盐水体积达到峰值，分别为均质情况的1.48倍和1.57倍（表4）。B1、B2中近陆侧裂隙的高度分别小于、等于近海侧裂隙，故盐水楔尾部的扩张程度随其位置更高而加剧。B3—B5中近陆侧裂隙不再与盐水楔接触，此时近海侧裂隙限制了盐水楔尾部的扩张高度。随着近陆侧裂隙的位置升高，含水层下部淡水逐渐从被近陆侧裂隙汇聚转变为被高度较低的近海侧裂隙汇聚，淡水对盐趾长度的抑制作用增加（图4（b））。此外，盐水楔质心高度对近陆侧裂隙高度变化的敏感性较小，这一方面是由于近海侧裂隙高度对盐水楔尾部垂直扩张程度的限制，另一方面是由于近海侧裂隙的位置决定了盐水楔中大部分盐水的分布。

受近海侧裂隙的限制，近陆侧裂隙对混合带质心的影响主要针对其水平位置。B1中最大混合带垂直厚度（26.08 cm）远高于其他工况（表4）。这是由于此时近陆侧裂隙靠下，其汇聚的水大部分又被近海侧裂隙汇聚，为近海侧裂隙上部咸淡水混合带的发展提供了条件。B2中2条裂隙连接，出现最大混合带水平宽度的峰值为1.631 5 m（表4），大于A1中的1.293 4 m（表3），这是由于B2中裂隙与淡水的接触面积更大;B3中近陆侧裂隙仍能吸引下方盐水，形成较宽的混合带;但随着近陆侧裂隙与盐水楔的距离更远（B4—B5），其对盐水楔的影响陡然减小，表现为最大混合带水平宽度的突降。5组工况中的混合体积、混合带平均厚度和平均宽度均高于均质情况。

2.2.3 含水层中部裂隙的水平位置变化对海水入侵特征的影响

图5展示了水平双裂隙含水层中含水层中部裂隙位于不同水平位置时的海水入侵，相关数值总结于表5。

图5和表5显示，随着含水层中部裂隙更靠近内陆，盐趾长度总体呈现增加的趋势、而盐水体积总体呈现减小的趋势。这是由于随着含水层中部裂隙更靠近内陆，其对淡水的汇聚逐步增强，排泄区淡水流速增大而对盐水楔上部的盐水形成更强的挤压和稀释作用，盐水体积从而减小。同时，含水层底部裂隙对淡水的汇聚因含水层中部裂隙汇聚更多淡水而减弱，其对盐趾长度的抑制作用减小。图5（b）中盐水楔质心的水平位置先增大后减小，对应了含水层中部裂隙对淡水的汇聚增强使盐水楔更伸向内陆（C1—C3）、又在远离盐水楔时排放淡水而对盐水楔上部形成挤压的过程（C4—C6）。盐水楔质心高度与均质情况接近，变化不大。

如图5和表5所示，混合体积、混合带厚度和宽度总体均呈现先增大后减小的变化趋势。这是由于随着含水层中部裂隙逐渐靠近盐水楔中部（C1—C4），吸引了底部裂隙中的混溶水，主混合区逐渐从含水层中部裂隙的上方和含水层底部裂隙附近转移至裂隙间，并随着含水层中部裂隙右移而扩张;随着含水层中部裂隙更靠近内陆（C5—C6），主混合区转移至含水层底部裂隙附近（图5（a））。因此，混合带质心位置先向海洋和含水层顶部变化，再向內陆和含水层底部变化（图5（b））。

2.2.4 含水层底部裂隙的水平位置变化对海水入侵特征的影响

图6展示了水平双裂隙含水层中含水层底部裂隙位于不同水平位置时的海水入侵，相关数值总结于表6。

图6（b）显示盐趾长度的变化趋势为先增大后减小，且总大于均质情况。这是由于含水层底部裂隙在逐步脱离盐水楔的过程中（D1—D4），其对盐水的汇聚和吸引的增强令盐趾长度变长;随着远离盐水楔（D5），其左端排放的淡水使盐水楔尾部被推向海洋，令盐趾长度陡降0.706 2 m（表6），显示两者的相对位置决定了裂隙对盐水楔的作用（吸引或抑制）。盐水楔质心高度变化不大，水平位置则先逐步减小然后突变增大再逐步减小，对应了底部裂隙逐渐脱离盐水楔再逐渐远离盐水楔的过程。同时，盐水体积也因含水层底部裂隙汇聚更多淡水而挤压变小（图6（b）、表6）。

结合图6和表6来看，D2中混合体积达到峰值同时盐水体积陡降，此时含水层底部裂隙右侧开始汇聚内陆淡水，又被中部裂隙下方的低流速区所影响，2条裂隙之间产生了较大的混合带;随着含水层底部裂隙远离盐水楔，其左侧开始排放大量淡水挤压盐水楔尾部，也令混合带减小。图6（b）显示混合体积、混合带厚度和宽度的变化趋势均为先增大后减小，对应于含水层底部裂隙水平位置向内陆变化时，2条裂隙对咸淡水混合的促进作用先增大后减小。混合带位置的变化趋势也表明2条裂隙之间和含水层底部裂隙附近的溶质混合先加强后减弱。

3 结 论

本文通过室内试验及数值模拟研究了不同水平双裂隙分布含水层中的海水入侵特征与机制，主要结论如下：

（1） 位于海水入侵范围的水平双裂隙能够增强盐水楔以及混合带的发展，但水平双裂隙对海水入侵的影响会彼此限制，位于近海侧和含水层中部的裂隙作用更明显。裂隙的汇聚淡水作用会在其上、下方产生低流速区，故裂隙在盐水楔附近时会增强咸淡水混合，且裂隙间的混合带尤为显著;裂隙远离盐水楔时则会促进内陆淡水向海洋排放，抑制海水入侵。

（2） 水平双裂隙含水层中近海侧裂隙高度越高，盐水楔更向内陆和含水层顶部发展，混合带逐渐从其上部转移到裂隙间且混合作用逐渐加剧，并始终高于均质情况。近陆侧裂隙的主要影响范围是盐水楔尾部，虽受到近海侧裂隙高度限制而对海水入侵影响相对较弱，但其高度变化对盐趾长度、盐水体积、混合带范围等仍有重要作用，且可决定盐水楔质心和混合带质心的水平位置。

（3） 水平双裂隙含水层中含水层中部裂隙水平位置越靠近内陆，盐趾长度更大、盐水体积减小;当其远离盐水楔时，含水层中部的混合带急剧变小。含水层底部裂隙水平位置越靠近内陆，盐趾长度更大、盐水体积减小，而当其远离盐水楔时，盐趾长度陡降;其对混合带参数的影响相对较小。

参考文献：

［1］WERNER A D，BAKKER M，POST V E A，et al.Seawater intrusion processes，investigation and management：recent advances and future challenges［J］.Advances in Water Resources，2013，51：3-26.

［2］杨蕴，吴剑锋，林锦，等.控制海水入侵的地下水多目标模拟优化管理模型［J］.水科学进展，2015，26（4）：579-588.（YANG Y，WU J F，LIN J，et al.A multi-objective simulation-optimization model constrained by the potential seawater intrusion［J］.Advances in Water Science，2015，26（4）：579-588.（in Chinese））

［3］陈飞，徐翔宇，羊艳，等.中国地下水资源演变趋势及影响因素分析［J］.水科学进展，2020，31（6）：811-819.（CHEN F，XU X Y，YANG Y，et al.Investigation on the evolution trends and influencing factors of groundwater resources in China［J］.Advances in Water Science，2020，31（6）：811-819.（in Chinese））

［4］王焰新，甘义群，邓娅敏，等.海岸带海陆交互作用过程及其生态环境效应研究进展［J］.地质科技通报，2020，39（1）：1-10.（WANG Y X，GAN Y Q，DENG Y M，et al.Land-ocean interactions and their eco-environmental effects in the coastal zone：current progress and future perspectives［J］.Bulletin of Geological Science and Technology，2020，39（1）：1-10.（in Chinese））

［5］李闯，倪金，宋庆春，等.大连典型海水入侵区地下水供水适宜性及合理利用对策［J］.节水灌溉，2021（7）：65-70.（LI C，NI J，SONG Q C，et al.Suitability and rational utilization strategy of groundwater supply in a typical seawater intrusion area of Dalian［J］.Water Saving Irrigation，2021（7）：65-70.（in Chinese））

［6］董健，曾献奎，吴吉春.不同类型海岸带海水入侵数值模拟研究进展［J］.高校地质学报，2018，24（3）：442-449.（DONG J，ZENG X K，WU J C.Advances in numerical simulation of seawater intrusion in different coastal zones［J］.Geological Journal of China Universities，2018，24（3）：442-449.（in Chinese））

［7］PARK H Y，JANG K，JU J W，et al.Hydrogeological characterization of seawater intrusion in tidally-forced coastal fractured bedrock aquifer［J］.Journal of Hydrology，2012，446/447：77-89.

［8］GIESE M，BARTHEL R.Review：saltwater intrusion in fractured crystalline bedrock［J］.Hydrogeology Journal，2021，29（7）：2313-2328.

［9］MACALLISTER D J，JACKSON M D，BUTLER A P，et al.Remote detection of saline intrusion in a coastal aquifer using borehole measurements of self-potential［J］.Water Resources Research，2018，54（3）：1669-1687.

［10］SUN N B，ZHUANG H B，WANG S Y，et al.Study on the causes of seawater intrusion in coastal tuff hill area［J］.Water Resources and Hydropower Engineering，2021，52（S1）：298-305.

［11］譚博，刘曙光，代朝猛，等.滨海地下水交互带中的胶体运移行为研究综述［J］.水科学进展，2017，28（5）：788-800.（TAN B，LIU S G，DAI C M，et al.A review of colloid transport behaviors based on coastal groundwater interaction［J］.Advances in Water Science，2017，28（5）：788-800.（in Chinese））

［12］LU C H，CHEN Y M，ZHANG C，et al.Steady-state freshwater-seawater mixing zone in stratified coastal aquifers［J］.Journal of Hydrology，2013，505：24-34.

［13］MAHMOODZADEH D，KARAMOUZ M.Seawater intrusion in heterogeneous coastal aquifers under flooding events［J］.Journal of Hydrology，2018，568：1118-1130.

［14］MEDICI G，SMERAGLIA L，TORABI A，et al.Review of modeling approaches to groundwater flow in deformed carbonate aquifers［J］.Ground Water，2021，59（3）：334-351.

［15］ZHU W W，KHIREVICH S，PATZEK T W.Impact of fracture geometry and topology on the connectivity and flow properties of stochastic fracture networks［J］.Water Resources Research，2021，57（7）：e2020WR028652.

［16］邓亚平，张烨，施小清，等.非均质裂隙介质中重非水相流体运移［J］.水科学进展，2015，26（5）：722-730.（DENG Y P，ZHANG Y，SHI X Q，et al.Study on the migration of dense non-aqueous phase liquids in heterogeneous fractured media［J］.Advances in Water Science，2015，26（5）：722-730.（in Chinese））

［17］罗明明，季怀松.岩溶管道与裂隙介质间溶质暂态存储机制［J］.水科学进展，2022，33（1）：145-152.（LUO M M，JI H S.Mechanism of solute transient storage between Karst conduit and fissures［J］.Advances in Water Science，2022，33（1）：145-152.（in Chinese））

［18］成建梅，罗一鸣.岩溶多重介质地下水模拟技术及应用进展［J］.地质科技通报，2022，41（5）：220-229.（CHENG J M，LUO Y M.Simulation technology and application progress of Karst multi-media groundwater［J］.Bulletin of Geological Science and Technology，2022，41（5）：220-229.（in Chinese））

［19］ETSIAS G，HAMILL G A，CAMPBELL D，et al.Laboratory and numerical investigation of saline intrusion in fractured coastal aquifers［J］.Advances in Water Resources，2021，149：103866.

［20］SEBBEN M L，WERNER A D，GRAF T.Seawater intrusion in fractured coastal aquifers：a preliminary numerical investigation using a fractured Henry problem［J］.Advances in Water Resources，2015，85：93-108.

［21］KOOHBOR B，FAHS M，ATAIE-ASHTIANI B，et al.Uncertainty analysis for seawater intrusion in fractured coastal aquifers：effects of fracture location，aperture，density and hydrodynamic parameters［J］.Journal of Hydrology，2019，571：159-177.

［22］ETSIAS G，HAMILL G A，THOMSON C，et al.Laboratory and numerical study of saltwater upconing in fractured coastal aquifers［J］.Water，2021，13（23）：3331.

［23］宋明春，林少一，楊立强，等.胶东金矿成矿模式［J］.矿床地质，2020，39（2）：215-236.（SONG M C，LIN S Y，YANG L Q，et al.Metallogenic model of Jiaodong Peninsula gold deposits［J］.Mineral Deposits，2020，39（2）：215-236.（in Chinese））

［24］MATTHI S K，BELAYNEH M.Fluid flow partitioning between fractures and a permeable rock matrix［J］.Geophysical Research Letters，2004，31（7）：L07602.

［25］SHEN C J，FAN Y，WANG X B，et al.Effects of land reclamation on a subterranean estuary［J］.Water Resources Research，2022，58（6）：e2022WR032164.

Abstract：Exploring the seawater intrusion mechanism in coastal aquifer with horizontal double fractures is critical for understanding the seawater intrusion mechanism in coastal region.Sand-tank experiments and numerical simulation based on discrete fracture matrix model were performed.The effects of horizontal double fractures with different spatial distribution on the seawater intrusion were analyzed in terms of the saltwater wedge and mixing zone.Results showed that：① The development of saltwater wedge and mixing zone can be enhanced significantly by the horizontal double fractures.The seawater-freshwater mixing between the fractures is more pronounced.The individual effects of the two fractures are limited to each other.The effects are more significant for the fracture that is horizontally close to the offshore and vertically more close to the middle of the aquifer.② A higher offshore-side fracture can exacerbate the development of the saltwater wedge and the mixing between fractures.However，the height of land-side fracture has a significant influence only on the tail of the saltwater wedge.③ When the upper or lower fracture is horizontally closer to the inland，the saltwater wedge becomes wider with the total volume being shrank.However，once they are distributed so far to the inland that being out of the range of the saltwater wedge，the mixing zone area and the width of the saltwater wedge can drop sharply.

Key words：seawater intrusion;groundwater dynamics;horizontal double fractures;laboratory experiment;numerical simulation