大慶市貼不貼泡區(qū)多層介質(zhì)滲透規(guī)律及中石油類(lèi)污染的數(shù)值模擬研究

多層水系統(tǒng)在自然界中廣泛分布。例如，松嫩平原和松遼平原分布著這樣的地下水系統(tǒng)，但一些大型油田，如大慶和遼河油田，只分布在同一個(gè)地方。這類(lèi)地下水系統(tǒng)往往遭受來(lái)自油田開(kāi)發(fā)區(qū)石油類(lèi)污染質(zhì)的污染。含水系統(tǒng)中不同層的含水介質(zhì)具有不同的水文地質(zhì)特性,特別是其中的含水層之間的弱透水層,往往既是下部承壓含水層中優(yōu)質(zhì)地下水的補(bǔ)給通道,又是污染通道。因此,本文針對(duì)這一實(shí)際問(wèn)題,選取大慶市的一個(gè)典型的納污湖泡——貼不貼泡區(qū)作為研究區(qū),對(duì)該區(qū)多層介質(zhì)含水系統(tǒng)中地下水的石油類(lèi)污染進(jìn)行預(yù)測(cè)評(píng)價(jià),從而為此類(lèi)地下水污染的綜合治理提供科學(xué)依據(jù)和有效措施,確保油田的可持續(xù)發(fā)展。1貼不貼泡在北北環(huán)境地區(qū)的分布研究區(qū)位于石油城市大慶市喇嘛甸鎮(zhèn)西北,北緯51°76′～51°79′,東經(jīng)216°29′～216°34′。研究區(qū)為一小型閉流盆地,面積約12.28km2,東面與向榮屯、良種場(chǎng)接壤,西面與沈家北泡毗鄰,南面與三棱泡交界,北側(cè)與北十里相連(圖1)。區(qū)內(nèi)貼不貼泡呈北西—南東向分布,縱向長(zhǎng)約3km,橫向?qū)捈s1km,四周被風(fēng)沙、崗地圍成盆狀。地層巖性主要有四層,自下而上依次是:第三系泰康組和第四系白土山組砂礫巖和砂礫石承壓含水層,厚約80m;第四系荒山組飽水粘性土層(弱透水層),厚約51m;最上部為第四系哈爾濱組的細(xì)砂、亞粘土、亞砂土層,厚13～19m。其中,第四系荒山組弱透水層連接著上部的潛水和下部的承壓水,成為地下水和污染質(zhì)運(yùn)移的通道。區(qū)內(nèi)貼不貼泡已成為石油化工業(yè)的廢水排放處,石油類(lèi)含量在排污口處為30.60mg/l,湖泡中為0.5mg/l。其下部潛水被污染的范圍是以湖泡為中心,向四周呈橢圓形狀擴(kuò)散,構(gòu)成污染暈。其中心的污染質(zhì)濃度為0.5mg/l,向外逐漸降至0.03mg/l。2土試樣滲壓實(shí)驗(yàn)第四系哈爾濱組的細(xì)砂潛水層和第三系泰康組和第四系白土山組砂礫巖和砂礫石承壓含水層符合通常條件下的達(dá)西定律,由亞粘土、亞砂土組成的包氣帶土層則符合非飽和流滲透規(guī)律,而第四系荒山組飽水粘性土組成的弱透水層,其滲透規(guī)律與其它各層有明顯的不同。將飽水的原狀粘性土試樣置于滲壓容器中,利用滲壓儀進(jìn)行不同水頭差下的滲透實(shí)驗(yàn),獲得V～I(xiàn)關(guān)系曲線,即OABC折線(圖2),及V～I(xiàn)關(guān)系方程:線段OA:V1=K1I,(0<I<I01)(1)線段AB:V2=V01+K2(I-I01)(I01<I<I02)(2)線段BC:V3=V02+K3(I-I02),(I>I02)(3)式中I——水力梯度;I01——第2階段初始水力梯度(0.291);I02——第3階段初始水力梯度(0.440);V01——第2階段初始滲透速度(0.938×10-4m/d);V02——第3階段初始滲透速度(1.938×10-4m/d);K1——第1階段滲透系數(shù)(3.223×10-4m/d);K2——第2階段滲透系數(shù)(6.711×10-4m/d);K3——第3階段滲透系數(shù)(19.954×10-4m/d)3水系統(tǒng)中污染轉(zhuǎn)移值的模擬3.1地下水邊界的組成研究區(qū)內(nèi)含水介質(zhì)自下而上依次為:砂礫巖和砂礫石含水層,飽水粘性土層(弱透水層),細(xì)砂、亞粘土、亞砂土層。全區(qū)所有側(cè)向邊界均為第二類(lèi)邊界,其中承壓含水層側(cè)向邊界為已知流量邊界,其余為隔水邊界。頂部邊界湖泡水體為第一類(lèi)邊界,其余為入滲補(bǔ)給邊界。承壓含水層底板為底部隔水邊界。地下水在多孔介質(zhì)中的運(yùn)動(dòng)符合Darcy定律,污染質(zhì)濃度變化符合Fick定律,無(wú)化學(xué)反應(yīng)。地下水運(yùn)動(dòng)為非穩(wěn)定流,含水介質(zhì)為非均質(zhì)各向異性。3.2數(shù)學(xué)模型根據(jù)上述含水系統(tǒng)的概化,可建立起描述地下水含水系統(tǒng)的數(shù)學(xué)模型。該模型由兩部分組成:水分運(yùn)移模型和污染質(zhì)運(yùn)移模型。(1)導(dǎo)水率及保水性的計(jì)算為了將飽和帶和非飽和帶進(jìn)行統(tǒng)一描述,采用壓力水頭為狀態(tài)變量。數(shù)學(xué)模型為:{A(ψ)οψοt=οοx(chóng)(Κ(ψ)x)οψοx(chóng))+οοy(Κ(ψ)yοψοy)+οοz(Κ(ψ)z)οψοz)+οοz(Κ(ψ)z)+G-W?(4)(x,y,z)∈Ω?t>0ψ(x,y,z,t)|t=0=ψ0(x,y,z)(x,y,z)∈Ω?t=0ψ(x,y,z,t)|Γ1=ψ1(x,y,z?t)(x,y,z)∈Γ1?t>0Κ(ψ)οψο→n|Γ2=f1(x,y,z?t)(x,y,z)∈Γ2?t>0???????????????????????????????????????????A(ψ)οψοt=οοx(chóng)(K(ψ)x)οψοx(chóng))+οοy(K(ψ)yοψοy)+οοz(K(ψ)z)οψοz)+οοz(K(ψ)z)+G?W?(4)(x,y,z)∈Ω?t>0ψ(x,y,z,t)|t=0=ψ0(x,y,z)(x,y,z)∈Ω?t=0ψ(x,y,z,t)|Γ1=ψ1(x,y,z?t)(x,y,z)∈Γ1?t>0K(ψ)οψοn?∣∣Γ2=f1(x,y,z?t)(x,y,z)∈Γ2?t>0式中ψ——壓力水頭;G——補(bǔ)給項(xiàng)(如降水入滲等);W——為排泄項(xiàng)(如蒸發(fā)排泄、抽水等);ψ0——初始?jí)核^;ψ1——第一類(lèi)邊界上的壓力水頭;f1——第二類(lèi)邊界上的水分通量;→nn?——邊界外法線方向;Γ1、Γ2——分別為第一、第二類(lèi)邊界;其中?A(ψ)={οθοψ(ψ<0)μs(ψ≥0)(5)其中?A(ψ)={οθοψμs(ψ<0)(ψ≥0)(5)式中θ——有效孔隙度(飽和)或體積含水率(非飽和);μ*s——貯水率;其它符號(hào)同前。導(dǎo)水率K(ψ)由下式給出。Κ(ψ)={Κs[1-(αψ)n1-1?(1+(αψ)n1)-m]2[1+(αψ)n1]12m?(ψ<0)Κs(ψ≥0))(6)K(ψ)={Ks[1?(αψ)n1?1?(1+(αψ)n1)?m]2[1+(αψ)n1]12m?(ψ<0)Ks(ψ≥0))(6)式中α、n1、m——經(jīng)驗(yàn)系數(shù);Ks——飽和滲透系數(shù);其它符號(hào)同前。在實(shí)際計(jì)算中,式(6)中的Ks在飽和弱透水粘性土層中的取值由(1～3)式獲得:Κs={Κ10<Ι≤Ι01Κ2Ι01<Ι≤Ι02Κ3Ι≥Ι02(7)Ks=?????K1K2K30<I≤I01I01<I≤I02I≥I02(7)(2)初始污染濃度cz-cvy三維飽和-非飽和污染質(zhì)運(yùn)移的定解問(wèn)題為:{Rdο(θC)οt=οοx(chóng)[θ(DxxοCοx(chóng)+DxyοCοy+DxzοCοz-CVx)]+οοy[θ(DxyοCοx(chóng)+DyyοCοy+DyzοCοz-CVy)]+οοz[θ(DxzοCοx(chóng)+DyzοCοy+DzzοCοz-CVz)]+GCe-WC?(x,y,z)∈Ω?t>0C(x,y,z,t)|t=0=C0(x,y,z)?(x,y,z)∈Ω?t=0C(x,y,z,t)|Γ1=C1(x,y,z,t)?(x,y,z)∈Γ1?t>0θDοCο→n|Γ2=f2(x,y,z,t),(x,y,z)∈Γ2,t>0(8)式中C——污染濃度;C0——初始污染濃度;C1——一類(lèi)邊界上的污染濃度;Vx、Vy、Vz——x、y、z方向的滲透速度;Rd——阻滯因子;D——彌散系數(shù)張量;Ce——污染源的濃度;其它符號(hào)同上。3.3初支模型的建立先將研究區(qū)進(jìn)行平面三角剖分(圖3),共剖分為137個(gè)結(jié)點(diǎn),226個(gè)三角單元。垂向上分6層,7個(gè)層面分別為地表面、地下5m、地下10m、弱透水層頂板、弱透水層頂板以下8m、弱透水層底板和承壓含水層底板。采用三棱柱單元輸入信息,然后計(jì)算機(jī)程序自動(dòng)將計(jì)算區(qū)剖分成959個(gè)結(jié)點(diǎn),4068個(gè)四面體單元,邊界結(jié)點(diǎn)除湖水面上的結(jié)點(diǎn)為第一類(lèi)邊界點(diǎn)外,其余的均為第二類(lèi)邊界點(diǎn)。然后采用特征線法與有限單元法相結(jié)合來(lái)求解上述數(shù)學(xué)模型。計(jì)算時(shí)的初始流場(chǎng)和初始污染質(zhì)濃度場(chǎng)由專門(mén)布設(shè)的20個(gè)測(cè)壓管式觀測(cè)孔的實(shí)測(cè)值,經(jīng)泛Kriging空間插值獲得。由于計(jì)算區(qū)面積比較小,地層巖性平面上變化不大,故將計(jì)算區(qū)平面上作為一個(gè)參數(shù)區(qū),而在垂向上分出4個(gè)參數(shù)區(qū),即地面至地下5m、地下5m至弱透水層頂板、弱透水層和承壓含水層,參數(shù)初值參考已有水文地質(zhì)資料和室內(nèi)外實(shí)驗(yàn)測(cè)定值給出。以1995年11月20日實(shí)測(cè)地下水水頭及石油類(lèi)濃度為初始水動(dòng)力場(chǎng)和初始石油類(lèi)污染濃度場(chǎng),選取1995年11月20日～1996年7月9日期間的實(shí)測(cè)資料,分8個(gè)時(shí)段對(duì)所建模型進(jìn)行識(shí)別。經(jīng)反復(fù)計(jì)算,計(jì)算水頭與實(shí)測(cè)水頭之差的絕對(duì)值小于0.5m的觀測(cè)孔數(shù)占觀測(cè)孔總數(shù)的90%,計(jì)算石油類(lèi)濃度與實(shí)測(cè)濃度之差的絕對(duì)值小于0.1mg/l的觀測(cè)孔數(shù)占觀測(cè)孔總數(shù)的85%。最后得到的模型參數(shù)見(jiàn)表1。用1996年7月9日～11月20日的實(shí)測(cè)資料對(duì)模型進(jìn)行驗(yàn)證,證實(shí)所建模型可靠。4地下水中石油類(lèi)污染質(zhì)運(yùn)移特征利用所建立的大慶市貼不貼泡區(qū)地下水污染數(shù)值模擬模型即可對(duì)該區(qū)地下水中的石油類(lèi)污染質(zhì)的運(yùn)移做出預(yù)測(cè)。預(yù)測(cè)時(shí),源、匯項(xiàng)和邊界均按多年平均來(lái)處理,弱透水層的滲透系數(shù)由(1)～(3)式確定。下面分三個(gè)方案對(duì)地下水中石油類(lèi)污染質(zhì)的運(yùn)移做出預(yù)測(cè)評(píng)價(jià)。方案1保持現(xiàn)有條件,地下水中石油類(lèi)污染質(zhì)的濃度變化不大,10年后污染質(zhì)鋒面(0.025mg/l)未穿透弱透水層,但已污染了弱透水層4m多;方案2加大開(kāi)采目的層開(kāi)量,造成承壓水頭下降10m,此時(shí),污染質(zhì)擴(kuò)散范圍有所加大,但速度也比較慢,10年后污染質(zhì)鋒面(0.025mg/l)進(jìn)入弱透水層近8m;方案3減小開(kāi)采目的層開(kāi)采量,并將潛水和承壓水之間的水頭差保持在14.55m,結(jié)果是污染質(zhì)的下移基本被控制?？梢?jiàn),要防止下部?jī)?yōu)質(zhì)地下水源被污染,一是減少污染質(zhì)的排放量,降低湖泡水位和湖水中石油類(lèi)污染質(zhì)的含量;二是節(jié)約用水,減少承壓水的開(kāi)采量,縮小劣質(zhì)潛水與優(yōu)質(zhì)承壓水之間的水頭差。5儲(chǔ)層出水污染預(yù)測(cè)模型分析(1)多層介質(zhì)含水系統(tǒng)是自然界中普遍存在的。對(duì)其中地下水污染的預(yù)測(cè)評(píng)價(jià)需建立在對(duì)各層介質(zhì)滲透規(guī)律和地下水流場(chǎng)的認(rèn)知基礎(chǔ)上。(2)應(yīng)用地質(zhì)統(tǒng)計(jì)方法和數(shù)值模擬技術(shù),建立起的包括包氣帶、潛水、弱透水層和承壓水在內(nèi)的多層介質(zhì)含水系統(tǒng)中石油類(lèi)污染質(zhì)運(yùn)移的數(shù)值模擬模型,為解決此類(lèi)復(fù)雜含水介質(zhì)地區(qū)地下水污