高放废物处置库北山预选区区域水文地质概念模型

王海龙，郭永海

（核工业北京地质研究院，中核集团高放废物地质处置评价技术重点实验室，北京100029）

地下水是核素运移的载体，因此，高放废物处置库场址的适宜性、安全性在很大程度上取决于水文地质条件，取决于地下水流的状态［1］。而场址水文地质评价不仅需要把握预选区区域水文地质特征，更需要把握未来地下水的演化规律和进程。如何才能正确认识区域地下水的循环交替规律和流场特征？怎样才能预测地下水的未来变化？答案是必须依赖地下水流数值模拟研究［2］。根据我国北山预选区目前的选址研究阶段和已经掌握的数据资料情况，拟对该预选区开展区域地下水流数值模拟研究。主要介绍区域水文地质概念模型的建立。通过对研究区地质、水文地质资料的分析研究，利用GMS软件的地质建模模块，建立了北山区域三维水文地质结构模型。结合边界条件、流场特征、地下水均衡项和水文地质参数的分析和描述，建立了北山地区区域水文地质概念模型，为地下水数值模拟模型的建立提供依据［3-4］。

1 模拟范围确定

从地下水系统的角度考虑，北山区域地下水流数值模型模拟范围应该包括完整的北山流域。由于地形高程差异是北山区域地下水流动的主要驱动力，地下水水位与地形之间存在明显的相关性，因此，可以认为地表分水岭便是地下水的分水岭。运用ArcGIS软件的水文分析模块，利用地表数字高程模型（DEM）数据对北山地区进行水系提取和流域划分，从而确定区域地下水流数值模型的范围［5］。

1.1 利用DEM提取数字河网

1.1.1 DEM数据的获取

DEM（地表数字高程模型）数据可利用全站仪或激光测距仪等工具测量获取，也可通过航空或者遥感影像及地形图进行提取。目前可以免费获取的精度较高的DEM数据是NASA（美国国家航空航天局）发布的SRTM3（Shuttle Radar Topography Mission）DEM 的“Finished”版本数据，分辨率为3弧秒（相当于90 m）。该版本数据在原始DEM数据基础上，剔除了与周围高程相差达100 m的数据异常点，将海平面高程视为0 m，并视湖面高程为一个常数，对于无数据区域进行插值填补，使得数据精度大大提高（垂直精度小于16 m，水平精度小于20 m），数据下载地址为http://seamless.usgs.gov。

1.1.2 流向确定

流向确定是DEM水文分析的基础。目前应用最广泛的计算方法为D8算法，即假设水流只能从一个网格流入与之相邻的某个网格中。在选取的3×3的窗口中，求取中心网格与相邻网格的距离权落差（即各网格与中心网格的高程差除以其对应距离），然后以取值最大的网格作为流出网格，取中心网格的流向为流出方向。

1.1.3 汇流累积量的确定

汇流累积量计算的基本原理是赋予DEM栅格一个单位的水量，按照水流从高往低流的规律，以区域水流方向数据为基础计算流过每个栅格的水量数值，得到该栅格的汇流累积量。汇流累积量反映栅格汇水能力的强弱程度，其值越大则表示汇流能力越强，反之则越弱，若汇流累计量值为零，则表示无水流流向目标栅格，即栅格所处位置为该区域内的高地。

1.1.4 河网的提取

通过汇流累积量的确定，每一个栅格都被赋予一个汇水能力的特征值，它表示能够流入该栅格的所有栅格的数量。当特征值大于给定的阈值时，则认为该栅格位于水道之上，赋值为1，反之，则认为该栅格为产流区，赋值为0，然后，按有效水流方向连接各河道产生流域河网。

按照上述步骤进行北山区域水系提取，得到北山地区地表水系分布图并对流域进行划分（图1），北山流域范围北起外蒙戈壁阿尔泰山，南至河西走廊，西部以马鬃山一带为边界，东部以黑河流域的额济纳旗盆地为排泄区。

1.2 模型范围确定

由图1可见，该区季节洪水形成的沟系十分发育，沟谷洪流主流向为自SW向NE方向流动，最终汇入黑河下游平原区。该子流域西部边界是黑河流域西部边界的一部分，沿地表分水岭延伸；北部边界以最新的黑河流域边界为准，该边界位于蒙古共和国境内EW向分布的阿塔斯山脉；南部边界是河西走廊平原北缘与北山地区山前斜坡接触地带；东部边界则为黑河下游平原西缘与北山地区东部山前斜坡接触带。

图1 北山预选区地表水系分布图Fig.1 Surface water system of candidate Beishan area

