地表水-地下水耦合模拟的计算方法

摘要

本发明公开了一种地表水-地下水耦合模拟的计算方法，其特征在于：包括产流的计算方法：产流即为径流量R，采用指数型蓄水容量曲线进行计算；还包括地下水蓄水库演算方法及基流量Qg的计算方法，地下水埋深D及地下水埋深D符合Gamma分布时基流量Qg的计算方法和流域汇流的计算方法。本发明提供的一种地表水-地下水耦合模拟的计算方法，改进了流域水文模型，建立考虑地下水埋深统计分布特征的降水入渗补给量、潜水蒸发量计算方法以及河川径流量与地下水埋深的计算方法，实现了不同灌溉开采量影响下的地下水位动态变化模拟，以及流域降雨-径流响应分析，并为研究地下水位变化对地表-地下水转化影响提供了研究方法。

说明书

技术领域

本发明涉及一种地表水-地下水耦合模拟的计算方法，尤其涉及一种改进 了流域水文模型的地表水-地下水耦合模拟的计算方法。

背景技术

在极易形成暴雨，降水量年际变化很大的地区，流域主要是雨洪径流，在 多雨年份，汛期雨量集中，洪水暴涨暴落，极易造成灾害；还有，不同地貌对 地下水的利用情况也不同。

本发明为研究地下水开采以及地下水位下降对蒸发和径流的影响，针对淮 河流域降雨-径流关系、径流和基流与地下水位关系，以及地下水位时空变化 的统计特征，本发明改进了流域水文模型，建立考虑地下水埋深统计分布特征 的降水入渗补给量、潜水蒸发量计算方法以及河川径流量与地下水埋深的计算 方法，实现了不同灌溉开采量影响下的地下水位动态变化模拟，以及流域降雨 -径流响应分析，并为研究地下水位变化对地表水-地下水转化影响提供了研究 方法。

发明内容

本发明所要解决的技术问题是，将产流概化为径流量，提供一种流域产流 的计算方法，揭示了降雨-径流关系；进一步地，本发明提供一种地下水蓄水 库演算方法及基流量Qg的计算方法，揭示了径流和基流与地下水位关系，建 立了考虑地下水埋深统计分布特征的降水入渗补给量、潜水蒸发量计算方法以 及河川径流量的计算方法，实现了不同灌溉开采量影响下的地下水位动态变化 模拟；更进一步地，本发明提供一种地下水埋深D及地下水埋深D符合Gamma 分布时基流量Qg的计算方法，揭示了地下水位时空变化的统计特征，建立了 地下水埋深的计算方法；更进一步地，本发明提供一种流域汇流的计算方法， 通过构建适用于各子流域以及嵌套子流域地表水与地下水耦合的水文模型，研 究不同土地利用条件下地下水动态变化及其对降雨-径流过程响应的影响。

为解决上述技术问题，本发明采用的技术方案为：

地表水-地下水耦合模拟的计算方法，其特征在于：包括产流的计算方法： 产流即为径流量R，采用指数型蓄水容量曲线进行计算：

PE_t＝P_t+Wg-E_t  (1)

式(1)中，PE为净雨量；PE_t为t时刻的净雨量；P为时段降雨量；P_t为t时刻的降雨量；Wg为灌溉开采量，单位为mm；E为土壤蒸发量，采用三 层蒸发模式；E_t为t时刻的土壤蒸发量；

一般灌溉开采量Wg不产流，但部分回归地下水；

当PE_t≤0，则R_t＝0  (2)

当PE_t＞0，若PE_t+A_t≥WMM，则R_t＝W_t+PE_t-WM  (3)

当PE_t＞0，若PE_t+A_t＜WMM，则R_t＝W_t+PE_t-W_t+1  (4)

W_t+1＝W_t+PE_t+Eg_t  (5)

式(2)-(5)中，Eg为潜水蒸发量；Eg_t为t时刻的潜水蒸发量；W_t和 W_t+1分为t和t+1时刻流域平均蓄水量；A为流域单点最大蓄水量，A_t为t时刻 的流域单点最大蓄水量；R为径流量；R_t为t时刻的径流量；WMM为流域单 点最大蓄水容量，WM为流域最大平均蓄水容量；

所述流域单点即为流域单元，即将所述流域划分为若干流域单点。

