第5章 华北平原非饱和带中地下水滞留时间

包气带中地下水运动时间(地下水自地表运动至潜水面所需的时间)是影响浅层地下水年龄,评价地下水是否受到地表污染影响的重要因素。

点上的包气带中的地下水运动时间可通过示踪剂的运移模拟估算,计算示踪剂进入包气带到达潜水位的运动时间作为地下水在包气带中的运动时间,

而区域尺度上非饱和带中地下水运动时间可通过地下水补给量和非饱和带的水文地质特性估算。

5.1 区域地下水补给计算方法

在干旱、半干旱地区,农业依赖于灌溉,因此地下水补给以及蒸发量的准确估算对实现地下水资源的可持续管理极为重要。

然而,由于土壤质地和植被类型的空间变化、以及气象因素,如降水和蒸发的时空变异,地下水补给量通常不易估算,并且在灌溉的条件下,地下水的补给过程会更为复杂(Gee和Hillel,1988)。

传统的测量地下水补给的物理方法(如蒸渗仪和渗透仪lysimeter,infiltrometer)和测量蒸散发的物理方法(如涡度相关和蒸渗仪),以及诸多经验方法通常要求大量的试验数据,因此往往局限在小尺度应用。

1. **蒸渗仪（Evapotranspiration）**：蒸渗仪用于测量土壤中的蒸发和植被蒸腾。它通常包括一个放置在土壤表面的容器，测量水分流失速率来确定蒸发和蒸腾的水量。然而，它需要对土壤特性和植被影响进行详细的测量，并且对于大范围的应用来说可能需要大量的设备和数据处理。

2. **渗透仪（Lysimeter）**：渗透仪是一种用于测量土壤水分运动的装置。它通常由一个埋在土壤中的容器组成，可以测量土壤中水分的变化量。虽然渗透仪可以提供土壤水分运动的重要信息，但它的使用通常受到设备成本和数据分析复杂性的限制。

3. **Infiltrometer（入渗仪）**：入渗仪用于测量土壤的入渗速率，即水分渗入土壤的速率。它通常涉及将水或其他液体施加在土壤表面，并测量水分渗透的速率。尽管入渗仪是一种常用的测量土壤渗透性的方法，但它在大范围应用时可能需要大量的实验和数据处理。

4. **涡度相关（Eddy Covariance）**：涡度相关是一种直接测量地表蒸发和植被蒸腾的方法。它通过测量地表和大气之间的气体交换来估算水分的流失。尽管涡度相关提供了一种直接的测量方法，但它通常需要复杂的仪器和数据处理，因此在小尺度应用上可能不太实用。

这些物理方法和经验方法在测量地下水补给和蒸散发时提供了重要的信息，但它们的应用通常受到设备成本、数据处理复杂性和实验量的限制，因此在大范围或小尺度应用中可能存在局限性。近年来，随着技术的发展，一些新的方法和技术（如遥感技术、数值模拟等）也被引入到水文学研究中，以弥补这些传统方法的局限性。

1. **蒸渗仪（Evapotranspiration）**：蒸渗仪用于测量土壤中的蒸发和植被蒸腾。它通常包括一个放置在土壤表面的容器，测量水分流失速率来确定蒸发和蒸腾的水量。然而，它需要对土壤特性和植被影响进行详细的测量，并且对于大范围的应用来说可能需要大量的设备和数据处理。

2. **渗透仪（Lysimeter）**：渗透仪是一种用于测量土壤水分运动的装置。它通常由一个埋在土壤中的容器组成，可以测量土壤中水分的变化量。虽然渗透仪可以提供土壤水分运动的重要信息，但它的使用通常受到设备成本和数据分析复杂性的限制。

3. **Infiltrometer（入渗仪）**：入渗仪用于测量土壤的入渗速率，即水分渗入土壤的速率。它通常涉及将水或其他液体施加在土壤表面，并测量水分渗透的速率。尽管入渗仪是一种常用的测量土壤渗透性的方法，但它在大范围应用时可能需要大量的实验和数据处理。

4. **涡度相关（Eddy Covariance）**：涡度相关是一种直接测量地表蒸发和植被蒸腾的方法。它通过测量地表和大气之间的气体交换来估算水分的流失。尽管涡度相关提供了一种直接的测量方法，但它通常需要复杂的仪器和数据处理，因此在小尺度应用上可能不太实用。

