CN109522655B - 基于变饱和水分运动系统的区域地下水补给量计算方法 - Google Patents

Abstract

本发明提供基于变饱和水分运动系统的区域地下水补给量计算方法，包括：步骤1.收集区域资料；步骤2.进行空间和时间的划分；步骤3.运行非饱和带水分运动模块，计算各个子区域的地下水补给量；步骤4.将该时间步内所有子区域的补给量输入饱和带地下水运动三维模型中，获得该时间步内地下水的水位；步骤5.根据推算出非饱和带厚度做为下一次计算的初始条件，并且进行收敛性判断，若判断结果为是，则获得该时间步内地下水补给量的真实值，然后进入下一个时间步的计算；否则返回步骤3开始该时间步的下一个迭代过程，直至满足收敛条件；步骤6.重复步骤3～5，获得所有时间步的地下水补给量真实值。本方法能够全面、准确地反映地下水补给量。

Description

基于变饱和水分运动系统的区域地下水补给量计算方法

技术领域

本发明属于地下水计算领域，具体涉及基于变饱和水分运动系统的区域地下水补给量计算方法。

技术背景

随着社会经济的发展，水资源供需矛盾日益突出。地下水作为我国北方地区的主要供水水源，对人民生活、工农业生产、城市建设以及社会经济发展都起着极其重要的作用。地下水补给量是地下水系统的重要组成部分，是区域地下水资源评价与开发利用的基础，只有准确定量地下水补给量，才能谈及地下水资源的合理开发和保护。

区域地下水补给量的计算方法有很多，根据非饱和带的土壤水分均衡分析获得地下水补给量是方法之一。非饱和带中的土壤水作为地下水的重要补给来源与消耗途径，通过对其水分运动过程的准确把握，可以较为准确的获取地下水的补给量。

目前，利用非饱和带土壤水分均衡计算地下水补给量的方法存在以下不足：(1)不考虑地下水的变动对于地下水补给量的影响。非饱和带的土壤水与饱和带的地下水是一个相互影响的系统，地下水的变化会影响非饱和带土壤水分的运动状况。在地下水位较深或者计算时间较短的情况下，这种影响较小，可以忽略。但是在其他情况下，仅通过非饱和带土壤水分运动计算地下水的补给量，难以从整体上把握土壤水-地下水运动系统，从而无法有效代表实际的地下水补给量。(2)非饱和带水分运动过程的处理过于简单，忽略非饱和带垂向向上的水分通量。且现有方法常忽略非饱和带的非均质性，并采用较多经验参数，因此导致现有方法难以准确代表非饱和带的土壤水分状况，从而影响地下水补给量的计算效果。

发明内容

本发明是为了解决上述问题而进行的，目的在于提供一种基于变饱和水分运动系统的区域地下水补给量计算方法，本方法能够全面地反映地下水补给量的时空变化，通过构造新的非饱和带水分运动计算模块，并且对饱和-非饱和水分运动系统进行整体求解，从而获得地下水补给量的准确计算结果。

本发明为了实现上述目的，采用了以下方案：

如图1所示，本发明提供了一种基于变饱和水分运动系统的区域地下水补给量计算方法，其特征在于，包括以下步骤：

步骤1.收集区域数字高程模型(DEM)数据、土地利用类型、气象资料以及水文地质资料，并进行数据清理和格式统一化处理；

步骤2.进行空间和时间的划分：

根据区域的尺度与研究的精度，结合不同的水文地质条件，将区域划分为多个子区域；并按照所需的时间精度，将计算时间划分为一个个时间步；

步骤3.运行非饱和带水分运动模块，计算地下水补给量：

对任一子区域，在一个时间步内，以地下水位为基准面，通过入渗水分配、重力作用下的土壤水分对流运动、源汇项影响、基质势作用下的土壤水分扩散运动来计算地下水的补给量：

R^p＝q_I+q_A-q_S-q_D，

