不同降深的稳定流抽水试验求参在梭罗树水源地地下水资源评价中的应用

Application of Steady Flow Pumping Tests with Different Drawdown Depths in Evaluation of Groundwater Resources in Thoreau Tree Water Source Area

摘要: 岚山区是日照市重要的工业基地，经济发展及人口增加迅速，工业项目较多，而岚山区主要供水水源地为绣针河水源地，该水源地下游汾水–荻水水源地由于海水入侵已报废。现在只有中游大朱曹水源地能够利用，但因其水量小且已不能保障，供水不足已严重制约了岚山区的经济发展。而梭罗树水源地由于距岚山城区较近，该水源地淡水资源丰富，可作为城市生活用水，应急供水水源地。因此，文章在查明梭罗树水源地的水文地质条件后，预测了地下水动态变化，并利用2次不同降深的稳定流抽水试验计算有关水文地质参数，最终算出地下水允许开采量，这对水资源紧缺导致的岚山区经济可持续发展、居民饮用水等问题提供了数值化的建议，对于我国其他区域水资源的开采利用也具有一定的参考价值。

Abstract: Lanshan District is an important industrial base of Rizhao City, with rapid economic development and population growth, and more industrial projects, while Lanshan District’s main water supply source is the Xiuzhen river source, which is the downstream Fenshui-Dishui water source, but has been abandoned due to seawater invasion. At present, only the middle reaches of the Dazhucao water source can be used, but because of its small water quantity and inability to guarantee, the shortage of water supply has seriously restricted the economic development of Lanshan District. As the water source of the Thoreau tree is close to Lanshan City, the water source is rich in freshwater resources, which can be used as urban living water and emergency water supply sources. So, after finding out the hydrogeological conditions of the water source of the Thoreau tree, this article initially predicts the variation of groundwater, and calculates the hydrogeology parameters by two steady flow pumping tests with different depths, finally gives the allowances production of groundwater. It not only provides the numerical value of advice to the shortage of water resources of Lanshan District about the problems of drinking water for residents and economic sustainable development, but also gives reference value to other areas of water resources exploitation in our country.

1. 引言

岚山区是日照市重要的工业基地，经济发展及人口增加发展迅速，工业项目较多，相应的工业、生活用水量逐年急剧增加，目前每天需水量达6万m³。而岚山区主要供水水源地为绣针河水源地。该水源地下游汾水–荻水水源地，由于海水入侵已报废，现在只有中游大朱曹水源地能够利用，但因其水量小，每天仅供1.5万m³，且已不能保障，需要从日照水库及莒南县大山水库买水，供水不足已严重制约了岚山区的经济发展，如遇连续干旱年份及突发事件，缺水现象更加严重，因此梭罗树水源地由于距岚山城区较近，可作为城市生活用水，应急供水水源地[1]-[3]。由于该水源地淡水资源丰富，水资源需求市场较大。因此对梭罗树水源地进行勘查与开发，提出科学合理的开发利用方案，对岚山区社会稳定、经济可持续发展，确保居民生活用水安全，将产生巨大的社会效益和经济效益[4]。

2. 勘查区地质概况

日照市岚山区属暖温带潮湿海洋季风气候区，降水多集中在七、八、九三个月份，降水量占全年的50%~70%左右。多年平均蒸发量为1146.2 mm，以五、六月份蒸发量最大，占全年蒸发量的51%。全区年平均干燥度为0.91~0.97，年平均湿度72%左右。岚山区地貌特征主要受梭罗树构造断裂的控制，地貌上呈梯状地形异常明显。梭罗树断裂东西两侧，从山地到河谷阶地过渡区间内，发育着山地和倾斜台地。构成了本区从山地构造剥蚀到河谷侵蚀堆积的完整地貌类型[5] [6]。根据本区地貌特征，可分为尖棱状山地和倾斜台地及河谷阶地三种地貌类型。

3. 水文地质条件和特征