这些物理方法和经验方法在测量地下水补给和蒸散发时提供了重要的信息，但它们的应用通常受到设备成本、数据处理复杂性和实验量的限制，因此在大范围或小尺度应用中可能存在局限性。近年来，随着技术的发展，一些新的方法和技术（如遥感技术、数值模拟等）也被引入到水文学研究中，以弥补这些传统方法的局限性。

区域地下水补给量的估算通常采用的方法有:水均衡法、潜水位动态法(WTF)、化学和同位素方法以及数值模拟方法(Simmers,1987;Healy和Cook,2002;Scanlon等,2002;Xu和Beekman,2003)。

估算区域地下水补给量是水资源管理和水文学研究的重要任务之一。常用的方法包括水均衡法、潜水位动态法（WTF）、化学和同位素方法，以及数值模拟方法。这些方法为我们提供了从不同角度评估地下水补给量的手段。

1. **水均衡法**：水均衡法是通过对流域或地下水流系统进行水量输入与输出的比较来估算地下水补给量。该方法考虑了地表水和降水对地下水补给的影响，以及地下水流动和蒸发蒸腾的过程。虽然水均衡法是一种常用的估算方法，但其精度受到数据质量和模型假设的影响。

2. **潜水位动态法（WTF）**：潜水位动态法是通过监测地下水位的变化来估算地下水补给量。该方法基于地下水位变化与降雨、蒸发蒸腾和地下水抽取之间的关系。通过分析潜水位的时空变化，可以推断地下水的补给量和流动过程。

3. **化学和同位素方法**：化学和同位素方法利用地下水中的化学成分和同位素组成来推断地下水补给量和来源。通过分析地下水中不同元素和同位素的含量，可以揭示地下水的补给过程和流动路径。这种方法在解释地下水系统中的水文地质过程方面具有重要意义。

4. **数值模拟方法**：数值模拟方法利用数学模型对地下水流动和水文过程进行模拟和预测。这种方法基于对流域地质、地下水位、降水和蒸发等数据的输入，通过求解水文方程来估算地下水补给量。数值模拟方法可以考虑复杂的地下水流动过程和水文地质条件，提供了一种全面而有效的估算手段。

这些方法在估算地下水补给量时各有优劣，通常需要结合使用以提高估算精度和可靠性。随着技术的进步和水文学研究的深入，这些方法也在不断发展和完善，为地下水资源的管理和保护提供了重要支持。

仅应用某种方法估算地下水补给量的准确度一般很难判定,因此建议同时采用多种方法以对不同方法的结果互相验证以降低结果的不确定性(Scanlon 等,2002)。

另外,由于气象因素、土地利用类型以及地表物理和生物特性的空间变异,将点状计算结果插值到区域上往往造成较大误差(Scanlon等,2002)。

因此,对地下水补给估算方法的选择取决于研究的空间和时间尺度,以及研究区的数据情况(Scanlon等,2002;Delin等,2007)。

估算区域地下水补给量的主要方法为土壤水量平衡法,其中,土壤带或根植层底界的地下水流出量作为净入渗量(net infiltration)或潜在补给量(potential recharge)。

基于该种方法的数值模型很多,例如美国地质调查局(USGS)的INFIL 模型(USGS,2008),BALAN 模型(Samper和García-Vera,1992),WAVES模型(Zhang和Dawes,1998),和美国环保署(U.S. EPA)的HELP模型(Schroeder  等,1994)。

虽然这些模型能够简便快速的估算地下水补给量的区域分布,前提条件是地下水潜水面要低于根植层的底界,并假设根植层底部的渗出量即为地下水补给量(潜在补给量=实际补给量),认为自根植层底界至潜水面的非饱和带内的水流过程是可以忽略的。

然而,当非饱和带,或包气带厚度较大时,非饱和带将成为控制补给时间和水分传输的主要因素,忽略其中的水流过程是不可取的 (Harter和Hopmans,2004;Hunt等,2008)。

鉴于以上各种方法的优缺点,非饱和带数值模拟方法逐渐成为估计区域地下水补给时空分布的通用方法(例如,Keese等,2005;Small,2005;Liggett等,2009)。

此类模型主要包括 DAISY(Hansen等,1990),UNSAT-H(Fayer,2000),SWAP(Kroes和 van Dam,2003)和HYDRUS -1D(Šimůnek等,2005)。

虽然这些模型能够模拟非饱和带中的水流并估算蒸发量和地下水补给量,其应用过程中通常需要指定模型下边界条件(如潜水面位置),没有耦合饱和带地下水流动过程。

因此其结果不能直接反映出气候变化、土地利用和土壤质地对地下水动态的影响。

