CN110851969A - 水文地质特征与区域地下水循环模拟方法 - Google Patents

Abstract

本发明公开了水文地质特征与区域地下水循环模拟方法，包括GIS系统、SRTM DEM数据与MODFLOW模型，GIS系统将大量区域尺度的数据进行水文地质信息数据的集成，SRTM DEM数据用于进行水文空间分析，MODFLOW模型用于建立区域地下水流动模型。该水文地质特征与区域地下水循环模拟方法，可分析得出各个模拟时期粒子运移的变化趋势为：1万年时三个区粒子运移变化都不大，主要随着地下水流动方向运动，运移长度在2～5km，综合分析不同时期4种不同情景下的3个区域内粒子运动规律，从水文地质角度出发可得以下结论：多种情景下其粒子运动变化不大，反应出其地下水流动较稳定，受外界条件干扰较小，地下水流速相对缓慢，粒子运移速度小。

Description

水文地质特征与区域地下水循环模拟方法

技术领域

本发明涉及地下水的循环数值模拟研究领域，具体是水文地质特征与区域地下水循环模拟方法。

背景技术

地下水是核素运移的载体，因此地下水流动与核素在地下复杂地质体中运移的定量评价是处置场址评价与选择的关键。对水文地质、地下水循环与核素运移评价的研究与探索是现今处置场址评价最热门的研究焦点之一。

近几十年来，已经广泛进行了数值模拟，以分析各个地区的地下水流。Kihm使用三维数值模拟预测了地下水流量与土地变形之间的关系。Blessent等人在核燃料现场离散裂缝的结晶基岩中进行了流体流动和溶质运移模型。Tatti等人基于地下水循环井建立了被污染的低渗透层的数值模型。Pétré等人使用数值模型分析了跨界Milk River含水层的区域地下水流系统。Mengistu等人使用数值模拟来评估Moab Khutsong深金矿的地下水监测网络。

为了实现对中国新疆库鲁克塔格区域危险化学品处置场址的评价，利用地下水流动数值模型模拟多种情景下的地下水流动的时空分布特征，研究和探索地质介质结构及参数的不确定性对地下水流动的影响。研究地下水流动的机理与模式，预测长时间尺度下区域地下水流动的变化规律，从而为中国新疆库鲁克塔格区域危险化学品处置场址的评价提供依据。

发明内容

本发明的目的在于提供一种水文地质特征与区域地下水循环模拟方法，以解决上述背景技术中提出的问题。

为实现上述目的，本发明提供如下技术方案：

一种水文地质特征与区域地下水循环模拟方法，包括GIS系统、SRTM DEM数据与MODFLOW模型，所述GIS系统将大量区域尺度的数据进行水文地质信息数据的集成，所述SRTM DEM数据用于进行水文空间分析，所述MODFLOW模型用于建立区域地下水流动模型。

作为本发明进一步的方案：所述GIS系统数据包括地质、地形、泉点位置与水井信息。

作为本发明再进一步的方案：所述SRTM DEM数据的水文空间分析，即数字河网的提取及流域划分，以此从流域整体对边界进行把握。

作为本发明再进一步的方案：所述MODFLOW模型中利用地下水实测水位对模型进行校正，使得区域地下水模型能够更加真实地反应真实的地下水流场。

与现有技术相比，本发明的有益效果是：

本发明通过分析新疆库鲁克塔格区域的水文地质条件，从地下水动力学的角度初步分析地下水流场特征、不同含水单元之间的水力联系、地下水与地表水之间的关系，以及地下水补给与排泄关系等方面的特征，采用地下水动力学与野外现场水文地质调查相结合的研究手段，并在此基础上进行区域地下水流动模拟研究，预测长时间尺度下区域地下水流动的变化规律，选择位于工作区东部和中部的三个花岗岩区进行粒子示踪模拟研究，使用MODPATH分别对六种情景进行模拟计算，得到了不同时间下核素运移轨迹，反映了地下水流动和核素运移的影响因素和变化规律，由以上叠加结果可分析得出各个模拟时期粒子运移的变化趋势为：1万年时三个区粒子运移变化都不大，主要随着地下水流动方向运动，运移长度在2～5km，综合分析不同时期4种不同情景下的3个区域内粒子运动规律，从水文地质角度出发可得以下结论：多种情景下其粒子运动变化不大，反应出其地下水流动较稳定，受外界条件干扰较小，地下水流速相对缓慢，粒子运移速度小。

附图说明

图1为发明的工作区地表高程图。

图2为发明中吐－哈盆地及其邻区的构造格局的图。

图3为发明中库鲁克塔格构造略图。

图4为发明中工作区DEM图。

图5为发明中工作区三维立体图。

图6为发明中工作区水文空间分析结果图。