式中，R为单位面积上该时间步内的地下水净补给量，上标p表示迭代次数，q_I为单位面积上通过入渗水分配计算得到的通过地下水位面的水分通量，q_A为单位面积上通过重力作用下的土壤水分对流作用计算得到的通过地下水位面的水分通量，q_S为单位面积上在源汇项影响下通过地下水位面的水分通量，q_D为单位面积上基质势作用下的土壤水分扩散运动通过地下水位面的水分通量；

基于上述地下水补给量计算方法，依次计算出各个子区域的R^p；

步骤4.将该时间步内所有子区域的R^p输入饱和带地下水运动三维模型中，计算获得在该时间步内地下水的水位

步骤5.根据

推算出非饱和带厚度du^p+1，并将该非饱和带厚度du^p+1设置为下一个计算的初始条件；

并根据以下公式进行收敛性判断：

式中，ε_H为容忍误差；

若判断结果为是，则本次迭代过程中步骤3计算出的各个子区域的R^p为该时间步内地下水补给量的真实值，然后进入下一个时间步的计算；否则返回至步骤3开始该时间步的下一个迭代过程计算非饱和带模块，直至满足该收敛条件；

步骤6.重复步骤3到步骤5，获得所有时间步的地下水补给量真实值。

本发明提供的基于变饱和水分运动系统的区域地下水补给量计算方法，还可以具有以下特征：在步骤2中，在水平方向上，将非饱和带划分为一个个子区域；在垂直方向上，将非饱和带划分为一个个相互串联的土层，以用于计算；饱和带则根据实际地质情况与饱和带模块的计算需求进行划分。

本发明提供的基于变饱和水分运动系统的区域地下水补给量计算方法，还可以具有以下特征：步骤3包括以下子步骤：

步骤3-1.计算入渗水分配：

q_I＝max(I-∑M×θ_de,0)，

其中I是单位面积的入渗水的总量；θ_de为某土壤层总的水分亏缺量，即每个土壤层饱和含水率与实际含水率之间的差值；M为该层土壤的厚度；即如果入渗水量大于非饱和带的水分亏缺，则其多余的水分通量即为q_I，q_I方向为垂直向下；

步骤3-2.计算重力作用下的土壤水分对流运动：

式中，K(θ)是非饱和水力传导度，采用经验公式描述，即

其中K_s为饱和水力传导度，θ_f为田间持水率，θ_s为饱和含水率；对于每层土壤而言，只有当土壤含水率大于田间持水率时，才会发生重力作用下的土壤水分对流运动，且对于该过程的计算为从上至下依次进行；对以上公式积分、化为递推公式并进行整理，得到通过地下水位面的通量q_A的计算公式为：

式中，θ为地下水位面以上的土壤层的含水率，M为地下水位面以上的土壤层的厚度。q_A方向为垂直向下；

步骤3-3.计算源汇项影响下通过地下水位面的水分通量：

q_S＝max(ET-ET_v,0)，

式中，ET为实际蒸散发量，ET_v为发生在非饱和带的实际蒸散发量；即如果非饱和带不能满足实际蒸发的需要，则需要消耗饱和带的地下水；q_S的方向为垂直向上；

步骤3-4.基质势作用下的土壤水分扩散运动的计算公式为，

式中，D(θ)为扩散度，

为描述土质垂向变异性的修正项，N代表描述土质变异性的参数个数，l表示第l个表示土质变异性的参数，参考van Genuchten模型，则描述土壤水分特征曲线的参数有4个，即θ_r，θ_s，α和n，且其中α的影响非常小，忽略不计，

除了非饱和带已有的土质分层，在地下水位面以下增设一个饱和的虚拟土壤层，加入离散方程中，将基质势作用下的土壤水分扩散运动方程采用有限差分法离散：

式中，

式中，下标i表示土壤分层，为计算开始前的步骤2中确定；上标t代表时间，t-1时刻为上一个计算时间的结果，为已知值；S为饱和度，S＝θ-θ_r/θ_s-θ_r；G为与n及θ有关的参数；i-1/2与i+1/2表示相邻两层(即i与i-1层，i与i+1层)土壤相关参数的均值，一般取其几何平均；采用追赶法计算方程组；