近年来开发的耦合非饱和带和饱和带水流过程的数值模型,如HYDRUS-3D(Šimůnek,2006),MODFLOW-VSF(Thom等,2006)和HydroGeoSphere(Therrien等,2005),是对上述非饱和带数值模型的巨大改进。

然而,这些模型的理论基础是刻画非饱和带水流过程的理查德(Richards’)方程,其应用主要受到计算机计算能力的限制(van Walsum和Groenendijk,2008)。

而对于流域尺度的模拟计算,将非饱和带内水分运动假设为一维流动通常是合理的(Harter和Hopmans,2004)。

因此,许多模型,如MIKE-SHE(Refsgaard等,1995)和MODFLOW-UZF(Niswonger等,2006)采用了非饱和带一维流动和饱和带三维流动的耦合方式以降低耦合模型的复杂程度和提高计算效率。

在这些模型中,一维的理查德(Richards’)方程被进一步近似成运动波(kinematic-wave)方程并采用特征值方法求解。

该求解方法的优点在于不需要对非饱和带进行网格剖分(Niswonger等,2006)。

基于非饱和带水流过程模拟估算区域地下水补给量的主要困难在于模型参数的获取和尺度转化,如土壤的渗透系数。

对于该过程,Harter和Hopman(s2004) 对多种地质统计方法提升土壤渗透系数的空间尺度做了细致的综述性讨论,但这些方法是基于大量的小尺度范围内的渗透系数估算结果(Vereecken等,2007)。

作为直接测量原状土壤渗透系数的替换方法,土壤转换函数方法(PTFs)(Wösten,2001)能有效的用于区域非饱和带渗透系数的估算。

利用土壤转换函数方法可以直接测量土壤特征指标(通常附带于土壤调查数据库)并估算非饱和的水力参数(Li等,2007)。

许多广泛应用的土壤转换函数方法是基于统计模型的方法(Pachepsky和Rawls,2004)估算土壤水特征函数和饱和渗透系数。

所利用的土壤特征指标通常为土壤颗粒分级,土壤干密度和有机成分含量(Nemes等,2003)。

绝大多数非饱和带数值模拟软件最初的开发目的是用于农业领域研究工作,需要大量反应植被生长特征的参数。

虽然土地利用和植被类型分布可以通过调查统计直接获得,遥感解译也有助于获取植被参数,但对于植被类型和土地利用不是地下水补给主控因素的地区,简化植被生长过程模拟,减少模型参数更有利于减少模型的不确定性和模型校正过程中的不唯一性。

蒸散发(ET),作为土壤水均衡的主要排泄项,对利用包气带数值模拟方法计算地下水补给有着重要影响(环境同位素方法估计地下水补给不需要蒸散发的计算)。

上述不同的模型中对蒸散发的计算过程也各不相同。

最常用的计算方式首先是基于常规气象观测数据计算潜在蒸发或参考蒸散发(ETp或ET0)(Lu等,2005),然后利用经验公式计算裸土蒸发和植物蒸腾(如Šimůnek等,2005)。

实际蒸发和植物蒸腾根据土壤含水量情况和植被信息确定,并加和得到实际蒸散发(ETa)。

但是,这种方法通常要求较小的时间尺度,如日气象数据。

因此,饱和带水流模型的时间步长也需要相应减小,但对于饱和带地下水流、特别是大区域地下水流模型来说,显然是不合适的。

作为一种折衷方案,可以利用子时间步长,即在饱和带水流模型时间步长的基础上再进行剖分,得到的更小时间步长作为非饱和带模型的时间步长,如耦合HYDRUS和MODFLOW的程序包(Seo等,2007)。

基于互补关系的蒸散发模型(Bouchet,1963)是另外一类用于估算区域蒸散发的常用模型,这类模型仅利用常规气象观测数据,不需要刻画陆面边界的详细信息。

因此,此类模型避免了复杂的土壤-植被-大气这一复杂系统,不需要土壤含水量数据和植被参数(Hobbins等,2001)。

不同的气候类型地区与区域水均衡的蒸发计算结果对比已经证实了该类方法的适用性(如,Hobbins等,2001;Xu和Singh,2005)。

此外,此类模型认为陆表大气湿度包括了所有陆面过程的影响,已经考虑了灌溉对蒸发的影响(例如,Ozdogan和Salvucci,2004)。

5.2 华北平原地下水补给研究现状

前人研究认为华北平原降水入渗补给占到地下水自然补给量的70-80 % (陈望和,1999;陈梦熊和马凤山,2002)。

