地下水開采引發(fā)的土體變形對(duì)污染物遷移影響研究

馬青山，賈軍元，田福金，雷 廷

(中國(guó)地質(zhì)調(diào)查局南京地質(zhì)調(diào)查中心，南京 210016)

0 前言

地下水開采引發(fā)的地下水污染問(wèn)題一直以來(lái)備受國(guó)內(nèi)外許多學(xué)者的高度重視[1-3]。當(dāng)前，地下水污染物遷移的研究多以剛性多孔介質(zhì)中溶質(zhì)運(yùn)移的對(duì)流-彌散理論為主[4]，而忽略了多孔介質(zhì)變形對(duì)污染物遷移產(chǎn)生的影響。事實(shí)上，在地下水開采過(guò)程中，地下水位下降，土體有效應(yīng)力增加，土體壓縮變形，含水層中孔隙水的滲流特征會(huì)發(fā)生變化，進(jìn)而影響污染物隨土體孔隙中水的運(yùn)動(dòng)。因此，研究變形介質(zhì)中污染物遷移特征對(duì)保護(hù)地下水環(huán)境具有重要的現(xiàn)實(shí)意義。

在過(guò)去二十年中，多孔介質(zhì)中污染物遷移的研究已經(jīng)開始考慮土體的固結(jié)耦合效應(yīng)。Loroy等人首先分析了溶質(zhì)遷移過(guò)程中孔隙介質(zhì)的固結(jié)效應(yīng)，然而卻沒(méi)有正確反映多孔介質(zhì)變形對(duì)對(duì)流-擴(kuò)散方程中對(duì)流項(xiàng)的影響[5]。此后Smith等人通過(guò)對(duì)粘土固結(jié)進(jìn)行應(yīng)變分析，提出了飽和變形多孔介質(zhì)的一維污染物運(yùn)移理論[6]。Alshawabkeh等研究了黏土層體積變化對(duì)一維污染物溶質(zhì)運(yùn)移的影響，提出了同時(shí)考慮固結(jié)作用和溶質(zhì)擴(kuò)散的方程[7]。lewis等考慮了大變形情況下的污染物遷移情況，探討了固結(jié)變形、對(duì)流、有效擴(kuò)散系數(shù)等因素對(duì)污染物遷移的影響[8]。王水林等考慮了二維情況下孔隙介質(zhì)變形對(duì)溶質(zhì)運(yùn)移的影響，提出了考慮孔隙介質(zhì)固結(jié)效應(yīng)的溶質(zhì)運(yùn)移方程，并給出了相應(yīng)的求解方法[9]。張志宏等提出了黏土防滲層中污染物遷移轉(zhuǎn)化的一維數(shù)學(xué)模型，并在合理簡(jiǎn)化的基礎(chǔ)上給出了模型的解析解[10]。

但上述成果都限于局部尺度，主要研究荷載引起襯墊層變形對(duì)污染物遷移的影響。本文考慮了大尺度松散沉積層固結(jié)過(guò)程中污染物的遷移，以地下水滲流理論、太沙基一維固結(jié)理論和溶質(zhì)運(yùn)移對(duì)流-彌散理論為基礎(chǔ)，同時(shí)考慮孔隙度、滲透系數(shù)、水動(dòng)力彌散系數(shù)的動(dòng)態(tài)變化，通過(guò)建立變形土體地下水污染物遷移三維耦合數(shù)值模型，研究了地下水開采引起的土體變形對(duì)污染物遷移造成的影響。

1 理論基礎(chǔ)

1.1 地下水三維非穩(wěn)定流

根據(jù)滲流連續(xù)性方程和達(dá)西定律，對(duì)于各向異性介質(zhì)，以主滲透方向?yàn)樽鴺?biāo)軸方向，則描述地下水運(yùn)動(dòng)的三維非穩(wěn)定流方程可表示為[11]：

(1)

式中：Kxx、Kyy分別為水平滲透系數(shù)，m/d；Kzz為垂向滲透系數(shù)，m/d；h為t時(shí)刻點(diǎn)(x,y,z)的水頭值m；W為即源匯項(xiàng)，1/d；μs為儲(chǔ)水率，1/m；t為時(shí)間，d；Ω為計(jì)算域。

1.2 污染物遷移對(duì)流-彌散方程

地下水中污染物的遷移和轉(zhuǎn)化是物理、化學(xué)及生物作用的一個(gè)綜合過(guò)程[12]。在不考慮污染物吸附、化學(xué)反應(yīng)以及生物降解的情況下，描述任意時(shí)刻多孔介質(zhì)污染物遷移規(guī)律的水動(dòng)力彌散方程可表示為(選取滲流主方向與坐標(biāo)軸方向一致)[13]：