图7为发明中模型剖分图。

图8为发明中工作区地质图。

图9为发明中水平渗透系数分区图。

图10为发明中工作区渗透系数表。

图11为发明中模型各渗透系数分区值图。

图12为发明中工作区年平均降水量统计图。

图13为发明中罗北凹地不同埋深大气降水入渗补给图。

图14为发明中各地区不同埋深潜水蒸发系数一览表。

图15为发明中各区潜水蒸发量一览表。

图16为发明中基本模型稳定流流场图。

图17为发明中实测与模拟地下水水位等值线对比图。

图18为发明中示踪粒子位置图。

图19为发明中基本情景粒子运移模拟结果图。

图20为发明中情景一渗透系数分区值图。

图21为发明中情景一粒子运移模拟结果图。

图22为发明中情景二渗透系数分区值图。

图23为发明中情景二粒子运移模拟结果的图。

图24为发明中情景三粒子运移模拟结果图。

图25为发明中情景四粒子运移模拟结果图。

图26为发明中四个情景的粒子运移模拟综合分析图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

请参阅图1～26，本发明实施例中，一种水文地质特征与区域地下水循环模拟方法，包括GIS系统、SRTM DEM数据与MODFLOW模型，GIS系统将大量区域尺度的数据进行水文地质信息数据的集成，GIS系统不仅精确地描述工作区的水文地质状况，提高数据管理和使用效率，同时为区域地下水模型的建立提供必要的数据支持，SRTM DEM数据用于进行水文空间分析，使用SRTM DEM数据用于进行水文空间分析，用于解决工作区面积大，水文地质边界较难以确定的难题，MODFLOW模型用于建立区域地下水流动模型，GIS系统数据包括地质、地形、泉点位置与水井信息，SRTM DEM数据的水文空间分析，即数字河网的提取及流域划分，以此从流域整体对边界进行把握，在水文分析的基础上，根据工作区的水文地质条件、流动状况和补排关系进行模型边界的确定，MODFLOW模型中利用地下水实测水位对模型进行校正，使得区域地下水模型能够更加真实地反应真实的地下水流场，按照地层时代和岩性对区域地下水渗透系数进行划分。同时为考虑工作区强烈的非均质性，在模型中对断层进行单独刻画，以分析断层对地下水流动的影响。

水文地质概况

工作区西部为海拔约1050m的博斯腾湖，向东则进入山区。东部北侧为觉罗塔格山区，南侧为库鲁克塔格山，中间夹杂着山间盆地。工作区东部主要为南湖戈壁，工作区中部北侧为吐鲁番盆地，南部为罗布泊湖盆。罗布泊属塔里木盆地的最低洼区，为现代干涸盐湖区，面积10350km2。其外围区，西部为广阔的塔克拉玛干大沙漠，东部为北山及阿其克谷地，南部为阿尔金山山前冲洪积扇和沙丘，北部为库鲁克塔格低山丘陵区。罗布泊北、西、东部分布有大面积的雅丹或风蚀残丘、台地。罗布泊呈一望无际的盐漠平原，坡度仅0.04％，地表呈现不同形态的盐壳，多呈现出翻耕地状、龟裂状等微地形，盐壳坚硬、锋利，起伏高度变化在10-60cm左右。工作区地表高程见图1。

工作区属于西北内陆干旱区，多年平均降水量小于60mm，且均集中在夏季。工作区西部的库米什、托克逊年均降水量分别为61.6mm、9.2mm；北部的鄯善、哈密年均降水量分别为26.3mm、56.9mm；南部的若羌、米兰地区为47.8mm、42.9mm。

工作区中南部的罗布泊，位于塔里木盆地东端的最低凹处，北、东、南为山区，属典型的大陆性干旱气候。该区气候干燥，如铁千里克年均降水量不足20mm，年均蒸发量高达2600mm，年均气温10.6℃，7月平均气温26.7℃，1月平均气温-9.5℃。区内多大风，主风向为NE。

工作区北部为吐鲁番－哈密盆地，吐鲁番一哈密盆地是天山褶皱山系东段中的一个山间盆地，是在海西期回返的复向斜基础上发生、发展起来的中新生代沉积盆地，盖层厚达5000～6000m，最大厚度达8700m。盆地南面觉罗塔格山海拔2000m左右。它们系由上古生界火山岩、火山碎屑岩和碎羼岩组成。了墩隆起将吐鲁番一哈密盆地分为东西两个第四纪沉积盆地，其西部就是吐鲁番盆地，面积约21000km2。吐－哈盆地构造格局如图2所示。