华北平原地下水资源评价工作中通常采用渗透系数方法估算地下水补给量,入渗系数一般采用均衡试验得到,主要受土壤岩性控制,与降水量无关。

计算得到的区域地下水补给量范围较大,为 70-180 mm/yr(Zhu和Zheng,1983;刘昌明和魏忠义,1989;陈望和,1999;段永侯和肖国强,2003;任宪韶,2007)。

该方法是对地下水补给过程的高度简化,相当于认为地下水补给和降水量或灌溉量之间是一种线性关系(例如,Kendy 等,2003),然而实际的补给系数受土壤特性、植被、包气带厚度和气象条件等多种因素控制,且所有这些控制因素都是随时间和空间变化的(Scanlon等,2002;Li,2009)。

应用环境同位素方法 (Wang等,2008;Liu等,2009;Song等,2009;von Rohden等,2010;Yuan等,2011)和土壤水均衡和非饱和带数值模拟方法(例如,Kendy等,2004;Lu等,2011)可以估算小尺度地下水补给量。

这些小尺度地下水补给量估算结果在量级上基本一致(100s mm/yr;见表 5-2),这些计算结果仅能代表研究场地一定范围内特征地下水补给模式(降水入渗补给或灌溉回归补给),并与所采用的补给方法有关(见表 5-2)。

华北平原缺乏对集中式补给的定量研究,但这些补给仅限于局部地区(如,Yuan 等,2011)。

鉴于此,本研究进行了潜水面的地下水补给量估算。

影响华北平原地下水补给的一个关键因素是非饱和带厚度的增大。

华北平原长期大量开采地下水造成包气带厚度从70年代的2-15 m 增大到8-30 m甚至局部地区达到30-56 m(张光辉等,2004)。

忽略非饱和带中的水流过程,以及非饱和带厚度对地下水补给的影响,特别是对地下水污染的评估和地下水保护策略的制定都是不可取的。

虽然已有研究利用土壤含水量和土水势观测评估非饱和带厚度对入渗过程的影响(张光辉,1992;张光辉等,2007),但很少涉及非饱和带对区域地下水补给的影响研究。 

另一个关键问题是如何将局部地下水补给估算结果和模型参数提升到区域尺度,即水资源方案制定层面必须考虑的尺度范围。

一般来说,包括非饱和带厚度和渗透系数,降水模型,灌溉类型及模式,土地利用以及蒸散发在内等大量数据是提升补给量计算尺度的必要条件。

然而,考虑到华北平原的空间尺度及现有数据和研究成果,利用局部地下水补给量估算区域地下水补给量显然是不现实的。

虽然最近有研究利用数据融合技术耦合模型计算结果和环境同位素信息计算结果来估算区域地下水分布(如,Ng等,2009),但现阶段来说对区域研究计算量仍很大。

本次研究中采用的方法对区域地下水补给估算是适用的,并能够考虑控制地下水补给过程的一些关键因素。

 本研究耦合土壤水均衡模型与MODFLOW的UZF程序包,建立区域地下水流数值模型,模型能够模拟非饱和带水流过程并估算地下水补给量。

该方法能够利用多种有效估算区域陆面蒸散发的方法,如遥感方法,并且适用于区域尺度的模型校正和应用。

本研究的主要目的是评价土壤质地类型,非饱和带厚度变化对华北平原区域地下水补给的影响。

华北平原的空间尺度,土壤转换函数和已有数据能够支持建立合理适用的非饱和带-饱和带耦合水流模型来模拟地下水的动态演变。

模型校正不是本章研究的主要目的,但通过对比模型计算的地下水补给量和与其他估算方法取得的结果,以及计算地下水位和观测地下水位来验证计算得到的地下水补给量(见表 5-2)。

 

5.3 华北平原非饱和带-饱和带耦合数值模型

5.3.1 土壤水均衡模型

当没有明显地表径流时,某特定时段内根植层的土壤水均衡可表示为:

式中

Rp为根植层底部地下水渗出量,或潜在地下水补给;

P为降水;

Ir为灌溉;

ETa为实际蒸发量;

为土壤水储存量变化量。 

5.3.2 非饱和带水流模型

描述均质非饱和带中垂向一维水流过程的理查德方程可表示为: 

式中

θ土壤水体积含水量;

S为单位长度的蒸散量;

K(θ)为渗透系数函数,与含水量有关;

z为垂向上高程;

t为时间。

MODFLOW的UZF非饱和带程序包(Niswonger 等,2006)对公式(5-2) 进行了简化,去掉了扩散项,假设垂向上的地下水流动仅有重力驱动。

UZF中的控制方程表示为:

