基于FEFLOW的广州金沙洲地区地下水流场数值模拟研究

陈小月，黄健民，卢 薇

(1.广东省广州市地质调查院，广东 广州 510440，2.广州地理研究所，广东 广州 510070)

自2007年以来，广州金沙洲陆续发生岩溶地面塌陷、地面沉降地质灾害[1-3]，据已有的相关研究，上述地质灾害的发生与地下水的异常波动密切相关，并主要发生在地下水位剧烈波动的下降及恢复期间。为了更好、更深入分析金沙洲地下水变化情况与地面塌陷、地面沉降地质灾害的关系，本文根据金沙洲地带的水文地质条件、钻孔、地下水的水位动态及抽水试验成果，应用 FEFLOW软件建立金沙洲地下水系统的数值计算模拟模型。模拟区范围以金沙洲界线为界，将整个金沙洲作为地下水系统的模拟计算范围。

FEFLOW是由德国柏林水资源规划与系统研究所(WASY)开发的一种对地下水水量及水质进行模拟的软件系统，利用它不仅可以计算出水位、溶液浓度和温度等标量数据，而且还可以模拟降水、地表水、地下水的流动与转换，计算出流速、流线及流径线等向量数据，很适合分析模拟区地质灾害与地下水流场的关系[4-6]。

1 模拟区水文地质概况

金沙洲四面环水，东临沙贝海—白沙河，西侧约1 km为水口涌，沙贝海—白沙河及水口涌均为流溪河支流，两支流在金沙洲南部汇入珠江。区域水文主要受珠江干流控制，局部受流溪河影响。在金沙洲模拟区内，其西侧低丘陵的地层为石炭系下统大赛坝组的砂岩、页岩、泥岩等，这些岩层隔水性好;东侧石炭系壶天群、石磴子组灰岩及白垩系大塱山组泥质粉砂岩与珠江白沙河及沙贝海相连，地下水与珠江水存在水力联系;南面石磴子组灰岩及大塱山组泥质粉砂岩与佛山市黄岐相接，与区外地下水存在水量交换。因此，金沙洲具有完整的补、径、排条件，形成了一个相对独立的一级水文地质单元。在自然状态下，金沙洲地下水顺势总体自北西向南东流动，最终排泄于珠江水系。

2 地下水系统概念模型

为建立模拟区地下水系统的概念模型，首先要对实际的水文地质条件加以概化，简化或忽略与系统无关的某些系统要素和状态，以便于数学描述，并建立地下水系统模拟模型。

2.1 模拟区含水层结构的概化

根据金沙洲区内的地下水系统含水介质的物质组成及水文地质特性，可将其内部结构概化为四层，由上至下分别为:

(1)弱含水层(弱透水层):由人工填土、淤泥质粘土、淤泥组成。

(2)第一含水层(第四系含水层):主要由砂砾层夹砂，中粗砂、中砂、粗砂组成。

(3)隔(弱透)水层:由残积砾质粘性土组成。

(4)第二含水层(基岩岩溶裂隙含水层):主要由壶天群和石磴子组灰岩组成。

2.2 模拟区边界条件概化

金沙洲模拟区东面紧邻珠江水系，并与模拟区地下水有水量较缓，可视为定水头边界;西面丘陵地区为金沙洲地下水含水系统的分水岭，其山脊线可视为零流量边界，即为隔水边界;在模拟区南侧，第一层含水层和第二含水层与外围地下水有水量交换，主要为地下径流侧向补给，可作为径流边界;人类工程活动抽排地下水导致的降水中心概化为模型的第四类边界。潜水含水层自由水面为系统的上边界，通过该边界的潜水与系统外发生垂向水量交换;下部以基岩中风化层底板为底部边界，按零通量边界处理。

3 地下水系统数值模拟模型

3.1 地下水系统数学模型

从金沙洲模拟区的空间上看，地下水流整体上以水平运动为主、垂向运动为辅，方向性明显，且地下水系统的输入、输出随时间和空间变化，为了准确模拟各含水层的相互影响，将模拟区内的地下水流作为非均质、各向异性介质中的三维非稳定流处理，建立地下水系统数学模型如下:

式中:Ω为渗流区域;h为含水层的水位标高(m);Kx、Ky、Kz为分别为 x、y、z方向的渗透系数(m/d);Kn为边界面法向方向的渗透系数(m/d);S为自由面以下含水层储水系数(1/m);μ为潜水含水层在潜水面上的重力给水度;ε为含水层的源汇项(1/d);p为人工开采和降水等(1/d);h0为含水层的初始水位分布(m);φ1为第一类边界上的已知函数;Γ1为渗流区域的上边界，即地下水的自由表面;Γ2为渗流区域的侧向边界;Γ4为渗流区域的下边界，即承压含水层底部的隔水边界;n→为边界面的法线方向;q(x，y，z，t)为定义为二类边界的单宽流量(m2/d·m)，流入为正，流出为负，隔水边界为0。

在所建立地下水系统概念模型中，所有边界均按照第二类边界条件来处理。

3.2 模拟区网格划分

将模拟区范围内平面(即每一个模拟层)剖分为9 350个结点，14 092个三角单元。模型垂向上共分4个模拟层，4层模拟层的总结点数为9 350×4=37 400个，总单元数为14 092×4=56 368个。金沙洲模拟区的三维网格剖分图见图1。

图1 金沙洲地下水数值模拟区网格剖分图

