基于三维水动力模型的海砂开采对海洋环境的影响分析

胡 庆 唐国乐 郑川贵

（中电建路桥集团有限公司，北京 100048）

0 引言

砂石资源在我国占据极其重要的地位，它不仅为基础设施的建设提供了强大的支持，也为经济社会的发展提供了坚实的保障。由于陆地砂资源逐渐减少，加上经济的快速增长，海洋砂资源变得越来越重要并且已成为当前勘探开发的热点[1]。但在海砂开采过程中，由于开采施工而影响海洋环境也是不可忽略的一大重点，因此采用技术手段模拟海砂开采后海洋环境变化并分析对海洋环境的影响，保证合理开采利用海砂，弥补陆地矿产资源的不足，发挥海砂资源的经济效益，保证国家重大项目的海砂供应具有积极意义。

1 工程概况

海砂开采工程项目位于广东省汕尾市管辖海域碣石湾施公寮岛东南侧，根据矿区钻孔30 个组合样碎屑矿物分析结果，海砂的矿物组成主要为石英，次要矿物为长石，其中石英含量为88.24%~95.54%。砂层呈松散结构，平行层理构造或块状构造，紧密密度1700kg/m³，堆积密度1530kg/m³。

根据储量核实报告，广东省汕尾市管辖海域JH21-09 矿区海砂矿圈定回填用海砂矿体2 个，上部为“滨海—三角洲砂”、下部为“滨浅海砂”。采用地质块段法估算，上部“滨海—三角洲砂”矿体控制资源量为578.24 万m³，平均含泥量为16.03%；矿体推断资源量488.55 万m³，平均含泥量为15.57%；下部“滨浅海砂”矿体控制资源量为1387.69 万m³，平均含泥量为16.07%；矿体推断资源量为1052.99 万m³，平均含泥量为16.26%。

2 三维水动力模型建立与验证

采砂区水深约15~20m，平均采砂深度为9.08m。工程所在海域以潮流作用为主，采砂施工使区域地形发生变化，潮流动力受到一定影响，产生的悬浮泥沙随海流输运扩散，对采砂区及附近海域会产生一定影响。该文通过建立附近海域三维水动力模型来分析采砂工程对海洋水文动力环境的影响。

2.1 三维水动力模型基本方程

三维水动力数学模型是基于三维不可压缩流，采用雷诺平均化N-S 方程，结合Boussinesq 假设和和静水压力假设，考虑湍流的影响、流体的密度变化，并考虑温度平衡的因素，以更好地模拟水动力学过程[2]。三维水动力模型建立的基本方程主要包括连续方程（如公式（1）所示）和动量方程（如公式（2）~公式（5）所示）。

2.1.1 连续方程

连续方程是流体力学的基本方程之一，反映的是流体运动和流体质量分布的关系，是质量守恒定律在流体力学中的应用，具体如公式（1）所示。

2.1.2 动量方程

动量方程是流体运动的最基本运动学原理，即找出流体运动和它所受作用力之间的关系的数学表达式，依据的理论原理是牛顿的运动定律或动量定理，具体如公式（2）~公式（5）所示。

式中：x、y、z为坐标系3 个分量；D为总水深（m），D=H+η；H为平均海平面下的水深（m）；η为平均海平面起算水位（m）；u为x方向（东方向）流速（m/s）；v为y方向（北方向）流速（m/s）；w为σ坐标系流速（m/s），；f为科氏参数；AM为垂直湍流黏滞系数；AH为水平湍流黏滞系数；ρ0为参考密度，取1025kg/m³；ρ为海水密度。

2.1.3 初始条件

初始速度场、水位场均为0。

2.2 模型搭建及地形处理

2.2.1 模型搭建

模型包括内、外2 个模态。进行计算时，外模态忽略垂直结构，考虑水平对流和扩散，计算二维变量；内模态考虑垂向分层使用Sigma 坐标，模型层数与实测海流对应，分为3 层计算三维变量[3]。按稳定性要求，水流数学模型的计算步长为30s。

模型求解采用非结构网格中心网格有限体积求解，其优点为计算速度较快。模型计算域涵盖汕尾市、陆丰市海域，西至惠州双月湾，东至陆丰石角尾山，海域东边界为东经114°50′19″，西边界为东经115°54′42″，南边界至北纬22°19′43″，北至汕尾陆丰各个岸线。从外海至近海岸区域网格尺寸逐渐变小，其中近海岸海域网格尺寸为100~300m，外海网格尺度为2km~3km，采砂区海域周边3km 范围内网格疏密根据需要确定。工程海域附近网格布置较密，非工程海域相对较疏，采砂区范围内网格尺寸为30m，采砂区周边10km 范围内网格尺寸为30~100m。模型网格节点数为275276，单元数为111145，模型如图1 所示。

