南通市地下水壓縮開采對水質咸化的控制效應

陸 瑋，李 兆，駱祖江

(河海大學地球科學與工程學院，江蘇 南京 211100)

南通市地下水資源較為豐富，一直以來以開采優質的第Ⅲ承壓水來滿足日常生產生活用水需求，但隨著地下水的開采和地下水位的下降，從20世紀80年代起出現了水質咸化現象[1-4]，進入90年代后水質咸化現象日趨嚴重[5]。王琦等[6-7]通過對比分析地下水動力場和水化學場的變化特性，認為過量開采第Ⅲ承壓水，導致上覆地層中的咸水入侵進入第Ⅲ承壓水是咸化的主要原因。目前地下水開采對第Ⅲ承壓水咸化影響的定量分析研究成果較少，難以為第Ⅲ承壓水咸化的防控提供科學依據。本文以地下水滲流理論和溶質運移理論[8]為基礎，建立了南通市地下水非穩定滲流與溶質運移三維耦合數值模型，通過模擬預測地下水現狀開采和壓縮開采條件下第Ⅲ承壓水的咸化趨勢，定量評價了南通市地下水壓縮開采對第Ⅲ承壓水咸化的控制效應，可為第Ⅲ承壓水咸化的科學防控提供參考。

1 地質概況

南通市東臨黃海，南望長江，與上海隔江相望，地理位置優越。研究區范圍包括崇川區、港閘區和通州區，如圖1所示。區域內除有少量基巖出露外，其余地區均為第四紀松散沉積物所覆蓋，厚度一般為200～360 m，垂向上相互交錯分布著多層砂層，構成了一個錯綜復雜的含水層系統[9-11]，從上至下依次分為：全新統的潛水含水層(組)、上更新統的第Ⅰ承壓含水層(組)、中更新統的第Ⅱ承壓含水層(組)、下更新統的第Ⅲ承壓含水層(組)。圖2為研究區水文地質剖面，各含水層之間以弱含水的黏性土層分隔或直接接觸，存在著較為強烈的水力聯系，貯存著豐富的地下水資源，其中潛水、第Ⅰ承壓水和第Ⅱ承壓水均為咸水，第Ⅲ承壓水為淡水，水質好，水量大。

圖2 研究區A—A′斷面水文地質剖面示意圖Fig.2 A-A′ hydrogeological profile of the study area

2 概念模型

本次模擬計算平面上以行政區域為界，垂向上包括整個第四紀地層，四周側向邊界均概化為通用水頭邊界和溶質運移通量邊界，底部以第Ⅲ承壓含水層底板為界，是一層厚度30～50 m的致密亞黏土，隔水性能良好，概化為隔水邊界和溶質運移零通量邊界。頂部在接受大氣降水補給的同時，地下水又通過其蒸發排泄，在水力學上既是補給邊界又是排泄邊界，在溶質運移方面概化為通量邊界；將大氣降水入滲補給和地下水蒸發排泄兩者合并考慮，概化為綜合有效入滲補給強度。整個計算區含水層無論是平面上還是剖面上，富水性變化較大，概化為非均質層，垂直方向和水平方向的滲透性相差較大，概化為各向異性。各含水層中的地下水除水平運動外還通過層間的黏性土弱含水層或“天窗”發生水力聯系，地下水流態為三維非穩定流，地下水開采按點井給出。

3 數學模型

3.1 地下水滲流模型

根據南通市的水文地質概念模型，將主滲透方向設為與坐標軸方向一致，建立了南通市地下水系統三維非穩定滲流數學模型[12-13]，模型的微分方程、初始條件和流量邊界分別為

(1)

H(x,y,z,t)|t=0H0(x,y,z,t0) (x,y,z∈Ω)

(2)

(3)

(4)

式中：Kxx、Kyy、Kzz分別為含水層各向異性主滲透方向的滲透系數，m/d；H為點(x,y,z)在t時刻的水頭值，m；W為源匯項，L/d；μs為儲水率，m-1；t為時間，d；Ω為計算域；H(x,y,z,t)為滲流場水頭，m；H0(x,y,z,t0)為初始水頭值，m；q(x,y,z,t)為已知流量；cos(n,x)、cos(n,y)、cos(n,z)分別為流量邊界外法線方向與坐標軸方向夾角的余弦；μ為飽和差或給水度；qw為自由面單位面積上的大氣降雨入滲補給量。式(3)和式(4)分別為第二和第三類流量邊界。