(2)

式中：c為流體中溶質(zhì)的質(zhì)量分?jǐn)?shù)，mg/L；n為孔隙度；I為源匯項(xiàng)；Dxx、Dyy、Dzz為沿坐標(biāo)軸方向的水動(dòng)力彌散系數(shù)分量，m2/d；vx、vy、vz為介質(zhì)中流體沿坐標(biāo)軸方向的平均線性速度，m/d。

1.3 土體變形計(jì)算方法

針對(duì)地下水開采引起的土體固結(jié)變形問(wèn)題，本文采用太沙基一維固結(jié)理論進(jìn)行闡述。根據(jù)Terzaghi有效應(yīng)力原理，假定土體總應(yīng)力保持不變，并且認(rèn)為土體只在垂向上發(fā)生變形，則土體一維變形計(jì)算方程如下[14]：

Δb=-Δhμskb0

(3)

式中：Δb為含水層的變形量，m；μsk為含水層土體的骨架儲(chǔ)水率，1/m；Δh為水頭變化值，m；b0為含水層的初始厚度，m。

根據(jù)含水層有效應(yīng)力與前期固結(jié)應(yīng)力的相關(guān)關(guān)系，將含水層變形分為彈性變形階段和非彈性變形階段，對(duì)應(yīng)土體變形的不同階段，土體骨架儲(chǔ)水率的變化分別為[15]：

(4)

1.4 參數(shù)動(dòng)態(tài)變化

1.4.1 孔隙度與滲透系數(shù)的非線性

地下水開采過(guò)程中，根據(jù)有效應(yīng)力原理，隨著地下水的不斷開采，水位不斷降低，孔隙水壓力逐漸減小，含水層土體的有效應(yīng)力逐漸增加，孔隙度逐漸減小，含水層滲透系數(shù)隨之不斷變化。這些變化又影響孔隙流體的流動(dòng)和壓力分布。因此，有必要考慮多孔介質(zhì)中應(yīng)力場(chǎng)與滲流場(chǎng)之間的相互作用。

假設(shè)土層厚度為b0，由于降水引起的沉降量為Δb，不考慮土層的側(cè)向變形，則對(duì)于固結(jié)過(guò)程中孔隙度n有[15]：

(5)

一般地，沉降量與土層厚度相比非常小，因此，式(5)可簡(jiǎn)化為：

(6)

根據(jù)Kozeny-Carman方程，滲透系數(shù)的變化可用表示為：

(7)

1.4.2 水動(dòng)力彌散系數(shù)的非線性

土體固結(jié)變形下的污染物遷移問(wèn)題屬于多場(chǎng)耦合問(wèn)題，即應(yīng)力場(chǎng)、滲流場(chǎng)和濃度場(chǎng)之間的耦合。當(dāng)其中一個(gè)物理場(chǎng)發(fā)生變化時(shí)，必將會(huì)引起其他物理場(chǎng)變化。土體在壓密過(guò)程中，孔隙度、滲透系數(shù)將發(fā)生變化，進(jìn)而引起水動(dòng)力彌散系數(shù)的改變。水動(dòng)力彌散系數(shù)的動(dòng)態(tài)表達(dá)式如下[16]：

(8)

式中：Dm為分子擴(kuò)散系數(shù)，m2/d；Dxx、Dyy、Dzz為水動(dòng)力彌散系數(shù)張量對(duì)角線元素；αL為縱向彌散度，m；αT為橫向彌散度，m。

1.5 初始條件

地下水滲流初始條件為：

h(x,y,z,t)|t=0=h0(x,y,z)，(x,y,z)∈Ω

(9)

式中：h(x,y,z,t)為滲流場(chǎng)水頭，m；h0(x,y,z)為初始水頭值，m。

溶質(zhì)運(yùn)移初始條件為：

c(x,y,z,t)|t=0=c0(x,y,z)，(x,y,z)∈Ω(10)

式中：c(x,y,z,t)為濃度場(chǎng)質(zhì)量分?jǐn)?shù)，mg/L；c0(x,y,z)為初始質(zhì)量分?jǐn)?shù)值，mg/L。

地面沉降初始條件：

Δb(x,y,z,t)|t=0=0，(x,y,z)∈Ω

(11)

1.6 邊界條件

地下水滲流邊界條件包括：

h(x,y,z,t)|Γ1=h1(x,y,z,t)，t>0，(x,y,z)∈Γ1

