氢氧同位素（D和¹⁸O）作为天然稳定示踪剂，被广泛用于研究历史气候变化与水文过程^[1]。地下水的氢氧同位素特征也可用来指示地下水循环过程^[2]，例如识别地下水的补给来源、补给速率和遭遇的地质事件等^[3-5]。在大陆区域，大气降水的氢氧同位素含量具有高程效应，即一般随着高程的增加，气温逐渐降低，加速了水汽凝结成雨和同位素分馏，云团上升到高海拔地区后重同位素含量越来越少，导致D和¹⁸O含量下降^[1]。当大气降水入渗补给形成地下水时，其高程效应的信息会伴随地下水循环过程而保留在地下水的氢氧同位素特征中。利用这个性质，可以评估泉水的补给高程^[2,6]。

区域尺度的地下水循环可以发育不同级次的地下水流系统^[7-9]并影响地下水化学、年龄和同位素信息的空间分布^[10-14]。氢氧同位素在含水层的空间分布同样会受到地下水流系统的控制。不过，受取样条件的限制，大多数区域尺度的地下水同位素调查研究得到的是D和¹⁸O含量随平面位置的变化，而对其深度效应的研究相对较少。张之淦等^[10]通过分析不同深度的井孔水样测试数据，发现华北平原浅层地下水总体上比深层地下水更加富集D和¹⁸O，并初步认为这是浅层地下水蒸发作用导致的。王恒^[15]收集了鄂尔多斯高原都思兔河沿岸地区83口不同深度的水井取样测试数据，表明在埋深1000 m范围内，D和¹⁸O含量随深度的增加而减小，并推测这种现象可能与地下水流系统发育特征有关^[16]。王振等^[17]在青藏高原黄河源的玛曲盆地进行地下水调查，获得25组不同深度（40 m以内）的潜水样品，发现地下水D和¹⁸O含量随深度增加而减小。这种地下水D和¹⁸O含量随埋深增加而衰减是不是普遍规律、具有什么样成因机理？目前还缺乏系统的研究。

如果把含有D和¹⁸O的水分子视为地下水中的溶质，有望通过溶质运移数值模拟的方法从理论上研究地下水流系统对氢氧同位素空间分布的影响。Krabbenhoft等^[18]对威斯康星州北部的一个湖泊及周边构建三维地下水流和溶质运移模型，并使用观测的水位数据和¹⁸O数据对模型进行了校正，成功对湖泊和地下水交换量进行了估算。Shurbaji等^[19]开发了一个耦合水热气和同位素运移的数值模型，可以用来重塑土壤中D和¹⁸O同位素在剖面上的变化。Braud等^[20]开发了适用于裸土的耦合水热和稳定同位素运移的SiSPAT-Isotope模型，并对该模型进行了检验，分析认为同位素浓度对土壤表面和大气之间同位素迁移阻力的形成非常敏感。Caschetto等^[21]使用SEAWAT软件对意大利北部波河低地典型区域构建了数值模型，并使用水化学、D和¹⁸O等同位素对模型进行了校正，估算了沿海含盐含水层的地下水滞留时间和补给模式。Jiang等^[22]开发了一个同位素运移模块并将其耦合到TOUGH2软件中，使其能够模拟由于水的相变和密度变化引起的对流、弥散、水岩相互作用和分馏影响下地下水中D和¹⁸O的演化过程。同时，将该模型应用到贵德盆地，发现温度变化引起的密度驱动水流可导致低密度水和高密度水之间发生同位素分馏。那金等^[23]进一步使用TOUGH-Isotope 程序对康定—老榆林地区地热系统的D和¹⁸O同位素的迁移过程进行数值模拟，发现在高温地热系统中D和¹⁸O同位素的运移过程会严重受到对流弥散作用的影响。这些研究注重于模拟方法，没有探讨大气降水氢氧同位素高程效应向地下水氢氧同位素深度效应的转化机理。

笔者认为，大气降水中氢氧同位素分布的高程效应，对补给区的地下水具有直接的影响，并能够通过地下水循环过程转化为排泄区氢氧同位素的深度效应。从理论上解释这个转化机理并不一定需要十分复杂的溶质运移模拟方法。本文以没有地下水开采等人类活动影响的天然状态为背景，构建带有普遍特征的区域尺度地下水循环模型，尝试利用常规的对流-弥散过程数值模拟方法研究地下水氢氧同位素的运移规律，从一般意义上模拟D和¹⁸O同位素含量的空间分布特征，探讨氢氧同位素高程在地下水流系统转化为深度效应的机理。研究结果具有抛砖引玉的作用。

