基于數值模擬的水文地質分析

吳潔河

(廣東省水文局 佛山水文分局，廣東 佛山 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，與蒙特卡洛模擬的結果具有一致性，進一步驗證了本研究數值模擬的可行性。

5 結 論

本文以某研究區為研究對象，基于數值模擬，建立該研究區的數學模型，對其地下水情況進行分析，并對其結果進行敏感性分析和不確定性分析，結論如下：

通過對比實測值與數值模擬結果可得，通過數值模擬對地下水位的預測結果最為準確，雖然實測值與計算值存在一定的差異，但其差值均小于1%，說明采用數值模擬計算弱透水層的地下水位準確性較高。

對承壓含水層1和承壓含水層2敏感性最高的參數均為μ/Ss，對弱透水層敏感性最高的參數為垂直滲透系數。

該研究區的某一測點的地下水位主要集中在13～16 m，與蒙特卡洛模擬的結果具有一致性，進一步驗證了本研究數值模擬的可行性。