还包括地下水蓄水库演算方法及基流量Qg的计算方法：

a).含水层水平衡及补排项计算方法：

当浅层含水层概化为地下水蓄水库，其水量平衡公式为：

SG_t+1＝SG_t-Eg_t+Pr_t-Wg_t-Qg_t  (6)

式(6)中，SG_t、SG_t+1分别为t、t+1时刻流域平均地下水库蓄水量，Eg_t为t时刻的潜水蒸发量；Pr_t为t时刻的入渗补给量；Wg_t为t时刻的灌溉开采 量；Qg_t为地下水蓄水库t时刻的出流量，即基流量；

b).入渗补给量Pr的计算方法：

流域内某一点地下水埋深为D时，径流量R进入地下水蓄水库的入渗补 给量Pr用下列经验公式计算。

若P+Wg-E＞0，假设单元i入渗补给量与土壤含水量成正比，

即Pr_i＝α(P+Wg-E)(W_i/WMM)  (7)

则在流域面上积分得到：

$P_r={\int }_0^1{P}_{\mathrm{ri}}\mathrm{da}=\alpha (P+\mathrm{Wg}-E)(1-{(1-\frac{A_0}{\mathrm{WMM}})}^b+(1-\frac{W_0}{\mathrm{WM}})+\frac{A_0}{\mathrm{WMM}}{(1-\frac{A_0}{\mathrm{WMM}})}^b)---\left(8\right)$

若P+Wg-E＜0，则Pr＝0(9)

式(7)-(9)中，W_i为fi/F部分流域的蓄水量，其中fi为土壤蓄水饱和 时的面积，F为流域总面积；A₀为初始时刻单点最大蓄水量；W₀为初始时刻流 域平均蓄水量；WM为流域最大平均蓄水容量；α，b均为参数；

fi为土壤蓄水饱和时的面积，再产生降水就会流走形成径流；fi/F为已经 蓄水满的面积与总面积之比；

将每一个流域单点再划分为若干单元，单元i为若干单元中的任意一个；

c).潜水蒸发量Eg的计算方法：

Eg＝eta×EP×(1-D/D_max)ⁿ  (10)

式(10)中，eta为蒸发折算系数；D_max为潜水蒸发极限埋深；EP为同气 温条件下水面蒸发量；D为流域内某一点地下水埋深；

d).地下水灌溉开采量Wg的估算方法：

由于待研究区地下水开采量主要用于农业灌溉，故根据灌溉需水量W_guan乘以一个开采系数q_cai来计算灌溉开采量Wg及其变化过程；

Wg＝q_cai×W_guan  (11)

e).地下水蓄水库出流量Qg的计算方法：

Qg＝k×SG  (12)

SG＝Sy×(Dz-D)  (13)

式(12)-(13)中，Sy为给水度；Dz为基流量接近于0时的埋深，相当 于河流切割深度，Dz的取值范围在2～4m之间；SG为流域平均地下水库蓄水 量；k为渗透系数。

还包括地下水埋深D及地下水埋深D符合Gamma分布时基流量Qg的计 算方法：

地下水埋深D空间分布服从Gamma分布：

$f_D\left(d\right)=\frac{{\alpha }^{\gamma }}{\Gamma \left(\gamma \right)}{d}^{\gamma -1}{e}^{-\mathrm{\alpha d}},\alpha >0,\gamma >0---\left(14\right)$

式(14)中，Γ(γ)为Gamma函数；α和γ为参数，γ为形状系数；d代表 含水层厚度；

当地下水埋深D服从Gamma分布时，流域平均地下水埋深为Gamma分 布的数学期望；

$\bar{D}=E\left(d\right)=\frac{\gamma }{\alpha }---\left(15\right)$

$\alpha =\frac{\gamma }{\bar{D}}---\left(16\right)$

假设地下水埋深分布函数中形状系数γ不变，则α随平均地下水埋深变 化而变化；

假设最大平均地下水埋深D_m时的参数为α，则某一地下水平均埋深为时，参数α₀：

${\alpha }_0=\alpha \frac{\overline{D_m}}{\bar{D}}---\left(17\right)$

地下水埋深服从Gamma分布时，基流量Qg、潜水蒸发量Eg的计算方法 如下：