1.   概念模型、控制方程与公式

1.1   概念模型框架

本文采用如图1所示的剖面二维概念模型从理论上研究地下水氢氧同位素从补给区到排泄区的再分布机理。该模型表现了一个流域或盆地的典型特征，即地下水在地势高处获得补给，经过一定纵深的水循环路径，在地势低处排泄返回地表。图1中A-E和C-D是两条垂线表示区域尺度的地下分水线，E-D为含水层底部隔水边界（假设为水平延展）。大气降水在补给区（A-B区间）的氢氧同位素具有高程效应，并传递到补给区的潜水面，导致不同地点补给的地下水具有不同的氢氧同位素含量。地下水在循环过程中，携带氢氧同位素发生对流运移，同时发生弥散混合作用，使氢氧同位素在含水层空间重分布。地下水到达排泄区（B-C区间）时，氢氧同位素呈现随深度的有规律变化。而泉水及排泄区的地表水混合了不同来源的地下水。

图  1  大气降水同位素高程效应与地下水同位素深度效应示意图

Figure  1.  Schematic diagram of the isotopic altitude effect in precipitation and isotopic depth effect in groundwater

下载: 全尺寸图片 幻灯片

在概念模型框架下，本文尝试通过地下水流和溶质运移数值模拟的方法研究地下水氢氧同位素重分布的宏观规律。模型的构建包含一些简化的假设条件，然后在此基础上选择具体的控制方程和经验公式。

1.2   地形变化的描述

假设地面高程的变化可以表示为以下三角函数组合形式：

式中：z_top——地面高程（海拔高度）/m；

z₀——参考高程/m；

L——盆地含水层长度/m；

x——相对E点的距离/m；

Δz_m——区域尺度的地形变化幅度/m；

Δz_f——局部尺度的地形变化幅度/m。

当Δz_f =0时，式（1）描述的是一个单向倾斜的盆地，只发育区域地下水流系统。当Δz_f ≠0时，地面会存在局部的波动，导致补给区A-B之间存在一个洼地，一定条件下可发育局部流动系统。

1.3   地下水循环的动力学条件

模型的侧面（A-E、C-D）和底部（E-D）均假设为隔水边界，地面接受大气降水的均匀入渗补给。含水层中形成地下水的稳定流场，且潜水面下方饱和带服从孔隙介质地下水渗流的Darcy定律，可用以下控制方程进行描述：

式中：H——水头/m；

K_h、K_v——含水层水平方向、垂向的渗透系数/（m·d⁻¹）， 渗透各向异性的主轴与坐标轴一致。

地下水渗流的上边界，在补给区以潜水面为边界，在排泄区以地表为边界，简化描述为：

式中：z_w——潜水面的高度/m，是距离x的函数；

w——入渗补给强度/（m·d⁻¹）；

C_d——等效排水因子/d⁻¹，表示地下水向泉水或地表水排泄的能力。

式（3）描述补给区上边界，根据地下水动力学中潜水面的边界条件^[24]推导而来。式（4）是地下水溢出带的经验描述。地下水的上边界与饱和带渗流场处于耦合状态，模型需要根据这种耦合关系联合求解出H（x, z）和z_w（x）。

本文考虑含水层渗透性随深度衰减的普遍规律。假设水平渗透系数随埋深的增大服从以下指数衰减趋势：

式中：K₀——地表处渗透系数/（m·d⁻¹）；

z_top−z——埋深/m；

β_k——渗透系数的衰减系数/m⁻¹。

含水层垂向渗透系数根据各向异性比K_h/K_v来计算，并假设各向异性比不随埋深变化。假设含水层的有效孔隙度（ϕ）随着埋深增大而呈现指数衰减，即：

式中：ϕ₀——地表处有效孔隙度；

β_p——孔隙度的衰减系数/m⁻¹。

1.4   氢氧同位素定量

含有重同位素的水分子一般类型为HDO和H₂¹⁸O，称为e类型分子。在纯水中的浓度可量化为：

式中：C_e——e类型分子的质量浓度/（g·m⁻³）；