从地下水系统角度来看，该区构成了一个完整的、相对独立的基岩裂隙地下水系统。西部边界主要由地下分水岭及与地下水流方向平行的零通量边界构成；北部边界为沿山脊构成的地下水分水岭构成，南部、东部边界则主要为地下水排泄边界。区内地下水主要为当地大气降水入渗补给形成，西部地形较高的地区为区域地下水补给区，广大的中部、东部和南部地区为地下水径流区，这些地区同时也接受大气降水入渗补给，黑河下游平原及河西走廊则为地下水的排泄区，具备了一个完整地下水系统应有的特征。模型的模拟区面积为8.24×104km2。

2 三维水文地质结构模型

三维水文地质结构模型能够清晰、准确地展示出整个研究区内地层及含水层系统的结构、各含水层之间的相互关系，为水文地质概念模型的建立提供直观、可靠的依据。采用GMS7.1.10（地下水模型系统）建立了研究区三维水文地质结构模型。GMS（Groundwater Modeling Systems） 是 美 国 Brigham Young University环境模型研究实验室和美国陆军排水工程实验工作站开发的一个综合性的用于地下水模拟的图形界面软件。在其所包含的模块中，Borehole Data、TIN（Triangulated Irregular Networks）和Solid实体模块是建立水文地质结构模型的主要功能模块［6］。

建立三维水文地质结构模型，需要了解研究区的地貌、地质条件，构造发育情况。研究中收集并整理了北山地区1∶20万、1∶50万区域水文地质报告、区域构造分布图、区域地质图、地表数字高程（DEM）以及钻孔数据资料等，由于研究区区域面积大，水文地质条件复杂，建模时根据模型需要并结合实际水文地质条件对资料进行如下整理和概化［7-8］。

（1）合并含水岩组，将洪积、湖积松散沉积物定义为第1含水岩组；将灰岩、页岩、砂岩等沉积岩类定义为第2含水岩组；将岩浆岩、变质岩定义为第3含水岩组。

（2）钻孔数据相对欠缺，为了更准确建立水文地质结构模型，通过利用研究区地质剖面图、地质图等资料的相关内容，解决钻孔资料不足的问题。

（3）将调整后的钻孔资料进行整理，建立适用于GMS录入的文本格式。

2.1 钻孔数据整理

本次共收集的地质、水文地质钻孔柱状资料276个，根据地质资料虚拟钻孔83个。将钻孔数据按岩性进行分组，按地层沉积顺序进行分层。

在Borehole Date模块中将钻孔数据进行处理，输入钻孔坐标、高程、岩性和分层等数据，生成borehole文件（图2）。

图2 钻孔位置分布图Fig.2 Distribution of boreholes

2.2 地表高程处理

根据需要将研究区剖分成有限个不规则的三角形，并根据钻孔高程数据插值计算出每个节点上的高程值，从而确定了研究区地表高程空间展布。

2.3 三维水文地质结构模型

GMS（地下水模型系统）提供了三维水文地质结构模型建立的插值方法，包括反距离加权法、自然临近法、Kriging法和对数法等。研究中选用自然临近法，软件将输入的钻孔数据根据不同分层进行插值计算，生成各含水岩组的底板高程线，并连接各钻孔的同一含水层底板高程点，形成最终的三维水文地质结构模型（图3）。

图3 研究区三维水文地质构造模型示意图Fig.3 Sketch map showing three-dimensional hydrogeological structure model of candidate Beishan area

在上述研究区三维水文地质结构模型基础上，结合区域水文地质条件，如地下水系统的边界条件、初始地下水流畅、含水层结构、均衡要素以及水文地质参数分布等，建立水文地质概念模型，旨在为地下水流数值模型建立提供依据。

3 水文地质概念模型

3.1 边界条件确定

模型北部边界为戈壁阿尔泰山山脉，是北山流域北部的天然分水岭，为模型零通量边界（隔水边界）；西部边界亦为北山流域的天然分水岭，为模型零通量边界 （隔水边界）；南部以河西走廊疏勒河、北石河一带为边界，为模型一类水头边界条件（给定水头边界）；东部边界分两部分，南东边界为黑河，为模型一类水头边界（给定水头边界），北东边界为北山流域排泄区，设为模型三类边界（通用水头边界）；在垂向上，将潜水含水层自由水面作为上边界；模型底部边界为隔水边界。

3.2 垂向含水层结构的概化

研究区岩性分布以花岗岩等基岩为主，局部分布灰岩、砂砾岩、页岩、变质岩以及松散岩类。为便于处理和计算，将研究区地下水系统按照类同等厚度处理方法概化为11层含水层：根据实际调查结果，第1层为潜水含水层，厚度50 m；第2层为承压含水层，厚度50 m；第3至第11层均为承压含水层，每层厚度100 m；第11含水层以下为隔水底板。

3.3 初始流场确定