区内地下水赋存条件与分布规律，严格受地层岩性、地形地貌及构造等自然因素的控制。岩浆岩体由于长期出露地表，遭受风化剥蚀，裂隙多被填充。而风化裂隙仅在浅部发育，且受地形坡度大的影响，补给源不足，蓄水能力较弱、富水性差。但是，在地形地貌有利部位，地下水富集，水量较大，对供水有一定意义[7]。第四系分布区，岩性由残坡积，坡积的亚砂土及亚粘土夹碎石透镜体组成。富水性差，而在中南部，第四系则比较发育，主要分布在龙王河和绣针河河谷中。第四系厚度4~12 m，最厚19.25 m。底部普遍存在粗砂砾石层，直接覆盖在变质岩之上。含水层一般3~8 m，岩性主要为中细砂、粗砂及砾石。地下水埋藏浅，上覆弱隔水亚粘土层较薄，易于接受补给。因此，地下水比较丰富，具有较为理想的供水意义[8] [9]。

本次勘查的水源地地下水类型为构造裂隙水，主要构造为梭罗树断裂，该断裂总体走向330˚~340˚，倾向南西，倾角70˚~80˚。勘查区内长度大于3.0公里，宽度50~80 m，隐伏于第四系之下。梭罗树断裂岩石破碎强烈，富水性强，并多发生绿泥石化，碳酸盐化等，根据带内岩石特征，其活动性质以左行压扭为主，晚期张性活动强烈。根据调查，日照市石棉矿自1972年~2002年在梭罗树–水车沟蛇纹岩体内开采石棉矿，开采水平以−50 m、−100 m水平为主，两个中段的日排水量达5000 m³以上，长期疏干排水，没有引起周围村庄地下水位下降，及泉水断流现象，枯水期与丰水期排水量变化不大。该矿于2002年8月停采，目前−50 m、−100 m水平坑道内已储存约400万m³地下水，形成地下水库。

4. 水文地质调查与钻探

为查明水源地水文地质条件，在南起疏港大道北至巨峰河，东起大旺山，西至平山，即东径119˚14'15''~119˚17'15''，北纬35˚10'20''~35˚15'00''范围内进行水文地质调查。查明了地层产状、岩性、含水层分布、构造发育状况及规律、泉水与井点的水文地质条件。在区域内布置了6条激电联剖测线。将激电测深剖面布置在激电联剖曲线有低电阻率异常的地段[10]。激电联合剖面装置采用AO = 110 m，MN = 20 m，点距 = 20 m。激电测深装置，采用不等比对称垂向测深装置，最大AB/2 = 340 m。测量所用设备为重庆地质仪器厂生产的DZD-4多功能电法仪，以2.0 kw发电机配以整流器为供电电源。供电方式为：5 s × 5 s，延时0.2 s，观测1.5个周期。并在充分分析激电联剖资料的基础上，结合测区地质条件，分析了激电异常的形成原因。最终，依据激电联剖和激电测深资料，结合测区地质条件，建议在1线125点、2线218点及4线212点布置水文地质孔查证。

在前面的基础上钻探16个观测孔进行抽水试验，观测孔分布图见图1。本次布设了两条观测线，并对水位、水温，冲洗液消耗量，风化带深度，含水厚度，起止深度、漏水位置等均进行详细观察记录，观测孔的孔深、孔径等符合有关规范及设计要求，抽水试验采用15JQS-120型潜水泵，电机功率55千瓦，流量120 m³/h，杨程60 m，抽水井及观测孔水位观测采用万用电表，流量观测利用大型三角堰，水温用玻璃棒温度计，排水利用直径1.0 m的筒状农用塑料布沿河排泄，长度450 m，已排出采空区外围，防渗效果较好，抽水过程中，资料现场整理，计算，随时指导抽水工作。

Figure 1. Distribution of observation holes in pumping test

图1. 抽水试验观测孔分布图

5. 水文地质参数计算

本次稳定流抽水试验主要分为两个阶段，第一阶段：以4910.64 m³/d的定流量抽水，延续抽水36 h即停止抽水，之后水位开始恢复；第二阶段：以3364.56 m³/d的定流量抽水延续抽水20 h即停止抽水，之后水位开始恢复达90 h。

