高放废物深地质处置地下水流数值模拟方法研究进展

李露露，周志超，邵景力，崔亚莉，赵敬波

（1.中国地质大学(北京)水资源与环境学院，北京 100083；2.核工业北京地质研究院，北京 100029）

核能已成为人类社会中必不可少的能源之一。然而，核科技与核能技术的发展伴随产生了一定的放射性废物，特别是高水平放射性废物（简称高放废物），具有放射性极强、毒性极大、持续时间长等特点[1]，一旦进入人类生存环境，危害极大且难以消除，因此，如何对其安全处置仍然是世界性难题。目前国际上公认的最主要的处置方式是深地质处置，是指将固体形式的高放废物埋藏在距地表500～1 000 m 的地质体中，使之永久与人类生存环境隔离[2]。然而，在放射性废物最终处置过程中，有害的放射性核素迁移到人类环境中最有可能的方式就是地下水的搬运。处置场建成运转后，填埋废物体的核素释放到地下水中，并随之迁移扩散，造成地下水甚至地表水体的污染。因此，要有效地进行高放废物的深地质处置，就要了解区域深部地下水形成演化机制，利用数值模拟技术预测未来相当长时间系统内地下水的流动特征已逐渐成为高放废物处置库选址和场址评价中必不可少的研究内容。

放射性核素在地质体中的迁移机理及过程极为复杂，有对流、弥散、吸附、衰变和微生物等物理、化学反应过程，目前准确模拟核素在介质迁移的全过程仍然十分困难。以上过程中，核素随地下水流迁移是最主要、最直接的方式，也是分析和了解其他过程的基础。由于处置场岩石空隙类型、构造以及对地质和水文地质资料掌握程度的不同，建立的水文地质概念模型和数学模型亦不相同。此外，由于深部资料相对匮乏，加之基岩含水介质高度的非均质和各向异性，地下水流数值模拟结果会出现较大的不确定性，因而不确定性分析工作必不可少。现阶段已有很多可用于地下水流数值模拟的软件，如何选择适用的软件也是值得关注的问题。

本文就以上问题展开综述，以了解世界各国高放废物深地质处置地下水流数值模拟方法研究进展及其应用，可为我国高放废物处置场选址、安全评价和后续放射性核素迁移转化的预测模拟提供基础支撑。

1 深地质处置地下水流数值模拟方法

深地质处置中涉及的围岩类型大多为低渗透性基岩，常见的有花岗岩、黏土岩、凝灰岩和岩盐等。岩体中的含水空间主要有裂隙介质或裂隙-孔隙双重介质。孔隙发育较为均匀，而裂隙的隙宽、展布方向、延伸长度等均具强烈的非均质性和各向异性。因此，深入了解裂隙岩体渗流理论及过程是深地质地下水流数值模拟的基础。目前用于描述裂隙岩体渗流的模型可分为等效连续介质模型、离散裂隙网络模型、双重介质模型、等效-离散耦合模型[3]。

1.1 等效连续介质（EC）模型

EC 模型是将裂隙中的水流等效平均到整个岩体中，将裂隙岩体看作具有对称渗透张量的各向异性连续体，利用经典的、成熟的连续介质理论进行分析。实际上就是求解地下水在多孔介质中流动的控制方程：

式中：H—水头/m；

K—渗透系数张量/（m·d−1）；

µs—贮水率/m−1；

W—源汇项/d−1。

这类模型适用于大区域、长序列、裂隙发育程度较高或较均匀的地区，其优点是技术方法成熟、需要数据获取较为容易、操作性好。但等效连续介质方法存在着明显的缺陷：不能很好地刻画裂隙的特殊导水作用，同时由于裂隙的非连续特征，不是所有的裂隙岩体都可以等效成连续介质，典型单元体的大小和等效水力参数较难确定。因而适用范围有限，不适合场地和储罐区的裂隙水流的精细模拟。Long 等[4]认为应用EC 模型首先要确定裂隙岩体的等效渗透张量，指出渗透张量必须要无条件地适用于动力场相似的水流系统。Long 等[5]进一步考虑采用几何形态法确定渗透张量，探讨了裂隙间连通程度对裂隙渗透率的大小和性质的影响。Snow[6]通过假定裂隙面是无限延伸的，推导了单个裂隙、一组平行裂隙的渗透张量公式，并提出多组裂隙岩体的渗透张量可以叠加。Oda[7]将包含大量地质不连续面的岩体视为各向异性、弹性多孔介质，建立了求解应力和流体流动耦合的控制方程组，并利用现场实例验证了该模型的适用性。田开铭等[8]建立渗透张量的改进模型，并将统计测量与现场试验方法求得的渗透张量进行对照和校正，提高了渗透张量计算的准确性和可靠性。