作区中西部为库鲁克塔格地区，该地区构造见图3。震旦纪和古生代盖层大多以地台型碳酸盐和碎屑岩建造为主。总厚度较大，达8，000m，最厚达14，000-15，000m；当第三纪陆相碎屑岩沉积发育时，盖层厚度则大大减少，仅100-200m。由于差异性的新构造运动，使该区地貌形态处于较年青的发展阶段，特别是在西北部形成了许多高峻陡峭的山脉。该区与北面的卡泽尔塔格丘陵地带(地槽区)和南面的塔里木盆地(地块区)，在地形上有明显的差别。

库鲁克塔格山地广泛分布着古生代及前古生代的碎屑岩、岩浆岩和碳酸岩，山前狭长地带和局部山间盆地中分布有厚度不大的第四纪松散层。第三纪碎屑岩类主要分布在罗布泊湖盆地东北部的剥蚀台地和山间洼地中，地层产状一般较为平缓，岩性主要为泥岩、砂质泥岩夹砂岩和砂砾岩。第四纪松散沉积物主要分布在湖盆地及其周边低洼地带，其中南部阿尔金山山前倾斜平原和东部阿其克谷地中分布有较厚的第四纪冲洪积砂和砂砾石层。

本工作区内主要有岩石圈断裂3条，壳断裂2条，一般断裂更多。根据其展布方向，可以大致分为三组：北北东向断裂组、近东西向断裂组和北东东向断裂组，其中又以北东东向断裂最发育。近东西向断裂组发育于阿尔金山区，属基底式断裂，多呈舒波状展布，规模较大，以压扭性为主，具阻水性，控制基岩山区和盆地的界线；北东和北西向断裂组成共轭断裂带，发育于库鲁克塔格、北山和阿其克谷地一带，其延伸较远，规模较大，具左旋压扭性，控制着罗布泊干盐湖的北界，并穿过罗由泊干盐湖，将其分为几个单元，沿断裂有泉点分布。

主要发育于库鲁克塔格和北山一带，数量众多，形态多样。其中以北山褶皱带中的规模较大，以复式褶皱为主，轴向一般呈北东东向和北东向。库鲁克塔格一带规模较小，轴向具多向性。在罗北的红土堡，有新生代背斜发育。

在库鲁克塔格和北山一带，岩石中普遍有节理裂隙发育，其密集带多分布于断裂两侧，灰岩、花岗岩等脆性岩石中节理裂隙的发育程度要明显高于软弱岩石。

运用地下水系统论的观点，把具有相对独立补给、径流、排泄的完整地下水运移空间在区域上划分为6个地下水一级系统，再根据含水层特征、补排条件、水动力特征和地下水水质特征，在这6个一级系统内划分若干个二级系统(亚区)。

大气降水在山区转变为阿尔金山和库鲁克塔格及北山的基岩裂隙水和地表水，基岩裂隙水一部分向深部运移转变为深循环水，另一部分沿沟谷或山前以泉的形式出露地表补给地表水，地表水流出山口后进入阿尔金山和库鲁克塔格及北山山前倾斜砾质平原以及山间洼地前缘沿断层出露成泉后又转变成地表水，并以侧向补给形式补给阿其克谷地、孔雀河三角洲和塔里木河三角洲冲洪积平原，然后再补给罗布泊湖积平原。在这种运移过程中，通过水面蒸发、潜水蒸发、植物蒸腾和人工开采等方式排泄，在垂向上地下水系统相对稳定由下向上运动，水交替呈混合型，中深层承压水水头相对较高，通过断层或相对隔水层补给浅层承压水，浅层承压水再补给潜水，在此过程地下水形成了动态均衡，在水平上水质也由HCO3－Ca·Mg型水转变为SO4·Cl－Ca·Na型，最终变为Cl－Na型水；在垂向上水质也由SO4－Ca·Mg转变为SO4·Cl－Ca·Na型水；最终变为Cl－Na型水，而后经过长期的蒸发浓缩形成了高矿化的富钾卤水。

区域地下水流动模式模拟分析

由于当今地下水模型技术的发展以及计算机技术的迅速提高，地下水数值模型已经成为定量研究地下水运动规律的重要工具手段，并在解决具体的水文地质问题过程中起着日益重要的作用。到目前为止，唯有通过模型，可以把不同手段方法(如物探、钻探、调查、监测、遥感等)得到的数据放进一个统一的系统，进行数据的集成与管理，分析彼此之间是否存在矛盾。如果对水文地质已有了一些预先的理解，还可以分析出这些数据和我们的理解的是否一致，模型是重要手段。利用模型进行水文地质数据综合处理的过程中，可以发现哪些水文地质信息已经足够，哪些信息还缺乏，缺乏信息的地方就需要布置勘探工作量。所以，模型也是帮助进一步部署野外调查工作量的重要手段，可以把有限的投资放在真正需要得到数据的地方。我们对地下水的流动、排泄、补给有很多假定，检验这些假定是否正确的有效方法是把假设放进模型，比较模型运行结果与可靠的实测结果是否一致。通过这个方式可以加深和协助我们对水文地质条件的认识。