利用稳定抽水试验形成的水位观测资料以及Aquifer Test软件，进行参数计算。本次抽水试验均采用Theis模型求参，下面为绘制的部分图件(见图2~6)及最终的结果统计表(见表1)：

① 第一次稳定流抽水试验，根据抽水试验基础资料统计试验测得的水位及埋深数据。

Figure 2. Observation hole ZK1 using Theis model to fit the curve

图2. 观测孔ZK1利用Theis模型拟合曲线

Figure 3. Observation hole ZK4 using Theis model to fit the curve

图3. 观测孔ZK4利用Theis模型拟合曲线

Figure 4. Observation hole ZK10 using Theis model to fit the curve

图4. 观测孔ZK10利用Theis模型拟合曲线

② 第二次稳定流抽水试验

Figure 5. Observation hole ZK3 using Theis model to fit the curve

图5. 观测孔ZK3利用Theis model拟合曲线

井孔布置时一条平行地下水流向，另一条垂直地下水流向。设定平行地下水流向的渗透系数为K₁，垂直地下水流向的渗透系数为K₂。根据计算可得K₁ = 25.81 m/d，K₂ = 17.51 m/d。

在此基础上，根据稳定流抽水试验第一落程初始时间段的降深与时间的关系曲线，利用泰斯公式求解含水层的弹性释水系数μ* = 0.000378以及给水度μ = 0.145，计算结果见表2、表3。可知弹性释水系数远小于给水度，因此，水资源评价时，可忽略含水层的弹性释水量。

Figure 6. Observation hole ZK4 using Theis model to fit the curve

图6. 观测孔ZK4利用Theis模型拟合曲线

Table 1. The final result of K

表1. K的最终结果表

井号	第一次稳定流抽水	第二次稳定流抽水	选取值
Pw	21.3	19.4	20.35
ZK1	21.9	22.2	22.05
ZK2	22.7	21.05	21.88
ZK3	18.9	17.45	18.18
ZK4	17.85	11.7	14.78
ZK5	16.4	10.74	13.57
ZK6	18.7	21	19.85
ZK7	12.05	12.97	12.51
ZK8	19.13	13.01	16.07
ZK9	15.7	11.68	13.69
ZK10	17.8	11.53	14.67
ZK13	15.4	16.2	15.8
ZK18	22.6	15.1	18.85
ZK20	25.4	27	26.2

注：K的单位为m/d。

Table 2. Result statistics of μ

表2. μ的结果统计

井号	ZK7	ZK8	ZK9	ZK11	ZK20
参数	0.000973	0.000153	0．000258	0.000227	0.000279

对μ值取均值可得μ = 0.145。

Table 3. Result statistics of μ*

表3. μ*的结果统计

井号	ZK7	ZK8	ZK9	ZK11	ZK20
参数	0.000973	0.000153	0．000258	0.000227	0.000279

对μ*值取均值可得μ* = 0.000378。

6. 计算区数值模型的建立与预测

数值法的目的是通过对研究区水文地质条件的分析和已获取的地下水流场建立计算区地下水系统的数值模拟模型，并通过对已知地下水动态水位的拟合与验证，确定模型的可靠性，预测拟建水源地开采条件下的地下水流场的动态变化趋势，为本地区地下水资源的可持续开发利用提供决策依据。

为了简化野外实际问题，便于对该地下水系统进行分析，需对水文地质条件概化，然后结合计算区实际条件和已有资料，本次模拟采用等效连续介质模型。地下水系统渗流运动要素随时空变化，故地下水为非稳定流。最终将计算区地下水流系统概化为非均质、各向异性、二维非稳定地下水流系统。用下述定解问题来描述：