3.2 溶質運移模型

在不考慮溶質吸附、化學反應以及生物降解的情況下，將坐標軸方向設為與溶質運移方向一致，建立了南通市地下水系統三維非穩定流溶質運移數學模型，模型的微分方程[14]、初始條件和邊界條件分別為

(5)

ρ(x,y,z,t)|t=0ρ0(x,y,z,t0) (x,y,z∈Ω)

(6)

x,y,z∈Γ2)

(7)

g(x,y,z,t) (t>0;x,y,z∈Γ3)

(8)

式中：ρ為流體中溶質的質量濃度，mg/L；n為孔隙度；I為源匯項；Dxx、Dyy、Dzz為沿坐標軸方向的水動力彌散系數分量，m2/d；De為分子有效擴散系數，m2/d；vx、vy、vz為介質中流體沿坐標軸方向的平均線性速度，m/d；αL、αT分別為縱向和橫向彌散度，m；ρ(x,y,z,t)為溶質的質量濃度分布，mg/L；ρ0(x,y,z，t0)為已知溶質的質量濃度分布，mg/L；Ω代表整個模型區域；f(x,y,z,t)為Γ2邊界上已知的彌散通量；cos(n,x)、cos(n,y)、cos(n,z)分別為邊界Γ2外法線方向與坐標軸方向夾角的余弦；g(x,y,z,t)為已知函數，表示Γ3邊界法線方向上的對流彌散總通量。式(7)和式(8)分別為第二和第三類邊界條件。

4 模型的識別驗證

采用地下水三維模擬軟件Visual Modflow[15-17]進行模擬計算，在平面上剖分為159×100的網格單元，如圖3(a)所示，其中有效單元為8 332個，無效單元為7 568個。剖面上從上至下分別將潛水、第Ⅰ承壓、第Ⅱ承壓、第Ⅲ承壓含水層以及各含水層之間的黏性土弱含水層剖分成獨立的層位，共分7層，圖3(b)(c)分別為第35行和第43列網格剖分結果。選取2014年8月31日至2016年8月31日和2016年8月31日至2017年12月31日分別作為模型的識別和驗證階段，把每個月作為一個應力期，識別階段共分為24個應力期，驗證階段共分為17個應力期，每個應力期為一個時間步長。

(a) 平面剖分

(b) 第35行剖分

圖3 研究區網格剖分
Fig.3Gridmeshingofthesimulationfield

研究區分布有85口開采井，10口水位觀測井，17口水質監測井，基本覆蓋全區。各含水層的初始等水位線圖由實測資料得到，各含水層之間黏性土弱含水層的初始等水位線圖由其上下含水層插值得到。選取地下水中分布廣泛且較為穩定的氯離子作為地下水咸化的模擬因子[18]，各含水層的氯離子初始質量濃度場由實際監測資料給出，各含水層之間黏性土弱含水層的氯離子初始質量濃度場由其上下含水層插值得到。參數的初始值根據前人資料結合經驗給出[19]，經過識別、驗證，模型共分為59個參數分區，其中潛水含水層4個，第Ⅰ黏性土弱含水層4個，第Ⅰ承壓含水層8個，第Ⅱ黏性土弱含水層4個，第Ⅱ承壓含水層5個，第Ⅲ黏性土弱含水層14個，第Ⅲ承壓含水層20個。綜合有效入滲補給系數為2×10-6。圖4為第Ⅲ承壓含水層水文地質參數分區，表1為第Ⅲ承壓含水層各分區水文地質參數。圖5和圖6分別為識別階段和驗證階段第Ⅲ承壓含水層水位及水質監測井的擬合情況，可見實測值與計算值吻合良好，表明反演計算所得的參數準確、可靠，模型模擬結果符合實際情況。