图1 三维潮流模型

2.3 模型验证

利用工程地区已有的4 个潮位站（C1~C4）的潮位资料和12 个流速测点（S1~S12）的流速、流向资料进行模型的验证，并将模型计算所得的潮位、流速、流向数据与当地各站记录的实测数据进行对比，各个站点数据偏差/误差见表1。

表1 各个站点数据偏差/误差

2.2.2 地形处理

进行数学模型计算时，根据项目区实测水下地形确定水深，以海平面为标准进行转换计算，并将各个高程点赋值于模型节点中[4]。

2.2.3 边界处理

模型共设一个潮位开边界，外海开边界潮位以9 个调和常数的形式给出，由中国海域潮汐预报软件China tide 计算获得，主要考虑4 个半日分潮（M、S2、N 和Kz）、4 个全日分潮（Ki、Oi、Pr 和QI）及一个长周期分潮（Sa）[5]。

由表1 可以看出，模型计算得出的数据与当地记录的实测数据基本吻合，工程海域12 个潮流点的计算流速、流向和实测值较吻合，流速值的相对误差都在8%以下，表明所建模型可以用来模拟研究工程实施造成的水动力变化情况。

3 模型计算结果与分析

3.1 现状潮流场分析

采用经过验证的潮流数学模型，分别对该工程近水海域的大潮涨急、落急的表层、中层、底层潮流场进行计算，结果如图2 所示（以中层大潮涨急为例进行展示）。

图2 现状工程海域大潮涨急流场（中层）

根据计算结果，涨潮时，海水从西边界进入碣石湾海域，并向西北向传播。在碣石湾湾口和湾内，受地形影响，潮流呈往复流。落潮时碣石湾内海水朝东南方向流向外海。在施公寮岛外围形成回流，涨潮潮流呈逆时针方向流动，落潮潮流呈顺时针方向流动。

采砂区位于施公寮岛东南侧约3.2km 处，大潮涨急流向大致在0°~45°，大潮落急流向大致在180°~225°。涨潮流流经采砂区后传播进入施公寮岛周边，落潮流流经该采砂区后流向外海。受岛屿地形影响，涨落潮在岬湾附近形成流速高值区。

工程区域水深较大，大潮涨潮平均流速在0.1~0.44m/s，大潮落潮平均流速在0.1~0.61m/s。工程海域潮流流速总体不大，大潮涨落潮平均流速在0.12m/s 左右，项目海区潮流流速整体较弱。

3.2 施工后潮流场分析

根据海域海砂资源储量分析，按最不利情况考虑，即海砂开采区被平均挖深到24.08m。通过改变该采砂区海底地形，在采砂区周边水域选取43 个代表点对施工后该工程近水海域的大潮涨急、落急的表层、中层、底层潮流场进行计算，计算结果如图3所示（以中层大潮涨急为例进行展示），对比结果见表2（以海砂开采区中层大潮涨潮为例，选取JH21-09 区块附近10 个点位进行展示）。

表2 施工前、后中层大潮涨潮平均流速流向变化对比表

图3 施工后工程海域大潮涨急流场（中层）

根据计算结果，采砂工程施工后，施工范围海域的潮流影响相对较大并且表现为随距离施工区域越远，影响越小的规律。采砂施工完成后，对施工海域南北方向的影响主要为潮流减速，对东西方向的影响主要为潮流加速，影响范围为南北向和东西向约6km，工程前、后流速变化幅度均在0.022m/s 以下。

从工程区域大范围海域流速变化来看，涨落潮时采砂工程对周边海域产生的流速变幅除了采砂区内及采砂区附近的t1、t18、t23、t27、t28、t31、t32、t36、t40 点位外，基本都在0.002m/s 及以下，流向变幅基本在1°以内。t23、t28 位于工程区内，流速变幅率在-15.9%~-21.4%，流向变化幅度在-1.1°~1.5°。t18、t32 靠近工程区，其流速变幅率不超过11%，流向改变幅度较小，不超过5°。因此可认为采砂工程对采砂区内流场的影响主要集中于采砂区内，采砂区2.5km以外海域潮流场基本无变化，对大范围海域潮流场的影响较小。

4 结论

该文通过采用数学模型对JH21-09 区块所在海域的采砂工程进行海洋环境影响模拟分析，得出如下结论：1）采砂工程施工完成后对海域水动力条件的影响主要集中在海砂开采区2.5km 范围内，在2.5km 范围外流速变化小于0.01m/s，流向变化小于1°。2）采砂工程对大范围海域流场影响较小，对红海湾开发区周边海域不会造成明显影响，不会造成海岸线坍塌等危害并且采砂工期短，结束后可恢复正常。