一种同时确定弱透水层水文地质参数和非达西渗流控制参数的方法与流程

1.本发明属于水文地质参数确定技术领域，尤其涉及一种同时确定弱透水层水文地质参数和非达西渗流控制参数的方法。


背景技术：

2.弱透水层是含水层系统的重要组成部分，黏土、亚黏土等是弱透水层的典型组成成分，其在冲积平原地区广泛分布。在中国东部地区，由含水层和弱透水层组成的地下水系统十分常见。研究表明，弱透水层在地面沉降当中起到重要作用。除此之外，弱透水层在水文地球化学领域也有重要的作用。最近几年的研究中发现，弱透水层在解决地下水渗漏和地下水水质污染的相关问题中存在着不可替代的作用，因此，确定弱透水层的水文地质参数，如渗透系数和贮水率，将为地下水资源评价、地面沉降、地下水污染等多个领域提供不可或缺的帮助。
3.虽然对于弱透水层水文地质参数确定方法已有不少研究，但其基本都认为弱透水层渗流遵循达西定律。然而，实际情况却是弱透水层渗流并不完全为达西渗流，故以达西渗流为前提进行参数确定得到的结果可能具有较大误差。因此，亟需开发一种考虑非达西渗流情况下的弱透水层水文地质参数确定方法。


技术实现要素：