假设计算得到的虚拟层的含水率为

则基质势作用下的土壤水分扩散运动通过地下水位面的水分通量q_D为：

其中，M_x为虚拟层的厚度，q_D方向为垂直向上。

发明的作用与效果

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

1.本发明将土壤水-地下水流动系统作为一个整体进行分析与计算，进而通过在系统层面上对水分运动的精确描述，得到地下水的补给量。该方法弥补了地下水补给量计算中仅从局部考虑的不足，为地下水补给量的计算提供了新思路，并可以准确全面地反映下水补给量的时空变化情况。

2.本发明在对非饱和带水分运动过程的描述中，全面考虑了入渗水分配、重力作用下的土壤水分对流运动、源汇项影响、基质势作用下的土壤水分扩散运动，使得对非饱和带水分运动的描述更加完善准确，提高了地下水补给量的计算精度与可靠性，为地下水补给量的计算提供了新方法、新思路。

附图说明

图1为本发明实施例中涉及的基于变饱和水分运动系统的区域地下水补给量计算方法的流程图；

图2为本发明实施例中涉及的永联灌区地理位置(a)、土地类型划分(b)和DEM(c)的示意图；

图3为本发明实施例中涉及的变饱和水分运动系统空间示意图；

图4为本发明实施例步骤3中地下水补给量计算方法示意图；

图5为本发明实施例中涉及的永联灌区农地地下水补给量随时间的变化情况示意图。

具体实施方式

以下结合附图对本发明涉及的基于变饱和水分运动系统的区域地下水补给量计算方法进行详细地说明。

<实施例>

本实施例中，以中国内蒙古河套灌区永联灌区地下水补给量计算为例进行说明，如图1所示，本实施例所提供的基于变饱和水分运动系统的区域地下水补给量计算方法包括以下步骤：

步骤1.数据收集与整理：

收集区域DEM数据、土地利用类型、气象资料以及水文地质资料，并进行数据清理和格式统一化处理；图2中显示了内蒙古河套灌区永联灌域的图示、以及其土地利用类型、DEM等数据。

步骤2.进行空间和时间的划分：

根据区域的尺度与研究的精度，结合不同的水文地质条件，在水平方向上将非饱和带划分为一个个子区域；在垂直方向上，将非饱和带划分为一个个相互串联的土层，以用于计算。饱和带则根据实际地质情况与饱和带模块的计算需求进行划分。并按照所需的时间精度，将计算时间划分为一个个时间步。

本实施例中，如图2所示，根据土地利用类型将该区域在水平方向划分为3个子区域，分别为农地、村庄和荒地。垂向上将各子区域的非饱和带划分为一个个相互串联的土壤层，每层厚0.1m。饱和带根据实际地质情况进行划分。本算例中变饱和土壤水分模型空间示意情况见图3所示。将计算时间设置为5天一次。时间精度太小会导致计算量极大的增加，太大会使得计算结果过于粗糙。

步骤3.初始条件设置：

设置用于该时间步计算所需的条件，包括空间信息和时间信息。空间信息如步骤1中的DEM数据和水文地质参数，在第一次设置后保持不变。而同属于空间信息的非饱和带厚度则是随着计算过程而变动，其变动由饱和带与非饱和带的计算结果共同决定。时间信息如属于气象资料的降雨量，是根据实际观测结果而不断变化(随着时间的变化而变化)。

假设前一时刻为t-1，计算第t时刻的地下水补给量，时间步长为Δt。在一个计算时间步中，采用迭代计算地下水补给量，其目的是保证总体的质量平衡，即以下步骤4～步骤7为该计算时间步中的迭代计算，将迭代次数记为p。

假设本实施例中共划分出了K个子区域，对其中任一个子区域，设置其第p次迭代的非饱和带厚度du^p，则在每个时间步的迭代计算中，仅更新非饱和带的厚度du^p，详见后续步骤。

在步骤2中已将其划分为一个个垂向相互串联的土层，如图4所示，对其从上至下编号，第i个土层的厚度为M_i(i＝1,…,J)。