(12)

(13)

(14)

式中：h1(x,y,z,t)表示邊界Γ1上點(diǎn)(x,y,z)在時(shí)刻t的水頭值，m；n為邊界Γ2的外法線方向；q(x,y,z,t)為已知函數(shù)，表示單位面積上流入(流出)的量，m/d；cos(n,x)、cos(n,y)、cos(n,z)分別為流量邊界外法線方向與坐標(biāo)軸方向夾角的余弦；μ為給水度；z(x,y,t)為各計(jì)算節(jié)點(diǎn)標(biāo)高，m；cosθ為自由面外法線方向與鉛垂線之間的夾角。

污染物運(yùn)移邊界條件包括：

c(x,y,z,t)|Γ1=cb(x,y,z,t)，(x,y,z)∈Γ1

(15)

(16)

(x,y,z)∈Γ3

(17)

式中：Γ1表示定濃度邊界；cb(x,y,z,t)為任意時(shí)刻滲流域內(nèi)第一類邊界(Γ1)上流體中已知污染物的濃度分布；q(x,y,z,t)表示Γ2邊界上已知的彌散通量；cos(x,n)、cos(y,n)、cos(z,n)分別為邊界Γ2外法線方向與坐標(biāo)軸方向夾角余弦。若邊界為不透水邊界時(shí)，邊界上沒(méi)有物質(zhì)交換，即邊界上的彌散通量為零時(shí)，則q(x,y,z,t)=0。g(x,y,z,t)為已知函數(shù)，表示Γ3邊界法線方向上的對(duì)流彌散總通量。

采用伽遼金有限元法對(duì)上述控制方程進(jìn)行離散，結(jié)合初始條件，邊界條件和參數(shù)的動(dòng)態(tài)變化模型即可運(yùn)用Fortran語(yǔ)言編制相應(yīng)的有限元程序進(jìn)行求解。

2 應(yīng)用算例

2.1 研究區(qū)概況

研究區(qū)位于長(zhǎng)江三角洲平原，沉積了巨厚的第四紀(jì)松散沉積物，地下水以第四系松散孔隙水為主。由于水文地質(zhì)條件的不同加之成因歷史的差異，“上咸下淡”的地下水水質(zhì)格局普遍存于相鄰含水層中，即水質(zhì)較差的咸水貯存在淺部的含水層中，而水質(zhì)較好的淡水則分布在下部較深的含水層中。根據(jù)水文地質(zhì)條件，將研究區(qū)含水層細(xì)分為7層。垂向從上往下依次為：潛水含水層、第I承壓含水層，第II承壓含水層、第III承壓含水層以及各含水層之間的粘性土弱含水層，其中潛水含水層巖性以粉砂為主，底板埋深20～30m，靜水位埋深0.9～1.2m；第I承壓含水層巖性以粉砂亞砂土為主，底板埋深62～70m，靜水位埋深3.2～3.8m；第II承壓含水層巖性亞砂土、細(xì)中砂為主，底板埋深93～110m，靜水位埋深7.6～8.18m；第III承壓含水層巖性以細(xì)中砂、中粗砂為主，底板埋深125～150m，靜水位埋深11.5～13.4m。

2.2 模型離散

模型的平面尺度為1 000m×1 000m，垂向深度為150m。模型采用六面體八節(jié)點(diǎn)單元進(jìn)行空間離散，平面上剖分為2 400個(gè)矩形網(wǎng)格單元，一共分成16 800個(gè)網(wǎng)格單元(圖1)。垂向上從上往下依次剖分為：潛水含水層、第I、第II、第III承壓含水層及其之間的粘性土弱透水層共7層。各層均概化為均質(zhì)各向異性，并且均按獨(dú)立的層位參與計(jì)算，另外各層之間均發(fā)生水力聯(lián)系，地下水流態(tài)為三維非穩(wěn)定流。潛水含水層、第I、第II、第III承壓含水層底板標(biāo)高分別為-25、-67、-107、150m。每層劃分為一個(gè)參數(shù)分區(qū)，共劃分為7個(gè)參數(shù)分區(qū)。模型四周概化為第一類邊界條件-已知水頭和已知濃度邊界，頂部為自由面邊界。底部與下覆承壓含水層之間有較厚的粘土層，將其概化為隔水邊界和零通量邊界。

圖1 模型剖分示意Figure 1 A schematic diagram of model subdivision

2.3 地下水開采條件下污染物運(yùn)移模擬研究