$\{\begin{array}{l}\frac{\partial }{\partial x}\left(k\left(H-Z_b\right)\frac{\partial H}{\partial x}\right)+\frac{\partial }{\partial y}\left(k\left(H-Z_b\right)\frac{\partial H}{\partial y}\right)+W=\mu \frac{\partial H}{\partial t}\left(x,y\right)\in D,t\ge 0\\ H\left(x,y,t\right)|{}_{t=0}=H_0\left(x,y\right)\left(x,y\right)\in D,t=0\\ k\left(H-Z_b\right)\frac{\partial H}{\partial \overrightarrow{n}}|{}_{{\Gamma }_2}=q\left(x,y,t\right)\left(x,y\right)\in {\Gamma }_2,t>0\end{array}$ (1)

式(1)中：$H$ ——含水层水头(m)；$H_0\left(x,y\right)$ ——含水层初始水头(m)；$Z_b$ ——含水层底板高程(m)；$k$ ——含水层渗透系数(m/d)；$\mu$ ——含水层给水度；$W$ ——含水层垂向补给强度(m/d)；$E$ ——地下水蒸发排泄强度(m/d)；$P$ ——含水层开采强度(m/d)；${\Gamma }_2$ ——已知流量边界；$q\left(x,y,t\right)$ ——含水层侧向单宽补排量(m²/d)；$\overrightarrow{n}$ ——边界上的外法线方向；$D$ ——计算区范围。

求解方法是在计算区域内采用矩形剖分和线性插值，将上述数学模型离散为有限单元方程组，然后利用有限差分法求解。本次运用Visual Modflow 4.0软件，建立了合理的Windows菜单界面与可视化功能，增强了模型数值模拟能力、简化了三维建模的复杂性。该系统还具有强大的图形显示功能，可以很好地可视化显示区域地形地貌概况，地下水流场及其动态过程线。

利用MODFLOW对地下水流进行模拟，通过把研究区在空间和时间上的离散，建立研究区每个网格的水均衡方程式，所有网格方程联立成为一组大型的线性方程组，迭代求解方程组可以得到每个网格的水头值。在不考虑水的密度变化的条件下，含水介质中地下水在三维空间的流动可以用下面的偏微分方程来表示：

$\frac{\partial }{\partial x}\left(K_{xx}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{yy}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial z}\left(K_{zz}\frac{\partial h}{\partial z}\right)-W=S_s\frac{\partial h}{\partial t}$ (2)

式(2)中：$K_{xx}$ ，$K_{yy}$ 和$K_{zz}$ ——分别为渗透系数在x、y和z方向上的分量，假定渗透系数的主轴方向与坐标轴的方向一致；h——水头；W——单位体积流量，用以代表流进汇或来自源的水量；S——孔隙介质的贮水率；t——时间。

根据能量守恒和质量守恒的原则，建立连续性方程，结合达西公式，并用有限差分的方法进行空间上的离散后，得到该计算单元的方程：

$\begin{array}{l}C{R}_{i,j-\frac{1}{2},k}\left(h_{i,j-1,k}-h_{i,j,k}\right)+C{R}_{i,j+\frac{1}{2},k}\left(h_{i,j+1,k}-h_{i,j,k}\right)+C{C}_{i-\frac{1}{2},j,k}\left(h_{i-1,j,k}-h_{i,j,k}\right)\\ +C{C}_{i+\frac{1}{2},j,k}\left(h_{i+1,j,k}-h_{i,j,k}\right)+C{V}_{i,j,k-\frac{1}{2}}\left(h_{i,j,k+1}-h_{i,j,k}\right)+C{V}_{i,j,k+\frac{1}{2}}\left(h_{i,j,k+1}-h_{i,j,k}\right)\\ +Q{S}_{i,j,k}=S{S}_{i,j,k}\left(\Delta r_j\Delta c_i\Delta v_k\right)\frac{\Delta h_{i,j,k}}{\Delta t}\end{array}$ (3)

式(3)中：$\Delta r_j,\Delta c_i,\Delta v_k$ ——分别为单元格沿行、列、层方向上的距离；CR、CC、CV——分别为沿行、列、层方向上的水力传导系数；$S{S}_{i,j,k}$ ——给水度或释水系数；$Q{S}_{i,j,k}$ ——各类源汇项；$h_{i,j,k}$ 、$h_{i,j-1,k}$ 、$h_{i,j+1,k}$ 、$h_{i+1,j,k}$ 、$h_{i-1,j,k}$ 、$h_{i,j,k+1}$ ——分别为计算单元(i, j, k)以及与之相邻的计算单元的水位值。

Figure 7. Water table prediction map of the study area

图7. 研究区地下水位预报图

进一步对时间离散化，根据有限差分计算的方法，可以将水头对时间的偏导数用差商来近似表示：${\left(\frac{\Delta h_{i,j,k}}{\Delta t}\right)}_m\approx \frac{h_{i,j,k}^m-h_{i,j,k}^{m-1}}{t_m-t_{m-1}}$ 。将所有包括未知水头的项移到方程的左侧，将所有的已知项移动到方程的右侧，则有：

$\begin{array}{l}C{V}_{i,j,k-\frac{1}{2}}{h}_{i,j,k}^m+C{C}_{i-\frac{1}{2},j,k}{h}_{i-1,j,k}^m+C{R}_{i,j-\frac{1}{2},k}{h}_{i,j-1,k}^m+(-C{V}_{i,j,k-\frac{1}{2}}-C{C}_{i-\frac{1}{2},j,k}-C{R}_{i,j-\frac{1}{2},k}-C{V}_{i,j,k+\frac{1}{2}}-C{C}_{i+\frac{1}{2},j,k}\\ -C{R}_{i,j+\frac{1}{2},k}+HCO{F}_{i,j,k})h_{i,j,k+1}^m=RH{S}_{i,j,k}\end{array}$ (4)

上式(4)就是地下水流有限差分方程式。对模型所包含的计算单元逐个写出差分方程，可以得到一个线性方程组。根据这两个矩阵，联立方程组，通过迭代法对所求的水头矩阵求解，得到整个研究区的水位分布情况。

模型中水文地质参数初值的确定主要根据前人工作成果和本次野外抽水试验的数据计算所得。结合第四部分水文地质参数计算，根据所建立的水文地质概念模型和数学模型进行预报，预报时段为2009年11月20日至2010年11月19日，共计365天。以2009年11月20日的地下水流场作为预报时段的初始流场。预报时段降水量采用多年平均降水量814.9972 mm。地下水排泄量根据解析法进行估算，由于距离抽水井较远，计算结果显示水源地建成后，预报时段内南侧排泄量不变。将拟建水源地计算开采量(8500 m³/d)作为开采井，输入模型中。经模型运行后，预报得到2010年11月地下水水位见图7。

根据预报结果可知，拟建水源地以8500 m³/d开采量开采一年后，水源地中心水位为12.25 m，水位降深为19.54 m。由地下水位预报图分析，水源地降落漏斗范围不大。根据预报结果显示，水源地中心水位趋于稳定。因此，拟建水源地选取8500 m³/d的开采量，水资源量完全能够保证，且水源地有扩大开采量的潜力。

7. 结论

本文根据岚山区的经济发展受水资源制约，如遇连续干旱年份及突发事件，缺水现象更加严重，梭罗树水源地由于距岚山城区较近，可作为城市生活用水、应急供水水源地。因此，本文对梭罗树水源地进行勘查与开发，建立地下水动态监测系统，对开采主井及水源地内观测孔进行长期动态观测，通过水文地质调查与钻探、水文地质参数计算，进而完成计算区数值模型的建立与预测。根据预报结果可知，拟建水源地以8500 m³/d开采量开采一年后，水源地中心水位为12.25 m，水位降深为19.54 m。由地下水位预报图分析，水源地降落漏斗范围不大，水资源量完全能够保证，且水源地有扩大开采量的潜力。本文的研究对岚山区社会稳定、经济可持续发展，确保居民生活用水安全，将产生巨大的社会效益和经济效益，并对我国其他区域水资源的开采利用也具有一定的参考价值。