步骤4.运行非饱和带水分运动模块，计算地下水补给量：

对任一子区域，在一个时间步内，以地下水位为基准面，通过入渗水分配、重力作用下的土壤水分对流运动、源汇项影响、基质势作用下的土壤水分扩散运动来计算地下水的补给量：

R^p＝q_I+q_A-q_S-q_D，

式中，R为单位面积上该时间步内的地下水净补给量，上标p表示迭代次数，q_I为单位面积上通过入渗水分配计算得到的通过地下水位面的水分通量，q_A为单位面积上通过重力作用下的土壤水分对流作用计算得到的通过地下水位面的水分通量，q_S为单位面积上在源汇项影响下通过地下水位面的水分通量，q_D为单位面积上基质势作用下的土壤水分扩散运动通过地下水位面的水分通量；具体地，步骤4包括以下子步骤(采用“tipping-bucket”方法处理子步骤4-1和4-2)：

步骤4-1.计算入渗水分配：

q_I＝max(I-∑M×θ_de,0)，

其中I是单位面积的入渗水的总量；θ_de为某土壤层总的水分亏缺量，即每个土壤层饱和含水率与实际含水率之间的差值；M为该层土壤的厚度；即如果入渗水量大于非饱和带的水分亏缺，则其多余的水分通量即为q_I；

步骤4-2.计算重力作用下的土壤水分对流运动：

式中，K(θ)是非饱和水力传导度，

其中K_s为饱和水力传导度，θ_f为田间持水率，θ_s为饱和含水率；对于每层土壤而言，只有当土壤含水率大于田间持水率时，才会发生重力作用下的土壤水分对流运动，且对于该过程的计算为从上至下依次进行；对以上公式积分、化为递推公式并进行整理，得到通过地下水位面的通量q_A的计算公式为：

式中，θ为地下水位面以上的土壤层的含水率，M为地下水位面以上的土壤层的厚度；

步骤4-3.计算源汇项影响下通过地下水位面的水分通量：

q_S＝max(ET-ET_v,0)，

式中，ET为实际蒸散发量，ET_v为发生在非饱和带的实际蒸散发量；即如果非饱和带不能满足实际蒸发的需要，则需要消耗饱和带的地下水；q_S的方向为垂直向上；

步骤4-4.基质势作用下的土壤水分扩散运动的计算公式为，

式中，D(θ)为扩散度，

为描述土质垂向变异性的修正项，N代表描述土质变异性的参数个数，l表示第l个表示土质变异性的参数；参考van Genuchten模型，则描述土壤水分特征曲线的参数有4个，即θ_r，θ_s，α和n，且其中α的影响非常小，忽略不计；n采用如下经验方程描述，即

除了非饱和带已有的土质分层，在地下水位面以下增设一个饱和层单元，加入离散方程中，将基质势作用下的土壤水分扩散运动方程采用有限差分法离散：

式中，

式中，下标i表示土壤分层，为计算开始前的步骤2中确定；上标t代表时间，t-1时刻为上一个计算时间的结果，为已知值；S为饱和度，S＝θ-θ_r/θ_s-θ_r；G为与n及θ有关的参数；i-1/2与i+1/2表示相邻两层(即i与i-1层，i与i+1层)土壤相关参数的均值；采用追赶法计算方程组。

假设计算得到的虚拟单元的含水率为

则基质势作用下的土壤水分扩散运动通过地下水位面的水分通量q_D为：

其中，M_x为虚拟单元的厚度，q_D方向为垂直向上。

基于上述地下水补给量计算方法，依次计算出各个子区域的R^p，本实施例中子区域个数为K个，因此得到K个R^p。

如图3所示，q_I与q_A的方向是垂直向下的，而q_S与q_D的方向是垂直向上的。如图4所示，根据地下水位调整非饱和带的垂向分层，使之与实际的情况相适应。在本实施例中，共有3个子区域(K＝3)，所以本步骤应该重复3次，依次计算出各子区域在该时间步的地下水补给量。

