陸 瑋,李 兆,駱祖江
(河海大學(xué)地球科學(xué)與工程學(xué)院,江蘇 南京 211100)
南通市地下水資源較為豐富,一直以來以開采優(yōu)質(zhì)的第Ⅲ承壓水來滿足日常生產(chǎn)生活用水需求,但隨著地下水的開采和地下水位的下降,從20世紀(jì)80年代起出現(xiàn)了水質(zhì)咸化現(xiàn)象[1-4],進(jìn)入90年代后水質(zhì)咸化現(xiàn)象日趨嚴(yán)重[5]。王琦等[6-7]通過對(duì)比分析地下水動(dòng)力場(chǎng)和水化學(xué)場(chǎng)的變化特性,認(rèn)為過量開采第Ⅲ承壓水,導(dǎo)致上覆地層中的咸水入侵進(jìn)入第Ⅲ承壓水是咸化的主要原因。目前地下水開采對(duì)第Ⅲ承壓水咸化影響的定量分析研究成果較少,難以為第Ⅲ承壓水咸化的防控提供科學(xué)依據(jù)。本文以地下水滲流理論和溶質(zhì)運(yùn)移理論[8]為基礎(chǔ),建立了南通市地下水非穩(wěn)定滲流與溶質(zhì)運(yùn)移三維耦合數(shù)值模型,通過模擬預(yù)測(cè)地下水現(xiàn)狀開采和壓縮開采條件下第Ⅲ承壓水的咸化趨勢(shì),定量評(píng)價(jià)了南通市地下水壓縮開采對(duì)第Ⅲ承壓水咸化的控制效應(yīng),可為第Ⅲ承壓水咸化的科學(xué)防控提供參考。
南通市東臨黃海,南望長(zhǎng)江,與上海隔江相望,地理位置優(yōu)越。研究區(qū)范圍包括崇川區(qū)、港閘區(qū)和通州區(qū),如圖1所示。區(qū)域內(nèi)除有少量基巖出露外,其余地區(qū)均為第四紀(jì)松散沉積物所覆蓋,厚度一般為200~360 m,垂向上相互交錯(cuò)分布著多層砂層,構(gòu)成了一個(gè)錯(cuò)綜復(fù)雜的含水層系統(tǒng)[9-11],從上至下依次分為:全新統(tǒng)的潛水含水層(組)、上更新統(tǒng)的第Ⅰ承壓含水層(組)、中更新統(tǒng)的第Ⅱ承壓含水層(組)、下更新統(tǒng)的第Ⅲ承壓含水層(組)。圖2為研究區(qū)水文地質(zhì)剖面,各含水層之間以弱含水的黏性土層分隔或直接接觸,存在著較為強(qiáng)烈的水力聯(lián)系,貯存著豐富的地下水資源,其中潛水、第Ⅰ承壓水和第Ⅱ承壓水均為咸水,第Ⅲ承壓水為淡水,水質(zhì)好,水量大。