m_e——体积为V（m³）的纯水中e类型分子的总 质量/g；

M_e——e类型分子的摩尔质量/（g·mol⁻¹）；

N_e——e类型分子的摩尔浓度/（mol·m⁻³）。

将e类型分子摩尔浓度与普通水分子（H₂O）摩尔浓度的比值记为R_e，则式（7）可转化为：

式中：M_w——普通水分子的摩尔质量/（g·mol⁻¹）；

ρ_w——普通水分子的密度/（g·m⁻³）。

对于自然形成的水体，ρ_w可以近似用纯水的密度代替。

以国际原子能机构发布的标准化海水（Standard Mean Ocean Water，SMOW）同位素比值$ {R}_{{\rm{e}}}^{*} $为参考，重同位素的相对浓度δ_e可以表示为以下千分差值（‰）：

对于HDO和H₂¹⁸O，标准化海水同位素比值$ {R}_{{\rm{e}}}^{*} $分别为0.15576‰和2.0052‰^[25]。其它类型的重同位素分子含量相对可以忽略不计，因此水体中D和¹⁸O的测试结果基本反映的是HDO和H₂¹⁸O的相对浓度。对比式（8）和式（9），可以根据千分差值得到重同位素分子的实际浓度：

由式（10）可知，重同位素的实际浓度与相对浓度满足确定的线性关系。

大气降水中δD的高程效应一般表示为：

式中：δ^*D——等效海平面处大气水的D含量/‰；

η_D——D含量随高程线性衰减的系数/m⁻¹。

相应地，大气降水中δ¹⁸O的高程效应可表示为：

式中：δ^*18O——等效海平面处大气水的¹⁸O含量/‰；

η_O——δ^*18O随高程线性衰减的系数/m⁻¹。

在自然状态中，大气降水的氢氧同位素不仅仅受到高程的影响，也可能受到随纬度或其它地理因素的影响，还可能存在季节性变化乃至千年、万年尺度的趋势性变化。本文模型暂且只考虑高程效应，即单纯研究高程效应所引发的地下水同位素重分布特征。

1.5   地下水中氢氧同位素的运移

完成同位素分子浓度的定义之后，可以把HDO和H₂¹⁸O视为地下水的溶质，采用对流-弥散方程描述两者在地下水流系统的运移过程^[26]：

式中：s_i、s_j——空间坐标/m；

i、j——三维正交坐标系x, y, z任意坐标轴方向 （本文模型忽略y方向），且微分算子遵守 爱因斯坦求和约定；

[D_ij]——弥散系数张量/（m²·d⁻¹）；

V_i——Darcy流速/（m·d⁻¹）；

r——补给源的强度（单位体积含水层获得的补给流量）/d⁻¹；

${C}_{{\rm{e}}}^{{\rm{s}}}$——e类型分子补给源的分子质量浓度/（g·m⁻³）。

弥散系数张量的系数与弥散度有关，一般可以表述为^[27]：

式中：D₀——分子扩散系数/（m²·d⁻¹）;

α_L——纵向弥散度/m；

α_T——横向弥散度/m。

弥散度具有尺度效应。在本文模型中，假设V_y=0且$\partial {C}_{{\rm{e}}}/\partial {y} $=0，同时假设弥散度、分子扩散系数不随埋深变化。

地下水中氢氧同位素的来源是大气降水。从大气降水入渗到形成地下水补给，包气带水分运移过程也可能引起氢氧同位素的变化，而且补给氢氧同位素也可能随时间变化。本文模型假设这种变化可以忽略，从而使潜水面处地下水的氢氧同位素与相同x坐标处地面降水的氢氧同位素浓度近似相同，仅仅研究完全稳定状态（t→∞）的氢氧同位素分布特征。水岩相互作用导致的氧同位素漂移等现象，本文暂不考虑。另外，在干旱、半干旱区，潜水蒸发作用也可能导致排泄区浅部的地下水氢氧同位素发生变化。本文模型忽略这种潜水蒸发作用。关于诱发氢氧同位素变化的其它因素的影响，会在给出模拟结果之后加以讨论。

2.   情景设计和模拟方法

2.1   算例情景及参数设置