$\left(\begin{array}{c}\mathrm{Qg}={\int }_0^{\mathrm{Dz}}{f}_D\left(D\right)\mathrm{Kg}*\mathrm{Sy}*(\mathrm{Dz}-D)\mathrm{dD}\\ \approx \underset{i=1}{\overset{m}{\Sigma }}\frac{{(\alpha *\frac{D_m}{\bar{D}})}^{\gamma }}{\Gamma \left(\gamma \right)}{\left(\frac{\mathrm{Dz}*i}{m}\right)}^{\gamma -1}{e}^{-\alpha *\frac{D_m}{\bar{D}}\left(\frac{\mathrm{Dz}*i}{m}\right)}*\mathrm{Kg}*\mathrm{Sy}*(\mathrm{Dz}-\frac{\mathrm{Dz}*i}{m})*\frac{\mathrm{Dz}}{m}\end{array}\right)---\left(18\right)$

$\mathrm{Eg}\approx \underset{i=1}{\overset{m}{\Sigma }}\frac{{(\alpha *\frac{D_m}{\bar{D}})}^{\gamma }}{\Gamma \left(\gamma \right)}{\left(\frac{D_{\max }*i}{m}\right)}^{\gamma -1}{e}^{-\alpha *\frac{D_m}{\bar{D}}*\left(\frac{D_{\max }*i}{m}\right)}*\mathrm{eta}*\mathrm{EP}*{(1-\frac{D_{\max }*i}{m*D_{\max }})}^n*\frac{D_{\max }}{m}---\left(19\right)$

m代表把地下含水层分为m层，其中任意一层为用i表示。

流域汇流的计算方法为：流域径流量R一部分进入地下水蓄水库，再按式 (18)排泄至河道；另一部分按自由蓄水库曲线划分为地表径流Rs和壤中流 Ri；地表径流Rs直接进入河道，壤中流Ri经过线性水库调蓄进入河道；三者 汇入河道后经过马斯京根汇流流到出口断面。

待研究区的流域中，下游流域多为与上游流域相嵌套的流域，上游流域与 下游流域衔接的出口即为出口断面；因此这些河道水流演算需要考虑上游出口 断面入流以及傍测入流的演算过程；上游流域模拟的出流量作为下游流域的入 流量，再考虑下游流域产汇流进行连续演算；所述出流量即为径流量R。

傍测入流即河流与地下水交换量，包括潜水蒸发量Eg、灌溉开采量Wg 和入渗补给量Pr。

水文气象实测资料的测量时间间隔为逐日；在模拟的逐日径流量基础上， 统计逐月径流量作为目标函数；模拟了各流域每5日平均地下水埋深。

本发明提供的一种地表水-地下水耦合模拟的计算方法，改进了流域水文 模型，建立考虑地下水埋深统计分布特征的降水入渗补给量、潜水蒸发量计算 方法以及河川径流量与地下水埋深的计算方法，实现了不同灌溉开采量影响下 的地下水位动态变化模拟，以及流域降雨-径流响应分析，并为研究地下水位 变化对地表-地下水转化影响提供了研究方法。

附图说明

图1为本发明产流计算示意图；图中ΔW为时段降雨量P后W₀的增加量。

具体实施方式

下面结合附图对本发明作更进一步的说明。

如图1所示，本发明以颖涡地区为例：

沙颍河是淮河的最大支流，地处河南省腹地，发源于河南省伏牛山脉摩天 岭(又名没大岭)东麓，东南流经鲁山、平顶山、叶县、漯河、周口、项城、 沈丘等县市，至界首县城关镇附近进入安徽省，往下经太和、阜阳，于颍上县 沫河口入淮河。全长619km，流域面积近39880km2，其中山区面积9070km2， 丘陵区面积5370km2，平原面积20000km2，干流长418km。较大的支流有北 汝河、澧河、颍河、贾鲁河、新运河、汾泉河和黑茨河等8条。流域内有耕地 3180万亩，2400万人，有丰富的煤炭资源，是我国重要的能源基地，工农业 生产发展前景广阔。阜阳闸水文站为颍河中下游的闸坝控制站。阜阳站以上有 沙河、贾鲁河、茨河、泉河等众多支流汇入颍河，流域面积38240km2。流域 呈大陆性季风气候，气象变化受季风影响，为南北气候的过渡地带。在汛期(6～ 9月)，由于东南暖湿气流内移，加之西部地形影响，极易形成暴雨，为河南省 暴雨中心地区之一。年平均降雨量西部山区为800～1000mm，东部平原为 700～900mm，降雨量一般集中在汛期，占年降水量的60％以上。降水量年际 变化很大，最大值和最小值可达5倍。因此，流域主要是雨洪径流，在多雨年 份，汛期雨量集中，洪水暴涨暴落，极易造成灾害。