圖2 研究區(qū)A—A′斷面水文地質(zhì)剖面示意圖Fig.2 A-A′ hydrogeological profile of the study area
本次模擬計(jì)算平面上以行政區(qū)域?yàn)榻?,垂向上包括整個(gè)第四紀(jì)地層,四周側(cè)向邊界均概化為通用水頭邊界和溶質(zhì)運(yùn)移通量邊界,底部以第Ⅲ承壓含水層底板為界,是一層厚度30~50 m的致密亞黏土,隔水性能良好,概化為隔水邊界和溶質(zhì)運(yùn)移零通量邊界。頂部在接受大氣降水補(bǔ)給的同時(shí),地下水又通過其蒸發(fā)排泄,在水力學(xué)上既是補(bǔ)給邊界又是排泄邊界,在溶質(zhì)運(yùn)移方面概化為通量邊界;將大氣降水入滲補(bǔ)給和地下水蒸發(fā)排泄兩者合并考慮,概化為綜合有效入滲補(bǔ)給強(qiáng)度。整個(gè)計(jì)算區(qū)含水層無論是平面上還是剖面上,富水性變化較大,概化為非均質(zhì)層,垂直方向和水平方向的滲透性相差較大,概化為各向異性。各含水層中的地下水除水平運(yùn)動(dòng)外還通過層間的黏性土弱含水層或“天窗”發(fā)生水力聯(lián)系,地下水流態(tài)為三維非穩(wěn)定流,地下水開采按點(diǎn)井給出。
根據(jù)南通市的水文地質(zhì)概念模型,將主滲透方向設(shè)為與坐標(biāo)軸方向一致,建立了南通市地下水系統(tǒng)三維非穩(wěn)定滲流數(shù)學(xué)模型[12-13],模型的微分方程、初始條件和流量邊界分別為
(1)
H(x,y,z,t)|t=0H0(x,y,z,t0) (x,y,z∈Ω)
(2)
(3)
(4)
式中:Kxx、Kyy、Kzz分別為含水層各向異性主滲透方向的滲透系數(shù),m/d;H為點(diǎn)(x,y,z)在t時(shí)刻的水頭值,m;W為源匯項(xiàng),L/d;μs為儲(chǔ)水率,m-1;t為時(shí)間,d;Ω為計(jì)算域;H(x,y,z,t)為滲流場(chǎng)水頭,m;H0(x,y,z,t0)為初始水頭值,m;q(x,y,z,t)為已知流量;cos(n,x)、cos(n,y)、cos(n,z)分別為流量邊界外法線方向與坐標(biāo)軸方向夾角的余弦;μ為飽和差或給水度;qw為自由面單位面積上的大氣降雨入滲補(bǔ)給量。式(3)和式(4)分別為第二和第三類流量邊界。
在不考慮溶質(zhì)吸附、化學(xué)反應(yīng)以及生物降解的情況下,將坐標(biāo)軸方向設(shè)為與溶質(zhì)運(yùn)移方向一致,建立了南通市地下水系統(tǒng)三維非穩(wěn)定流溶質(zhì)運(yùn)移數(shù)學(xué)模型,模型的微分方程[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)

式中:ρ為流體中溶質(zhì)的質(zhì)量濃度,mg/L;n為孔隙度;I為源匯項(xiàng);Dxx、Dyy、Dzz為沿坐標(biāo)軸方向的水動(dòng)力彌散系數(shù)分量,m2/d;De為分子有效擴(kuò)散系數(shù),m2/d;vx、vy、vz為介質(zhì)中流體沿坐標(biāo)軸方向的平均線性速度,m/d;αL、αT分別為縱向和橫向彌散度,m;ρ(x,y,z,t)為溶質(zhì)的質(zhì)量濃度分布,mg/L;ρ0(x,y,z,t0)為已知溶質(zhì)的質(zhì)量濃度分布,mg/L;Ω代表整個(gè)模型區(qū)域;f(x,y,z,t)為Γ2邊界上已知的彌散通量;cos(n,x)、cos(n,y)、cos(n,z)分別為邊界Γ2外法線方向與坐標(biāo)軸方向夾角的余弦;g(x,y,z,t)為已知函數(shù),表示Γ3邊界法線方向上的對(duì)流彌散總通量。式(7)和式(8)分別為第二和第三類邊界條件。
采用地下水三維模擬軟件Visual Modflow[15-17]進(jìn)行模擬計(jì)算,在平面上剖分為159×100的網(wǎng)格單元,如圖3(a)所示,其中有效單元為8 332個(gè),無效單元為7 568個(gè)。剖面上從上至下分別將潛水、第Ⅰ承壓、第Ⅱ承壓、第Ⅲ承壓含水層以及各含水層之間的黏性土弱含水層剖分成獨(dú)立的層位,共分7層,圖3(b)(c)分別為第35行和第43列網(wǎng)格剖分結(jié)果。選取2014年8月31日至2016年8月31日和2016年8月31日至2017年12月31日分別作為模型的識(shí)別和驗(yàn)證階段,把每個(gè)月作為一個(gè)應(yīng)力期,識(shí)別階段共分為24個(gè)應(yīng)力期,驗(yàn)證階段共分為17個(gè)應(yīng)力期,每個(gè)應(yīng)力期為一個(gè)時(shí)間步長(zhǎng)。