1.2 离散裂隙网络（DFN）模型

DFN 模型是以裂隙网络中的单个裂隙为研究对象，认为基质岩块本身不透水，地下水只在裂隙中流动，将裂隙网络视为非连续体。以单个裂隙内水流基本公式为基础，利用流入和流出各裂隙交叉点流量相等的原则建立方程，通过求解方程组获得各裂隙交叉点的水头值。按照Wittke 线素模型[9]的基本原理，把裂隙交叉处作为节点，节点之间的裂隙称为线单元，每个线单元流向共同节点的流量等于0（稳定流）或等于储存量的变化量（非稳定流），因此表征裂隙岩体渗流网络中i节点单元域内地下水均衡方程式为：

式中：qj(j=1,2,···,N′)— 表征单元域内某一时刻流进流出各衔接线元的流量/（m3·d−1）；

wj(j=1,2,···,N′)—表征单元域中每个线元上的垂向补给量/（m3·d−1）；

Qi—i节点上的源汇项/（m3·d−1）；

Hfi—节点上的水头/m；

N′—线单元个数；

Si—裂隙以i点为中心的表征单元域内弹性贮水（释水）系数；

ej—裂隙的隙宽/m；

lj—线元的长度/m。

DFN 模型比较真实地描述了裂隙网络中地下水流特性，比较适合解决中小尺度（处置场地、储罐尺度）、需要精细刻画的地下水流问题。Wilson 等[10]提出了两种模拟二维裂隙网络水流的有限元技术，指出当裂隙交叉点处的水流干扰小到可以忽略不计时，可以使用线性单元程序。三角形单元程序可用于确定基质-裂隙耦合流。Wittke 等[9]提出了类似电路分析回路中的网络线素法来模拟裂隙网络。王明玉等[11]利用有限元方法建立三维多边形DFN 模型，模拟裂隙岩体非均质性和各向异性等渗流特征，进而确定裂隙岩体典型单元体及渗透张量的大小。Long 等[12]提出了三维圆盘随机DFN 模型，假定裂隙为不可渗透基质中的圆盘状不连续面，使用混合解析-数值技术计算通过该裂隙网络模型的稳定流量。Nordqvist 等[13]提出了可变隙宽的DFN 模型，用该模型研究了裂隙岩体中水流和溶质运移规律。但大多数情况下，尽管在处置场区、储罐区投入大量的勘探、测量、试验等相关工作，但很难准确获得每条裂隙的裂隙特征（裂隙的产状、张开度、间距、迹长等）及其连通性，同时存在工作量大、耗时多的缺点。

到目前为止，复杂的三维裂隙网络渗流，尚无成熟的算法和模拟软件，通常利用随机模拟生成裂隙网络并进行渗流模拟。该方法首先要形成DFN，但由于裂隙系统的非均质各向异性，裂隙的分布和特性存在着很大的不确定性，因而发展出统计学方法生成DFN。目前常用的统计学方法是Monte-Carlo 法，通过实测地表裂隙数据的统计规律反推地下空间岩体裂隙几何参数近似解。DFN 的建立首先需要测量和收集裂隙特征空间分布数据，形成裂隙面几何参数的概率统计模型，应用Monte-Carlo 随机方法生成DFN，最终对模型构建结果进行验证，调整模型参数确定更接近实际情形的DFN[14]。Andersson 等[15]运用三维裂隙网络模型模拟了不同水力条件下的渗流过程。Dverstorp 等[16]建立了稀疏裂隙岩石中的DFN 模型，并模拟了该岩石中的示踪试验结果。何杨等[17]针对稀疏的三维主干岩体裂隙网络，以渗流区内的水量平衡原理为依据推导了三维岩体裂隙网络非稳定渗流的基本方程，并利用Monte-Carlo 模拟了岩体中的裂隙网络。宗自华等[18]基于已有地表裂隙统计数据，对北山BS03 钻孔周边岩体进行了裂隙网络建模，并分析了其连通性。黄帆等[19]基于Monte-Carlo 法，引入裂隙迹线长度和开度的相关性函数，提出改进的随机裂隙网络生成算法，验证了算法的有效性和程序的正确性，最终借助DFN渗流模型对生成的裂隙网络进行渗流分析。

