基于OpenGeoSys的北京市南口地区地下水三维模型

孙　峰1　杨忠山1　黄振芳1　Olaf Kolditz2

（1.北京市水文总站　北京　100089；2.德国亥姆霍兹环境研究中心　莱比锡）

【摘　要】　20世纪90年代以来，水资源短缺和地下水污染已经成为北京主要关注的环境问题之一。近几年来，由于城市化的快速发展和供水量需求日益增加，作为城市供水和农业灌溉的主要水源地，当地的地下水含水层水位在日趋下降。因此，对当地水文地质系统的研究是水资源可持续管理的关键基础所在。本文以收集到的钻孔资料为基础，基于有限元方法的OpenGeoSys（OGS）建立了南口地区三维地下水模型。应用模型独立参数估计PEST进行模型的校准和灵敏度分析。校准后的模型结果与观测值比较一致。模型的瞬态模拟结果反映了过去9年里地下水水位下降和由于集中开采在某一地区形成的地下水漏斗。

【关键词】　地下水模型　OpenGeoSys　PEST　南口

1　引言

南口地区位于北京平原西北部（图1），总研究面积197.55km2。20世纪90年代以来，该地区的城市化发展导致用水量急剧增加。根据北京市水务局统计数据，该地区70%的饮用水源来自地下水。但是该地区1998～2008年的年平均降雨量仅为428mm，比多年（1956～2000年）平均低28%。一方面地下水的开采导致地下水位下降（图2），另一方面地下水中总硬度也在随之升高。当地政府也在采取相应的措施解决这些问题。从2003～2008年，大约有6亿m3的地表水从山西和河北等地的水库引入北京以缓解北京水资源短缺的问题。南水北调工程实施后，每年将有12亿m3的水引入北京地区［1］。2006年6月，在南口地区的南部（图1）建立了一个北京市地下水应急水源地，每年可供3亿～4亿m3的水［2］。应急水源地的运行导致该地区局部水位下降。根据北京市水文总站对该地区长期的水质监测结果，在南口地区的北部发现一块硝酸盐浓度较高的地区。因此，南口地区被选作一个典型污染区，进行地下水多因素污染治理和水资源管理的研究。为了更好地理解该地区的水文地质状况，建立一个地下水模型是很必要的。一个校准好的地下水模型可以为后续的污染物输移模型和污染治理情景模拟研究提供水文地质参数和边界条件。

2　研究区水文地质条件

南口地区在行政区域上属于北京市昌平区主管，在地理位置上属于温榆河流域。地面高程在40.00～197.00m之间。研究区属大陆季风性气候，降雨量年内变化较大，75%的降雨量集中在6～8月。多年（1960～2000年）平均月降雨量为551.9mm，平均年蒸发量为2080 mm，平均温度为11.6℃。

研究区西面和北面均与山脉相邻，北部为军都山山脉，西面为西山山脉，这些山脉控制着南口地区第四纪沉积物和岩相特征的形成。该地区第四纪沉积层最后可达600 m，位于研究区东南部，最浅仅有15m，位于西北部。从西北部到东南部，沉积物由砾石逐渐过渡为细砂，水位也变得越来越浅；含水层也由单层潜水过渡到多层非承压－承压水［图3（a）］。

研究区的地下水排泄最主要的是抽水，其次是南部边界地下水出流。东沙河由于十三陵水库蓄水很少泄流，因此模型中认为东部边界没有水量交换。模型的西北部边界由平原与山区的分界线定义。模型上边界地表高程数据由地形图插值获得，下边界根据第四纪沉积层厚度和钻孔资料确定。

3　岩性概化

研究区共有38个钻孔资料可以利用，大部分井以抽水为目的。钻孔资料揭示，研究区地层情况较为复杂，沉积物由近20种岩性类型构成，见表1。基于颗粒尺寸和水文地质特征，可以把岩性概化为4组，即砾石层、砂层、黏土层和岩石层。沿着地质沉积顺序，从底层开始，对不同土层之间的联系进行数字化。用以上概化和数字化的钻孔数据，利用GMS软件建立三维水文地质模型。它共有6个黏土层和5个砂砾石层。图3（b）展示了南口地区实体模型中的三个剖面。

4　三维水文地质结构模型