(a) 平面剖分

(b) 第35行剖分

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

圖4 第Ⅲ承壓含水層水文地質(zhì)參數(shù)分區(qū)Fig.4 Hydrogeological parameter division for aquifer No.Ⅲ
表1 第Ⅲ承壓含水層各分區(qū)水文地質(zhì)參數(shù)
Table1HydrogeologicalparametersofeachzonesforaquiferNo.Ⅲ

分區(qū)滲透系數(shù)/(m·d-1)KxxKyyKzz有效孔隙度儲(chǔ)水率/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) 識(shí)別階段

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

(a) 識(shí)別階段

(b) 驗(yàn)證階段圖6 地下水氯離子質(zhì)量濃度擬合結(jié)果Fig.6 Fitting graph of the groundwater chloride ion concentration
按2017年現(xiàn)狀開采布局及開采量對(duì)2018—2034年的地下水水位變化和地下水咸化趨勢(shì)進(jìn)行模擬預(yù)測(cè)?,F(xiàn)狀共有85口開采井,其中第Ⅰ承壓含水層分布有59口開采井,總開采量為444.21萬m3/a;第Ⅲ承壓含水層分布有26口開采井,總開采量為179.51萬m3/a,其中港閘、崇川兩區(qū)開采量為89.51萬m3/a,通州區(qū)開采量為90.00萬m3/a。降水量采用2007—2017年的年平均值,同時(shí)考慮區(qū)內(nèi)蒸發(fā)作用,確定綜合降雨入滲補(bǔ)給系數(shù)為2×10-6。圖7和圖8分別為2018年1月1日第Ⅲ承壓含水層等水位線和氯離子質(zhì)量濃度場(chǎng)。

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

圖8 第Ⅲ承壓含水層氯離子質(zhì)量濃度場(chǎng)(單位:mg/L)Fig.8 Chloride ion concentration for aquifer No.Ⅲ(unit: mg/L)
a. 地下水水位預(yù)測(cè)。圖9為2024年和2034年12月31日第Ⅲ承壓含水層預(yù)測(cè)等水位線,可以看出第Ⅲ承壓含水層地下水水位整體上出現(xiàn)了回升現(xiàn)象。2018—2034年,最靠近降落漏斗的觀山音街道附近水位回升最快,水位上升7.86 m,狼山鎮(zhèn)街道水位回升4.12 m,其余地區(qū)水位也呈緩慢上升趨勢(shì)。
b. 地下水咸化預(yù)測(cè)。圖10為2024年和2034年12月31日第Ⅲ承壓含水層預(yù)測(cè)氯離子質(zhì)量濃度場(chǎng),可以看出第Ⅲ承壓含水層中氯離子質(zhì)量濃度大于 250 mg/L 的面積從2018年初的147.35 km2增加到2034年底的 355.17 km2,其中質(zhì)量濃度大于270 mg/L的面積達(dá)到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日?qǐng)D9 現(xiàn)狀開采條件下第Ⅲ承壓含水層預(yù)測(cè)等水位線(單位:m)Fig.9 Predicted flow field for aquifer No.Ⅲ under current mining conditions(unit: m)


(a) 2024年12月31日