一般来说,基于理查德方程的非饱和带水流模型需要较细的网格剖分,因此对于计算时段较长的模型来说,要求较高的计算能力(Hunt 等,2008)。

再者,建立非饱和带渗透系数与土壤含水量关系函数的参数在区域范围内不能直接获得。

对于区域尺度上的非饱和带水流模拟,能够获得控制方程的闭合形式解,并保证数值计算的稳定性,因此更有利于大空间和长时间尺度的模型建立。

鉴于此,UZF采用运动波(Kinematic Wave)近似表达均质非饱和带中的一维理查德方程,这种描述方程更适用与大空间尺度模型(Niswonger 等,2006)。 

Brooks-Corey 模型(Brooks和Corey,1964)用来描述非饱和带渗透系数和土壤含水量之间的关系:

式中

Ks为饱和渗透系数;

θr为残余含水量;

θs为饱和含水量;

ε为Brooks-Corey 常数。  

本次研究的主要目的是估计区域地下水补给,此过程中将非饱和带在垂向上假设成均质条件一般认为是合理的(Niswonger和Prudic,2009)。

UZF的主要输入参数包括入渗强度(FINF),垂向渗透系数(VKS),饱和含水量(THTS) 和Brooks-Corey系数(EPS)(见表5-1).

上述土壤水均衡计算得到的正值Rp作为入渗量输入,而负值Rp表示水分从根植层下部非饱和带向根植层运动,作为蒸散发输入模型。 

5.3.3 模型数据输入

蒸发量采用第5章中CRAE模型计算的实际蒸发量。

土壤质地(包括砂,壤土和粘土含量)及土壤容重从中国1比100万土壤数据库获得(图 5-1)。

该数据库为分辨率1km的栅格结构,按照模型网格内主要质地类型提升到模型适用2km网格。

van Genuchten模型(van Genuchten,1980)中的七个土壤水力参数利用Schaap(2001)利用美国土壤数据开发的土壤转换函数计算得到。

θs和Ks(图 5-2)对应与UZF中的VKS和THTS。

Brooks-Corey 常数根据Morel-Seytoux 等(1996)导出公式计算: 

土壤含水量数据从中国农业气象数据库获得,该数据库在华北平原范围内包括104个观测站,观测数据为10,20和50 cm处的土壤重量含水量。

每月8号、18号和28号进行观测。

数据为相对含水量(土壤含水量与田间持水量之比),利用土壤容重和田间持水量数据将相对重量含水量转换为体积含水量后输入模型。  

华北平原没有详细精确灌溉量分布可供利用,小麦-玉米轮作灌溉需水量约为4,200 m3/ha/yr(柳长顺等,2004)。

华北平原小麦一般灌溉4次或5次,玉米一般灌溉3或4次。

模型将灌溉处理成7次,每年3月至7月各一次,9月和10月各一次,不区分作物类型,每次灌溉量为60 mm。

5.3.4 与饱和带水流模型耦合

根据上述讨论,本次研究中非饱和带-饱和带耦合水流模型实际上可认为是一个三层模型。

第一层模型为根植层的土壤水均衡模型,该模型将土壤水储变量和蒸散发从降水量中扣除,剩余量作为根植层底界和下部非饱和带的水分交换量;

第二层模型为非饱和水流数值模型,根植层底界渗出水分通过该模型传输到下部饱和带地下水;

第三层为饱和带水流模型,其表达与第5章中区域地下水模型一致。

模型模拟时段为1993年1月至2008年12月,每月作为1个应力期。

模型网格剖分与区域地下水模型一致。

模型初始条件采用区域非稳定水流模型计算得到的1992年底的水位分布。

非稳定流模型中开采量为年地下水总开采量,因耦合模型认为农业灌溉开采量占总开采的70%,并平均分配到3月至7月,9月及10月。

其余30%开采量作为工业开采和生活开采量平均分配到各月。

地下水开采的空间分布与非稳定流模型一致,而灌溉量均匀分布输入模型。

非饱和带中水流过程模拟中,初始含水量的确定通常采用非饱和带稳定流模拟,或者反复运行模型较长时间获得(Keese等,2005;Wang等,2009)。

然而,由于华北地下水已经处于超采状态,不可能建立现状开采条件下的稳定流模型。

本次研究中,采用0.2为模型初始含水量,从1960开始,运行模型至1992年,得到的含水量分布作为研究模拟期的初始含水量。

因为开采初期华北平原地下水埋深范围1-2 m,1960年-1992年模拟期间由于初始含水量带来的计算误差将会很快被修正。