涡河为淮河第二大支流，淮北平原区河道。发源于河南省尉氏县，东南流 经开封、通许、扶沟、太康、鹿邑和安徽省亳州、涡阳、蒙城，于怀远县城附 近注入淮河。长380km，流域面积1.59万km2。涡河流域属温暖带半湿润大 陆性季风气候区。冬春干旱少雨，夏秋季太平洋副热带高压增强，降水量集中， 易造成洪涝。非汛期(1～6月及10～12月)9个月份涡河月平均最大流量 9.1m3/s，汛期(7～9)三个月涡河月平均最大流量为149m3/s，年最大径硫总 量为8.25亿m3。多年平均降雨量600～900mm。6～9月多年平均降雨量占全 年降水量的70％，多年平均无雨天数为260天，多年平均气温14.5℃，各月平 均气温以1月份最低，为-0.1℃，7月份最高，为38.0℃。涡河流域内地势由 西北向东南倾斜，地面高程为55～40m，地面坡降1/4500，河床宽一般40～ 100m，局部较窄，水深一般1～3m，河床比降为1/6000～1/9000。堤内近河堤 处由于筑堤取土等原因，多分布有坑塘、洼地。涡河河道弯曲，河流沿岸为河 谷地貌形态。

本发明主要选取颍河阜阳闸以上流域和涡河玄武闸以上流域进行研究。将 阜阳闸以上颍河流域划分为12个子流域(区间)，加上涡河玄武闸以上子流域 共计13个流域，通过构建适用于各子流域以及嵌套子流域地表水与地下水耦 合的水文模型，研究不同土地利用条件下地下水动态变化及其对降雨-径流过 程响应的影响。

颍涡地区降雨-径流模拟：本发明针对颍河阜阳闸以上区域和涡河玄武闸 以上区域，将阜阳闸以上控制区域划分为12个子流域，各个子流域从上游到 下游互相嵌套，上游模拟的出流量作为下游流域入流量，再考虑下游流域产汇 流进行连续演算，对每个子流域用本发明建立的流域水文模型进行产汇流模 拟。选取1997-2011年逐日水文气象实测资料进行模型计算，其中1997-2006 年作为模型参数率定期，2007-2011年作为模型验证期。首先通过模型计算的 河流下断面流量与实测流量数据对比，优化模型参数，验证模型模拟结果的可 靠性。再通过对比模拟与实测流域平均地下水埋深，以及模拟与分割的基流量， 进一步验证模型模拟结果的可靠性。

河川径流量模拟：考虑到研究区河流多有闸坝控制，且本研究目的是地下 水位变化对设计径流量影响，本发明在模拟的逐日流量基础上，统计逐月流量 作为目标函数。如表1所示，率定期实测与模拟流量效率系数NSC在0.32～0.91， 大多数流域率定期和验证期NSC在0.7以上，表明模型能较好地模拟流量过程。

地下水埋深模拟：为了对模型进一步验证，模型同时模拟了各流域每5日 面平均地下水埋深，并与根据每5日实测地下水位资料计算的流域面平均地下 水埋深进行对比。如表1所示，有地下水位观测资料的流域模拟与实测多年平 均埋深接近，相关性系数R²大都在0.7以上，表明模型能较好地模拟地下水位 动态变化。

河道基流量模拟：模型模拟了12个子流域基流量，如表1所示，模拟的 多年平均基流量与分割的基流量较为接近。

模型通过在颍涡研究区各子流域的应用，分别模拟了径流量、地下水埋深 及基流量，取得了较好的模拟效果，说明该模型在颍涡地区具有较好的适用性。

表1.颍涡地区模型模拟统计结果

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通 技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰， 这些改进和润饰也应视为本发明的保护范围。