基于数值模拟的水文地质分析

2022-10-18吴洁河

水利科技与经济 2022年10期

吴洁河

(广东省水文局佛山水文分局，广东佛山 528100)

1 概述

随意开发地下水资源会危害当地的水文地质条件，合理开发和利用地下水资源已成为当前的热点研究课题。

李果等[1]基于Gibbs模型，对某矿区地下水的变化特征开展研究，分析影响水质的因素。林若静等[2]以临汾盆地为研究对象，分析地下水位与各因素之间的关系。贾建伟等[3]基于补给法和排泄法，对长江流域的地下水资源情况开展研究，并分析其空间特征。李明等[4]利用GMS软件，对某污染场地地下水中污染物质的迁移规律开展研究。张彦[5]等利用年代波动性分析，对某灌区地下水的响应特征进行分析，得出影响地下水位的主要因素。

本文以某研究区为研究对象，基于数值模拟，建立该研究区的数学模型，对其地下水情况进行分析，并对其结果进行敏感性分析和不确定性分析。

2 工程概况

本研究对某市研究区地下水情况进行分析，该地区的土质以黏土和砂土为主，主要供水来源于地下水，故对其地下水的相关参数进行分析，有利于优化和涉及地下水的使用方案。该市的相关气候及地质情况见表1。

表1 气候及地质情况

由表1可知，该市的年均降雨量小于年均蒸发量，说明该市有轻微的干旱现象，由降雨带来的供水较少，主要需要依靠地下水及湖泊等方式进行供水。为优化该市的地下水供给结构，对该市的水文地质情况进行数值模拟，并分析地下水参数对其敏感性的影响。

3 模型的建立

根据前人对于该市的勘察资料，可根据地质结构将该研究区分为3个部分，根据该市的水文地质概念模型，建立该试验区的数学模型，公式如下：

(1)

式中：Γ2为二类通量边界；Γ0为地下水的上边界；Ω为研究区域；ε为源汇项；μ为给水度；K为渗透系数；X、Y、Z分别为不同的方向；h(x，y，z，t)为时间为t时的标高；h0为初始标高；S为单位储水系数；q为水流通量。

利用有限元软件，对该研究区进行模拟，根据地质勘察结果，将研究区分为3个部分，分别为承压含水层1，承压含水层2和弱透水层。其相关的水文地质情况见表2-表4。

表2 承压含水层1水文地质参数

表3 承压含水层2水文地质参数

表4 弱透水层水文地质参数

为验证以上模型的可行性，将数值模拟的结果与实际观测的结果进行对比，分别对比承压含水层1和弱透水层的模拟结果与观测结果，其时间-水位曲线见图1、图2。

由图1可知，随着时间的增大，测点205和测点1422的水位整体呈下降趋势，而测点Q32z的水位随时间的变化具有一定的波动性。

测点205的实测值与计算值的差距较小，计算值曲线的数值与变化趋势均与实测值保持一致。

图1 承压含水层1计算实测对比图

测点1422的计算值与实测值的差距较小，随着时间的增大，测点1422的实测值与计算值逐渐出现差异性，且差距越来越明显，但是其计算曲线变化趋势与实测值保持一致。

测点Q32z的实测值与计算值的差异大于上述两个测点，随着时间的增大，计算值与实测值的差距显著。

综合以上分析可得，当地下水位较高时，数值模拟对承压含水层地下水位的预测结果较为准确。随着地下水位的增大，数值模拟结果与实测值会产生一定的差距，但其差值均在1%以内，说明采用数值模拟对承压含水层的地下水位进行预测是可行的。

图2 弱透水层计算实测对比图

由图2可知，随着时间的增大，测点H102与测点H162的水位呈先下降后上升的趋势，而测点186的水位随时间的增大逐渐下降。

测点186计算值与实测值的水位变化趋势一致，但是其水位的数值具有一定的差异性。当时间为50～100 d和200～250 d时，二者之间的差距较为明显，在其他时间下，实测值与计算值的差距较小。

综合以上分析可得，通过数值模拟对弱透水层的地下水位进行预测是可行的。其中，对测点H162的预测结果最为准确，虽然测点H102和测点186的实测值与计算值存在一定的差异，但其差值均小于1 %，说明采用数值模拟计算弱透水层的地下水位准确性较高。

4 水文地质参数不确定性分析

为研究各因素对该研究区地下水的影响，对各因素进行敏感性分析，计算其灵敏度系数(β)，以体现各因素对地下水位的影响程度。敏感性分析计算公式如下：

(2)

其中：H′为各参数增加或减少10%后的地下水位；H为初始地下水位。

根据式(2)对各参数的灵敏度系数进行计算，计算结果见表5-表7。

表5 承压含水层1敏感性系数

表6 承压含水层2敏感性系数

表7 弱透水层敏感性系数

表5为承压含水层1的灵敏度分析结果。由表5可知，当层号为1-4时，各因素的灵敏度系数最大，其中μ/Ss增大10%的灵敏度最大，垂直渗透系数增大10%的灵敏度最小。说明μ/Ss对承压含水层1的影响最大，且在层号为1-4时最为突出。而垂直渗透系数对承压含水层的影响最小，其数值的改变对承压含水层的地下水位的改变影响较小。

表6为承压含水层2的灵敏度分析结果。由表6可知，当层号为1、11、14时，各因素变化对应的灵敏度系数较大。其中当层号为1、μ/Ss减小10%时，其对应的灵敏度系数最大。垂直渗透系数的灵敏度系数最小，且随着层号的变化，其灵敏度系数变化趋势较为平缓，说明垂直渗透系数对承压含水层2的地下水位影响较小，当垂直渗透系数发生变化时，承压含水层2的地下水位发生的变化较小。μ/Ss的灵敏度系数最大，说明μ/Ss对承压含水层2的地下水影响程度较大，当μ/Ss发生变化时，地下水发生变化的程度较大。

表7为弱透水层的灵敏度分析结果。由表7可知，当层号为2、4，且垂直渗透系数增加10%时，弱透水层的灵敏度系数最大，说明垂直渗透系数对弱透水层地下水位的影响较大。当水平渗透系数和μ/Ss发生变化时，弱透水层各层号的灵敏度系数较小，且其灵敏度系数变化趋势较为平缓。说明以上两种因素对弱透水层的地下水位的影响较小，其数值上的变化不会影响地下水位发生变化。

为分析该研究区水文地质的不确定性情况，对该研究区的地下水位采用蒙特卡洛模拟，以测点H102为例，利用Matlab生成100组随机数。当地下水位为15.42 m时，其频率最大，其值为23.5%；当地下水位为18.24 m时，其频率最小，其值为1.21%。对比前文对于地下水位的分析可得，该测点的地下水位主要集中在13～16 m，与蒙特卡洛模拟的结果具有一致性，进一步验证了本研究数值模拟的可行性。