1.3 双重介质模型

双重介质模型认为裂隙岩石是由大裂隙系统和岩块孔隙或细小裂隙系统共同构成的连续介质的统一体，大裂隙导水性强，主要起到地下水流的通道作用，称为主干裂隙介质。岩块孔隙或细小裂隙分布均匀储水空间大，主要起到储存和释放地下水的作用，称为孔隙介质或分枝裂隙介质。根据实际问题的特点以及模拟时的要求，又分为：（1）双重孔隙度模型。该模型认为很多情况下，岩块的渗透率远小于裂隙的渗透率，在模拟计算时可认为岩块的渗透系数为0，使得控制方程组中有两个孔隙度，而只有一个渗透系数。（2）双重渗透性模型。该模型认为若岩块中的渗透率不能忽略，则岩块系统中不仅存在分子弥散作用，而且水流速度不为0。该模型对裂隙和基质分别建立水流方程，主干裂隙建立离散裂隙网络模型，分枝裂隙连同基质块建立等效连续介质模型[20]，然后根据边界通量相同的条件将二者联系起来。控制方程为：

式中：Ha—主干裂隙水头/m；

Hr—分枝裂隙水头/m；

Ka—主干裂隙介质的渗透张量/（m·d−1）；

r—迁移系数；

c—比例系数；

t—时间/d。

显然，这类模型能较为全面地反映裂隙岩体的渗流特征，比等效连续介质模型更为合理，参数获取比裂隙网络模型更为容易，可用来解决区域尺度裂隙水流问题。Barenblatt 等[21]首先提出了“水力双重介质”的概念，认为裂隙岩石是由空隙性差、导水性强的裂隙系统和空隙性好、导水性差的岩块孔隙系统共同构成连续介质的统一体，但该模型忽略了裂隙岩体渗流的各向异性特征。Warren 等[22]对裂隙岩体渗透特性的各向异性做了新的假设和描述，但该模型只能应用于均质的正交裂隙网络。Zimmerman 等[23]针对裂隙/多孔介质的单相流体流动，开发了一种新的双孔隙度模型，以集中参数方式处理矩阵块，与Warren-Root 模型对比，节省了近90%的计算时间。杨栋等[24]根据裂隙发育规模与工程尺度的关系，将裂隙岩体看作是由离散介质和拟连续介质组成的广义双重介质岩体，提出了广义双重介质岩体水力学模型，并对有限元解法进行了较为详细的研究。虽然该模型在一定程度上刻画出优先流的现象，但丝毫没有涉及介质中导水裂隙的具体展布位置，并不能表现出裂隙介质的各向异性、不连续性等特征，模型仍建立在多孔介质渗流理论的基础上。因此模型应用的缺点有：（1）模型中溶质交换项很难确定，物质交换系数对模型精度影响太大；（2）裂隙网络不一定能等效为连续介质，适用范围受限制；（3）求解需迭代计算，较复杂。

1.4 等效-离散耦合模型

该模型是双重连续介质模型和离散裂隙网络模型的发展，利用区域分解方法将研究区域按裂隙发育的情况进行分区，不同的分区使用不同的数学模型[25]，即对于裂隙密度大的区域采用等效连续介质模型，对于裂隙密度较小的地区采用离散裂隙网络模型。该类模型的优点是更符合一般地质条件下裂隙渗流的特征，但也存在水交换量难以确定、两类模型耦合技术上的问题。Cacas 等[26]运用混合模型对法国某个花岗岩铀矿不同尺度裂隙中的地下水流过程进行了模拟。黄勇等[27]将基于区域分解法的地下水耦合模型应用于锦屏水电站坝址区三维渗流场的模拟中，该方法有效避免了计算量过大的问题。

