基于地面沉降的地下水可開采資源量規(guī)劃評價

基于地面沉降的地下水可開采資源量規(guī)劃評價

0地下水開采量與地面沉降耦合模型隨著我國社會經(jīng)濟和城市建設的快速發(fā)展，地下水需求急劇增加，地下水資源過度開采，導致土壤沉降的地面沉降災害日益嚴重。特別是中國的長江三角洲、華北平原和蘇北平原經(jīng)歷了非常厚的第四紀慢沉層，巖石結構復雜，分布著多種承受水層。以往的地下水可采資源規(guī)劃評價由于受計算機技術和計算方法等的限制,一般按單層平面二維流或準三維模型單獨進行評價,且不考慮地面沉降的影響。結果一方面使地下水的開發(fā)本質(zhì)上和地面沉降控制相脫離,另一方面單層單獨進行評價又容易造成地下水可采資源量的重復計算,再加上不考慮地面沉降導致的含水層孔隙度和滲透系數(shù)等的變化,往往使計算所得的地下水可采資源量偏大,這是主觀上造成地下水過量開采,誘發(fā)地面沉降的原因。筆者針對我國區(qū)域第四紀松散沉積層地區(qū)地下水開采引起的地面沉降地質(zhì)災害問題,根據(jù)滲流理論和土力學理論,建立起地下水三維非穩(wěn)定滲流與地面沉降耦合數(shù)學模型;并考慮了含水層孔隙度、滲透系數(shù)、儲水率隨含水層發(fā)生固結沉降的變化特征,采用三維有限元數(shù)值分析方法,以江蘇省南通市地下水開采為例,基于地面沉降的控制要求,規(guī)劃評價出了各鄉(xiāng)鎮(zhèn)各含水層的地下水可采資源量,達到了在獲得最大地下水可采量的同時,將地面沉降控制到最小的目的。1地下水流假說與地面沉降數(shù)學模型1.1潛水面營養(yǎng)品質(zhì)指標t,有以下幾種基本概念地下水三維非穩(wěn)定滲流數(shù)學模型的一般形式可概括為{??x(Κxx?Η?x)+??y(Κyy?Η?y)+??z(Κzz?Η?z)+W=μs?Η?t,(x,y,z)∈Ω,t>0;Η(x,y,z,t)|t=0=Η0(x,y,z,t0),(x,y,z)∈Ω,t=0;Η(x,y,z,t)|Γ1=Η1(x,y,z,t),(x,y,z)∈Γ1,t>0;Κxx?Η?xcos(n,x)+Κyy?Η?ycos(n,x)+Κzz?Η?zcos(n,z)|Γ2=q(x,y,z,t),(x,y,z)∈Γ2,t>0;Η(x,y,z,t)|Γ3=z(x,y,t),(x,y,z)∈Γ3,t>0;Κxx(?Η?x)2+Κyy(?Η?y)2+Κzz(?Η?z)2-(Κzz+qw)?Η?z+qw|Γ3=μ?Η?t,(x,y,z)∈Γ3,t>0。(1)?????????????????????????????x(Kxx?H?x)+??y(Kyy?H?y)+??z(Kzz?H?z)+W=μs?H?t,(x,y,z)∈Ω,t>0;H(x,y,z,t)|t=0=H0(x,y,z,t0),(x,y,z)∈Ω,t=0;H(x,y,z,t)|Γ1=H1(x,y,z,t),(x,y,z)∈Γ1,t>0;Kxx?H?xcos(n,x)+Kyy?H?ycos(n,x)+Kzz?H?zcos(n,z)|Γ2=q(x,y,z,t),(x,y,z)∈Γ2,t>0;H(x,y,z,t)|Γ3=z(x,y,t),(x,y,z)∈Γ3,t>0;Kxx(?H?x)2+Kyy(?H?y)2+Kzz(?H?z)2?(Kzz+qw)?H?z+qw|Γ3=μ?H?t,(x,y,z)∈Γ3,t>0。(1)式中:μs為儲水率(1/m);Kxx,Kyy,Kzz為含水層各向異性主方向滲透系數(shù)(m/d);H為點(x,y,z)在t時刻的水頭值(m);W為源匯項(1/d);t為時間(d);Ω為計算域;H0為點(x,y,z)處的初始水頭值(m);H1(x,y,z,t)為第一類邊界上的已知水頭值(m);q(x,y,z,t)為第二類邊界上單位面積的已知流量(m/d);cos(n,x)、cos(n,y)、cos(n,z)分別為流量邊界外法線方向與坐標軸方向夾角的余弦;μ為潛水含水層在潛水面上的重力給水度;qw為潛水面大氣降雨入滲補給強度與潛水面蒸發(fā)強度的代數(shù)和(m/d);Γ1、Γ2、Γ3分別為第一類邊界、第二類邊界和自由面邊界。該模型考慮了地下水在x、y、z這3個方向上的流動狀態(tài),將含水層和弱透水層中的地下水流均作三維滲流處理,即將各含水砂層之間的黏性土層作為弱含水層和含水層一起直接參與計算。1.2土體骨架結構參數(shù)的變化表征飽和土層的固結沉降主要是土體在附加壓力和自重的作用下,由于孔隙水的排出,土體總體積壓縮而引起的。附加壓力與土體的自重壓力、土層的厚度與壓縮性、天然孔隙度與滲透性決定著沉降量值的大小。因此,由地下水位下降引起的含水層壓縮量的計算模型為{Δb1=-ΔΗ(1-n+nw)μskeb0=-ΔΗμfe,Δb*1=-ΔΗ(1-n+nw)μskvb0=-ΔΗμfv,Δb2=-ΔΗμskeb0=-ΔΗμfe,Δb*2=-ΔΗμskvb0=-ΔΗμfv。(2)其中∶{μske=ρwg[3(1-2)v2G(1+v)],Δσ′≤Δσ′max;μskv=0.434CcρwgΔσ′0(1+e0),Δσ′zz≥Δσ′max。(3)式中:Δb1、Δb*1分別為潛水含水層的彈性變形量和非彈性變形量;Δb2、Δb*2分別為承壓含水層的彈性變形量和非彈性變形量;n為含水層的孔隙率;nw為水位以上作為多孔介質(zhì)總體積一部分的濕氣容量;b0為可壓縮含水層的初始厚度(m);ΔH為孔隙水頭變化值(m);μske、μfe分別為土體骨架成分的彈性儲水率(1/m)和彈性儲水因子(μfe=μskeb0,無量綱);μskv和μfv分別為土體骨架成分的非彈性儲水率(1/m)和非彈性儲水因子(μfv=μfkvb0,無量綱);G為剪切模量(MPa);v為泊松比;g為重力常數(shù);Cc為土體的壓縮指數(shù);σ′0為初始有效應力(kPa);e0為初始孔隙比;ρw為水的密度(kg/m3);Δσ′zz為土層垂向有效應力的變化量(kPa);σ′max為土層的前期固結應力(kPa)。在沉降模型中,μske、μskv為土體力學參數(shù)的函數(shù),它考慮了土體的非線性固結問題,反映了應力狀態(tài)對土體骨架儲水率的影響。假設整個含水層的骨架都可發(fā)生壓縮并產(chǎn)生沉降作用,參數(shù)μfe、μfv是在整個含水層規(guī)模上的等效參數(shù)。當各含水層的壓縮量確定后,總地面沉降量也由此而求出。2地下水運行方程與土壤變形方程的結合2.1固結過程中的滲透系數(shù)k根據(jù)柯晉-卡門(Kozeny-Carman)經(jīng)驗方程可得出Κ≈Κ0(nn0)3?(1-n01-n)2。(4)式中:K0和K分別為初始滲透系數(shù)和固結過程中的滲透系數(shù);n0和n分別為含水層的初始孔隙度和固結過程中的孔隙度。設某土層厚度為bi,由于開采地下水引起的垂直沉降量為di,假若不考慮土層的側向變形,則對于固結過程中的孔隙度n有n=n0-dibi-di≈n0-dibi。(5)根據(jù)滲透系數(shù)K和孔隙度n的關系,通過沉降模型中得出的不斷變化的孔隙度來不斷修正滲透系數(shù)的值,使水流模型反映土體的實際壓縮情況;水流模型中的參數(shù)如滲透系數(shù)、儲水率在沉降過程中是變化的,是有效應力和孔隙比的函數(shù)。2.2地下水流數(shù)學模型的耦合為計算由土體壓縮帶來的釋水量變化,將模型(1)中的地下水運動連續(xù)性方程的右端項分為2項,由傳統(tǒng)定義的(未考慮含水層土體變形耦合效應)單位時間內(nèi)單元體的釋水量μs?Η?t和單位時間內(nèi)單元體由土體壓縮帶來的釋水量μsk?Η?t兩部分組成。即??x(Κxx?Η?x)+??y(Κyy?Η?y)+??z(Κzz?Η?z)+W=μ′s?Η?t。(6)式中:μ′s為考慮含水層土體變形耦合效應條件下的含水層儲水率(1/m),μ′s=μs+μsk;μs為未考慮含水層土體變形耦合效應條件下的含水層儲水率(1/m);μsk(等于μske或μskv)為含水層土體骨架部分的儲水率(1/m);μsk?Η?x為土體骨架壓縮帶來的釋水量(1/d)。通過水頭項和土體儲水項二者之間的內(nèi)在聯(lián)系相互影響、相互制約,使地下水運動連續(xù)性方程與含水層土體變形方程形成耦合,達到地下水流數(shù)學模型與地面沉降數(shù)學模型的耦合效果。上述模型采用有限單元法進行求解,并采用預處理共軛梯度法迭代求解代數(shù)方程組,其中的自由面求解采用改進復合單元滲透矩陣調(diào)整法,自由面流量補給矩陣采用高斯公式將求自由面邊界的面積分轉化為求單元體積分與各表面積分之差來計算,以避免計算自由面的具體位置。將整個數(shù)學模型的求解過程,用VisualFortran95語言編制成計算機程序,在PIV3.0微機上進行計算。3使用示例3.1含水層的沉降南通市地處江蘇省東南部長江口北岸,瀕江臨海,總面積8500km2,區(qū)內(nèi)沉積了一套巨厚的第四系松散地層,其厚度由南向北逐漸增大,垂向上多層砂層相互疊置,構成錯綜復雜的含水層系統(tǒng)。根據(jù)巖性和時代,區(qū)內(nèi)孔隙含水層自上而下分布有全新統(tǒng)潛水含水層、上更新統(tǒng)第Ⅰ承壓含水層、中更新統(tǒng)第Ⅱ承壓含水層、下更新統(tǒng)第Ⅲ承壓含水層,各含水層之間均以黏性土弱含水層相分隔,并發(fā)生強烈的水力聯(lián)系。20世紀南通市累積地面沉降達70mm,市區(qū)最大累計沉降量約150mm,沉降中心與市區(qū)各承壓水水位降落漏斗一致,目前地面沉降仍在發(fā)展,市區(qū)最大累計沉降量達200mm左右。這是由于系統(tǒng)內(nèi)從上往下各含水層之間均有一定厚度的黏土、亞黏土弱含水層相分隔,局部地段缺失又使相鄰含水層直接接觸,發(fā)生強烈的水力聯(lián)系;系統(tǒng)內(nèi)主采層第Ⅲ承壓含水層的開采將導致各含水層和黏土、亞黏土弱含水層孔隙水壓力的降低,有效應力的增加,使各層均發(fā)生壓縮,導致地面沉降。3.2地下水邊界深刻性分析將南通市從上往下自潛水含水層至第Ⅲ承壓含水層作為一個統(tǒng)一的水文地質(zhì)系統(tǒng),計算目的層為第Ⅰ、Ⅱ、Ⅲ承壓含水層,各含水層之間的黏性土、亞黏性土弱含水層也按獨立的層位參與計算,共分7層,各層均概化為非均質(zhì)各向異性。系統(tǒng)的頂部邊界為潛水含水層的自由面,既是一補給邊界,又是一排泄邊界;側向邊界上,系統(tǒng)內(nèi)多年的地下水開采已與北部鹽城、西部泰州之間形成了地下水分水嶺,以行政區(qū)劃為界概化為隔水邊界。東北部含水層延伸至黃海底部,以黃海為界概化為透水邊界。系統(tǒng)南部的長江邊界,從上往下切割至第Ⅰ承壓含水層,故將潛水含水層、第Ⅰ承壓含水層及之間的黏性土弱含水層的長江邊界概化為定水頭邊界,下部各層由于多年的地下水開采,與蘇州、上海之間已形成了分水嶺,以行政區(qū)劃為界概化為隔水邊界;系統(tǒng)底部是一隔水邊界。各含水層之間的垂向水力聯(lián)系極為復雜,地下水位均受地下水開采的影響,地下水流態(tài)呈現(xiàn)出強烈的三維非穩(wěn)定狀態(tài)。3.3地下水滲透與土壤沉降的結合值的模擬3.3.1地下水的水文地質(zhì)圖將南通市整個第四紀地下水系統(tǒng)垂向上從上往下,按潛水含水層、第Ⅰ、第Ⅱ、第Ⅲ承壓含水層及各含水層之間的黏性土弱含水層,共剖分成7層。平面上采用1∶2.5萬的水文地質(zhì)圖,每層剖分成7276個矩形有限元網(wǎng)格單元,共剖分成50932個矩形有限元網(wǎng)格單元,詳見圖1。用單元邊作折線近似模擬區(qū)域邊界及水文地質(zhì)參數(shù)分區(qū)界線,界線交接點及水位觀測孔均以節(jié)點表示,抽水孔組及供水點均置于單元內(nèi)部。3.3.2水文地質(zhì)參數(shù)分區(qū)以2004-12-01—2005-12-31作為模型識別時段,以2006-01-01—2006-12-01作為模型驗證時段,將2004年12月和2006年10、11月分別單獨作為一個地下水開采制度期,其余時段均以每個季度作為一個制度期,共分9個制度期,每個制度期又分10個計算時間步長;各含水層的初始流場采用2004-12-01的實測水位值,采用Kriging插值獲得同一層各單元中心點及各含水層之間的黏性土弱含水層的初始流場。各層的初始壓縮沉降量對模擬計算影響甚微且由于受實測資料的限制,統(tǒng)一取0;大氣降水入滲和地下水蒸發(fā)按模型要求綜合起來,概化成有效入滲補給強度,各含水層地下水的開采均以鄉(xiāng)鎮(zhèn)為單位,按大井處理,其開采量由實際調(diào)查獲得。由于南通市第四紀松散沉積層無明顯地質(zhì)構造分布,因此將各含水層和黏性土弱含水層的滲透系數(shù)在水平方向上作各向同性考慮,即Kxx=Kyy,垂向滲透系數(shù)按經(jīng)驗取水平方向滲透系數(shù)的1/10給出初值。潛水含水層的初始孔隙度,按3%計算。潛水含水層水位以上孔隙中的濕氣容量,按1.5%計算。潛水含水層按非彈性變形考慮;系統(tǒng)內(nèi)承壓含水層的水位均高于含水層頂板,均按彈性變形考慮;松散土層中,承壓含水層的水體膨脹釋水要比含水層骨架壓縮釋水小得多,幾乎可以忽略不計,所以將承壓含水層的彈性釋水認為主要由含水層骨架壓縮引起,承壓含水層的彈性釋水率即為含水層的骨架壓縮釋水率。各含水層均有一定數(shù)量的地下水位觀測井用來進行水位擬合,共計34個觀測井,基本能控制全區(qū)地下水流場。經(jīng)過多次調(diào)試、優(yōu)選和識別、驗證,獲得了最優(yōu)參數(shù)分區(qū)和參數(shù)數(shù)組,其中:潛水含水層、第Ⅰ黏性土弱含水層、第Ⅰ承壓含水層各分為8個參數(shù)分區(qū);第Ⅱ黏性土弱含水層、第Ⅱ承壓含水層、第Ⅲ黏性土弱含水層,第Ⅲ承壓含水層分別分為10,11,32,36個參數(shù)分區(qū),共分113個參數(shù)分區(qū)。圖2為第Ⅲ承壓含水層水文地質(zhì)參數(shù)分區(qū)特征,表1給出了第Ⅲ承壓含水層的水文地質(zhì)參數(shù)值。以2005-10-01為例,圖3給出了第Ⅲ承壓含水層各觀測井的計算水位與實測水位擬合情況,其水位誤差均在1m以下,與實際流場相吻合。同時,模型識別計算所得的各含水層及各黏性土弱含水層地下水的總補給量與儲存量消耗之和基本等于總排泄量,均符合水均衡的要求,說明模型的收斂性和穩(wěn)定性較好。因此,模型計算所得結果正確、可靠,可以用來預測地下水系統(tǒng)的動態(tài)變化和地面沉降預警。3.4地下水開采計劃3.4.1地下水位下降南通市地下水主要開采第Ⅲ承壓水含水層,年開采量約為9856.45×104m3,但區(qū)內(nèi)開采井的分布不均勻,一般集中分布于鄉(xiāng)鎮(zhèn)區(qū)所在地,開采量相差懸殊,中東部地區(qū)及沿海一帶深層地下水開采強度較之西部如皋、海安大,地下水位下降幅度明顯大于西部,目前已形成了以如東和海門一帶為中心的區(qū)域性地下水位降落漏斗,誘發(fā)了嚴重的地面沉降地質(zhì)災害。根據(jù)南通市有關部門的要求,到2021年12月底,第Ⅰ、第Ⅱ承壓含水層的最低地下水位分別不低于-10m和-30m,且地下水位降落漏斗中心大致與行政區(qū)域中心相一致,第Ⅲ承壓含水層最低地下水位海安、如皋不低于-20m,通州、南通市區(qū)、海門、如東不低于-40m,啟東不低于-35m,各縣市的最大地面沉降量從2006年12月初至2021年12月底累計不超過50mm。本著這一原則,利用上述識別驗證過的數(shù)學模型,在盡可能滿足各鄉(xiāng)鎮(zhèn)對各含水層現(xiàn)有地下水需求量的同時,根據(jù)現(xiàn)有掌握的地下水系統(tǒng)動態(tài)資料,以2006年地下水開采布局為基礎,將2006-01-01作為預測的初始時刻,初始地面沉降量取值為0,預測各含水層2007年至2021年逐年年底的地下水流場分布特征及地面沉降的發(fā)展趨勢,以每年為1個應力期,共分15個應力期,通過調(diào)整現(xiàn)有的地下水開采格局,優(yōu)化確定出規(guī)劃到每個鄉(xiāng)鎮(zhèn)的地下水可開采資源量。3.4.2地下水開采量預測模型運行結果表明:按上述規(guī)劃方案進行開采,至2021年底,第Ⅰ、Ⅱ、Ⅲ承壓含水層中的最低地下水位和地面沉降量均能達到控制要求,地下水位降落漏斗日趨平緩,各區(qū)域地面沉降趨于穩(wěn)定,詳見表2。全市第Ⅰ、Ⅱ、Ⅲ承壓含水層地下水總的年規(guī)劃可開采量為17870.56×104m3,其中第Ⅰ、Ⅱ、Ⅲ承壓含水層的年規(guī)劃可開采量分別為:2882.11×104,2418.48×104,12569.97×104m3,分別占總可采量的16.1%,13.5%和70.4%。各縣市第Ⅰ、Ⅱ、Ⅲ承壓含水層地下水的年規(guī)劃可采量詳見表3,并將地下水可采量具體規(guī)劃到了各鄉(xiāng)鎮(zhèn)(略)。以第Ⅲ承壓含水層2015年和2021年為例,其地下水流場變化見圖4,各地地面沉降發(fā)展趨勢見圖5。從上述地下水可采量