步骤5.将该时间步内所有子区域的R^p输入饱和带地下水运动三维模型中，计算获得在该时间步内地下水的水位

本实施例中，将3个不同子区域的地下水补给量作为一个整体的数组R，传递给饱和带模块。R将作为饱和带模块的上边界条件，加入到饱和带模块进行计算。饱和带模块采用三维动力学模型，并不像非饱和带进行水平分区，而是作为整体统一计算。

步骤6.根据

推算出非饱和带厚度du^p+1，并将该非饱和带厚度du^p+1设置为下一次计算的初始条件；本实施例中，在饱和带模型计算完成后，得到的计算结果

即为在t时刻，第p次迭代的饱和带地下水水头分布情况。此时，重新计算非饱和带的厚度。非饱和带的厚度与地下水水头的和等于计算区域的总厚度。得到的新的非饱和带的厚度用于下次计算。

接着，根据以下公式进行收敛性判断：

式中，ε_H为容忍误差。

若判断结果为是，则将本次迭代过程中步骤4计算出的各个子区域的R^p为该时间步内地下水补给量的真实值，然后进入下一个时间步的计算；否则返回至步骤4开始该时间步的下一个迭代过程计算非饱和带模块，下一个迭代过程开始前应将非饱和带厚度更新为du^p+1；重复计算过程，直至满足该收敛条件。

值得强调的是，本方法中存在两个迭代循环，一个是与时间T有关的循环，每一次时间循环会进入新的时间步(T+ΔT)，然后计算该时间步中的地下水补给量。在时间循环中，还有一个迭代循环。进行迭代的目的是为了保证整个系统水量均衡。采用饱和带模型计算得到的地下水位作为迭代循环的判断指标，如果前后两次计算得到的地下水位差值的最大值在误差允许范围之内，那么就认为迭代收敛，此时进入新的时间步进行计算。

步骤7.重复上述步骤4到步骤6，获得所有时间步的地下水补给量真实值，进而获得地下水补给量随时间的变化情况。

如图5所示，为本实施例中通过上述方法计算得到的2014年10月至2016年10月农地的地下水补给量随时间的变化情况。每年的10月至次年的4月为非种植期，5月至9月为种植期。在每年的10月份，灌区存在大规模的秋浇，所以导致了此阶段会产生非常大的地下水补给量。而在非种植期的其他阶段，则变饱和带水分运动相对稳定。在作物种植季节，由于土壤蒸发作用和作物蒸腾作用的存在，会导致农地地下水的持续消耗，因而此阶段地下水补给量为负值，即消耗地下水。实际计算得永联灌区农地地下水年补给量为28mm/year。

综上，本实施例所提供的求解区域地下水补给量的方法，从土壤水-地下水运动系统的整体上考虑地下水的补给量的计算，弥补了以往方法中仅从局部考虑的不足，且更加完善的考虑非饱和带土壤水分的运动过程，使得该方法可以较为准确的获取地下水的补给量，为区域上地下水资源的评价与开发利用提供参考。

以上实施例仅仅是对本发明技术方案所做的举例说明。本发明所涉及的基于变饱和水分运动系统的区域地下水补给量计算方法，包括：变饱和水分运动系统的构建、拟三维迭代计算方法、非饱和带均衡计算方法、地下水补给量计算方法和其各分量的描述，并不仅仅限定于在以上实施例中所描述的内容，而是以权利要求所限定的范围为准。本发明所属领域技术人员在该实施例的基础上所做的任何修改或补充或等效替换，都在本发明的权利要求所要求保护的范围内。

Claims (3)

1.基于变饱和水分运动系统的区域地下水补给量计算方法，其特征在于，包括以下步骤：

步骤1.收集区域数字高程模型数据、土地利用类型、气象资料以及水文地质资料；

步骤2.进行空间和时间的划分：

根据区域的尺度与研究的精度，结合不同的水文地质条件，将区域划分为多个子区域；并按照所需的时间精度，将计算时间划分为一个个时间步；

步骤3.运行非饱和带水分运动模块，计算地下水补给量：