模型選擇氯離子為模擬因子，潛水含水層、第I、第II、第III承壓含水層的初始水位分別為-1、-3.5、-7.8、-12.2m，初始氯離子的質(zhì)量分?jǐn)?shù)分別為2 485、1 775 、1 065 和35.5 mg/L。模型模擬時(shí)間為30d。每天為一個(gè)時(shí)間步長(zhǎng)，共30個(gè)時(shí)間步長(zhǎng)。抽水井位于X=500m，Y=500m處，抽取第III承壓含水層地下水，抽水量為500m3/d。觀測(cè)點(diǎn)位于距離抽水井5m處的第III承壓含水層。

模型各參數(shù)分區(qū)參數(shù)值如表1所示。假設(shè)模型中各參數(shù)分區(qū)的縱向彌散度相同，值為10m，且彌散度的縱橫比采用常用的經(jīng)驗(yàn)值，即橫向彌散度與縱向彌散度的比值均為0.1。

表1 模型參數(shù)

3 結(jié)果分析

3.1 參數(shù)動(dòng)態(tài)變化特征

參數(shù)準(zhǔn)確性對(duì)模型計(jì)算結(jié)果的精度有著深刻的影響。考慮土體變形效應(yīng)，污染物遷移問(wèn)題實(shí)際上是應(yīng)力場(chǎng)、滲流場(chǎng)和濃度場(chǎng)之間相互影響、相互作用的多場(chǎng)耦合問(wèn)題。地下水開采引發(fā)土體變形、土體水力學(xué)參數(shù)和土力學(xué)參數(shù)以及水動(dòng)力彌散系數(shù)均會(huì)發(fā)生變化。因此，有必要正確模擬參數(shù)的變化過(guò)程。選取觀測(cè)點(diǎn)(495，495，-130)所在的單元進(jìn)行孔隙度、滲透系數(shù)和水動(dòng)力彌散系數(shù)的變化分析。

圖2 孔隙度沿抽水井徑向方向變化趨勢(shì)Figure 2 Variation trend of porosity along pumping well radial

圖3 孔隙度和地下水位動(dòng)態(tài)變化趨勢(shì)Figure 3 Dynamic variation trends of porosity and groundwater level

不同時(shí)刻孔隙度隨抽水井徑向距離的變化趨勢(shì)如圖2所示。從圖中可以看出，孔隙度的最大變化發(fā)生在抽水井附近，孔隙度減少至0.399 3。

圖3反映了抽水期間孔隙度和地下水位隨時(shí)間的變化趨勢(shì)。從圖中可以看出，孔隙度和地下水位的變化趨勢(shì)一致。地下水開采、地下水位降低、土體有效應(yīng)力增大，導(dǎo)致土體壓縮變形，從而引起孔隙度的變化。地下水位達(dá)到穩(wěn)定，孔隙度不再發(fā)生改變，此時(shí)地下水位下降11m，孔隙度從0.4減小至0.399 6。

圖4 滲透系數(shù)沿抽水井徑向方向變化趨勢(shì)Figure 4 Variation trend of permeability coefficient along pumping well radiation

圖5 滲透系數(shù)和地下水位動(dòng)態(tài)變化趨勢(shì)Figure 5 Dynamic variation trends of permeability and groundwater level

圖4反映了水平滲透系數(shù)在不同時(shí)刻隨抽水井徑向距離的變化趨勢(shì)。從圖中可以看出，在井附近，水平滲透系數(shù)的值在1天內(nèi)下降至0.993m/d。

滲透系數(shù)隨地下水位的變化關(guān)系如圖5所示?？梢钥闯觯瑵B透系數(shù)和地下水位變化趨勢(shì)一致，地下水位下降時(shí)，土體被壓密，孔隙度減小，土體出水能力變?nèi)?，從而?dǎo)致滲透系數(shù)減小。第15d時(shí)，地下水位達(dá)到穩(wěn)定，滲透系數(shù)從1m/d減小至0.995 7m/d。

水動(dòng)力彌散系數(shù)和地下水位之間的關(guān)系如圖6所示。圖中水動(dòng)力彌散系數(shù)與地下水的變化呈現(xiàn)相反的趨勢(shì)，開采地下水時(shí)，水位下降較快，水力梯度較大，滲流速度較快，對(duì)流作用明顯，因此，水動(dòng)力彌散系數(shù)變大。當(dāng)?shù)叵滤贿_(dá)到穩(wěn)定時(shí)，滲流場(chǎng)和應(yīng)力場(chǎng)狀態(tài)不再發(fā)生變化，孔隙度、滲透系數(shù)和滲流速度也不再變化，因此水動(dòng)力彌散系數(shù)也隨之達(dá)到穩(wěn)定狀態(tài)。