模型初始流场主要采用多年平均地下水水位资料，通过绘制地下水等水位线图获取。从绘制的等水位线图可以看出，研究区地下水位在山区较高，而地形低洼处较低，地下水主流向为自西向东，南部的局部地区地下水自北向南流动。

3.4 主要均衡要素确定

3.4.1 降雨入渗补给量

北山地区地下水主要由大气降雨入渗补给，其中包括地表径流（如洪水）入渗补给。计算降雨入渗补给量的公式为：

式中：Qr—大气降水入渗量，104m3·a-1；α—大气降水入渗系数；F—接受降水入渗面积，m2；Qj—多年平均降雨量，m。

根据北山研究区各个气象站的资料，利用泰森多边形法将北山研究区不同降水量的控制面积区域进行分区，降水入渗系数值主要是根据包气带岩性，通过类比确定的。

3.4.2 地下水侧向径流量

根据区域地下水稳定流流场，可以得到流域边界上的水力梯度值，再根据岩性确定渗透系数，利用达西定律公式进行计算，便可确定侧向径流量。

3.4.3 蒸发量

蒸发是研究区地下水系统的主要排泄方式之一。蒸发量的大小与潜水位埋深和气候条件等因素相关，参考区域水文地质报告，计算时认为水位埋深大于5 m的地区潜水蒸发很小，不同地下水埋深潜水蒸发强度通过计算获得。

3.5 水文地质参数确定

水文地质参数主要包括潜水含水层的渗透系数（k）、给水度（μ），承压水含水层的渗透系数（k）及储水率（Ss）。考虑到北山地区地下水基本处于天然状态，采用稳定流模型来刻画，因此给水度和储水率在模型中不予考虑。

裂隙介质的渗透性受裂隙度、裂隙连接度、裂隙发育程度、裂隙宽度、裂隙密度和裂隙倾向等因素的影响。但是相对于大尺度模型来说，可将裂隙介质近似为多孔介质来模拟，视渗透系数主要受控于岩性。确定各区的渗透系数（k）初值主要依据在区内所做的水文地质试验结果及相关水文地质研究成果。垂向上的渗透系数分布满足指数递减方程：

式中：Kz—渗透系数，m·d-1；K0—第1含水层渗透系数，m·d-1；λ—衰减指数；Z—深度，m。

4 结 语

世界上许多有核国家都在开展高放废物处置库选址研究，其中地下水流数值模拟研究是其中的重要内容之一。一般来说，高放废物处置库选址中的地下水流数值模拟研究是分阶段分层次开展的，对应区域水文地质调查，开展区域地下水流数值模拟研究；对应场地详细勘查，则开展场址乃至岩块规模的地下水流数值模拟研究。根据我国北山预选区目前的选址研究阶段和已经掌握的数据资料情况，对该预选区开展了区域地下水流数值模拟研究。地下水流数值模拟有多种方法及其适用条件，但都需要建立研究区水文地质概念模型。正确的水文地质概念模型应综合地下水系统的基本特征，反映地下水系统整体性和运动性的特点。笔者在建立北山地区区域水文地质概念模型过程中，收集、整理和分析了大量已有的水文地质资料，通过区内地表水系的生成，把握了地下水系统的整体性，并为正确确定模型边界位置和属性提供了先决条件；在分析了所有钻孔资料水文地质信息基础上，建立了三维水文地质结构模型；同时，对区域地下水流场特征、地下水均衡要素和水文地质参数等进行了综合分析和描述，由此，建立了研究区区域水文地质概念模型。该模型的建立，将对后续的地下水流数值模拟研究构建一个良好的基础。

［1］ Barenblatt G I，et al.Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks［J］.PMM，Sov.Appl.Math.Mech.，1960，24:852-864.

［2］ Long J C S，Witherspoon P A.The relationship of the degree of interconnection to permeability in fractured networks［J］. Journal of Geophysical Research，1985，90（B4）:3087-3098.

［3］ Sagar B，Runchal A.Permeability of fractured rock:Effect of fracture size and data concretion［J］.Water Resources Research，1982，18（2）:266-274.

［4］ Hsieh P A，Neumann S P.Field determination of the three dimensional hydraulic conductivity tensors of anisotropic media:1，theory［J］.Water Resour.Res.，1985，21（11）:1655-1665.

［5］郭永海，杨天笑，刘淑芬.高放废物处置库甘肃北山预选区水文地质特征研究［J］.铀矿地质，2001，17（3）:184-189.

［6］田开铭，万 力.各向异性裂隙介质渗透性的研究与评价［M］.北京：学苑出版社，1989:3-21.

［7］王恩志，王洪涛.各向异性裂隙岩体渗透系数计算方法探讨［J］.武汉水利电力大学学报，1997，30（2）:49-53.

［8］万 力，李定方，李吉庆.三维裂隙网络的多边形单元渗流模型［J］.水利水运科学研究，1993，（4）:347-353.