本研究中利用地下水模型软件GMS 6.0构建三维水文地质结构模型，GMS（Groundw⁃ater Modeling System）是一个由美国Aquaveo公司开发的地下水数值模拟软件的集成系统。在三维水文地质结构模型中，顶层根据地形等高线获得，底层由钻孔深度插值获取。在概化钻孔柱状图中，各种岩性沉积物分布在各层中，软件难以实现自动识别，因此各钻孔的沉积物顺序号人为分析给定。在每个钻孔分布配置一个水平序号后，三维水文地质模型利用GMS自动生成，并在OGS［3］中开发接口用以输入GMS三维结构模型。在本研究中，GMS的三维单元有三种类型，分别为四面体、金字塔锥形体和三棱柱。但是由于OGS软件目前仅识别四面体和三棱柱，因此当三维网格文件输入到OGS时，金字塔锥形体单元被转变成两个四面体。

5　地下水数值模型

5.1　控制方程

数学概念模型可以用多孔介质地下水水流方程表示，如下式［4］。

W——单位体积的流量，L/T；

Sy——孔隙介质的特定储水能力，L/T；

t——时间，T。

在模型中，W是空间和时间的函数，主要由补给量和排泄量决定。在上述给定初始条件和边界条件的方程中，每一步都可由OGS源代码解决。

5.2　空间和时间的离散

整个研究区域基于实体模型利用GMS空间离散成三种单元类型，即四面体、三棱柱和金字塔锥形体。金字塔锥形体单元的非多项式形状函数在一些应用中，经常产生一种不熟悉的数值属性。既然所有的GMS单元网格都是定义好的，把每个金字塔锥形体分成两个四面体时避免了网格总结点数的变化。网格化之后，整个研究区分成73324个单元格，每一水平层包括4597个结点，单元格密度约为200 m。

在本研究中，主要模拟1998～2007年期间的地下水水流情况。在假定1998年地下水水流为稳定流的条件下，进行参数识别。在稳态下，识别的参数主要为三种不同介质在水平方向和垂直方向的水力学传导系数。然后利用得到的参数进行1999～2007年的非稳定流模拟。根据获取的资料情况，非稳定流模拟应力期可以定义为1年。从2006年6月起，马池口应急水源地开始使用，作为方程的一个输出项，可以得到更可靠的校准结果。

采用固定时间步长进行时间离散，根据以前研究者定义的数值稳定性，设置时间步长为1天。

5.3　边界条件

根据研究区内流场特征及对地层结构分析，将研究区边界定义成三种不同的边界类型。在研究区，西部和北部边界是平原区和山区的地形分界。山前补给区，大部分西部和北部边界定义为特定流量边界，少部分平行于水流的方向，模型中定义为无流量。基于该区域现状，地下水水位等值线图（图1），沿着东部边界，区域地下水流也几乎平行于模型边界，于是在模型中，东部边界处理做无流量边界。至于研究区的南部边界，为区域地下水水流流出边界，如果过量地抽取研究区内的地下水，将可能和它的现状流向相反，因此定义为通用水头边界。上部边界是模型与外界发生降水、蒸发和开采等流量交换的界面，定义为流量流入流出边界。其中降雨入渗补给系数取经验值0.25。由于潜水表面蒸发随地下水位的降低而减少，根据以前的研究，当地下水位埋深低于2.5m时，北京平原潜水蒸发量几乎为0［5］。因此在本研究中没有考虑潜水蒸发量，因为研究区内最小的地下水埋深也为4m。由于地下水井抽水数据缺乏，因此对于整个研究区只用一个平均的取水量。1998～2007年地下水取水量和降雨补给量详见图4。模型下边界多数为黏土层和石灰岩，因此设定为不透水边界。

5.4　模型的识别与验证

利用参数反演识别校验的目的是使模拟值和观测值之间的差值最小。在本研究中，不同多孔介质的水力传导系数和山前地区地下水侧向补给为校验参数。在输入OGS模型前，用PEST来定义每组参数的初始值和范围。在采用PEST演绎时，可以通过控制模型计算水头和观测水头尽可能精确地拟合来调整参数值。因此，用PEST计算出一组参数值直到数值解与观测值的方差和达到最小值。