圖6 水動(dòng)力彌散系數(shù)和地下水位動(dòng)態(tài)變化趨勢(shì)Figure 6 Dynamic variation trends of hydrodynamic dispersion coefficient and groundwater level

3.2 污染物遷移規(guī)律研究

圖7～圖8反映了考慮土體固結(jié)效應(yīng)和未考慮土體固結(jié)效應(yīng)兩種條件下，水平水動(dòng)力彌散系數(shù)和垂向水動(dòng)力彌散系數(shù)隨時(shí)間的變化趨勢(shì)圖。從圖中可以看出，土體固結(jié)效應(yīng)條件下，水平水動(dòng)力彌散系數(shù)和垂向水動(dòng)力彌散系數(shù)的值均小于未考慮土體固結(jié)效應(yīng)下的值。這是因?yàn)橥馏w固結(jié)變形、孔隙度和滲透系數(shù)減小，地下水平均線性速度減小，進(jìn)而導(dǎo)致水動(dòng)力彌散系數(shù)的減小。

圖7 考慮土體固結(jié)和未考慮土體固結(jié)條件下水平方向水動(dòng)力彌散系數(shù)對(duì)比Figure 7 Contrast of horizontal hydrodynamic dispersion coefficients both considered and not considered soil mass consolidation

圖8 考慮土體固結(jié)和未考慮土體固結(jié)條件下垂向方向水動(dòng)力彌散系數(shù)對(duì)比Figure 8 Contrast of vertical hydrodynamic dispersion coefficients both considered and not considered soil mass consolidation

圖9反應(yīng)了考慮土體固結(jié)效應(yīng)和未考慮土體固結(jié)效應(yīng)兩種情況下觀測(cè)點(diǎn)處Cl-的質(zhì)量分?jǐn)?shù)變化趨勢(shì)，從圖中可以看出：考慮土體固結(jié)效應(yīng)條件下Cl-的質(zhì)量分?jǐn)?shù)變化曲線較未考慮土體固結(jié)效應(yīng)的曲線平緩，濃度值比未考慮土體固結(jié)效應(yīng)模擬時(shí)的偏小，到模擬最后時(shí)刻第30天，考慮土體固結(jié)效應(yīng)觀測(cè)質(zhì)量分?jǐn)?shù)為401.5mg/L，小于未考慮土體固結(jié)效應(yīng)時(shí)的418.3mg/L。其原因主要是地下水開采，土體有效應(yīng)力增大，土體骨架壓縮變形，地下水平均線性速度減小，導(dǎo)致水動(dòng)力彌散系數(shù)增幅減小(圖7～圖8)，進(jìn)而延緩了溶質(zhì)運(yùn)移的過(guò)程。由此可見(jiàn)，土體壓縮變形對(duì)污染物傳輸起抑制作用。

圖9 第III承壓含水層觀測(cè)點(diǎn)處氯離子變化趨勢(shì)Figure 9 Variation trends of chloridions on confined aquifer III observation points

4 結(jié)論

地下水開采過(guò)程中引發(fā)的土體變形對(duì)污染物遷移的影響研究是地下水環(huán)境保護(hù)工作中面臨的一個(gè)新的重要難題，通過(guò)本次模擬研究可以得出如下結(jié)論：

1)基于地下水滲流理論、太沙基一維固結(jié)理論和溶質(zhì)運(yùn)移理論，同時(shí)考慮孔隙度、滲透系數(shù)、水動(dòng)力彌散參數(shù)的非線性變化，建立了變形土體污染物運(yùn)移三維耦合數(shù)值模型。在此基礎(chǔ)上，通過(guò)對(duì)實(shí)際案例的模擬分析，研究了土體固結(jié)變形與污染物遷移三維耦合作用下污染物遷移的規(guī)律，模擬計(jì)算了污染物濃度的變化趨勢(shì)，計(jì)算結(jié)果更加科學(xué)真實(shí)。

2)通過(guò)對(duì)實(shí)際案例模擬分析可知，地下水開采，土體中的孔隙水壓力減小，土體有效應(yīng)力增大，土體發(fā)生固結(jié)變形，孔隙度減小，導(dǎo)致滲流速度減慢，從而減緩了污染物的對(duì)流擴(kuò)散機(jī)制，也就是說(shuō)，土體固結(jié)作用對(duì)污染物在多孔介質(zhì)中的運(yùn)移起抑制作用。