本文模拟3种典型情况下的氢氧同位素重分布特征。算例Case-I 模拟地面单向倾斜盆地区域地下水流系统的情形。算例Case-II和Case-III均模拟发育局部洼地的盆地，但前者只发育1个区域地下水流系统，而后者潜水面“切过”洼地形成并在洼地两侧发育局部地下水流系统。这些算例的控制参数如表1所示。模型中盆地的总长度为L，即图1中E点与D点之间的距离。含水层的最小厚度，取决于C点高程z_top(L)与D点的高程z_D的差。算例均以A点为最高点，高程1300 m，以C点为地形最低点，高程100 m，含水层最小厚度均为 500 m。模型物理参数在一般水文地质条件和大气降水同位素的经验取值范围内。

表  1  不同模拟情景的控制参数

Table  1.  Control parameters in different simulation scenarios

参数		模拟情景
		Case-I	Case-II	Case-III
地貌 形态	z0/m	100	500	500
	Δzm/m	1200	800	800
	Δzf /m	0	400	400
	其它	L=12 km，zD= −400 m
渗透性	K0 / （m·d−1）	1.00	1.00	0.25
	其它	Kh/Kv=10，βk =0.002 m−1
补给与排泄		w=0.001 m/d，Cd =0.1 d−1
孔隙度		ϕ0=0.2，βp =0.001 m−1
弥散参数		αL=20 m，αT=2 m，D0=0.0063 m2/d
同位素D		δ*D= −50.0‰, ηD=0.0300 m−1
同位素18O		δ*18O= −7.5‰, ηO=0.0038 m−1

下载: 导出CSV 

| 显示表格

根据给定的大气降水氢氧同位素高程效应参数，利用HDO和H₂¹⁸O分子量以及纯水在标准状态下的密度，可以将大气降水H₂¹⁸O质量浓度计算为：

可以将大气降水HDO质量浓度计算为：

2.2   模拟方法

以上算例情景的模拟是利用GMS软件完成的，调用MODFLOW进行地下水流模拟，调用MT3DMS进行重同位素分子的运移模拟。按照最大高差范围形成的矩形区域划分有限差分网格，网格单元在剖面上为矩形，横向单元长度均为100 m，垂向单元厚度均为5 m，并把高于z_top的单元设置为无效单元。

MODFLOW水流模型，设置为稳定流模拟，利用矩阵计算工具设置不同埋深单元的渗透系数和孔隙度。降水入渗调用Recharge模块，施加在最高的有效单元上。潜在的排泄区顶部单元（图1中B-C段）采用Drain模块处理，排水高程为z_top，排水系数设置为参数C_d与单元面积的乘积。使用Drain模块是处理地下水渗出面或溢出带的一种方法。有限差分方程组的求解采用共轭梯度法，收敛精度为0.001 m。潜水面的计算形成是一个迭代过程：将水头初始状态设置为地面高程，MODFLOW在进行渗流模拟时，会把计算水头低于底部高程的单元自动处理为无效单元，从而逐渐疏干一部分形成包气带。如果迭代计算过程中某个疏干单元周边的水头超过其底部高程一定幅度，也可以让这个疏干单元重新变为含水的有效单元。通过调用MODFLOW的干湿单元切换功能，能够最终模拟形成稳定状态的潜水面和溢出带宽度。

在使用MT3DMS模拟HDO和H₂¹⁸O运移过程时，将初始浓度设置为式（16）和式（17）计算结果的最小值。溶质运移唯一补给源是降水入渗，根据单元中心x坐标对应的z_top计算入渗水的氢氧同位素浓度，输入到MT3DMS的补给源项中。溶质对流-弥散计算采用非稳态的运移过程模拟方法，最初的若干模拟阶段（Stress Period）运移时间长度为1000 d，时间步长为1 d，如此经过5个阶段的模拟之后把运移时间长度增加到5000～20000 d，时间步长数目保持为1000个，直到变化最慢的单元浓度值接近一个稳定值。本文模拟没有采取稳定态对流-弥散方程的求解，因为MT3DMS限定稳定态溶质运移只能使用上游加权的有限差分法，受到数值弥散的严重影响，误差很大。采用非稳态运移过程模拟方案，可以调用更准确、有效的修正特征值法进行求解，逐步逼近t→∞时的稳定态。

得到地下水的H₂¹⁸O浓度模拟结果之后，可以根据式（10）反算其千分差值：