圖4 第Ⅲ承壓含水層水文地質參數分區Fig.4 Hydrogeological parameter division for aquifer No.Ⅲ

表1 第Ⅲ承壓含水層各分區水文地質參數
Table1HydrogeologicalparametersofeachzonesforaquiferNo.Ⅲ

分區滲透系數/(m·d-1)KxxKyyKzz有效孔隙度儲水率/m-1彌散度/m縱向橫向150.050.05.00.4400.00095070.7220.020.02.00.4430.00065080.8317.517.51.60.2960.000960212.1422.022.02.20.4320.00040080.8550.050.05.00.4380.00095080.8660.060.06.00.4420.00065080.8725.025.02.50.4440.00065080.888.08.00.80.4280.000400151.597.07.00.70.4420.00009590.9101.01.00.10.2060.000500323.21115.015.01.50.2850.000850252.51218.018.01.80.4330.00160090.91310.010.01.00.4400.00160090.91430.030.03.00.4410.00040090.91512.012.01.20.4410.00045090.91625.025.02.50.4480.000910121.21720.020.02.00.4420.00040090.9181.51.50.20.4410.00160090.9198.08.00.80.2000.000500353.52030.030.03.00.2350.000500292.9

(a) 識別階段

(b) 驗證階段圖5 地下水水位擬合結果Fig.5 Fitting graph of the groundwater level

(a) 識別階段

(b) 驗證階段圖6 地下水氯離子質量濃度擬合結果Fig.6 Fitting graph of the groundwater chloride ion concentration

5 地下水水位與咸化趨勢預測

5.1 現狀開采條件下地下水水位與咸化趨勢

按2017年現狀開采布局及開采量對2018—2034年的地下水水位變化和地下水咸化趨勢進行模擬預測?，F狀共有85口開采井，其中第Ⅰ承壓含水層分布有59口開采井，總開采量為444.21萬m3/a；第Ⅲ承壓含水層分布有26口開采井，總開采量為179.51萬m3/a，其中港閘、崇川兩區開采量為89.51萬m3/a，通州區開采量為90.00萬m3/a。降水量采用2007—2017年的年平均值，同時考慮區內蒸發作用，確定綜合降雨入滲補給系數為2×10-6。圖7和圖8分別為2018年1月1日第Ⅲ承壓含水層等水位線和氯離子質量濃度場。

圖7 第Ⅲ承壓含水層等水位線(單位：m)Fig.7 Flow field for aquifer No.Ⅲ(unit: m)

圖8 第Ⅲ承壓含水層氯離子質量濃度場(單位：mg/L)Fig.8 Chloride ion concentration for aquifer No.Ⅲ(unit: mg/L)

a. 地下水水位預測。圖9為2024年和2034年12月31日第Ⅲ承壓含水層預測等水位線，可以看出第Ⅲ承壓含水層地下水水位整體上出現了回升現象。2018—2034年，最靠近降落漏斗的觀山音街道附近水位回升最快，水位上升7.86 m，狼山鎮街道水位回升4.12 m，其余地區水位也呈緩慢上升趨勢。

b. 地下水咸化預測。圖10為2024年和2034年12月31日第Ⅲ承壓含水層預測氯離子質量濃度場，可以看出第Ⅲ承壓含水層中氯離子質量濃度大于 250 mg/L 的面積從2018年初的147.35 km2增加到2034年底的 355.17 km2，其中質量濃度大于270 mg/L的面積達到30.67 km2。2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承壓水咸化速率分別為19.31 km2/a、12.76 km2/a、11.49 km2/a和9.59 km2/a。

(a) 2024年12月31日

(b) 2034年12月31日圖9 現狀開采條件下第Ⅲ承壓含水層預測等水位線(單位：m)Fig.9 Predicted flow field for aquifer No.Ⅲ under current mining conditions(unit: m)

(a) 2024年12月31日

(b) 2034年12月31日圖10 現狀開采條件下第Ⅲ承壓含水層預測氯離子質量濃度場(單位：mg/L)Fig.10 Predicted chloride ion concentration for aquifer No.Ⅲ under current mining conditions(unit: mg/L)

5.2 壓縮開采條件下地下水水位與咸化趨勢

