三江平原撓力河流域地下水三維數(shù)值模擬研究

三江平原撓力河流域地下水三維數(shù)值模擬研究

1地下水資源生態(tài)保護的意義三江平原是中國的一個重要糧食和能源基地，為中國的社會發(fā)展和經(jīng)濟發(fā)展做出了重要貢獻。但隨著三江平原區(qū)域社會和經(jīng)濟的迅猛發(fā)展,國民經(jīng)濟各部門各行業(yè)水利的保障要求越來越高,對水資源的需求量越來越大,對水環(huán)境的保護要求也越來越迫切。由于多年水利發(fā)展的滯后和水利工程建設(shè)的薄弱,給三江平原的水資帶來一系列問題:水利的安全保證能力低;水資源利用缺乏統(tǒng)一規(guī)劃,水資源供需矛盾突出;水環(huán)境惡化,生態(tài)環(huán)境質(zhì)量下降;平原內(nèi)地下水水位急劇下降,導(dǎo)致地面沉降。因此對三江平原地下水資源進行分析與評價具有重要的現(xiàn)實意義和指導(dǎo)意義。由于三江平原部分地區(qū)觀測井和水文地質(zhì)資料缺乏,本文主要以三江平原撓力河流域作為研究對象,對三江平原撓力河流域的地下水系統(tǒng)進行了系統(tǒng)分析,建立了三江平原撓力河流域水文地質(zhì)概念模型,并在水文地質(zhì)概念模型的基礎(chǔ)上,建立了三江平原撓力河流域地下水?dāng)?shù)學(xué)模型,還采用了國際通用的地下水?dāng)?shù)值模擬軟件VisualModflow,對撓力河流域地下水進行模擬仿真。2撓力河流域概況三江平原撓力河流域位于黑龍江省東部地區(qū),是三江平原的四大流域之一。地理坐標(biāo)為東經(jīng)131°31′~134°10′,北緯45°43′~47°35′。在行政區(qū)上包括寶清縣全部,友誼、富錦、集賢、雙鴨山和饒河五市縣部分以及國營農(nóng)場總局的建三江和紅興隆兩個國營農(nóng)場管理局所屬五九七、八五二、八五三、八五九、紅旗嶺、紅衛(wèi)、大興和勝利等共17個大型農(nóng)場,總面積為23219km2,其概況圖如圖1所示。撓力河流域?qū)儆谥袦貛Т箨懠撅L(fēng)氣候的半濕潤地帶,夏季高溫多雨,冬季干冷而漫長。根據(jù)本區(qū)歷年氣象資料統(tǒng)計,其年平均氣溫為2.7℃,其年平均降水量為505mm。該流域在大地構(gòu)造單元上屬合江內(nèi)陸中斷陷,自第三紀以來一直處于間歇性的沉陷,因而沉積了巨厚的第三系鶴立組、寶泉嶺組與富錦組泥巖、砂巖及砂礫巖和第四系砂、砂礫石和礫卵石層,其中第四系砂、砂礫石和礫卵石的孔隙中賦存有豐富的松散巖類孔隙水,第三系砂巖和砂礫巖的孔隙裂隙中賦存有孔隙裂隙水。另外在區(qū)內(nèi)殘山殘丘及山前臺地區(qū)分布的前第四系基巖裂隙中也賦存基巖裂隙水。撓力河流域的經(jīng)濟支柱主要是煤炭開采、火力發(fā)電和以農(nóng)墾總局為龍頭的商品糧生產(chǎn)及精加工業(yè)。3重力河流域地下水流模型的建立3.1地下水含水層結(jié)構(gòu)撓力河流域從平面上分析,沿山前—主流帶—兩側(cè)—泛流帶,其含水層由薄變厚,顆粒由粗變細,但河口附近含水層較薄;流域內(nèi)七星河、撓力河和烏蘇里江與淺層地下水水力聯(lián)系密切,地下水呈放射狀補給河流。從垂向分析,根據(jù)已有的研究結(jié)果,撓力河在垂向上的影響深度在350m以內(nèi),淺層和深層之間有層粘土組成的較穩(wěn)定的隔水層,厚度一般在15~50m,使得淺層和深層地下水存在水力關(guān)系,因此在垂直方向上需將含水層分成淺層含水層、深層含水層和弱透水層。其中淺層含水層為第四系沖積、沖洪積和河流及湖積相沉積物,分布在出山口河谷、沖積扇及沖積平原區(qū),屬孔隙含水層;深層含水層為下古生界寒武紀和奧陶紀碳酸鹽地層,屬于裂隙巖溶含水層。在覆蓋及隱伏區(qū)為承壓水,在裸露區(qū)為微承壓或無壓水。由于目前還沒有分層承壓含水層分布、水位和流量等資料,所以將研究區(qū)內(nèi)深層含水層概化為一個承壓含水層,模型所描述的承壓含水層的水力特征和參數(shù)等均為研究區(qū)內(nèi)所有深層含水層的等效值。弱透水層在淺層含水層和深層含水層之間,沉積有不連續(xù)弱透水夾層。本文將弱透水夾層概化為一個相對弱透水層,即越流層。3.2邊境條件的總結(jié)3.2.1發(fā)生水量交換模擬區(qū)的頂部以潛水面為界,通過該邊界淺層地下水與外部環(huán)境發(fā)生水量交換,如地下水接受大氣降水的入滲、農(nóng)田灌溉回滲和河渠入滲補給,以及潛水蒸發(fā)排泄等。承壓含水層底部分布著粘土和亞粘土弱透水層,其滲透性較差,可將它概化為隔水邊界。3.2.2層的邊界總結(jié)3.2.2.河流邊界性質(zhì)撓力河流域東南、南和西南邊界為完達山山脈,淺層地下水在該邊界與外界沒有水流交換,可視為隔水邊界;東北、北部及西北邊界均為撓力河和七星河以及蓮花河和別拉紅河地表水分水嶺,可概化為給定水頭邊界。通過系統(tǒng)分析表明,撓力河和七星河為模擬區(qū)內(nèi)的優(yōu)勢河流,主河道與淺層地下水有大量的水量交換,可概化為河流邊界。另外,東部邊界烏蘇里江也作為河流邊界考慮。小撓力河水位在每個應(yīng)力期內(nèi)變化不大,可視其為給定水頭邊界。淺層含水層邊界概化結(jié)果如圖2所示。3.2.2.流量邊界的確定撓力河流域深層地下水與外界存在一定的緊密水力關(guān)系,其東南、南部和西南邊界接受外界徑流補給,為徑流補給邊界,東北、北部和西北邊界為地下水排泄邊界。這些都可視為給定流量邊界。根據(jù)達西定律計算各邊界流量,以井流形式分配到相應(yīng)計算單元。3.3淺層地下水排放撓力河流域淺層地下水主要接受大氣降水入滲補給、農(nóng)田灌溉回滲補給及撓力河、七星河和烏蘇里江滲漏補給。工農(nóng)業(yè)及生活用水開采、淺層含水層的蒸發(fā)及越流排泄是淺層地下水的主要排泄方式。淺層含水層的越流補給和側(cè)向徑流補給是深層地下水主要補給來源,工農(nóng)業(yè)及生活用水開采和側(cè)向徑流排泄是深層地下水的主要排泄方式。3.4制度內(nèi)邊界面根據(jù)以上綜述,將撓力河流域地下水流系統(tǒng)概化為非均質(zhì)各向同性三維非穩(wěn)定地下水流系統(tǒng),其數(shù)學(xué)方程為:式中:?表示模擬滲流區(qū)域;(x,y,z)表示空間位置坐標(biāo);t表示時間;S0表示潛水面;S1表示給定水頭邊界;S2表示定流量邊界面;S3表示河流邊界面;H(x,y,z,t)為模擬滲流區(qū)內(nèi)的水頭分布;H0(x,y,z)表示初始時刻滲流區(qū)內(nèi)及邊界上的水頭分布;H1(x,y,z,t)表示滲流區(qū)第一類邊界的水頭函數(shù);HRIV為第三類邊界條件的河水位;q表示滲流區(qū)流量邊界上的單位面積流量,其中隔水邊界流量為零;n表示邊界的外法線方向;Kn表示邊界法線方向的滲透系數(shù);Kx、Ky和Kz表示在x、y和z方向含水層的滲透系數(shù);Ss表示自由面以下含水層的儲水率;μ表示潛水含水層在潛水面上的重力給水度;σ表示為河床堆積物的阻尼系數(shù),σ=M/Kz,其中M為河床堆積物的厚度,P為潛水面單位時間單位面積補入或排泄的水體積,包括降水入滲和蒸發(fā)等;W為單位時間單位體積含水層得到或失去的水量,用以代表匯源項。3.5數(shù)值模型的構(gòu)建3.5.1剖面分區(qū)和一般意義的分離本文采用國際通用的地下水模擬軟件VisualModflow對撓力河流域地下水進行數(shù)值模擬。全區(qū)按照長方形網(wǎng)格進行剖分,共剖分40行和45列,垂向上剖分三層,共剖分出48762個單元,如圖2所示。潛水含水層滲透系數(shù)和給水度有4個分區(qū),承壓含水層滲透系數(shù)和貯水率有3個分區(qū),各分區(qū)參數(shù)的初值如圖3和圖4所示。淺層含水層和深層含水層的初始流場是依據(jù)2005年1月1日的觀測井?dāng)?shù)據(jù)及模型邊界條件的水頭數(shù)據(jù),按照克里克插值獲得的,其初始流場如圖5和圖6所示。3.5.2地下水開采期3月采用間接法識別和驗證模型中的參數(shù)。以2005年1月1日至2005年12月31日作為模型識別時段,以2006年1月1日至2006年12月31日作為模型驗證時段,將每個月作為一個地下水開采期。3.5.3水文地質(zhì)特征根據(jù)水文地質(zhì)模型所建立的數(shù)學(xué)模型應(yīng)該反映實際流場的特點,因此在進行模擬預(yù)報前,必須對數(shù)學(xué)模型進行識別,即利用建立的數(shù)學(xué)模型,計算觀測孔所在單元的水頭,并和實測水頭進行對比,從而反求有關(guān)的水文地質(zhì)參數(shù)。本次計算選選擇2005年1月1日~12月31日的地下水動態(tài)監(jiān)測資料進行模型的識別,該時段歷經(jīng)枯水期、平水期和豐水期,地下水流場的特征能較好地反映含水層系統(tǒng)的特征。識別過程中,計算水位和觀測水位的擬合平均誤差為0.305m,絕對誤差的平均值為2.015,標(biāo)準均方根為1.869%。圖7和圖8分別為淺水含水層和深層含水層在5月31日的識別水位和實測水文擬合而成。從圖中可以看出,模擬期等水位線與實測水位線模擬效果較好,符合實際情況。識別后的淺層含水層、深層含水層水文地質(zhì)參數(shù)如表1和表2所示,弱透水層滲透系數(shù)kx=ky=0.01,kz=0.001,釋水系數(shù)s=0.001。為了進一步檢驗所建立的數(shù)學(xué)模型和模型識別后確定的水文地質(zhì)參數(shù)的可靠性,利用2006年1月1日~12月31日的地下水位動態(tài)監(jiān)測資料對模型進行檢驗,檢驗過程中模型平均誤差為0.052m,誤差絕對值的平均值為4.057,均方根為2.121%,表示該模型計算水位與觀測水位之間有較好的適應(yīng)性。淺水含水層和深層含水層等水位線擬合過程顯示,計算水位與實測水位等值線的整體擬合程度良好,計算水位基本上反映了實際水位空間的變化規(guī)律。模型識別和檢驗結(jié)果證明所建立的數(shù)學(xué)模型、邊界條件、水文地質(zhì)參數(shù)和源匯項的確定都是符合實際情況的,可以應(yīng)用該模型進行地下水動態(tài)的預(yù)測。4建立模型及計算邊界流量根據(jù)已掌握的地下水動態(tài)觀測資料,選擇了距本次模擬計算時間較近的2006年12月31日的地下水位作為預(yù)報時段的初始水位,根據(jù)觀測孔的實測水位資料采用克里克方法進行插值,然后導(dǎo)入模型作為預(yù)報時段的初始流場。模型預(yù)報前要先設(shè)計開采方案以及預(yù)報邊界條件和源匯項。二類邊界的預(yù)報采用迭代逼近方法,即先給出水力梯度的初估值,由達西定律計算相應(yīng)的邊界流量,運行預(yù)報模型,然后根據(jù)預(yù)報結(jié)果計算出二類邊界附近的水力梯度,運用達西定律計算邊界流量,再進行預(yù)報。重復(fù)進行上述迭代過程,直到相鄰兩次迭代中的邊界流量之差的絕對值小于允許誤差為止。最后一次迭代過程中的流量即為所求邊界流量。研究區(qū)源匯項主要是降水入滲和蒸發(fā)量。本文采用改進的BP神經(jīng)網(wǎng)絡(luò)對2009~2020年撓力河流域的降水量和蒸發(fā)量進行預(yù)測,同時,采用Elanm小波神經(jīng)網(wǎng)絡(luò)模型對地下水開采量進行預(yù)測。把用隨機方法預(yù)測出的降水量、蒸發(fā)量和地下水開采量代入到確定性模型中進行預(yù)報,實現(xiàn)了隨機模型與確定性模型相耦合,建立了確定性-隨機耦合模型。以下采用該模型進行撓力河流域地下水進行預(yù)報。其2020年12月31日預(yù)報的地下水流場圖見圖9。從預(yù)報的等水位線圖可以看出,按照預(yù)測的開采量進行開采,研究區(qū)的地下水位在預(yù)報期間下降較大,所算得的2007~2020年撓力河流域地下水系統(tǒng)水位平均降深1.48m,平均年降深0.1057m,因此,應(yīng)減少該流域的地下水的開采量,增加地表水的開采量,以免引發(fā)地下水位下降、地面下沉和水環(huán)境惡化等現(xiàn)象的發(fā)生。5建立模型及預(yù)測結(jié)果本文以三江平原撓力河流域地下水為研究對象,根據(jù)撓力河流域的地質(zhì)和水文地質(zhì)資料,建立了撓力河水文地質(zhì)概念模型,在水文地質(zhì)概念模型的基礎(chǔ)上建立了地下水三維數(shù)學(xué)模型。其次,對撓力河流域地下水三維數(shù)學(xué)模型進行分析,采用國際通用的地下水?dāng)?shù)值模擬軟件VisualModflow建立了數(shù)值模型,并對該模型進行識別和檢驗,確定了水文地質(zhì)參數(shù)。再其次,本文采用改