类似的，模拟得到地下水的HDO浓度之后，可以反算以下千分差值：

3.   模拟结果

3.1   Case-I情景

在算例情景Case-I 中，盆地具有单向的倾斜地貌，因此越靠近补给区的A点，大气降水的氢氧同位素含量越低。在补给区与排泄区之间，形成单一的区域地下水流系统，图2（a），但实际上属于闭合型地下水循环单元。模拟出的潜水面埋深在700 m以内，溢出带的最高点（B点）与模型右侧边界（C点）相距约为1.9 km。在距离C点1.0 km处取一个溢出点作为泉点，观察地下水中¹⁸O含量随埋深的变化。整体上地下水中¹⁸O的空间分布如图2（b）所示，即在B点附近的上游区δ¹⁸O最高，而距离这个峰值区越远的空间点处δ¹⁸O越小，呈环绕形分布特征。在泉点处，δ¹⁸O随埋深呈非线性衰减趋势，图2（c），且浅层衰减快、深层衰减慢，可以近似拟合为：

图  2  算例情景Case-I模拟结果

Figure  2.  Simulation results in the Case-I scenario

下载: 全尺寸图片 幻灯片

式中：d——埋深/m。

同样的方法也可用于获取δD的浓度变化，它与δ¹⁸O保持大气降水所决定的线性关系，同步变化。δ¹⁸O的变化范围明显小于δD的变化范围，如果地下水δ¹⁸O随深度具有显著趋势变化特征，那么δD的特征将会更加明显。因此，本文在情景Case-I 中展示δ¹⁸O模拟结果代表重同位素的深度效应，而在其他情景中展示δD的模拟结果，起到对照作用。

在浅层，地下水δ¹⁸O随埋深增大的衰减速率约为0.03 m⁻¹，显著大于降水δ¹⁸O随高程的变化率（η_O=0.0038 m⁻¹），这说明补给区氢氧同位素的高程效应在排泄区深度效应中得到了放大。在补给区，潜水面下部δ¹⁸O随深度增大也是呈现近似指数衰减趋势。

3.2   Case-II情景

算例情景Case-II 中，盆地在倾斜地貌的基础上叠加了波状的地形起伏，产生了一个上游洼地，如图3（a）中的G点所示。由于地形在排泄区附近比在Case-I中显著变陡，导致溢出带变窄，B点更加靠近C点（相距约0.5 km），因而流线在排泄区更加集中，图3（a），流速更大。相对而言，潜水面在补给区显得更加平缓，由于没有高出G点，仍然只发育单一的地下水流动系统。

图  3  算例情景Case-II模拟结果

Figure  3.  Simulation results in the Case-II scenario

下载: 全尺寸图片 幻灯片

地形的波状起伏使Case-II 中补给区地下水的氢氧同位素含量发生了比Case-I中更加复杂的变化。取D含量的模拟结果进行观察，可以发现补给区G点下部的地下水δD高于A点和F点下部，随流线的延伸而形成一种舌状的高δD值条带（图3b），并一直影响到排泄区下部δD的分布。在溢出带取距离C点350 m处作为泉点，观察δD随深度的变化。如图3（c）所示，δD变化呈S形曲线，在埋深约180 m处存在一个峰值（δD≈−77.7‰），其δD值略低于补给区G点下部地下水的δD值（−73.9‰），这是重同位素分子弥散运移并混合的结果。在埋深50 m以内，δD随深度增加的衰减速率约为0.16/m，也显著大于降水δD随高程的变化率（η_D=0.03/m）。在埋深200～500 m区间，δD随深度增加的平均衰减速率与降水δD随高程的变化率很接近，略偏大。这说明在Case-II 中排泄区深层地下水氢氧同位素的深度效应能够重现补给区氢氧同位素的高程效应。

3.3   Case-III情景

算例情景Case-III与Case-II的不同之处在于渗透系数偏小，因此补给区潜水面更加壅高，局部高于G点而溢出。模拟过程中利用MODFLOW的Drain模块在G点周边设置排水单元，地下水流场模拟结果如图4（a）所示。盆地发育3个局部流动系统和1个区域流动系统，G点下方发育一个地下水流的驻点作为地下水流系统的分界点。

图  4  算例情景Case-III模拟结果

Figure  4.  Simulation results in the Case-III scenario

下载: 全尺寸图片 幻灯片