地下水流动是一个复杂的地质过程，也就是说我们不能够直接观察到地下水在哪儿流，水以多快的速度流动，以及含水层水文地质参数的空间变化。我们最想知道的是地下水在地下的空间分布，地下水流动特征以及其随着时间的变化动态。为了达到这样的目标，科学工作者进行了大量的野外勘察和研究，以期掌握地下水流动在时间和空间上的分布特征与规律。而长期以来这一科学问题仍然无法到达满意的解决。其中的主要原因之一是含水层结构的非均质性及其参数的不确定性。因为人们只能根据有限的地质和水文地质勘探资料来分析和研究含水层系统。水文地质学界一直在寻求一种能够更合理描述含水层系统参数分布的理论与方法。所以利用地下水三维流动数值模型来描述含水层结构的空间变化特征与地下水流动规律是当前国际上水文地质领域的研究趋势。

本技术方案着重对工作区域地下水的流动模式进行了模拟分析，从而掌握区域地下水循环的模式和流动特征。

数值方法具有可以刻画复杂水文地质条件、人为活动条件、不规则边界条件、一系列时间变量等诸多因素影响的地下水流动规律。在20世纪70年代之后，随着计算机的发展与普及，数值方法已被成功地用于模拟地下水运动及地下水资源评价。以往的地下水流数值模拟研究与应用具有两个明显的特征：一是多数模型讨论数值算法与计算格式占很大篇幅，对边界条件的处理过于简单并通常论证不足。不少模型选用人为边界，人为边界选取是否合理，会影响到地下水的流动模式，人为边界在模拟过程中问题不大，可以通过调整局部参数使计算值与实测值达到预期设定目标，但在预测时，原来人为边界处的水头一般无法给出，只能作已知边界处理。二是平面二维模型为主，近年来使用三维模型有增加趋势，但实际应用中三维模型很少。在三维模型中垂向剖分层很少。较多模型采用理想化研究域及单元(矩形域正六面体)，视含水层为水平层状。计算单元过于理想化以及视含水层为水平层状，离实际距离较远，难以合理刻画含水层的实际空间结构变化。

本技术方案主要采用三维有限差分模型进行模拟地下水流动，从而分析区域地下水流动的时空分布特征，并对其影响因素进行了研究分析。

根据三维地下水流动水文地质概念模型，区内含水系统地下水流动的三维非稳定流数学模型可描述如下：

初始条件：

H(x,y,z,0)＝H₀(x,y,z) (2)

边界条件：

H(x,y,z,t)|_Γ1＝H(x,y,z)

H(x,y,z,t)|_Γ2-1＝Z

式中：

x,y,z为笛卡尔坐标轴；

t为时间；

H为已知水头；

Kxx,Kyy,Kzz为坐标轴方向的主渗透系数；

μs为比弹性给水度；

μ_d为重力给水度；

W为单位体积井流量，抽水时取负号；

Γ1为第一类边界；

Γ2-1为潜水面边界；

Γ2-2为零流量边界；

ε′为降水入渗补给量。

有限差分方法的基本思想是：用渗流区内的有限个离散点的集合代替连续的渗流区，在离散点上用差商近似代替微商，将微分方程及其定解条件化为未知函数在离散点上的近似值为未知量的代数方程，然后求解差分方程，进而得到微分方程的解在离散点上的近似值。

该方法根据地下水流动的连续性方程进行。按照连续性方程，流入和流出某个计算单元的水流之差的等于该单元储水量的变化。

地下水连续性方程可表示为：

∑Q_i＝SS·△h·△v (4)

式中：

Q_i为单位时间内流进或流出该计算单元的水量；

SS为含水层的贮水率；

△h为单位时间内水头的变化；

△v为计算单元的体积。

将渗流区进行剖分离散后，可以确定计算单元，根据连续性方程及达西公式：

∑Q_i＝SS·△h·△v

可以得到在行方向上由计算单元(i，j-1，k)流入单元(i，j，k)的流量，表达为：

式中：

q_i,j-1/2,k为通过格点(i，j，k)和格点(i，j-1，k)之间界面的流量(L³T^-1)；

KR_i,j-1/2,k为格点(i，j，k)和格点(i，j-1，k)之间的渗透系数(LT^-1)；

△c_i△v_k为横断面面积(L²)；

h_i,j-1/2,k为水头在格点(i，j-1，k)的值；

h_i,j,k为水头在格点(i，j，k)处的值；

△r_j-1/2为格点(i，j，k)和格点(i，j-1，k)之间的距离(L)。