4.发明目的：针对以上问题，本发明提出一种同时确定弱透水层水文地质参数和非达西渗流控制参数的方法，构建了考虑幂函数形式非达西渗流的一维弱透水层水动力响应模型，获得了弱透水层一侧水头瞬时降低某一常量并保持不变条件下弱透水层累计释水变形量的数值解，而后利用双对数坐标系中变形-时间曲线的前段斜率α提出了一种全新的配线法，可用于同时确定弱透水层渗透系数k和贮水率ss以及α的值。本发明理论严密，方法简单，容易掌握，有很好的推广应用价值。
5.技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种同时确定弱透水层水文地质参数和非达西渗流控制参数的方法，包括以下步骤：
6.步骤一，构建考虑幂函数形式非达西渗流的一维弱透水层水动力响应模型；
7.步骤二，通过有限差分方法计算获得弱透水层一侧水头瞬时降低某一常量并保持不变条件下弱透水层累计释水变形量的数值解；
8.步骤三，获取双对数坐标系中不同非达西渗流控制参数即非达西指数n和临界水力梯度i
l
组合条件下的弱透水层累计变形量标准曲线；
9.步骤四，通过标准曲线前段斜率α来表征弱透水层非达西渗流状态，并给出α与非达西指数n和临界水力梯度i
l
的数据参照表；
10.步骤五，根据参照表确定观测数据的理论变形曲线，利用前段斜率α的配线法，同时确定观测资料中弱透水层渗透系数k和贮水率ss的值。
11.进一步的，步骤一所述一维弱透水层水动力响应模型从下到上依次为下覆含水层、弱透水层、上覆含水层，满足以下条件：
12.1)弱透水层是均质的，且渗透系数k、贮水率ss、非达西指数n、临界水力梯度i
l
不随时间变化；
13.2)初始时刻弱透水层中各点水头相同，即弱透水层中任一点位置z在t＝0时刻的降深s(z，t)＝0；
14.3)弱透水层中的水流为垂向一维流；
15.4)在初始时刻，弱透水层上覆含水层水位始终保持恒定，下覆含水层水位瞬时降低后保持恒定；
16.5)弱透水层在水力响应过程中始终饱和。
17.进一步的，取垂直向上的一维坐标系，且零点位于弱透水层的下边界，一维弱透水层中任一点位置z在t时刻的降深s服从的水动力响应控制方程为
[0018][0019]
式中，v为t时刻位置z处的地下水流速，ss为贮水率；在水力梯度i小于其阈值即临界水力梯度i
l
时发生非达西渗流，则流速v用i的指数函数表示为
[0020][0021]
其中n为非达西指数，iu为单位水力梯度；当n＞1时，弱透水层中的渗流为非达西流，当n＝1时，弱透水层中的渗流为达西渗流；
[0022]
将式(2)代入式(1)，得到一维弱透水层水动力响应模型控制方程：
[0023][0024]
其中，
[0025][0026]
初始条件为
[0027]
s(z，0)＝0
ꢀꢀꢀ
(5)
[0028]
边界条件为
[0029][0030]
s(b，t)＝0
ꢀꢀꢀ
(7)
[0031]
式中，为初始时刻下覆含水层水位瞬时降低量，b为弱透水层厚度。
[0032]
进一步的，步骤二中，一维弱透水层水动力响应模型数值解计算方法如下：
[0033]
弱透水层累计释水变形量为
[0034][0035]
引入无量纲变量
[0036]
[0037]
td，zd，sd，和fd分别为无量纲形式的时间、位置、弱透水层水位降深、下覆含水层的瞬时降深以及弱透水层累计变形量；
[0038]
根据无量纲化的模型表示为
[0039][0040]
sd(zd，0)＝0
ꢀꢀꢀ
(11)
[0041][0042]
sd(1，td)＝0
ꢀꢀꢀ
(13)
[0043][0044]
若n＝1，则θ＝1且与i
l
不相关，弱透水层中的渗流为达西渗流，则式(14)简化为
[0045][0046]
此后，使用一维有限差分方法，将无量纲空间域[0，1]离散为n个元素，n 1个节点，将无量纲时间域[0，t
fd
]离散为m个网格区域，t
fd
为最终模拟时间；
[0047]
通过该方法，得出每一时空节点的水力梯度并通过式(4)求得θ，最终求出离散的弱透水层水位降深数据；然后通过式(14)求出弱透水层的累计变形fd的数值解。
[0048]
进一步的，步骤五中，根据参照表确定观测数据的理论变形曲线，利用前段斜率α的配线法，同时确定观测资料中弱透水层渗透系数k和贮水率ss的值，具体为：
[0049]
第一，获取累计变形量f随时间t的观测数据，并获取其前段直线部分，即双对数坐标系下f-t曲线初始阶段呈线性增长的部分；
[0050]
第二，从前段直线部分中随机读取两个坐标点位(t1，f1)和(t2，f2)并计算其斜率
[0051][0052]
第三，根据数据参照表以及计算出的斜率α，得出其对应的非达西指数n和临界水力梯度i
l
；
[0053]
第四，通过式(14)计算得出该n和i
l
下的理论变形曲线；在相同模数的双对数坐标系内，获得的实测f～t曲线与fd～td标准曲线在形状上相同；
[0054]
利用配线法确定k和ss的具体步骤如下：
[0055]
1)在双对数坐标系上绘制已确定α和的fd～td标准曲线；
[0056]
2)在另一个模数相同的坐标系中绘制实测f～t曲线；
[0057]
3)将实测曲线置于标准曲线上，在保持对应坐标轴彼此平行的条件下相对平移，直至实测f～t曲线与标准曲线重合；
[0058]
4)任意选取一匹配点，记录匹配点相应坐标值fd，td以及f，t；
[0059]
5)将相应坐标值代入式(19)和(20)，
[0060][0061]
[0062]
求取参数k和ss。
[0063]
有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：
[0064]
本发明所提供的一种同时弱透水层水文地质参数和非达西渗流控制参数的方法，考虑弱透水层释水变形过程中的非达西渗流，相较传统只考虑达西渗流的模型更接近实际，理论严密，结果更加精确。采用对比前段斜率α的新配线法，操作简单，易于上手。一次实验即可获得渗透系数、贮水率、前段斜率α、指数n和临界水力梯度il，获得的参数多，使用方便。
附图说明
[0065]
图1为本发明的实施流程图；
[0066]
图2为本发明一维弱透水层水动力响应模型概念图；
[0067]
图3为本发明双对数坐标系下不同n和i
l
取值下的fd～td标准曲线图；
[0068]
图4为本发明前段斜率α的新配线法求取水文地质参数原理图。
具体实施方式
[0069]
下面结合具体实例，进一步阐述本发明，应理解实施例仅应用于阐述本发明并不用于限制本发明的范围，在阅读了本发明后，本领域技术人员对本发明各种等价形式的修改都应落于本技术所附权利要求所限定范围内。
[0070]
本发明所述的同时确定弱透水层水文地质参数和非达西渗流控制参数的方法，如图1所示，包括以下步骤：
[0071]
步骤一，构建考虑幂函数形式非达西渗流的一维弱透水层水动力响应模型；如图2所示，所述一维弱透水层水动力响应模型从下到上依次为下覆含水层、弱透水层、上覆含水层，满足以下条件：
[0072]
1)弱透水层是均质的，且渗透系数k、贮水率ss、非达西指数n、临界水力梯度i
l
不随时间变化；
[0073]
2)初始时刻弱透水层中各点水头相同，即弱透水层中任一点位置z在t＝0时刻的降深s(z，t)＝0；
[0074]
3)弱透水层中的水流为垂向一维流；
[0075]
4)在初始时刻，弱透水层上覆含水层水位始终保持恒定，下覆含水层水位瞬时降低后保持恒定；
[0076]
5)弱透水层在水力响应过程中始终饱和。
[0077]
取垂直向上的一维坐标系，且零点位于弱透水层的下边界，一维弱透水层中任一点位置z在t时刻的降深s服从的水动力响应控制方程为
[0078][0079]
式中，v为t时刻位置z处的地下水流速，ss为贮水率；在水力梯度i小于其阈值即临界水力梯度i
l
时发生非达西渗流，则流速v用i的指数函数表示为
[0080][0081]
其中n为非达西指数，iu为单位水力梯度；应当注意n表征渗流的非线性程度，n越大，非线性程度越强，非达西流影响越明显，因此，当n＞1时，弱透水层中的渗流为非达西流，当n＝1时，弱透水层中的渗流为达西渗流；
[0082]
将式(2)代入式(1)，得到一维弱透水层水动力响应模型控制方程：
[0083][0084]
其中，
[0085][0086]
初始条件为
[0087]
s(z，0)＝0
ꢀꢀꢀ
(5)
[0088]
边界条件为
[0089][0090]
s(b，t)＝0
ꢀꢀꢀ
(7)
[0091]
式中，为初始时刻下覆含水层水位瞬时降低量，b为弱透水层厚度。
[0092]
步骤二，通过有限差分方法计算获得弱透水层一侧水头瞬时降低某一常量并保持不变条件下弱透水层累计释水变形量的数值解；计算方法如下：
[0093]
弱透水层累计释水变形量为
[0094][0095]
引入无量纲变量
[0096][0097]
td，zd，sd，和fd分别为无量纲形式的时间、位置、弱透水层水位降深、下覆含水层的瞬时降深以及弱透水层累计变形量；
[0098]
根据无量纲化的模型表示为
[0099][0100]
sd(zd，0)＝0
ꢀꢀꢀ
(11)
[0101][0102]
sd(1，td)＝0
ꢀꢀꢀ
(13)
[0103][0104]
若n＝1，则θ＝1且与i
l
不相关，弱透水层中的渗流为达西渗流，则式(14)简化为
[0105]
[0106]
此后，使用一维有限差分方法，将无量纲空间域[0，1]离散为n个元素，n 1个节点，将无量纲时间域[0，t
fd
]离散为m个网格区域，t
fd
为最终模拟时间；
[0107]
通过该方法，得出每一时空节点的水力梯度并通过式(4)求得θ，最终求出离散的弱透水层水位降深数据；然后通过式(14)求出弱透水层的累计变形fd的数值解。
[0108]
步骤三，获取双对数坐标系中不同非达西渗流控制参数即非达西指数n和临界水力梯度i
l
组合条件下的弱透水层累计变形量标准曲线。
[0109]
步骤四，通过标准曲线前段斜率α来表征弱透水层非达西渗流状态，并给出α与非达西指数n和临界水力梯度i
l
的数据参照表；本实施例所得数据参照见表1：
[0110]
表1前段斜率α与指数n和临界水力梯度i
l
的数据参照表
[0111][0112]
步骤五，根据参照表确定观测数据的理论变形曲线，利用前段斜率α的配线法，同时确定观测资料中弱透水层渗透系数k和贮水率ss的值；具体为：
[0113]
第一，获取场地或室内试验下累计变形量f随时间t的观测数据，并获取其前段直线部分，即双对数坐标系下f-t曲线初始阶段呈线性增长的部分；
[0114]
第二，从前段直线部分中随机读取两个坐标点位(t1，f1)和(t2，f2)并计算其斜率
[0115]
第三，根据数据参照表以及计算出的斜率α，得出其对应的非达西指数n和临界水力梯度i
l
；
[0116]
第四，通过式(14)计算得出该n和i
l
下的理论变形曲线；
[0117]
对td和fd同时取对数，可得
[0118]
log
10
td＝log
10
t log
10
[k/(ssb2)]
ꢀꢀꢀ
(17)
[0119][0120]
式(17)和(18)右边第二项都是常数，因此在相同模数的双对数坐标系内，试验获得的实测f～t曲线与fd～td标准曲线在形状上是相同的，只是横纵坐标分别平移了k/(ssb2)和因此，可以根据配线法确定待定常数k和ss。
[0121]
利用配线法确定k和ss的具体步骤如下：
[0122]
1)在双对数坐标系上绘制已确定α和的fd～td标准曲线；
[0123]
2)在另一个模数相同的坐标系中绘制实测f～t曲线；
[0124]
3)将实测曲线置于标准曲线上，在保持对应坐标轴彼此平行的条件下相对平移，直至实测f～t曲线与标准曲线重合；
[0125]
4)任意选取一匹配点，记录匹配点相应坐标值fd，td以及f,t；
[0126]
5)将相应坐标值代入式(19)和(20)，
[0127][0128][0129]
求取参数k和ss。
[0130]
本实施例具体实验如下：
[0131]
(1)实验模型
[0132]
其实验在一个高b＝2.05m，内径为38cm的有机玻璃缸内进行。实验层为粉质黏土，厚度为0.2m。上覆含水层由细砂组成，其中85％的颗粒大于0.075mm。下含水层为中砂层，50％的颗粒粒径大于0.25mm。含水层系统最初为静态，总水头1.56m。上覆含水层水头在整个试验过程中保持不变，而下覆含水层水头在固结试验开始时迅速下降个试验过程中保持不变，而下覆含水层水头在固结试验开始时迅速下降在上、下含水层界面设置两个表盘，记录两个界面的垂直位移。根据刻度盘指示器的读数计算出累积含水层固结程度。该实验试共持续4天。表2为实验弱透水层(黏土层)变形的原始数据。
[0133]
表2室内实验弱透水层变形原始数据
[0134][0135]
(2)参数确定
[0136]
由双对数坐标系下的f～t实测曲线(参见图3)计算出表征标准非达西渗流程度的前段斜率
[0137]
α＝log
10
(5.68/1.94)/log
10
(4705/495)＝0.4771。
[0138]
查表1，其n和i
l
的值分别选取为1.2和16。此后，通过n和i
l
的值获取fd～td标准曲线。其匹配点的坐标为：td＝0.12，fd＝0.69，t＝4115s以及f＝5.42mm.即图4中实心点。通过式(17)和(18)，可得
[0139]
贮水率为
[0140][0141]
渗透系数为
[0142]
k＝tdssb2/t＝(0.12
×
6.5
×
10-2
×
0.22)/4115＝7.6
×
10-8
m/s。
[0143]
(3)验证
[0144]
验证配线法求参与实测数据的拟合程度，可使用均方误差(mean square error，mse)对拟合程度进行检验，其公式如下
[0145][0146]
其中，其中m为观测压实数据总数，f
iobs
为第i次观测压实数据，f
isim
为与f
iobs
同时模拟的第i次压实数据。均方误差越小，参数确定的可靠性越高，反之亦然。
[0147]
本实例计算所得的mse为0.087mm2，而基于达西渗流所得的mse为0.099mm2。因此，考虑非达西渗流下情况下所求的参数值结果更为可靠。可以得出，对于黏土等细粒土组成的弱透水层，考虑非达西渗流作用所得弱透水层水文地质参数确定结果更能科学地反映弱透水层的水动力响应和释水变形规律。