Daniel J. Goode美国地质调查局宾夕法尼亚州马尔文普林斯顿大学新泽西州普林斯顿摘要本文提出了一种直接模拟地下水年龄的新方法该方法使用对流-弥散输运方程其中包含一个单位强度1的分布零级源项对应于老化速率。控制方程中的因变量是平均年龄即质量加权平均年龄。该控制方程是根据应用于概念性“年龄质量”的停留时间分布概念推导出来的。对于更一般的瞬态流动情况年龄的瞬态控制方程是根据应用于“年龄质量”的质量守恒原理推导出来的。年龄质量是水质量与其年龄的乘积并假设年龄质量在混合过程中是守恒的。边界条件包括所有无流和流入边界上的零年龄质量通量以及跨流出边界的无年龄质量弥散通量。对于瞬态流动条件必须已知年龄的初始分布。控制输运方程的解可以产生平均地下水年龄的空间分布并包括扩散、弥散、混合和交换过程而这些过程通常仅通过特定示踪剂的溶质输运模拟来考虑。传统方法依赖于平流输运来预测地下水运移时间和年龄的点值。所提出的方法保留了仅对流模型的简单性和与示踪剂无关的特性但结合了弥散和混合对体积平均年龄的影响。在两个理想化的区域含水层系统一个是均质的另一个是分层的中的年龄模拟实例证明了所提方法与传统粒子追踪方法之间的一致性并说明了使用所提方法确定扩散、弥散和混合对地下水年龄影响的用途。引言地下水系统中水的年龄对于定量分析非常有用。各种环境示踪剂可用于估计从各个地点采集的地下水样本自补给以来的时间。这些年龄数据反过来又可以用来约束流动和输运模型的参数。同位素信息通常通过使用运移时间方法来解释将地下水运动模拟为“活塞”流。Reilly et al. [1994]表明通过数值粒子追踪模拟的地下水年龄的平流模型与在浅层砂石含水层中观察到的氯氟烃和氚的分布一致。在平流模型中年龄由水的运移时间决定该运移时间根据达西定律计算然而地下水中的同位素输运通常不仅受平流影响。多项研究表明如果忽略弥散和混合对同位素浓度的影响地下水运移时间的估计可能会出现误差。Plummer et al. [1993, p. 287]对当前的测年方法进行了综述并指出正如许多其他学者所言“因为环境示踪剂是随地下水一起输运的溶解溶质通常有必要考虑流体动力弥散对模拟年龄的影响。”Walker and Cook [1991, p. 41]“展示了在补给率同样较低的非承压含水层中使用同位素碳 14 时忽略扩散如何导致严重低估地下水年龄。”Maloszewski and Zuber [1991]证明了基质扩散和交换反应对裂隙岩石中碳 14 运动及地下水年龄的影响。Geyh and Backhaus [1978]考察了抽水过程中扩散和混合对碳 14 分布和年龄的影响。从承压含水层进入相邻隔水层的碳 14 扩散以及由此产生的对解释年龄的影响由Sudicky and Frind [1981]进行了定量研究。在他们的综述中Mazor and Nativ [1992, p. 211]确定了地下水年龄解释中的“问题区域”包括“缺乏单一的补给和排放区……地下水被截留在‘死’体积中……以及不同年龄地下水的混合。”结合平流以外的输运过程的影响可以从示踪剂分布中提取更多信息。例如Torgersen et al. [1978]根据观测到的湖泊中氚/氦-3 分布估计了垂直扩散率。Weeks et al. [1982]使用氟碳化合物分布在非饱和带中估计土壤扩散系数。结合弥散的年龄模型有助于识别地下水系统的弥散特性[Robinson and Tester, 1984]以及平均流动特性。Egboka et al. [1983]通过拟合氚输运的一维模型从观测到的氚分布中估计了纵向弥散度。Musgrove and Banner [1993]使用同位素信息来帮助量化区域规模流动系统中不同咸水的混合。扩散、弥散和混合的影响可以纳入感兴趣的示踪剂或多种示踪剂的输运模拟中例如氚[Nir, 1964; Simpkins and Bradbury, 1992; Solomon et al., 1993]。氚是氢的一种同位素并入水分子的结构中。通过使用考虑水分子的弥散和扩散的模型来模拟氚的输运反映了潜在的弥散和扩散而这些水分子具有不同的年龄。因此水本身的年龄受到扩散、弥散和其他流体动力过程的影响。正如粒子追踪被用于生成仅平流情况下的地下水年龄一样直接模拟受弥散和混合影响的地下水年龄而无需诉诸单独示踪剂的单独模型也是有用的。这种分析是对感兴趣示踪剂输运模拟的补充但不能取代它示踪剂的输运可能会受到吸附和化学交换等过程的影响[Fritz et al., 1979]而这些过程不会影响水。停留时间分布的概念多年来一直被用于描述流动反应器中组分寿命的统计数据[MacMullin and Weber, 1935; Danckwerts, 1953]。选择这种方法是因为分子表现并不完全相同相反不同路径的混合、扩散和可变流动条件导致了不同分子的不同停留时间。停留时间分布已广泛应用于化学反应系统的分析[Levenspiel, 1972]。Eriksson [1958]、Bolin and Rodhe [1973]以及Nir and Lewis [1975]讨论了停留时间分布理论在稳态和瞬态地球物理系统中的应用。Robinson and Tester [1984]使用从示踪剂测试确定的停留时间分布来分析裂隙地热储层中的弥散。Campana and Simpson [1984]将其中的一些概念应用于离散状态而非连续状态室模型中的地下水同位素测年其中地下水系统的独立区域被视为混合水池。停留时间分布与理想示踪剂输运之间的内在联系可用于开发对应于平流-弥散输运方程解的停留时间分布[Danckwerts, 1953; Wen and Fan, 1975; Nauman, 1981; Zuber, 1986]。Danckwerts [1953]描述了如何通过实验从非反应性示踪剂的流出浓度中确定停留时间分布。在这种情况下一点处的平均年龄由浓度的瞬时积分决定如下所述受溶质输运方程中考虑的所有过程的影响包括弥散和混合。驻留浓度和概率运移时间方法与模拟地下水输运之间的联系由Shapiro and Cvetkovic [1988]以及Dagan and Nguyen [1989]描述他们讨论了运移时间方法对某些分析的优势。最近Harvey and Gorelick [1995]提出了一个将时间矩生成方程应用于反应性输运的通用框架。在本文中我展示了地下水年龄的分布遵循溶质输运方程的一种特殊形式。平均年龄可以直接在解析或数值输运模型中模拟模拟的结果即预测的“浓度”就是平均年龄。这种空间年龄分布是直接获得的无需在输运模拟之外进行进一步处理。对于稳态流动条件年龄输运方程是从由平流-弥散方程控制的系统中停留时间分布的先前结果推导出来的。这种形式以前曾被提出用于分析房间通风过程中的平均或“局部”空气年龄[Sandberg, 1981]。此处针对瞬态流动条件的更一般推导是基于想象的“年龄质量”守恒假设。该方法通过两个区域含水层系统的数值模拟进行了说明一个是均质的另一个是分层的。平均年龄输运方程的推导停留时间分布Danckwerts [1953]定义了一个函数来表征化学反应器中分子的停留时间分布。该函数对应于在零时刻以脉冲单位时间内的单位质量注入反应器入口的溶质在反应器出口处的浓度。对于一维活塞流仅平流条件函数除了在等于系统平流运移时间的时间点外均为零此时该函数是一个狄拉克 delta 函数。在完全混合的情况下是一个指数衰减函数流出浓度等于储层内的均匀浓度由于在入口处添加了不含示踪剂的流体该浓度随时间降低。对于多孔介质中的输运函数类似于流出质量通量每单位时间的质量除以一列水的质量即以通量测得的浓度[Kreft and Zuber, 1978]。在多维系统中注入所有流入边界的质量与跨边界的流体通量成正比。Levenspiel [1972]以及Wen and Fan [1975]展示了流动反应器的几种模型及其各自的函数。稳态域中的平均停留时间可以根据在零时刻以脉冲形式注入的示踪剂浓度确定计算公式为1其中是平均停留时间或反应器中分子的平均年龄[Spalding, 1958]。(1) 式的分子是浓度加权平均时间我将其表示为。这种形式类似于随机变量的数学期望其概率密度函数为。分母将分子归一化使得除以的时间积分具有概率密度函数的特性。这个积分在稳态流动下是恒定且均匀的[Spalding, 1958]记为。如果注入流入边界的质量与跨边界的流体通量成正比则该项对于多维系统也是均匀的[见 Harvey and Gorelick, 1995]。对于恒定密度流体中的平流和弥散输运浓度满足2其中 $\theta$ 是水分含量饱和流动的孔隙度$\mathbf{D}$ 是弥散张量$\mathbf{q}$ 是比流量向量渗透速度。地下水中弥散的标准模型假设水分含量与弥散张量 $\mathbf{D}$ 的乘积由下式给出3其中是扩散系数$\delta_{ij}$ 是克罗内克 delta如果 $ij$ 则 $\delta_{ij} 1$如果 $i \neq j$ 则 $\delta_{ij} 0$$\alpha_L$ 和 $\alpha_T$ 分别是纵向和横向弥散度$q_i$ 是比流量向量的一个分量$q$ 是比流量的大小。将 (2) 式乘以时间并对所有时间进行积分得出[Spalding, 1958]4(4) 式的左侧是通过对时间导数项进行分部积分获得的。除以 $M$它是均匀的且可以带入空间导数内并假设孔隙度也是均匀的得出5其中已代入定义的平均年龄 $A E/M$[参见 Sandberg, 1981 进行比较]。因此点处的平均年龄满足具有单位强度 (1) 的内部源项的稳态平流-弥散输运方程。 (5) 的一般形式将在下一节中通过假设年龄在混合过程中守恒来推导。该推导还为 (5) 提供了自然边界条件的方案从而完成了数学框架。通过粒子追踪计算年龄对应于不带弥散项的 (5) 式的解。在这种情况下粒子路径由控制方程的特性定义增加的运移时间对应于 (5) 式中的单位源项。方程 (1) 可用于根据粒子追踪分析计算小体积内的平均年龄以包括例如采样引起的混合。本文介绍的年龄模拟方法对于核算弥散和混合对其他涉及地下水和运移时间的分析的影响也可能很有用例如来自垃圾填埋场的污染[Lee and Kitanidis, 1993]和时间依赖性捕捉区[Bair et al., 1990]。年龄质量守恒在本节中年龄输运方程的更一般形式是根据质量守恒原理推导出来的。假设混合水的平均年龄是质量加权平均值那么平均年龄类似于保守溶质浓度。虽然年龄不是一个直接可测量的物理性质因此这一假设无法通过实验验证但对于我们的概念模型来说它似乎是合适的即当两个水质量混合时混合物的平均年龄是各组成部分的加权平均值。给定平均年龄为 $A$ 的一定质量的水可以假设其特征在于其“年龄质量”即平均年龄与水质量的乘积 $A \rho V$其中 $\rho$ 是水密度$V$ 是水体积。假设水的密度是恒定的两组分混合物的平均年龄是体积加权平均值$$6其中下标区分组分。这完全类似于保守溶质在恒定密度水中混合时的浓度。这种类比可以用来推导平均年龄的控制输运方程。年龄质量输运的控制方程可以从类似于Konikow and Grove [1977]推导质量输运方程的简单框平衡box balance中推导出来另见Bear [1979]。考虑含水层材料控制体积维度为 $\Delta x, \Delta y, \Delta z$内的年龄质量守恒。控制体积内的年龄质量是平均年龄 ($A$) 与水的质量 $\theta \rho \Delta x \Delta y \Delta z$ 的乘积。跨越控制体积边界的年龄质量通量每单位面积记为 $\mathbf{J}$分量为 $J_x, J_y, J_z$包括平流和水速以及弥散通量。在时间步长 $\Delta t$ 内最初在控制体积内的年龄质量由于 $\Delta t$ 与水质量的乘积而增加。此外还包括速率为 $F$ 的年龄质量内部净源项例如用于解释跨相的净交换。$F$ 中包含的物理过程描述如下。假设年龄质量守恒差分形式的守恒方程为7除以体积和时间步长并使控制体积的大小和时间步长在极限情况下趋于零给出平均年龄输运的控制偏微分方程8这种形式的平均年龄输运方程比 (5) 式更通用因为它允许水密度变化且方程是瞬态的。即使对于具有非稳态流动的含水层系统也可以确定平均年龄空间分布只要已知流动历史。模拟地下水年龄的瞬态演变能力可能很有用例如在评估气候变化对大型含水层系统的影响时。为了完成控制方程需要描述年龄质量通量 $\mathbf{J}$以年龄或其梯度的形式。在这里我采用了质量通量的标准模型即 $\mathbf{J}$ 由平流、扩散和建模为菲克扩散Fickian diffusion的弥散组成。也就是说(8) 中的 $\mathbf{J}$ 被替换为$$\tag{9}$$其中 $\mathbf{D}$ 是弥散张量并包括扩散项。通过代入控制方程变为$$\tag{10}$$如果孔隙度和密度在时间上是恒定的且在空间上是均匀的(10) 变为$$\tag{11}$$在稳态流动条件下(11)或 (10)中的时间导数可以设为零以推导稳态平均年龄空间分布的控制方程。稳态年龄分布不存在的唯一情况是 $\mathbf{q}$ 和 $\mathbf{D}$ 均为零。从 (11) 导出的稳态方程与从停留时间分布理论推导出的 (5) 相同带有一个额外的源项 $F$这在下文讨论。平均年龄输运方程的初始和边界条件可以从控制体积平衡中指定。边界和初始条件地下水的年龄是相对于水进入系统的时间而言的。也就是说到系统的补给被假设为年龄为零。从年龄质量的定义来看补给水的年龄质量通量也是零。此外跨越所有无流边界的年龄质量通量也为零。这些条件可以写成$$\tag{12}$$其中 $\Gamma_1$ 是无流或流入边界$\mathbf{n}$ 是其单位外法线。注意 (12) 使用了总年龄质量通量 $\mathbf{J}$ 的先前表示法其中包括弥散以及平流。因此假设弥散在边界上逆着流入流的方向发生。流出边界上的边界条件取决于物理情况。我将采用一个通用的假设即跨流出边界的质量通量仅通过平流发生。这种无弥散跨边界的条件可以写成$$\tag{13}$$其中 $\Gamma_2$ 是流出边界。这种条件对于排放到地表水体的情况最为合适。对于其他物理情况也可以制定替代的边界条件例如将平流和弥散合并到模拟域之外的单独含水层中但这里不予追求。平均年龄输运方程的一般形式是瞬态的。如果流场是稳态的则平均年龄的稳态分布由控制方程的解决定类似于 (5)。在这种情况下不需要初始条件。如果流场不是稳态的则必须指定含水层系统内每个点的年龄输运的初始条件。对于模拟相对于通过含水层系统的平流速率而言非常长的时间地下水年龄对初始年龄分布不敏感尽管在数学上需要初始年龄来产生解。年龄质量内部源年龄质量的内部净源 $F$ 可以解释几个过程。许多含水层系统仅在水平维度上建模因为垂直水头梯度和流速被假设为很小。这种情况下的控制方程可以通过在饱和厚度上进行垂直积分从 (5) 或 (11) 获得。在得到的一维或二维模型中$F$ 可以包括由于来自底层水文地质单元的流入而产生的年龄质量通量。例如对于通过隔水层渗漏的情况$F$ 将是渗漏水的质量通量率与渗漏水年龄的乘积对于流入感兴趣含水层的情况为正。类似地从准二维水平模型中进行的蒸散将被视为年龄质量的净汇且 $F$ 为负。对于通用三维模型$F$ 方案可以代表由于相变或多重连续体交换而产生的年龄质量源和汇。例如裂隙岩石中的流动可以建模为双结构域系统其中水在裂隙和岩石基质之间交换。裂隙和岩石基质之间的水交换将包括年龄质量的交换。对于裂隙域$F$ 将是水的质量从基质流入裂隙的速率与基质中水的年龄的乘积。裂隙方程中的 $F$ 项的大小相等但符号相反。内部年龄质量交换的最后一个例子是非饱和流流经部分冻结的土壤。在这种情况下流出的水可能会冻结充当年龄质量的汇或者系统原有的冰可能会融化充当系统流动的年龄质量源。当然与其他示例一样对于源水年龄增加的情况必须从单独的、可能是耦合的数学模型中指定或确定源水的年龄。低流量或停滞区不一定包含在这个内部源项中但可以直接在控制方程中处理。在这些区域平流项很小地下水年龄主要由扩散决定扩散包含在弥散张量 $\mathbf{D}$ 中。在没有扩散的情况下稳态流场中停滞点的水的年龄按定义是无限的。然而这个无限年龄仅适用于具有无限小流体体积的一个点。模拟示例前一节开发的模拟地下水年龄的数学理论被应用于区域含水层流动的剖面模型。此应用演示了使用年龄输运方程包含及不包含弥散对地下水年龄分布的实际模拟。考虑了两个假设的含水层系统一个是均质含水层另一个是在一定深度存在高渗透层的含水层系统。这些配置类似于Freeze and Witherspoon [1967]在一系列里程碑式论文中所分析的配置他们在论文中使用了数值流模型来研究区域地下水流的特征。使用块中心有限差分流模型 MODFLOW[McDonald and Harbaugh, 1988]求解流动方程。一种三维特征类溶质输运模型 MOC3D[Goode and Konikow, 1991]被修改为求解 (10) 中的年龄输运方程包括年龄的零级源项。与输运方程的有限差分或有限元数值解相反MOC3D 非常适合仅平流$\mathbf{D}0$的情况。所考虑的两个假设区域含水层的域几何形状和边界条件是相同的。该域长 1 公里厚 100 米在垂直方向上离散化为 10 行每行 10 米厚在水平方向上离散化为 50 列每列 20 米长。左侧、右侧和底部边界指定为无流条件。顶边界代表潜水面由指定通量和指定水头条件的组合建模。由于潜水面移动导致的域边界位置变化被认为是一个次要影响在此予以忽略。这里模拟的区域流动系统与Freeze and Witherspoon [1967]考虑的相似但使用了不同的潜水面边界条件。Freeze and Witherspoon [1967]使用数值流动模型来研究不均匀水力传导率对区域地下水流的影响。由于潜水面高度和净地下水补给对地下渗透性敏感因此使用了Freeze and Witherspoon [1967]讨论的指定水头条件。因此潜水面高度在响应各种水力传导率配置时会发生变化但补给分布在两种情况下是相同的。这些模拟比Freeze and Witherspoon [1967]的模拟产生了更大的水头变化但流量变化较小。模拟了两个不同的含水层系统均质含水层和分层含水层系统。均质含水层的等向水力传导率为 $10^{-6}$ m/s其孔隙度为 0.2。分层含水层系统的水力传导率在系统的上部 70 米高十倍为 $10^{-5}$ m/s在下部 30 米为 $10^{-6}$ m/s。两个层的孔隙度均为 0.2。图 2a 和 3a 显示了稳态条件下的水头分布和地下水流速图 2b 和 3b 显示了相应的流线。施加的边界条件和均匀属性导致了水头的平滑分布以及速度随均匀含水层向下的逐渐变化图 2a。在分层含水层系统中大部分流动发生在上层深处那里的水力传导率较大图 3a。均质含水层中的最大年龄结果高度依赖于含水层结构和边界条件配置。对分层含水层系统的模拟除了将下层的水力传导率降低到 $10^{-7}$ m/s 之外其余配置与此处所示相同产生了非常不同的结果。在这种情况下系统的上部比下部更具导水性上层的年龄与均质情况相似。然而由于低渗透性下层的通量和速度显著降低出现了非常陈旧的水特别是向域的右侧。仅考虑平流这种情况下的最大体积平均年龄超过 1000 年。弥散对地下水年龄分布的影响可以通过求解输运方程并包括m 和m 的弥散度以及m²/s 的扩散系数来轻松获得。该扩散系数比实际预期值高约 10 倍用于说明在这些假设含水层系统中扩散的最大可能影响。仅此扩散的影响如下所述。不同年龄水的弥散和扩散混合往往会限制最大年龄图 2d 和 3d。根据普遍接受的菲克定律模型溶质弥散通量的发生方向是浓度降低的方向。在完全类比的情况下年龄质量的弥散通量发生在年龄减少的方向即远离最大年龄区域。在平流通量相对较小的地方即速度较小的地方这种弥散通量对年龄分布的影响最大。纵向弥散的影响在年龄输运的情况下是减轻的因为沿着流线年龄由于老化而顺着流动方向平滑增加。因此通常与推进的溶质前缘相关的陡峭纵向梯度并不存在。然而年龄的陡峭梯度存在于垂直于流动方向的方向特别是在分层含水层系统中图 3c 和 3d。在这些陡峭梯度的区域弥散对地下水年龄的影响最大。图 4 是分层情况下地下水年龄百分比变化的等值线图范围从 -37 到 64%这归因于扩散和弥散。为了进一步检查扩散和弥散相对于先前结果的相对贡献并将在本文中开发的理论与停留时间分布粒子追踪方法进行比较在包含平流和扩散但不包含弥散的情况下进行了分层含水层系统的模拟。图 5 显示了通过本文提出的理论获得的结果这些结果与仅平流的结果相似但在最大年龄区域表现出轻微的减少。在速度很高的底层扩散对年龄分布没有可观察到的影响。图 5 还显示了应用于同一问题的粒子追踪模型的结果。粒子路径和运移时间通过线性速度插值[Goode, 1990]和随机游走[Kinzelbach, 1988]计算以模拟扩散。地下水年龄分布通过对 (1) 式进行数值积分计算得出。使用这些替代方法获得的结果之间的高度一致性进一步支持了此处提出的理论论据。总结提出了一种模拟地下水年龄的新方法。地下水年龄的空间分布受包含单位强度 (1) 内部源项对应于老化速率的输运方程控制。该控制方程既是从停留时间分布概念推导出来的也是从应用于概念性年龄质量的质量守恒原理推导出来的。这种地下水年龄模拟方法一方面优于现有的方法即仅受平流控制的地下水年龄模拟通过使用粒子追踪模型另一方面它优于溶质输运模型中感兴趣的同位素或化学标记物的模拟。使用此处介绍的理论允许进行一般的年龄模拟而无需特定于示踪剂的建模并通过整合扩散、弥散、混合和交换过程对年龄的影响。这些方法可以很容易地合并到现有的含水层输运模型中。致谢我感谢我的审稿人 Michael A. Celia 和 Allen M. Shapiro以及对本文初稿提出建设性意见的 WRR 审稿人 Charles F. Harvey、Eileen Poeter 和 D. Kip Solomon。这项工作得到了美国地质调查局水资源部研究生学院培训计划和普林斯顿大学的支持。