通过其他5个界面的地下水流量，均可类推，例如，沿行方向的格点(i，j+1，k)至格点(i，j，k)的地下水流量为：

同理，可依次得出：

以上公式表示了通过计算单元(i，j，k)的六个界面的地下水流量，将格点间距和渗透系数合并为一个变量，成为水力传导系数：

CR_i,j-1/2,k＝KR_i,j-1/2,k△c_i△v_k/△r_j-1/2 (12)

将水力传导系数应用于公式(6)～(11)得到：

q_i,j-1/2,k＝CR_i,j-1/2,k(h_i,j-1,k-h_i,j,k) (13)

q_i,j+1/2,k＝CR_i,j+1/2,k(h_i,j+1,k-h_i,j,k) (14)

q_i-1/2,j,k＝CC_i-1/2,j,k(h_i-1,j,k-h_i,j,k) (15)

q_i+1/2,j,k＝CC_i+1/2,j,k(h_i+1,j,k-h_i,j,k) (16)

q_i,j,k-1/2＝CV_i,j,k-1/2(h_i,j,k-1-h_i,j,k) (17)

q_i,j,k+1/2＝CV_i,j,k+1/2(h_i,j,k+1-h_i,j,k) (18)

这些公式用来计算单元(i，j，k)的六个边界面流入的地下水流量，此外，还需要考虑其他各种外部源和汇对计算单元的影响，例如河流、沟渠、生产井、注水井、蒸发蒸腾等，这些源汇流入单元的量可用一个通式表达:

a_i,j,k,n＝p_i,j,k,nh_i,j,k+q_i,j,k,n (19)

式中：

a_i,j,k,n为第n项外部源对计算单元(i，j，k)的补给量(L³T^-1)；

q_i,j,k,n，p_i,j,k,n为常数，单位为(L²T^-1)，(L³T^-1)。

考虑计算单元的六个相邻的格点以及该单元中所包含的所有源汇项，可将连续性方程公式(4)表示为：

式中：

SS_i,j,k为计算单元的贮水率(L^-1)；

为水头对于时间的偏导数之差分近似表达式；

△r_i△c_j△v_k为计算单元的体积(L3)。

将公式(13)～(18)以及(19)代入(20)，得：

将水头对时间的偏导数用差商近似表示：

则所有流量项均以某一时间段的结束时间t_m为准，则(21)变为：

对以上方程进行迭代求解，开始时，对每个水头未知的计算单元赋给初始水头或者估计水头，每次迭代的结果，用于下一次的计算。

根据差分方程，可以写出方程组的矩阵形式：

[A]{h}＝[q]

式中：

[A]为水头的系数矩阵；

{h}为所求水头矩阵；

[q]为各个方程的中包含的所有常数项和已知项。

在MODFLOW中，系数矩阵和右侧项是通过各个软件包来逐步建立起来的，最后MODFLOW根据这两个矩阵，通过迭代对{h}进行求解。

Visual MODFLOW 2000软件被用于建立本工作区地下水流动和粒子运移的有限差分模型。Visual MODFLOW是实际应用中最完整、易用的三维地下水流动和污染物运移模拟软件。这个完整的集成软件将MODFLOW、MODPATH和MT3D同直观强大的图形用户界面结合在一起，可以方便地确定模拟区域大小和选择参数单位、设置模型参数和边界条件、运行模型模拟(MT3D、MODFLOW和MODPATH)、对模型进行校正以及用剖面图和平面图的形式将模型网格、输入参数和结果加以可视化显示。

由于工作区面积很大，水文地质边界较难确定，因此首先使用SRTM DEM数据进行水文空间分析，及数字河网的提取及流域划分，以此从流域整体对边界进行把握。工作区的DEM高程图如图4所示。

SRTM是航天飞机雷达地形测图任务(The Shuttle Radar Topography Mission)的缩写，此计划于2000年2月完成，其数据绝大部分位于北纬60°和南纬56°之间，覆盖80％的陆地面积。SRTM包括二种分辨率：SRTM1，覆盖范围仅包括美国大陆，其分辨率为1弧秒；SRTM3，数据覆盖全球，分辨率为3弧秒。由于SRTM DEM数据存在无数据区，也可以使用CIAT(国际热带农业中心)处理过的SRTM DEM产品，称为CGIAR-CSI SRTM3，最新发布了V3版本，有arc-formatted ASCII和GeoTIFF两种格式。本研究中使用的是从CIAT下载的GeoTIFF格式CGIAR-CSI SRTM3 DEM，以下均简称SRTM3 DEM。SRTM3的高程基准是EGM96的大地水准面，平面基准是WGS84。SRTM3精度为弧秒，即将1度的面积分成了1200×1200个小区域，每一个小区域的大小是3弧秒(90m×90m)，每一个小区域有一个数值代表该区域中心的高程。SRTM3 DEM数据可相当于1：25万比例尺地图数据使用。工作区的三维立体图如图5所示。