稳态模型考虑1998年地下水监测井观测水位数据。考虑所有8个控制点，优化目标函数，加权残差平方和是55.27m，非零残差序列的均值为－0.3139m。其中最小残差值为－0.198m，位于327号监测井，如图5（a）所示。衡量拟合度，即相关系数为0.9452，一般相关系数在0.9以上都是可以接受。参数识别结果见表2。北部边界的水流通量为6.0×10－2～1.5×10－1m/d，西部边界的水流通量为1.7×10－2m/d。8眼监测井地下水水位的观测值和模拟值比较如图5（b）所示。

基于稳定流的模型识别，利用初步得到的水力传导系数和边界条件，建立1999～2007年期间的非稳定流模型。非稳定流模型的识别和稳定流模型的识别相似，仅是要考虑不同岩性介质的给水度。识别的不同岩性介质的给水度结果见表2。非稳态模型拟合相关系数是0.9716。参数反演进一步表明，模型随水平渗透系数的变化比垂向渗透系数的变化更灵敏。地下水水位增减的模拟值和监测值之间的关系如图6所示。

从地下水流模型获取整个地下水含水层1998年的水量衡算见表3。

5.5　模型不确定性分析

在地下水模型用于水资源管理和污染物输移模拟之前，应充分考虑模型的不确定性。模型的不确定性通常可以分为几类：输入数据、模型结构、模型参数和模型开发技术等［6］。本研究中，较为粗糙的模型数据是模型不确定性的主要来源。降雨对地下水的补给以及地下水的抽取在整个研究区域内的分布，在模型中假定是均匀的，这将引起一些观察井的模拟值过高或是过低。整个研究区的钻孔资料的分布也远远低于国家《供水水文地质勘察规范》（GB 50027—2001）的要求，因此水文地质结构模型与现实的水文地质结构之间也必定有一定的偏差。另外，有一些数据例如河流入渗补给以及农业灌溉退水补给在该模型中由于数据的缺失也没有考虑。因此，该模型在取得更多的历史数据以及钻孔和观测井水位观测资料的基础上可以得到进一步的改善。

6　结论与展望

本研究在南口地区建立了三维地下水模型，并应用已有的观测资料对模型进行了校准。降雨补给等水文过程和地下水开采等人类活动的交互过程导致地下水水位的变化。以上水文过程和人类活动等相关信息均在模型中得以应用，用于分析地下水流场在南口地区的变化，但并不足以识别所有的水文过程特征。因此只有通过制定更为完善的监测计划，收集更为详细的钻孔资料才能够更好地使地下水模型得以改善，减少模型的不确定性。

在本研究中，应用GMS软件建立水文地质结构模型作为地下水流场模型建立的基础。PEST软件用以进行参数估计。并将GMS和PEST与OGS耦合在一起进行地下水模拟研究。GMS、OGS、PEST三者代码之间的有效交流是该研究的主要成果之一。

研究区地下水水位的下降表明该地区的地下水资源管理在目前干旱的条件下是不可持续发展的。该三维地下水模型的建立只是该地区水资源综合管理工作的第一步，今后的水资源综合管理工作应充分考虑到气候变化、土地使用变化以及社会经济发展的影响等。关于该地区地下水污染和防治的工作也将很快在该模型的基础上进行。

参考文献

［1］　DUANW.Beijing water resources and the south to north water diversion project［J］.Canadian Journal of Civil Engineering，2005，32（1）：159－163．

［2］　谢振华，张兆吉，邢国章，等.华北平原典型城市地下水供水安全保障分析［J］．资源科学，2009，31（3）：400－405．

［3］　KOLDITZ O，BAUER S，BILKE L，et al.OpenGeoSys：an open－source initiative for numerical simula⁃tion of thermo⁃hydro⁃mechanical/chemical（THM/C）processes in porousmedia［J］．Environ Earth Sci，2012，67（2）：589－599．

［4］　FREEZE R A，CHERRY JA.Groundwater［M］.Prentice－Hall，1977．

［5］　许向科.北京平原区地下水流动数值模拟研究［D］．北京：中国地质大学，2006．

［6］　REFSGAARD JC，VAN DER SLUIJS JP，H JBERG A L，et al.Uncertainty in the environmentalmodel⁃ling process⁃a framework and guidance［J］.Environmentalmodelling＆software，2007，22（11）：1543－1556．