2 模型的不确定性分析

由于高放废物深地质处置深部资料的相对匮乏，加之基岩含水介质高度的非均质和各向异性，进行深地质处置地下水流数值模拟时结果会出现较大的不确定性，因而选用合适的不确定性分析方法进行深地质处置数值模拟不确定性分析显得尤为重要。

2.1 模型的不确定性来源

地下水流数值模拟的不确定性划分为三类：

（1）模型参数的不确定性

由于水文地质资料的局限性（包括缺乏资料或参数获取方法的局限）以及参数在时空上的变异性、尺度效应等都会导致模型参数的不确定性。参数的不确定性分析是影响地下水数值模拟结果可靠性的重要因素。

（2）水文地质概念模型的不确定性

由于资料的缺乏或对水文地质条件认识的局限，建立的水文地质概念模型存在不确定性，包括对地下水流动特征、含水介质的类型、模型边界及源汇项的估计等。

（3）观测数据的不确定性

观测数据的不确定性来源非常广泛，包括观测变量随机分布特征导致的观测误差、观测变量的抽样误差、间接测量误差以及测量仪器本身的观测误差和人为记录误差等。

以上的不确定性因素，对地下水流模型的稳定性、可靠性和精度有很大影响，因而有必要通过模型的不确定分析，了解模型的不确定因素及其对模型模拟结果的影响程度。

2.2 不确定性分析方法

（1）灵敏度分析

灵敏度分析是一种研究系统各种输入变化对输出影响程度的方法。在地下水流数值模拟中，通常以特定的模拟结果为敏感指标（如某个特定点的水位、水质点迁移到特定位置的时间等），以模型的某些要素（如水文地质参数、地下水补排量、边界条件等）作为敏感因子，改变敏感因子数值，通过模型计算敏感指标的变化程度，分析不同敏感因子的灵敏度。

灵敏度分析方法主要包括全局分析法以及局部分析法。局部分析法是研究某一个敏感因子的变化对敏感指标的影响，该方法主要有：①因子变化法，将敏感因子数值增加或减少一个单位幅度，如5%；②偏差变化法，将敏感因子增加或减少一个标准偏差。全局分析法研究在不同敏感因子的共同作用下，对结果产生的影响，该方法主要有Morris 法、多元回归法以及Sobol 法等。全局分析法操作相对复杂，目前主要应用于水文模型，在地下水流数值模型灵敏度分析研究中，仍然主要以局部分析法为主。局部灵敏度系数为：

式中：n—分析组数；

yi—参数变化后的输出值；

y0—参数变化前的输出值；

Pi—变化后参数值；

P0—初始参数值。

通过灵敏度分析法可将不同敏感因子对地下水流数值模型敏感指标的影响程度进行排序，得到对模拟结果影响较大的敏感因子及其分布区域。此外，通过灵敏度分析获取更为准确的灵敏度高的敏感因子，可大大提高模型精度，而对于灵敏度不高的敏感因子，则可适当减少获取该参数的工作量，或者在不确定分析中忽略该敏感因子，达到降维的目的。Mclaren等[28]、Bodvarsson 等[29]在场地尺度上充分考虑裂隙介质的特性，分析了影响地下水运动的参数敏感性，研究了介质对地下水运动的影响和对污染物的阻滞效应。Schwartz[30]利用Toughreact/EOS7R 软件模拟了裂隙花岗岩中放射性核素迁移过程，就裂隙间距对模拟结果的影响进行了敏感性分析。Shahkarami 等[31]提出一种可以解释任意长度的衰变链以及模拟核素通过任意通道迁移的Sobol 分析方法，并将其应用于模型输出的全局敏感性分析中。彭志娟[32]分别采用等效连续和双重介质模型概化裂隙建立高放废物处置库开挖扰动区（EDZ）模型，利用Tough2 软件研究了不同概化模型对129I 和135Cs 在EDZ 中迁移的影响，并对129I 的扩散系数和分配系数进行了敏感性分析。