3.3 模拟区水文地质参数划分

金沙洲地下水数值模拟模型主要的水文地质参数为渗透系数、给水度、储水系数及弹性释水系数等。

为保证对计算精度不产生影响，并使水文地质参数分区简化，将每一个渗透系数分区同时视为给水度、储水系数分区。渗透系数的初始分区按地形地貌、沉积类型和地层岩性等进行划分，如表1、表2及图2所示。各区的参数初始值根据抽水试验、沉积类型和地层岩性、地层时代等特征进行估值，垂向渗透系数根据类似地区的计算经验估计给值。

表1 含水层给水度和释水系数分区表

表2 第四系砂层渗透系数分区表 m/d

图2 金沙洲地下水数值模拟计算参数分区图

基岩裂隙含水层考虑到断层对含水层导水性的影响，沿断层的方向生成缓冲区，加大缓冲区的渗透系数，见表3。

3.4 源汇项处理

源汇项主要考虑地下水开采量、垂向入渗补给量和蒸发量。现状条件下，模拟区地下水开采主要为工程活动抽排地下水，概化为井单元处理。垂向入渗主要为大气降水的入渗补给和地表水的入渗补给，分别取多年加权平均降雨量和河流补给系数法确定补给量，模拟过程中不考虑蒸发量。

表3 基岩岩溶裂隙含水层渗透系数分区表 m/d

3.5 模型检验与校正

为了验证所建立的数学模型和模型识别后确定的水文地质参数的可靠性，采用试估—校正法，取 K36、K37和 K723三个钻孔观测水位数据与模拟水位数据进行对比，对数值模拟模型进行检验。根据收集到的地下水观测资料，取模型识别的时段为2009年7月13日至2009年12月31日，取时间步长为1天，共计1 171个时间步长，以实际水位和计算水位的拟合情况作为调整参数的依据。从图3中可以看出，经过识别的参数计算的水位线与实测水位线趋势相同，误差都在容许范围之内。

图3 实测水位与模拟水位检验图

4 地下水流场数值模拟计算

模拟区人类工程活动抽排地下水造成地下水位波动大致分为如下五个阶段(见表 4)，可概化为 1 200 m3/d、2 300 m3/d和3 200 m3/d等三种典型情况，现以抽排水量3 200 m3/d(第四阶段)为例进行金沙洲地下水流场数值模拟计算，模拟结果见图4和图5。

表4 模拟区抽排地下水阶段划分及其抽排水量

图4 抽排水量为3 200 m3/d时第四系砂层含水层模拟区流场分布图

图5 抽排水量为3 200 m3/d时基岩岩溶裂隙含水层模拟区流场分布图

从地下水流场分布图中可以看出，模拟区地下水向施工降水中心(竖井)流动，形成以人类工程施工降水点为中心的降落漏斗地段。地面塌陷的分布范围位于横沙村—沙贝村—凤冈隔水条带以西，说明隔水条带阻挡了工程降水并接受来自于东侧珠江河水的补给，使金沙中学附近的地下水沿断层强径流带由东北侧向工程施工降水点方向流动，导致这些地段的地下水流场动态变化加剧。

野外调查发现的T9-T14六宗地面塌陷地质灾害均分布于工程施工降水点附近及其所控制地面沉降范围内(BX1)(图4、图5);地面塌陷根据地下水监测资料，在2008年8月15日至2009年5月8日时间段内，金沙洲地下水出现异常波动，总体处于快速下降过程，部分民井逐渐干涸，在此期间，金沙洲共发生地面塌陷地质灾害14宗，占塌陷总数的58.3%(图6)。因此，不论从空间还是时间角度来看，模拟区地面塌陷地质灾害的发生都和模型分析结果吻合，说明了该模型的建立是合理的。

图6 岩溶地面塌陷与地下水位动态变化关系图

5 数值模拟的结论与认识

金沙洲模拟区内断裂构造发育，断层及不同规模的节理裂隙密集，造成可溶岩发育多层裂隙和溶洞，洞隙之间连通性好，构成区域地下水的主要径流场。人类工程施工降水是引起金沙洲地下水流场变化的主要原因:抽排水期间，地下水的流动方向由四周向降水中心流动，地下水流场分布图内断层发育地段的地下水流线加长，并形成以施工抽水点为中心的降落漏斗，地下水位变化明显，易引发地面塌陷和地面沉降地质灾害。因此，基于上述模拟结果，应严格控制金沙洲区域的工程活动施工降水强度，以有效防止地面塌陷及地面沉降地质灾害的发生。

[1]黄健民，郑小战，陈小月，等.广州金沙洲广州大学附属实验学校地面沉降形成机理研究[J].热带地理.2012，32(04):338-343.

[2]郭宇，黄健民，陈建新，等.广州市白云区金沙洲地区地质灾害风险区划[J].热带地理.2013，33(6):659-665.

[3]郭宇，黄健民，周志远，等.广州市白云区金沙洲地区地质灾害现状及防治对策[J].中国地质灾害防治学报.2013，24(3):100-104.

[4]赵国红，宁立波，王现国.新郑市浅层地下水流数值模拟及评价[J].地下水.2007，29(6):43-46.

[5]郭晓东，田辉，张梅桂，等.我国地下水数值模拟软件应用进展[J].地下水.2010，32(4):5-7.

[6]张立志.基于 FEFLOW的大兴区地下水动态模拟研究[D].地质大学(北京).1-66，2009.