(b) 2034年12月31日?qǐng)D10 現(xiàn)狀開采條件下第Ⅲ承壓含水層預(yù)測(cè)氯離子質(zhì)量濃度場(chǎng)(單位:mg/L)Fig.10 Predicted chloride ion concentration for aquifer No.Ⅲ under current mining conditions(unit: mg/L)
南通市地下水壓縮開采方案中共分為兩個(gè)階段,第一階段為2018—2019年,研究區(qū)分布有55口開采井,其中第Ⅰ承壓含水層分布有31口開采井,總開采量為352.05萬m3/a,相比現(xiàn)狀開采減少92.16萬m3/a;第Ⅲ承壓含水層分布有24口開采井,總開采量為154.60萬m3/a,相比現(xiàn)狀開采減少24.91萬m3/a,其中港閘、崇川兩區(qū)開采量為 74.50萬m3/a,通州區(qū)開采量為80.10萬m3/a。第二階段為2020—2034年,在第一階段基礎(chǔ)上不改變開采井的分布,繼續(xù)減少每口井的開采量,其中第Ⅰ承壓含水層總開采量減少為301.35萬m3/a,相比現(xiàn)狀開采減少142.86萬m3/a;第Ⅲ承壓含水層總開采量減少為108.90萬m3/a,相比現(xiàn)狀開采減少70.61萬m3/a,其中港閘、崇川兩區(qū)開采量為 55.80萬m3/a,通州區(qū)開采量為53.10萬m3/a。降雨入滲量同現(xiàn)狀開采條件。
a. 地下水水位預(yù)測(cè)。圖11為2024和2034年12月31日第Ⅲ承壓含水層預(yù)測(cè)等水位線,可見第Ⅲ承壓含水層地下水水位整體上是回升的。從2018—2034年,最靠近降落漏斗的觀山音街道附近水位回升最快,水位上升8.00 m,相比現(xiàn)狀開采水位上升了0.14 m,狼山鎮(zhèn)街道水位回升4.31 m,相比現(xiàn)狀開采水位上升0.19 m。
b. 地下水咸化預(yù)測(cè)。圖12為壓縮開采條件下2024和2034年12月31日第Ⅲ承壓含水層預(yù)測(cè)氯離子質(zhì)量濃度場(chǎng),可以看出第Ⅲ承壓含水層中氯離子質(zhì)量濃度大于250 mg/L的面積從2018年初的147.35 km2增加到2034年底的329.21 km2,其中質(zhì)量濃度大于 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。
由以上分析可知,現(xiàn)狀和壓縮開采兩種條件下,第Ⅲ承壓含水層地下水位均逐漸上升,與上覆含水層之間的水頭差均逐漸減小,因此上覆含水層的咸水入侵速度均逐漸減緩。相比于現(xiàn)狀開采,壓縮開采條件下第Ⅲ承壓水咸化速率明顯減慢,且咸化面積明顯減少。由此可以判斷,采用壓縮開采方案能有效控制南通市第Ⅲ承壓水的咸化問題。

(a) 2024年12月31日

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


(a) 2024年12月31日

(b) 2034年12月31日?qǐng)D12 壓縮開采條件下第Ⅲ承壓含水層預(yù)測(cè)氯離子質(zhì)量濃度場(chǎng)(單位:mg/L)Fig.12 Predicted chloride ion concentration for aquifer No.Ⅲ under compression mining conditions(unit: mg/L)
a. 保持現(xiàn)狀開采條件下,2034年底第Ⅲ承壓含水層中氯離子質(zhì)量濃度大于250 mg/L和270 mg/L的面積達(dá)到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年底第Ⅲ承壓含水層中氯離子質(zhì)量濃度大于250 mg/L和270 mg/L的面積分別達(dá)到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. 壓縮開采條件下,相比于現(xiàn)狀開采條件2018—2034年第Ⅲ承壓含水層中氯離子質(zhì)量濃度大于250 mg/L和270 mg/L的面積分別減少了25.96 km2和28.91 km2,咸化速率明顯減慢,采用壓縮開采方案能有效控制南通市第Ⅲ承壓水的咸化問題。