（2）蒙特卡罗（Monte Carlo）法

蒙特卡罗法是一种被广泛采用的分析复杂数值模型不确定性分析的方法。该方法是Warren 等[33]最早应用于地下水流研究中，将介质参数作为随机场，假定已知随机变量的概率分布和协方差函数，用伪随机数的生成技术生成成千上万组输入变量，产生介质参数随机分布的大量实现，最终获得成千上万组模型计算结果（如地下水模拟中的水头或溶质浓度）的统计值。该方法回避了随机分析中的数学困难，只要模拟的次数足够多，就可得到一个比较精准的概率分布，并且具有收敛速度与问题的维数无关、程序结构简单等优点。但也存在一定的缺陷，如计算工作量大、误差难以估计或控制等[34]。因此，目前蒙特卡罗法一般用于比较简单的模型不确定性分析或者用来验证其他方法的分析结果。

蒙特卡罗法在国内外的地下水流数值模型不确定性分析中应用十分广泛。Vaitinen 等[35]利用蒙特卡罗法在芬兰Romuvaara 场址生成了30 种不同的地下水流模型，用于确定随机连续介质类型分析的非均匀渗透张量。Joyce 等[36]应用多尺度嵌套建模技术，利用Connectflow 软件生成了10 个处置库尺度DFN 模型，分别进行了地下水流动和溶质运移数值模拟。Pandey 等[37]模拟了放射性核素129I 从废物处理设施释放并通过地质处置库迁移到生物圈的过程，利用蒙特卡罗方法进行了阻滞因子、流速以及路径长度等参数范围的不确定性分析。苏锐等[38]以一个假想的CRPGEORC 高放废物地质处置系统为例，运用Goldsim 软件对放射性核素79Se 和129I 在该系统远场中的迁移过程进行了灵敏度分析，同时利用蒙特卡罗法进行了概率统计模拟。林达[39]以某铀尾矿库为例，利用GMS软件结合Monte-Carlo 方法进行了核素U（Ⅵ）在地下水中迁移随机数值模拟，研究了渗透系数的变异性对尾矿库地域地下水流动及U（Ⅵ）迁移的影响。

3 数值模拟软件

目前用于高放废物深地质处置地下水流数值模拟方面的软件较多（表1），按照模拟方法的不同，可将软件分成4 个大类。

表1 主要软件介绍Table 1 Introduction of main software

（1）等效连续介质模型

Modflow、GMS 以及Porflow 等主要用于孔隙介质和裂隙介质等效连续模型。Belcher[40]利用Modflow-2000 对尤卡山地区进行了比较全面的地下水流预测模拟和分析。Löfman 等[41]通过3D 有限元的模拟方法评估Olkiluoto 场址地下500 m 深部岩体的渗透性，发现施工中的深钻孔使深部基岩中大量裂隙间的连通性增大。董艳辉等[42]利用Modflow 建立北山地区的地下水流模型，初步划分了地下水流动系统。王海龙[43]利用GMS 构建了北山区域地下水数值模型，对北山地区地下水流场形态进行了进一步刻画。季瑞利[44]以北山花岗岩为参考岩性、内蒙古高庙子膨润土为参考缓冲材料，使用Porflow 对放射性核素135Cs、129I 在处置单元及其围岩中的迁移行为进行了数值模拟研究。

（2）离散裂隙网络模型

Connectflow 主要用于离散裂隙网络的水流、溶质及热运移模拟。Werner 等[45]采用Connecflow 对瑞典Laxemar 场址高放废物深地质处置库建立地下水流数值模型。Joyce 等[36]考虑气候条件变化，利用Connectflow软件构建了区域—场址—处置库多尺度地下水流及溶质运移数值模型。

（3）双重介质模型