南通市地下水壓縮開采方案中共分為兩個階段，第一階段為2018—2019年，研究區分布有55口開采井，其中第Ⅰ承壓含水層分布有31口開采井，總開采量為352.05萬m3/a，相比現狀開采減少92.16萬m3/a；第Ⅲ承壓含水層分布有24口開采井，總開采量為154.60萬m3/a，相比現狀開采減少24.91萬m3/a，其中港閘、崇川兩區開采量為 74.50萬m3/a，通州區開采量為80.10萬m3/a。第二階段為2020—2034年，在第一階段基礎上不改變開采井的分布，繼續減少每口井的開采量，其中第Ⅰ承壓含水層總開采量減少為301.35萬m3/a，相比現狀開采減少142.86萬m3/a；第Ⅲ承壓含水層總開采量減少為108.90萬m3/a，相比現狀開采減少70.61萬m3/a，其中港閘、崇川兩區開采量為 55.80萬m3/a，通州區開采量為53.10萬m3/a。降雨入滲量同現狀開采條件。

a. 地下水水位預測。圖11為2024和2034年12月31日第Ⅲ承壓含水層預測等水位線，可見第Ⅲ承壓含水層地下水水位整體上是回升的。從2018—2034年，最靠近降落漏斗的觀山音街道附近水位回升最快，水位上升8.00 m，相比現狀開采水位上升了0.14 m，狼山鎮街道水位回升4.31 m，相比現狀開采水位上升0.19 m。

b. 地下水咸化預測。圖12為壓縮開采條件下2024和2034年12月31日第Ⅲ承壓含水層預測氯離子質量濃度場，可以看出第Ⅲ承壓含水層中氯離子質量濃度大于250 mg/L的面積從2018年初的147.35 km2增加到2034年底的329.21 km2，其中質量濃度大于 270 mg/L 的面積為1.76 km2。2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承壓水咸化速率分別為17.50 km2/a、11.92 km2/a、9.93 km2/a和7.52 km2/a。

由以上分析可知，現狀和壓縮開采兩種條件下，第Ⅲ承壓含水層地下水位均逐漸上升，與上覆含水層之間的水頭差均逐漸減小，因此上覆含水層的咸水入侵速度均逐漸減緩。相比于現狀開采，壓縮開采條件下第Ⅲ承壓水咸化速率明顯減慢，且咸化面積明顯減少。由此可以判斷，采用壓縮開采方案能有效控制南通市第Ⅲ承壓水的咸化問題。

(a) 2024年12月31日

(b) 2034年12月31日圖11 壓縮開采條件下第Ⅲ承壓含水層預測等水位線(單位：m)Fig.11 Predicted flow field for aquifer No.Ⅲ under compression mining conditions(unit: m)

(a) 2024年12月31日

(b) 2034年12月31日圖12 壓縮開采條件下第Ⅲ承壓含水層預測氯離子質量濃度場(單位：mg/L)Fig.12 Predicted chloride ion concentration for aquifer No.Ⅲ under compression mining conditions(unit: mg/L)

6 結 論

a. 保持現狀開采條件下，2034年底第Ⅲ承壓含水層中氯離子質量濃度大于250 mg/L和270 mg/L的面積達到355.17 km2和30.67 km2，2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承壓水咸化速率分別為19.31 km2/a、12.76 km2/a、11.49 km2/a和9.59 km2/a。

b. 壓縮開采條件下，2034年底第Ⅲ承壓含水層中氯離子質量濃度大于250 mg/L和270 mg/L的面積分別達到329.21 km2和1.76 km2，2018—2019年、2020—2024年、2025—2029年、2030—2034年第Ⅲ承壓水咸化速率分別為17.50 km2/a、11.92 km2/a、9.93 km2/a和7.52 km2/a。

c. 壓縮開采條件下，相比于現狀開采條件2018—2034年第Ⅲ承壓含水層中氯離子質量濃度大于250 mg/L和270 mg/L的面積分別減少了25.96 km2和28.91 km2，咸化速率明顯減慢，采用壓縮開采方案能有效控制南通市第Ⅲ承壓水的咸化問題。