对任一子区域，在一个时间步内，以地下水位为基准面，通过入渗水分配、重力作用下的土壤水分对流运动、源汇项影响、基质势作用下的土壤水分扩散运动来计算地下水的补给量：

R^p＝q_I+q_A-q_S-q_D，

式中，R为单位面积上该时间步内的地下水净补给量，上标p表示迭代次数，q_I为单位面积上通过入渗水分配计算得到的通过地下水位面的水分通量，q_A为单位面积上通过重力作用下的土壤水分对流作用计算得到的通过地下水位面的水分通量，q_S为单位面积上在源汇项影响下通过地下水位面的水分通量，q_D为单位面积上基质势作用下的土壤水分扩散运动通过地下水位面的水分通量；

基于上述地下水补给量计算方法，依次计算出各个子区域的R^p；

步骤4.将该时间步内所有子区域的R^p输入饱和带地下水运动三维模型中，计算获得在该时间步内地下水的水位

步骤5.根据

推算出非饱和带厚度du^p+1，并将该非饱和带厚度du^p+1设置为下一个计算的初始条件；

并根据以下公式进行收敛性判断：

式中，ε_H为容忍误差；

若判断结果为是，则将本次迭代过程中步骤3计算出的各个子区域的R^p为该时间步内地下水补给量的真实值，然后进入下一个时间步的计算；否则返回至步骤3开始该时间步的下一个迭代过程计算非饱和带模块，直至满足该收敛条件；

步骤6.重复步骤3到步骤5，获得所有时间步的地下水补给量真实值。

2.根据权利要求1所述的基于变饱和水分运动系统的区域地下水补给量计算方法，其特征在于：

其中，在步骤2中，在水平方向上，将非饱和带划分为一个个子区域；在垂直方向上，将非饱和带划分为一个个相互串联的土层，以用于计算；饱和带则根据实际地质情况与饱和带模块的计算需求进行划分。

3.根据权利要求1所述的基于变饱和水分运动系统的区域地下水补给量计算方法，其特征在于：

其中，步骤3包括以下子步骤：

步骤3-1.计算入渗水分配：

q_I＝max(I-∑M×θ_de,0)，

其中I是单位面积的入渗水的总量，θ_de为土壤层总的水分亏缺量，M为该层土壤的厚度，q_I方向为垂直向下；

步骤3-2.计算重力作用下的土壤水分对流运动：

式中，K(θ)是非饱和水力传导度，

其中K_s为饱和水力传导度，θ_f为田间持水率，θ_s为饱和含水率；对以上公式积分、化为递推公式并进行整理，得到通过地下水位面的通量q_A的计算公式为：

式中，θ为地下水位面以上的土壤层的含水率，M为地下水位面以上的土壤层的厚度，q_A方向为垂直向下，Δt为时间步长；

步骤3-3.计算源汇项影响下通过地下水位面的水分通量：

q_S＝max(ET-ET_v,0)，

式中，ET为实际蒸散发量，ET_v为发生在非饱和带的实际蒸散发量，q_s的方向为垂直向上；

步骤3-4.基质势作用下的土壤水分扩散运动的计算公式为，

式中，D(θ)为扩散度，

为描述土质垂向变异性的修正项，N代表描述土质变异性的参数个数，l表示第l个表示土质变异性的参数，

除了非饱和带已有的土质分层，在地下水位面以下增设一个饱和的虚拟土壤层，加入离散方程中，将基质势作用下的土壤水分扩散运动方程采用有限差分法离散：

式中，

式中，下标i表示土壤分层，为计算开始前的步骤2中确定；上标t代表时间，t-1时刻为上一个计算时间的结果，为已知值；S为饱和度，S＝θ-θ_r/θ_s-θ_r；G为与n及θ有关的参数；i-1/2与i+1/2表示相邻两层土壤相关参数的均值；

假设计算得到的虚拟层的含水率为

则基质势作用下的土壤水分扩散运动通过地下水位面的水分通量q_D为：

其中，M_x为虚拟层的厚度，q_D方向为垂直向上。

Images