水文分析是DEM数据应用的一个重要方面。利用DEM生成集水流域和水流网络分析地表水流状况是重要的应用，同时，对地球表面形状的理解也具有十分重要的意义。基于DEM的地表水文分析的主要内容是利用水文分析工具提取地表水流径流模型的水流方向、汇流累积量、水流长度、河流网络(包括河流网络的分级等)以及对工作区的流域进行分割等。通过对这些基本水文因子的提取和基本水文分析，可以在DEM表面之上再现水流的流动过程，最终完成水文分析过程。

水文分析主要为分析区域上地表水流动的方向、分析确定各个流域可能存在的地表分水岭，因此重点是对河网的提取及流域的划分。

基于DEM数据提取数字河网及流域划分的基本流程如下：

(1)DEM数据预处理。在利用DEM数据提取流域的数字河网之前，首先需要对DEM数据进行预处理，以消除DEM中存在的无效数据区域。此外，由于DEM数据本身往往存在很多洼地和平地，所以，在获取有效的DEM数据以后，还必须对DEM数据进行洼地和平地的确定以及洼地的填充和平地的抬升等预处理，从而得到可直接用于河网水系提取的有效DEM数据。

(2)栅格流向确定。在得到预处理后的有效DEM数据以后，就可以对数据中各栅格的流向进行确定。流向的判断方法主要有多流向法和单流向法。单流向法简单方便而被广泛的应用，其中应用最多的是D8法。D8法假设单个网格中的水流只有8种可能的流向，即流入与之相邻的8个网格中。它用最陡坡度法来确定水流的方向，即在3×3的DEM网格上，计算中心网格与各相邻网格间的距离权落差(即网格中心点落差除以网格中心点之间的距离)，取距离权落差最大的网格为中心网格的流出网格，该方向即为中心网格的流向。

(3)水流累积矩阵确定。根据确定的水流方向数据，计算每一个栅格单元的上游给水区范围，即确定有哪些上游栅格的流向是累积指向该栅格单元的。

(4)提取水系。对流向累积栅格设置集水面积阈值，根据阈值，即可得出该地区的河流水系栅格网络图，生成数字水系。

(5)流域划分。流域盆地是由分水岭分割而成的汇水区域。它通过对水流方向数据的分析确定出所有相互连接并处于同一流域盆地的栅格。流域盆地的确定首先是要确定分析窗口边缘的出水口的位置，也就是说，在进行流域盆地的划分中，所有的流域盆地的出水口均处于分析窗口的边缘。当确定了出水口的位置之后，其流域盆地集水区的确定类同于洼地贡献区域的确定，也就是找出所有流入出水口的上游栅格的位置。

应用GIS水文空间分析后可以得到工作区的河网和流域模型，如图6所示。

区域地下水模型从整体的角度考虑，尽量选取自然的水文地质边界。假设地表流域盆地的分水岭处也是地下水分水岭。通过水文空间分析得到的流域地表水分水岭，并作为地下水分水岭，以零通量边界处理。而由于降水补给及地形的不对称会引起地下水与地表水的分水岭不一致，但是对于大区域地下水的模拟来说影响是可以接受的。

根据分析结果，模型边界概化如下：

西部博斯腾湖为一类水头边界。

东部为侧向补给边界，定义为二类边界。

北部艾丁湖定义为一类水头边界。

南部罗布泊为蒸发边界。

其它边界为流量边界。

垂向上顶部边界为潜水边界，模型剖分至地表，其直接接受大气降水的补给和蒸发排泄，所以在模型中以源汇项的形式定义为流量边界。

垂向上底部边界为二类隔水边界(零通量边界)。

工作区采用矩形六面体(上、下两平面不一定平行)剖分。

模型在平面上离散为160行和285列，每个网格的长宽均为2km(如图7)。垂向上模型剖分为三层，每一层均为恒定厚度。三层的厚度由上向下分别为潜水面下0-1000m、1000-1800m、1800-2600m，每一层都不是水平分布的。垂向上分为三层是进行三维数值模拟最简化的剖分方式，随着工作区地质调查的深入，将会不断改进，以求尽量准确地刻画含水层的空间分布。模型底部为潜水面下2600m，假设在此深度上地下水主要以水平流动为主。

整个模型共计剖分单元136800个，其中有效单元79356个，无效单元57444个。

针对大区域模型地形数据不足的问题，使用SRTM DEM数据建立地下水空间模型，高精度的地表高程刻画可以有效提高地下水模型潜水蒸发量计算准确性。