由于出现了2个地下水排泄区和更多的地下水流系统，Case-III 中地下水氢氧同位素的空间分布比在Case-II 中更加复杂，见图4（b）。首先，补给区舌状的高δD值条带只在G点的下游方向发育，而G点上游方向的地下水δD值总体受局部流动系统的控制。G点正下方是δD低值区，把洼地的高δD值隔开了，这是洼地2个局部流动系统的流线向G点附近集中排泄的结果。其次，潜水面的整体抬高使溢出带向上游移动，B点附近的高δD值区具有更高的海拔，也导致距离C点350 m处的泉点只能获得更高来源的地下水补给，使得氢氧同位素更加贫化。从图4（c）可以看出，Case-III情景中泉点处地下水的δD值与Case-II 中相比，在400 m深度范围内出现整体的减小现象（变化范围也缩小了），尽管S形曲线的形态是相似的。

在上游洼地的排泄区，即G点附近，地下水δD值随埋深变化基本上呈现单调的衰减现象，见图4（c）。在埋深100～800 m范围，δD值随深度增加的衰减率约为0.008 m⁻¹，显著小于降水δD随高程的变化率（η_D=0.03 m⁻¹）。因此，大气降水氢氧同位素的高程效应在上游排泄区的深度效应中被减弱了。

4.   讨论

4.1   重同位素分子的弥散作用

本文的模拟研究是把重同位素分子HDO和H₂¹⁸O视为地下水的溶质，求解对流-弥散方程得到其空间分布。由于这些分子是水质点的一部分，相比其他的盐分离子，其更加服从水质点在多孔介质中的运移规律，即按照Darcy流速进行整体的对流运移，在此基础上叠加分子热运动和微观质点流速空间变异产生的机械弥散。目前，还不清楚弥散运移是否产生氢氧同位素分馏效应，但有这种可能性，因为它们的分子量比普通水分子大，可能弥散行为更弱一些。如果没有弥散作用，那么氢氧同位素将完全按照流线的路径进行投射，保持与补给来源一致，而弥散作用可以把不同来源的重同位素分子混合在一起。弥散度具有尺度效应，以上算例模拟采用了符合区域尺度特征的纵向弥散度数值（α_L=20 m）。为了考察弥散度取值不确定性的影响，对算例Case-I和Case-II进行了其它弥散度取值（α_L=1 m或 200 m）的试算，得到泉点处氢氧同位素含量随埋深的变化结果，如图5所示。显然，纵向弥散度取1 m和20 m得到的结果十分相近，说明这种情况下弥散作用并没有显著破坏对流作用形成的氢氧同位素分布格局。但是，如果弥散度进一步增大，例如增加到α_L=200 m，则泉点处的δD或δ¹⁸O剖面将发生显著变化：在Case-I中发生右移（图5a），也就是高海拔补给区的水源与低海拔水源显著混合，但近地表的δD或δ¹⁸O则有所降低，意味着来自低海拔水源的浅部地下水也受到了高海拔补给水源的影响；在Case-II中，δD或δ¹⁸O随深度的变化成单调减小的趋势（图5b），强烈的混合作用几乎完全抹平了S型曲线变化。弥散作用总体会削弱氢氧同位素从大气降水高程效应向地下水深度效应的转化。因此，弥散作用较强时，直接用氢氧同位素高程效应判断地下水补给来源可能会产生较大的误差。

图  5  泉点处地下水氢氧同位素模拟结果随纵向弥散度的变化

Figure  5.  Variation of simulated hydrogen and oxygen isotopes in groundwater with vertical dispersion at the spring point

下载: 全尺寸图片 幻灯片

4.2   对跨流域地下水循环的潜在指示作用

根据区域地下水流系统理论，地下水的分界行为与地表水有很大的差异。地表水的流域划分主要利用地貌分水岭（地形最高点）来进行，而地下水循环单元的划分需要确定流动系统的三维分界面，在剖面上流动系统的边界点可能处于潜水面最高点附近、但不等于潜水面最高点。这就意味着流域分水线（分水岭、分水点）并不能完全约束地下水循环，可以发生跨流域地下水循环的现象，即高处流域的地下水穿过分水岭断面流向低处流域。