Tough2 系列软件主要应用于双重介质的水流、溶质及热运移模拟。Liu 等[46]、Bodvarsson 等[29]选用Tough2 模拟软件还原尤卡山在区域尺度上非饱和带的地下水系统流动规律。Rechard 等[47−48]通过建立溶质运移数值模型验证尤卡山地区选址的可靠性。Schwartz[30]建立双孔隙度模型，利用Toughreact 软件耦合了水流和运移过程，模拟了可能从处置场地中泄露出来的核废物迁移过程。王礼恒[49]围绕甘肃北山高放废物地质处置预选区，分别针对区域—盆地—岩体3 级尺度开展研究，利用Modflow、Tough2 软件进行3 级尺度地下水流动数值模拟。曹潇元等[50]采用Tough2-MP/EOS3 和GRACE 重力卫星等方法，建立了北山区域地下水饱和-非饱和流模型。

（4）整合模型

整合模型主要指等效连续介质、离散裂隙网络、双重介质以及连续-离散耦合模型的整合，Feflow、Hydrogeosphere 均可用于此类模型的水流、溶质及热运移模拟。Blessent 等[51]采用Hydrogeosphere 软件，对区域地下水流场进行了模拟研究，并模拟了不同情景下，裂隙对区域流场的影响以及溶质在裂隙中的运移距离。朱君等[52]应用Feflow 软件构建三维地下水流数值模型，模拟计算了丘陵山区地下水流动特征下氚的迁移规律。包敏等[53]应用Feflow 建立了熔岩玻璃体239Pu 的溶解释放和迁移模型，模拟了10 万年内溶解态239Pu 和胶体态239Pu 的污染羽分布。

4 结论与展望

在高放废物深地质处置地下水流数值模拟研究中，选取何种合适的模型与对模拟区域的认识程度以及模型需要解决的实际问题有密切的关系。等效连续介质模型适用于大区域、长序列、裂隙发育程度较高或较均匀的地区，方法成熟、所需的数据和参数易于获取，但是不能精确刻画裂隙介质中地下水的流动特征。离散裂隙网络模型适合解决处置场地、储罐尺度等需要精细刻画的地下水流问题，但需要大量裂隙及其连通性数据、相关参数等，工作量大、耗时多。双重介质模型主要用于解决区域尺度裂隙水流问题，但并不能表现出裂隙介质的各向异性、不连续性等特征，适用范围存在一定的限制。等效—离散耦合模型可以通过区域分解法对裂隙密度大的区域采用等效连续介质模型，对于裂隙密度较小的地区采用离散裂隙网络模型，更符合一般地质条件下裂隙渗流的特征，但存在交换量难以确定、模型耦合技术问题等。

由于深部构造特征等资料的匮乏、地质屏障的各向异性以及建模要考虑的时间尺度很长，空间尺度差异很大（罐区—处置场—区域水文地质单元），都会导致模拟结果存在较大的不确定性。通过灵敏度分析将不同敏感因子对模型敏感指标的影响程度进行排序，从而提高模型精度、减小参数不确定性分析的工作量。蒙特卡罗法是目前常用的一种模型不确定性分析方法，原理简单，易于实现。

高放废物深地质处置工作是一个长期艰巨的任务，目前虽已取得长足的进展，还需注重以下几方面研究：

（1）数值模型仿真性研究。进一步开展地质、水文地质、裂隙测量等相关的现场调查及监测工作，综合野外实际观测资料及地球物理、水文地质手段等获取裂隙几何特征及水力学参数，深入开展裂隙发育规律、裂隙渗透性分布特征研究，为裂隙渗流模拟提供准确的裂隙网络模型和参数，这是提高模型仿真性的基础。

（2）模型的不确定性分析研究。区域大深度地下水流数值模型存在着较大的不确定性。因而需要根据水文地质结构、参数等建模要素的统计特征和概率分布函数，开展模型的不确定分析，为核废物处置场选址和安全评价提供服从一定概率分布条件下的要素置信区间。

（3）预测分析模拟研究。预测工况模拟要结合构建场地可能发生的工况，如处置库长期演变、废物罐失效情景、极端降雨等情景，开展不同时期、不同场景下的数值模型进行模拟分析。

（4）多介质耦合模型的研究。大区域构建等效连续孔隙介质模型，而场地尺度可以建立裂隙网络模型，通过两个模型共同边界流量和水位相等将两者耦合，既可精确刻画场地裂隙流特征，又可宏观把握水文地质单元地下水流场。