水平渗透系数分区主要根据工作区1：250万地质图进行划分，如图8所示，主要考虑占主导地位的水文地质特征，将花岗岩、页岩、碎屑岩、第四系沉积等分为不同系数分区反映在模型中，得到非均质参数模型。非均质特性的考虑更加符合实际水文地质条件，可以较好地反映区内地下水流动状况。

将其简化为3个参数分区：高渗透性能分区(K1)、中等渗透性能分区(K2)、低渗透性能分区(K3)。渗透系数参数分区是不连续的，每个分区包含模型中的若干单元，如图9所示。

根据工作区前人对于渗透系数的研究成果(见图10)，结合区域地下水模型的特点，三个分区渗透系数基本值定义(见图11)。

各分区的渗透系数值将在多情景模拟时进行放大或缩小，以分析参数的不确定性。

垂向渗透系数分区与水平渗透系数分区一致，根据工作区地质情况和经验值，在模型中垂向渗透系数与水平渗透系数比值取0.25，并在模型中进行试算确定。

工作区多年平均降水量(见图12)，根据不同地区的可进行降水量参数分区。

根据“95”国家重点科技攻关项目确定河西走廊地区(气候条件与罗布泊地区相近)的有效降水入渗率为30－50％。根据前人研究成果，罗北凹地不同埋深大气降水入渗补给量(见图13)。

根据不同的降雨入渗系数以及降水量，即可得到各地区的降雨入渗量。

计算公式：

Q_降＝A·β

式中：

Q_降为大气降水入渗补给量(10⁴m³/a)；

A为年有效降雨量；

β为有效降水入渗系数。

各地区的有效降水入渗量可由模型识别确定。

工作区不同地区潜水蒸发系数见表(见图14)，不同地区潜水蒸发量(见图15)。

关于潜水蒸发规律公式，我们采用的是经验公式(指数公式)：

ε＝ε₀·e^-bD

式中：

E为潜水蒸发强度；

ε₀为水面蒸发强度；

D为潜水埋深；

b为经验系数(与土质有关)，由模型识别确定。

由于工作区范围大，又处于西北偏远区，水文地质工作区程度较低，因此缺乏进行地下水非稳定流模拟所需的动态资料，如长期水位观测资料、逐月降水量、逐月蒸发量等资料，因此在研究中采用多年平均资料进行稳定流模拟，重点关注工作区区域水循环特征。

使用多年平均数据，渗透系数取基本组合进行稳定流数值模拟，模拟结果如图16所示。

由于区域地下水流场可以看出区内存在多个地下水流动系统。在工作区北部的吐鲁番盆地附近，地下水主要从觉罗塔格山及盆地东侧向艾丁湖汇流；觉罗塔格山西南侧和库鲁克塔格西北侧地下水向博斯腾湖汇集；库鲁克塔格南侧地下水主要沿孔雀河向罗布泊流动；工作区东北部地下水从西、南向哈密盆地的沙湖流动；工作区东部接受疏勒河下游的侧向补给并向罗布泊流动。

由于缺少工作区水位观测资料，因此使用区内若干实测地下水数据插值形成的地下水水头等直线与模拟初步流场进行对比，对其差异进行分析再根据不同地点的水文地质条件对模型中的参数进行调整。校正后的模拟流场与实测值插值流场对比见图17。

由图17可以看出模拟流场与实测流场整体上是一致的，但在一些局部存在一定的差异，这主要是有以下几点原因：

(1)实测地下水水位数据过少，造成插值流场局部失真，不能完全代表真实的地下水流场。工作区面积超过十万平方千米，但是用于插值的地下水数据只有几十个，特别是在工作区东部，插值点更少。以东水泉北部为例，在实测插值流场中此处有一个高水位点，而实际上此处并不是周围标高的最高点，也没有其它水源补给，只是由于附近没有其它插值点造成的。

(2)工作区面积大，但相关地质、水文地质信息较少，根据水量信息建立的模型不能完全模拟出区域内局部和细节的地下水流动特征。工作区基础资料缺乏，对于地质结构的刻画以及降水、蒸发等源汇项的处理不能达到太高的程度，因此缺乏对于山间小盆地等局部地区的小循环刻画的能力。

虽然在局部细节和精度上有所欠缺，但是该模型从区域整体上初步刻画了地下水流动的特征，并不影响从区域整体角度分析地下水流动特征，可以为区域地下水循环的分析提供支持。

为分析工作区内重点区域粒子运移状况，使用MODPATH软件进行多情景条件下的粒子示踪模拟，充分考虑参数的不确定性。