如果没有发生跨流域地下水循环，那么地下水循环及氢氧同位素含量的空间结构总体上类似于Case-I情景（图2），排泄区地下水的δD和δ¹⁸O值随深度单调衰减。在Case-II和Case-III情景中（图3—4），A-G-F段可代表高处流域，而F-B-C段可代表低处流域，从地下水流系统的模拟结果可以判断发生了从高处流域到低处流域的地下水循环，但一般水力学要素（水头）的观测很难表征这种跨流域地下水循环。低处流域排泄区地下水δD和δ¹⁸O值随深度的变化曲线呈现S形，恰好指示地下水获得了分水岭外侧低海拔处的补给源，而这个补给源属于高处流域。因此，非单调曲线形态的地下水δD和δ¹⁸O深度效应具有指示跨流域地下水循环的潜力。当然，这是单纯考虑大气降水高程效应所得到的结论，实际地下水的氢氧同位素变化还受到其他因素影响，需要做更进一步的分析才能做出准确判断。

4.3   假设条件的限制

本文的模拟结果受到不少假设条件的限制。由于仅考虑了大气降水氢氧同位素的高程效应，虽然能够反映补给水源同位素含量的空间变化，但忽略了氢氧同位素随时间的变化。大气降水的氢氧同位素存在显著的季节性升降和气候变化引起的长周期振荡。季节性变化从大气降水到地下水的传递一般并不显著，因为地下水补给一般代表多年平均状态，包气带缓冲作用过滤掉了季节性短周期波动。然而，长期气候变化可能转化为补给区地下水δD和δ¹⁸O的深度效应，并影响排泄区地下水氢氧同位素的深度效应。这意味着S形曲线的δD或δ¹⁸O剖面也可能是气候变化导致的，需要联合使用地下水年龄和氢氧同位素信息才能还原补给过程。干旱区强烈的土壤蒸发也可能导致补给区潜水面处氢氧同位素含量与大气降水的差异，从而破坏地下水补给源δD和δ¹⁸O与大气降水一致的假设。不过，有相当多的观测结果表明补给区地下水的氢氧同位素数据点一般落在大气降水线附近，而不是显著偏离到蒸发线上。其主要原因，可能在于实际的地下水补给一般是由包气带优势流提供的，而非均匀平摊的活塞式补给。至于这种优势补给通道所具有的尺度及其真实的同位素效应，是值得进一步研究的问题。此外，本文没有考虑排泄区蒸发作用导致地下水氢氧同位素向蒸发线漂移的情况，也没有考虑水-岩相互作用可能导致的氧漂移过程，这些因素如何影响地下水流系统氢氧同位素的空间分布，也需要进一步探讨。

5.   结论

本文建立了含有降水入渗补给区和地下水溢出带的区域尺度地下水循环模型，模拟了单向倾斜盆地与双峰波状盆地的稳态地下水流系统及其氢氧同位素分布，初步探讨了补给区大气降水氢氧同位素高程效应转化为排泄区地下水δD和δ¹⁸O值随深度变化的机理，得到以下主要结论：

（1）在单斜盆地，大气降水D和¹⁸O含量随高程的增加而减小，这个高程效应可以通过单一的区域地下水流系统转化为排泄区地下水δD和δ¹⁸O值随埋深增大而减小的趋势，即深度效应。即使补给区氢氧同位素的高程效应是线性的，在排泄区的深度效应也是非线性的，δD和δ¹⁸O值随深度的衰减近似符合负指数函数。

（2）在双峰波状盆地，当含水层渗透性相对入渗强度较大（K₀/w=1000）时，潜水面较低，也仅发育一个区域地下水流系统，但是在区域地下水的排泄区δD和δ¹⁸O随埋深增大而呈现S形曲线分布，总体上含水层底部δD和δ¹⁸O仍然低于浅表地下水。

（3）在双峰波状盆地，当含水层渗透性相对入渗强度较小（K₀/w=250）时，潜水面较高，可以发育多个局部地下水流系统，对氢氧同位素的空间分布形成扰动。这种情况下，区域地下水的排泄区δD和δ¹⁸O随埋深增大呈现S形曲线，而局部排泄区的δD和δ¹⁸O随深度增加仍然表现为单调衰减趋势。

在自然环境中，地下水氢氧同位素的时空变化受到很多因素的影响，高程效应向深度效应的转化发挥了什么作用，需要根据具体情况进行综合判断。