《Journal of Hydrology》：Two new analytical models for constant rate pumping test in double-porosity multilayered confined systems: new robin flow equation for aquitards
编辑推荐：
研究人员针对双孔隙度多层承压系统中的恒定速率抽水试验提出了两个新解析模型。提出了一种反映弱透水层效应的新Robin流动方程。该方程可简化为教科书中的典型Robin边界条件。预测了由基质储水效应引起的三峰形时间降深分布。
研究人员针对双孔隙度多层承压系统中的恒定速率抽水试验提出了两个新解析模型。提出了一种反映弱透水层效应的新Robin流动方程。该方程可简化为教科书中的典型Robin边界条件。预测了由基质储水效应引起的三峰形时间降深分布。
**论文解读文章**
**研究背景与问题**
含水层水力参数是地下水流动模型的基础输入，广泛应用于可持续地下水管理、气候适应性供水、地热能源开发和碳封存等领域。恒定速率抽水试验（CRT）是估算这些参数的标准方法，尤其适用于含双孔隙度（double-porosity）的裂隙含水层。双孔隙度概念由Barenblatt等（1960）提出，将裂隙含水层分为高渗透率的裂隙和低渗透率的基质，地下水主要沿裂隙流动，而基质向裂隙供水。
多层承压系统（包含多个含水层和弱透水层）中的CRT建模已有研究，但现有模型存在局限：第一类模型忽略弱透水层储水效应，导致20%的降深高估；第二类模型假设准三维流动，但在弱透水层-含水层相互作用问题中产生30%-60%的相对误差；第三类模型考虑三维流动，但缺乏对双孔隙度效应的全面考虑，仅Sedghi等（2018）的模型在底层含水层中考虑了双孔隙度。因此，有必要开发一个针对双孔隙度三层承压系统（两个含水层夹一个弱透水层）中CRT诱导三维流动的新解析模型。此外，传统数值解需要在弱透水层内进行精细空间离散，导致计算成本高昂，需要一种无需弱透水层离散的流动方程。
**主要技术方法**
研究人员采用以下关键技术方法：
1. 通过高斯散度定理（divergence theorem of Gauss）推导新Robin流动方程，该方程基于有限厚度单元，假设弱透水层内仅垂直流动且降深线性分布，从而无需弱透水层空间离散。
2. 应用拉普拉斯变换（Laplace transform）和韦伯变换（Weber transform）推导标准半解析解和Robin-flow半解析解，并使用Stehfest算法进行数值反拉普拉斯变换。
3. 利用Mathematica的NDSolve函数求解有限元数值解，并比较不同离散方案下的精度和效率。
4. 采用归一化敏感性系数（Si）分析参数相关性，避免参数估计的不可靠性。
5. 将模型应用于Greene（1993）在南达科他州西部Rapid City地区进行的现场CRT，通过Levenberg-Marquardt算法进行曲线拟合以估算水力参数。
**研究结果**
**3.1 新Robin流动方程的有效性**
通过比较标准半解析解（基于传统弱透水层流动方程）和Robin-flow半解析解（基于新Robin流动方程），验证了新方程的有效性。当满足以下条件时，两者在时空降深分布上的相对误差≤5%：弱透水层垂向水力传导率与上含水层径向水力传导率之比≤10?3，弱透水层厚度与含水层厚度之比≤101·1，以及log??(b'/b?)+7×log??(Ss,f'/Ss,f,?)≤11.6（类似条件适用于下含水层）。这些条件覆盖了大多数弱透水层的实际参数范围（如储水率10??~10?3 m?1）。敏感性分析表明，弱透水层径向水力传导率（Kr,f'）为不敏感参数，支持了仅垂直流动的假设。此外，弱透水层内降深的垂向分布近似线性，验证了线性假设的合理性。新Robin流动方程可简化为双孔隙度滞后Robin边界条件、Huang等（2020）的滞后Robin边界条件或教科书中的典型Robin边界条件。
**3.2 三峰形时间降深分布**
研究首次预测了双孔隙度三层承压系统中含水层时间降深的三峰形分布。通过分析标准半解析解计算的降深斜率（k）的波谷数量，发现三峰形降深分布包含三个平缓段（Flat 1~Flat 3），分别对应：抽水含水层基质储水效应、弱透水层基质储水效应、非抽水含水层基质储水效应。当改变抽水含水层的两个集总参数（Ss,m,?/Ss,f,?和C?b?/Kr,f,?）时，波谷数量可在1至3之间变化（图5）。对于非抽水含水层，同样存在三峰或双峰分布（图7）。这些结果可作为识别基质储水效应来源的指南：通过曲线拟合估计集总参数，然后对照图5和图7确定区域，再根据表2辨别产生平缓段的基质层。
**3.3 新Robin流动方程实现粗网格含水层离散**
Robin-flow数值解允许在弱透水层内无离散，在含水层内采用粗正方形离散，而标准数值解需要全域精细正方形离散才能达到准确的质量平衡。以精度指标AI（接近1表示质量守恒）和效率指标EI（标准数值解与Robin-flow数值解的计算时间比）为比较标准。在实例中，Robin-flow数值解在2002个节点下AI收敛至1.003，而标准数值解在均匀精细正方形离散（1,184,026个节点）下AI=0.99，但计算时间比EI=1503，即Robin-flow数值解的计算时间仅为标准数值解的1/1503。若采用粗三角形离散或矩形离散，标准数值解的AI仅0.72~0.91，表明其精度不足。因此，新Robin流动方程在计算效率上具有显著优势，尤其当弱透水层远薄于含水层时。
**3.4 现场CRT应用**
将模型应用于Greene（1993）在南达科他州的现场CRT（四天抽水试验，抽水率0.0429 m3/s，含水层系统由砂岩和喀斯特灰岩构成）。标准半解析解和Robin-flow半解析解均与观测降深数据拟合良好，标准误差（SEE）均小于Neuman和Witherspoon（1969）忽略双孔隙度效应的解。参数估计值接近，其中抽水含水层基质储水率（Ss,m,?）远大于其他储水率，表明基质储水是主要水源。敏感性分析显示，除垂向水力传导率（Kz,f,?和Kz,f,?）相关、以及弱透水层径向水力传导率（Kr,f'）、非抽水含水层水交换系数（C?）和储水率（Ss,f,?）不敏感外，其他参数估计可靠。通过联合观测井数据可进一步降低参数相关性。
**讨论与结论**
本研究提出了两个新的解析模型，用于双孔隙度三层承压系统中的CRT。新Robin流动方程通过高斯散度定理推导，可避免弱透水层空间离散，在满足条件时与标准解误差≤5%。该方程可简化为多种Robin边界条件，适用于大多数弱透水层。研究首次预测了三峰形时间降深分布，由三个基质储水效应引起，为理解双孔隙度多层系统流动行为提供了新视角。与标准数值解相比，新Robin流动方程可节省三个数量级的计算时间，在多层系统中优势更为显著。现场应用证明了其可行性。然而，该方程假设弱透水层内仅垂直流动和线性降深分布，仅在特定条件下适用；否则应使用传统弱透水层流动方程。
**研究结论**
本文提出了两个用于双孔隙度三层承压系统（两个含水层夹一个弱透水层）中恒定速率抽水试验的新解析模型。通过高斯散度定理基于有限厚度单元推导了一个新的Robin流动方程，而传统弱透水层流动方程基于无限小立方体建立。所谓的标准半解析解依赖于传统弱透水层流动方程，而Robin-flow半解析解应用了新Robin流动方程。附加解是Robin-flow半解析解的一部分，通过释放弱透水层内垂向线性降深假设，描述了CRT早期弱透水层中的瞬态流动。还开发了模型的有限元解。
主要发现如下：
- 当满足条件（b'/b?≤101·1、b'/b?≤101·1、log??(b'/b?)+7×log??(Ss,f'/Ss,f,?)≤11.6、log??(b'/b?)+7×log??(Ss,f'/Ss,f,?)≤11.6、Kz,f'/Kr,f,?≤10?3、Kz,f'/Kr,f,?≤10?3）时，标准半解析解和Robin-flow半解析解在时空降深分布上的相对误差≤5%。这表明新Robin流动方程适用于大多数弱透水层。新Robin流动方程可简化为一个新的双孔隙度滞后Robin边界条件、Huang等（2020）的滞后Robin边界条件，或地下水文教科书中的典型Robin边界条件。预测的含水层时间降深可呈现三峰形分布，包含三个平缓段，分别由含水层和弱透水层的基质储水效应引起。图5和图7中平缓段数量的分布可能有助于理解双孔隙度三层承压系统的流动行为。此外，在Greene（1993）的现场CRT中，标准半解析解和Robin-flow半解析解给出了接近的参数估计值。
- 采用新Robin流动方程的Robin-flow数值解允许在弱透水层内无离散，在含水层内采用粗正方形离散。然而，如果在弱透水层内应用精细正方形离散但在含水层内采用粗矩形离散，标准数值解则不准确。实例表明，具有2002个节点的Robin-flow数值解的计算时间仅为采用传统弱透水层流动方程和均匀精细正方形离散（1,184,026个节点）的准确标准数值解计算时间的1/1503。当弱透水层厚度远小于含水层厚度时，标准数值解会产生大量节点。因此，新Robin流动方程是传统弱透水层流动方程的有用且高效的替代方案。
总之，本研究展示了由基质储水效应引起的三峰形降深分布（三个平缓段），而文献中仅报道双峰形降深分布（两个平缓段，如Sedghi等，2018；Wang等，2024）。与传统弱透水层流动方程相比，新Robin流动方程在模拟含一个弱透水层的三层承压系统中的CRT时，实现了弱透水层内无空间离散，节省了三个数量级的计算时间。对于具有多个弱透水层的多层系统，计算效率更为重要。此外，新Robin流动方程适用于Greene（1993）的现场CRT。然而，新Robin流动方程假设弱透水层内仅垂直流动和线性降深分布，因此仅在上述条件下适用。当条件不满足时，应使用传统弱透水层流动方程。