MODPATH是一个对流传输模拟软件，它可以从MODFLOW中获取地下水模型和模拟结果数据，并计算质点运移轨迹。MODPATH使用一个半解析解的质点追踪方法来计算地下水运移路径和流动时间。该方法假定每个网格单元的各方向的速度分量在它们自己的坐标方向上呈线性变化。该假定使得我们可得到一个描述网格单元内的流线的解析表达式。已知单元内质点的任意初始位置，那么该单元内流线上的其它点的坐标及地下水质点经过这两点之间的所需的时间间隔都可以计算出来。

通过MODPATH的图形交互界面，用户可以使用鼠标点击来设置质点并进行质点运移的演示。大多数质点运移模拟软件都需要进行后处理才能得到路径和时间的可视化结果，MODPATH则可以同时进行计算和显示。MODPATH还对选定的任意模型层和时间步长提供多种数据显示功能，如同时显示水头等值线、降深等值线及流速与方向。

假定核素为理想质点，即暂不考虑核素—水—岩的反应，则可使用该程序模块，便捷地模拟核素粒子在不同情景下的运移轨迹。

根据研究地质图，选择位于工作区东部和中部的三个较大面积的花岗岩区进行粒子示踪模拟，示踪粒子位置见图18。在每个花岗岩区的中心区域，在以5km为半径的圆周上放置20个粒子。模拟中的粒子使用绿色表示，运移轨迹使用红色表示。

本情景使用渗透系数和降雨量参数组合，作为基本情景进行示踪粒子的运移模拟。进行了1万年时间的模拟，粒子示踪的结果见图19。

情景一

情景一考虑将花岗岩的渗透系数增大一倍时分析粒子运移的变化趋势，使用的渗透系数分区值见图20，分别进行了1万年时间的模拟，粒子示踪的结果见图21。

情景二

情景二考虑将花岗岩的渗透系数减为一半时模拟分析粒子运移的变化趋势，使用的渗透系数分区值见图22，分别进行了1万年时间的模拟，粒子示踪的结果见图23。

情景三

情景三考虑将降水量增加一倍时分析粒子运移变化趋势，分别进行了1万年时间的模拟，粒子示踪的结果见图24。

情景四

情景四考虑将降水量减少为一半时分析粒子运移变化趋势，分别进行了1万年时间的模拟，粒子示踪的结果见图25。

本发明的工作原理是：

使用时，通过分析新疆库鲁克塔格区域的水文地质条件，包括含水介质系统结构、边界条件，从地下水动力学的角度初步分析地下水流场特征、不同含水单元之间的水力联系、地下水与地表水之间的关系，以及地下水补给与排泄关系等方面的特征；

采用地下水动力学与野外现场水文地质调查相结合的研究手段，分析工作区区域水文地质条件，初步圈定水文地质条件相对简单，弱含水、低渗透、慢流速，且相对封闭的区域，并在此基础上进行区域地下水流动模拟研究，预测长时间尺度下区域地下水流动的变化规律。

针对大尺度模型的基础资料不足、不确定性强的问题，在模拟时运用多情景分析方法对模型在多条件和多情景下的地下水流动进行模拟分析。研究主要考虑了渗透系数、降水量变化及断层对地下水流动的影响；选择位于工作区东部和中部的三个花岗岩区进行粒子示踪模拟研究，使用MODPATH分别对六种情景进行模拟计算，得到了不同时间下核素运移轨迹，反映了地下水流动和核素运移的影响因素和变化规律。

为了对多个情景下的粒子运移状况进行综合分析评价，分别进行1万年运移时间将4个情景的粒子运移轨迹进行叠加，粒子示踪的结果见图26。

由以上叠加结果可分析得出各个模拟时期粒子运移的变化趋势为：1万年时三个区粒子运移变化都不大，主要随着地下水流动方向运动，运移长度在2～5km；

综合分析不同时期4种不同情景下的3个区域内粒子运动规律，从水文地质角度出发可得以下结论：多种情景下其粒子运动变化不大，反应出其地下水流动较稳定，受外界条件干扰较小，地下水流速相对缓慢，粒子运移速度小。

尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种水文地质特征与区域地下水循环模拟方法，其特征在于：包括GIS系统、SRTMDEM数据与MODFLOW模型，所述GIS系统将大量区域尺度的数据进行水文地质信息数据的集成，所述SRTM DEM数据用于进行水文空间分析，所述MODFLOW模型用于建立区域地下水流动模型。

2.根据权利要求1所述的水文地质特征与区域地下水循环模拟方法，其特征在于：所述GIS系统数据包括地质、地形、泉点位置与水井信息。

3.根据权利要求1所述的水文地质特征与区域地下水循环模拟方法，其特征在于：所述SRTM DEM数据的水文空间分析，即数字河网的提取及流域划分。

4.根据权利要求1所述的水文地质特征与区域地下水循环模拟方法，其特征在于：所述MODFLOW模型中利用地下水实测水位对模型进行校正。

Images