羅布泊鹽湖羅北凹地礦區(qū)鹵水?dāng)?shù)值模擬研究

王凌芬，于詠梅，李博昀，張凡凱，鄧宇飛，李文學(xué)，欽賀

(1.中化地質(zhì)礦山總局地質(zhì)研究院， 北京 100101； 2.國(guó)投新疆羅布泊鉀鹽有限責(zé)任公司， 新疆 哈密 839000)

0 引言

羅布泊位于塔里木盆地東部，羅北凹地是從羅布泊分割出來(lái)的一個(gè)次級(jí)盆地(劉成林等，2002)，行政區(qū)劃隸屬新疆維吾爾自治區(qū)巴音郭勒蒙古自治州若羌縣，地理坐標(biāo)：東經(jīng) 90°40′00″～91°22′00″，北緯 40°33′00″～41°05′00″，南北長(zhǎng)約60 km，東西寬約25 km，面積為1534 km2(李文學(xué)等，2018)，是一個(gè)賦存多層含鉀鹵水的大型液體鉀鹽礦床(陳永志等，2001；王弭力等，2006；李浩等，2008；劉成林等，2009)。目前，羅北凹地液體鉀鹽礦以開(kāi)采潛鹵水為主，潛鹵水礦層平均厚度為17.54 m，礦層平均埋深1.89 m，礦層的平均孔隙度為18.58%，給水度11.36%，鹵水平均比重1.22，氯化鉀平均品位1.45%，氯化鉀孔隙度儲(chǔ)量8268.03萬(wàn)t、給水度儲(chǔ)量3894.97萬(wàn)t，占羅北凹地內(nèi)儲(chǔ)量的46.5%(李文學(xué)等，2018；劉成林等，2020)。隨著公司生產(chǎn)規(guī)模的擴(kuò)大，采鹵范圍逐漸擴(kuò)大到深部承壓鹵水，鹵水資源不斷消耗，地下鹵水水位也隨之下降(劉傳福，2010；韓積斌等，2015)，由于液體鉀鹽礦與固體鉀鹽礦的不同，鹵水的賦存特征和礦體的厚度及組分也隨之發(fā)生較大的變化(王新民等，2003；李文學(xué)等，2017；王凱等，2020)。實(shí)踐證明，含鹵層為較為復(fù)雜的層組，地下鹵水的非均質(zhì)性及富水性的不均一性較為明顯(焦鵬程等，2003)。近十年來(lái)，前人在礦區(qū)開(kāi)展了一系列勘查、研究工作。焦鵬程等(2003，2006)采用深部承壓含水層鹵水抽水試驗(yàn)、14C測(cè)年技術(shù)、同位素和化學(xué)示蹤技術(shù)，研究晶間鹵水循環(huán)速率、鹵水流向和流速及鉀富集機(jī)理等鹵水演化過(guò)程；李文學(xué)等(2021)利用大氣聯(lián)通法及輔助孔注入鹵水法提高采鹵井涌水量；作者所在團(tuán)隊(duì)于2019—2021年先后在研究區(qū)開(kāi)展疏干開(kāi)采抽水試驗(yàn)、承壓水監(jiān)測(cè)網(wǎng)建設(shè)、礦區(qū)資源儲(chǔ)量核實(shí)工作、儲(chǔ)鹵層鉀資源勘探和中深部鹵水鉀資源調(diào)查以及承壓鹵水開(kāi)采方法研究等工作，積累了豐富的數(shù)據(jù)資料和研究基礎(chǔ)，但目前對(duì)地下鹵水資源的管理、規(guī)劃和開(kāi)采仍處于較初級(jí)的人工計(jì)算階段，存在耗時(shí)長(zhǎng)、計(jì)算量大、效率低、精度不夠等問(wèn)題，本文開(kāi)展了礦區(qū)地下鹵水?dāng)?shù)值模擬研究，通過(guò)數(shù)值模型的建立，實(shí)現(xiàn)對(duì)區(qū)內(nèi)鉆孔、采鹵井及地質(zhì)剖面、水文地質(zhì)參數(shù)等原始數(shù)據(jù)的有效利用，細(xì)化對(duì)羅北凹地礦區(qū)地層結(jié)構(gòu)、各含水層空間分布及相互間水力聯(lián)系的認(rèn)識(shí)，進(jìn)行地下鹵水流場(chǎng)的模擬預(yù)測(cè)，指導(dǎo)下一步采鹵井的布置，實(shí)現(xiàn)地下鹵水的合理化開(kāi)采。

1 礦區(qū)地質(zhì)和水文地質(zhì)

礦區(qū)地表出露和鉆孔已揭露地層均為第四系化學(xué)鹽類(lèi)沉積和碎屑沉積，并呈互層狀循環(huán)交替產(chǎn)出，反映了本區(qū)咸化和淡化沉積環(huán)境的規(guī)律性變化，其中鹽類(lèi)沉積層是區(qū)內(nèi)地下富鉀鹵水的主要賦存層位(顧新魯?shù)龋?003，2009；袁文虎等，2020)。以往資料將羅北凹地含鹽系劃分為7個(gè)(編號(hào)為S1～S7)含鹽組(圖1，表1)，含鹽組分別對(duì)應(yīng)相應(yīng)的鹽層和碎屑層，富鉀鹵水主要賦存于各鹽層中，形成相應(yīng)的儲(chǔ)鹵層，其中羅北凹地淺部(一般小于90 m)劃分出W1～W4四個(gè)儲(chǔ)鹵層(1個(gè)潛鹵水層、3個(gè)承壓鹵水層)。

2 模型概化及建立

2.1 地下鹵水滲流模型的建立

研究區(qū)地下水流整體上以水平運(yùn)動(dòng)為主、垂向運(yùn)動(dòng)為輔，地下水系統(tǒng)符合質(zhì)量守恒定律和能量守恒定律；在常溫常壓下地下水運(yùn)動(dòng)符合達(dá)西定律；考慮含水層之間的流量交換，地下水運(yùn)動(dòng)可以概化為空間三維流；地下水系統(tǒng)的垂向運(yùn)動(dòng)是由層間水頭差異引起的；地下水系統(tǒng)的輸入輸出隨時(shí)間、空間變化，故地下水為非穩(wěn)定流；參數(shù)隨空間變化，體現(xiàn)了系統(tǒng)的非均質(zhì)性，所以含水介質(zhì)概化為非均質(zhì)各向同性介質(zhì)。綜上所述，研究區(qū)可概化成非均質(zhì)、水平方向各向同性、垂向存在變異、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)，即地下水系統(tǒng)的概念模型。

圖1 羅北凹地含鹽系剖面圖(據(jù)新疆地礦局第二水文工程地質(zhì)大隊(duì)，2006①修改)

表1 羅北凹地含鹽系劃分表

(1)

式(1)中：Ω—滲流區(qū)域；h—地下水系統(tǒng)的水位標(biāo)高(m)；K—含水介質(zhì)的水平滲透系數(shù)(m/d)；Kz—含水介質(zhì)垂向滲透系數(shù)(m/d)；ε—含水層的源匯項(xiàng)(1/d)；h0—初始水位(m)；S—自由面以下含水層儲(chǔ)水率(1/m)；Γ0—滲流區(qū)域的上邊界，即地下水的自由表面；μ—潛水含水層重力給水度；p—潛水面的蒸發(fā)和降水入滲強(qiáng)度等(m/d)；Γ1—已知水位邊界；h1—已知邊界水位值(m)；Γ2—滲流區(qū)域的流量邊界；Kn—邊界面法線方向的滲透系數(shù)(m/d)，n—邊界面的法線方向；q—Γ2邊界的單位面積上的流量(m/d)，流入為正，流出為負(fù)，隔水邊界為0。

采用基于有限差分法的GMS軟件對(duì)上述數(shù)學(xué)方程進(jìn)行求解，地下水模擬系統(tǒng)(Groundwater Modeling System)，簡(jiǎn)稱GMS，是由美國(guó)Brigham young University的環(huán)境模型研究實(shí)驗(yàn)室在綜合Modflow、Modpath等已有地下水模型基礎(chǔ)上研發(fā)而成的，是一個(gè)具有綜合性、用于地下水模擬的圖形界面軟件。基于GMS軟件的Solid模塊建立的水文地質(zhì)結(jié)構(gòu)模型(圖2、圖3)，為分析研究區(qū)的沉積特征提供了方法和思路，為建立地下水?dāng)?shù)值模擬模型奠定了基礎(chǔ)。

整體上說(shuō)，模型在垂向上分為7層，包含4層含水層和3層隔水層。將模擬區(qū)剖分為220行、200列規(guī)則網(wǎng)格，各層均采用300 m×300 m網(wǎng)格大小進(jìn)行剖分，共剖分網(wǎng)格數(shù)目為118201個(gè)，其中，第一層有效單元(活動(dòng)網(wǎng)格)17077個(gè)，第二至第七層的活動(dòng)單元格數(shù)目均為16854個(gè)。

圖2 羅北凹地水文地質(zhì)結(jié)構(gòu)橫向剖面示意圖W1—潛鹵水層；W2—第一承壓鹵水層；W3—第二承壓鹵水層；W4—第三承壓鹵水層

圖3 羅北凹地水文地質(zhì)結(jié)構(gòu)縱向剖面示意圖W1—潛鹵水層；W2—第一承壓鹵水層；W3—第二承壓鹵水層；W4—第三承壓鹵水層

2.2 邊界條件及水文地質(zhì)參數(shù)的確定

圖4 研究區(qū)邊界條件概化圖

根據(jù)區(qū)內(nèi)流場(chǎng)特征和對(duì)地層結(jié)構(gòu)的分析，將研究區(qū)側(cè)向邊界確定為三種類(lèi)型(圖4)。其一是流量邊界，為北部、南部和西部臺(tái)地與凹地之間的自然邊界。這類(lèi)邊界接受山區(qū)融雪和降雨入滲到出山口的側(cè)向補(bǔ)給、西臺(tái)地的側(cè)向補(bǔ)給以及鹽田滲漏的側(cè)向補(bǔ)給。補(bǔ)給量的大小按照各塊段水文地質(zhì)條件和滲透系數(shù)大小有所不同，具體在計(jì)算過(guò)程中進(jìn)行計(jì)算分析；其二是定水頭邊界，將模型的南部邊界根據(jù)等水位線分布圖處理成定水頭邊界，從等水位線圖上看，南部邊界處為地下水位的高值區(qū)，呈現(xiàn)出類(lèi)似地下水分水嶺的特征，將其處理為一類(lèi)；其三是通用水頭邊界，即第三類(lèi)邊界條件，羅北凹地與騰龍臺(tái)地的邊界處，受附近采鹵井的影響，該邊界定義為通用水頭邊界。研究區(qū)上邊界為潛水含水層自由水面，模型底界以鉆孔揭穿最大深度處W4的底板作為模型的底部邊界，處理為隔水邊界。

在建模工作中，首先根據(jù)水文地質(zhì)條件和野外工作獲得的水文地質(zhì)參數(shù)，按參數(shù)分區(qū)給定參數(shù)初值，通過(guò)水位擬合進(jìn)行參數(shù)識(shí)別，最后確定各參數(shù)分區(qū)值，此處以潛水含水層滲透系數(shù)、孔隙度、給水度、第一承壓含水層彈性釋水系數(shù)為例(圖5～圖8)，其他各層不一一羅列。

圖5 研究區(qū)潛水含水層滲透系數(shù)分區(qū)圖

圖6 研究區(qū)潛水含水層給水度分區(qū)圖

圖7 研究區(qū)潛水含水層孔隙度分區(qū)圖

圖8 研究區(qū)第一承壓含水層釋水系數(shù)分區(qū)圖

2.3 模擬期確定及定解條件

根據(jù)現(xiàn)有各均衡要素資料和水位動(dòng)態(tài)資料情況，確定2019年11月—2020年11月這一個(gè)完整的水文年為模擬期，以月為單位共劃分為12個(gè)應(yīng)力期。

初始條件：以2019年11月地下水水位資料，采用內(nèi)插法和外推法獲得潛水含水層和第一承壓含水層的初始水位。根據(jù)現(xiàn)有的第二和第三承壓含水層水位觀測(cè)數(shù)據(jù)，各承壓含水層間水力聯(lián)系密切，故第二承壓含水層和第三承壓含水層初始流場(chǎng)采用第一承壓含水層初始流場(chǎng)。

2.4 模型識(shí)別與檢驗(yàn)

對(duì)研究區(qū)內(nèi)長(zhǎng)觀孔的實(shí)測(cè)水位和模擬水位擬合誤差進(jìn)行了分析。其中潛水含水層擬合的28個(gè)鉆孔中，平均誤差絕對(duì)值小于0.5 m的鉆孔有20個(gè)，平均誤差絕對(duì)值在0.5～1 m之間的鉆孔有7個(gè)，兩者占總數(shù)的96.43%；平均誤差絕對(duì)值大于 1 m 的鉆孔有1個(gè)，占總數(shù)的3.57%。承壓含水層長(zhǎng)觀孔點(diǎn)位較少，僅有11個(gè)，平均誤差絕對(duì)值小于0.5 m的鉆孔有8個(gè)，平均誤差絕對(duì)值在0.5～1 m 之間的鉆孔有3個(gè)，兩者占總數(shù)的100%。擬合誤差偏大的原因主要有以下三點(diǎn)：(1)誤差較大的地區(qū)集中在開(kāi)采降落漏斗附近，資料誤差比較大；(2)部分控制點(diǎn)受生產(chǎn)用井影響，觀測(cè)水位變動(dòng)過(guò)大，造成誤差偏大；(3)實(shí)測(cè)數(shù)據(jù)使用的是月平均值，計(jì)算導(dǎo)致誤差。

圖9 模擬期末潛水含水層流場(chǎng)擬合圖

潛水含水層和承壓含水層的流場(chǎng)擬合見(jiàn)圖9和圖10，典型長(zhǎng)觀孔水位過(guò)程線擬合情況見(jiàn)圖11和圖12。計(jì)算流場(chǎng)與實(shí)際流場(chǎng)的擬合程度較好，計(jì)算流場(chǎng)基本上反映了地下水的流動(dòng)與趨勢(shì)，含水層識(shí)別模型所對(duì)應(yīng)的水文地質(zhì)模型與實(shí)際水文地質(zhì)模型是基本吻合的，邊界概化較為合理，源匯項(xiàng)處理較正確，參數(shù)分區(qū)較為合理，參數(shù)基本正確。

圖10 模擬期末承壓含水層流場(chǎng)擬合圖

圖11 潛水觀測(cè)孔水位擬合曲線圖

圖12 承壓水觀測(cè)孔水位擬合曲線圖

3 鹵水開(kāi)采預(yù)測(cè)

利用所建立的地下水?dāng)?shù)值模型，結(jié)合現(xiàn)有開(kāi)采情況，預(yù)測(cè)了地下水流場(chǎng)變化趨勢(shì)，對(duì)地下水開(kāi)采降深進(jìn)行了評(píng)價(jià)。本次預(yù)測(cè)將2020年11月的統(tǒng)測(cè)水位確定的流場(chǎng)作為預(yù)測(cè)模型的初始流場(chǎng)，預(yù)測(cè)年限為3年和5年，即從2020年11月到2025年11月。圖13和圖14為3年末潛水位和承壓水位等值線圖；圖15和圖16分別為5年后潛水位和承壓水位等值線圖；圖17和圖18為3年后潛水水位降深和承壓水水位降深等值線圖；圖19和圖20為5年后潛水水位降深和承壓水水位降深等值線圖。

圖13 預(yù)測(cè)2023年11月潛水等水位線分布

圖14 預(yù)測(cè)2023年11月承壓水等水位線分布

從圖13與圖14可見(jiàn)，研究區(qū)潛水地下水水流方向總體上和初始流場(chǎng)形態(tài)基本一致，而承壓含水層在開(kāi)采井附近流場(chǎng)形態(tài)發(fā)生較大變化，出現(xiàn)明顯的降落漏斗。

經(jīng)3年開(kāi)采后，研究區(qū)潛水總體降深2.00～4.00 m左右(正為水位下降)，其中個(gè)別開(kāi)采區(qū)中心的最大降深可達(dá)到6.00 m以上。降深面積分布較為廣泛，開(kāi)采區(qū)降深2.00 m以上區(qū)域面積達(dá)376.77 km2。而開(kāi)采5年后降深2.00m以上區(qū)域面積達(dá)718.60 km2。

圖16 預(yù)測(cè)2025年11月承壓水等水位線分布

圖15 預(yù)測(cè)2025年11月潛水等水位線分布

圖17 預(yù)測(cè)2023年11月潛水水位降深分布圖

圖18 預(yù)測(cè)2023年11承壓水水位降深分布圖

圖19 預(yù)測(cè)2025年11月潛水水位降深分布圖

圖20 預(yù)測(cè)2025年11月承壓水水位降深分布圖

經(jīng)3年開(kāi)采后，研究區(qū)承壓水總體降深0.00～4.00 m左右，其中個(gè)別開(kāi)采區(qū)中心的最大降深可達(dá)到10.00 m以上。降深面積分布較為廣泛，開(kāi)采區(qū)降深4.00 m以上區(qū)域面積達(dá)194.11 km2。而開(kāi)采5年后降深4.00 m以上區(qū)域面積達(dá)338.00 km2。

4 結(jié)論

(1)使用GMS軟件建立水文地質(zhì)結(jié)構(gòu)模型，將模型在垂向上概化為7層，給出了參數(shù)分區(qū)及參數(shù)初值，確定的水文地質(zhì)概念模型為非均質(zhì)、水平方向各向同性垂向存在變異、空間三維結(jié)構(gòu)、非穩(wěn)定地下水流系統(tǒng)。

(2)從模擬期末等值線擬合、典型長(zhǎng)觀孔水位過(guò)程線擬合、水文地質(zhì)參數(shù)以及研究區(qū)地下水系統(tǒng)均衡分析來(lái)看，研究區(qū)地下水流數(shù)值模擬模型反映了羅北凹地地下水流動(dòng)規(guī)律和特征，符合研究區(qū)實(shí)際的水文地質(zhì)條件，反映了研究區(qū)內(nèi)地下水的運(yùn)動(dòng)特點(diǎn)和動(dòng)態(tài)變化趨勢(shì)，基本達(dá)到模型精度要求，模型識(shí)別和均衡計(jì)算結(jié)果與實(shí)際水文地質(zhì)條件相符。

(3)通過(guò)模型預(yù)測(cè)，得到開(kāi)采3年和5年后的地下水流場(chǎng)與降深圖。研究區(qū)潛水地下水水流方向總體上和初始流場(chǎng)形態(tài)基本一致，而承壓含水層在開(kāi)采井附近流場(chǎng)形態(tài)發(fā)生較大變化，出現(xiàn)明顯的降落漏斗。后期在開(kāi)采井部署時(shí)，可考慮遠(yuǎn)離降落漏斗分布區(qū)，在水位降深幅度較小地區(qū)加密采鹵井設(shè)計(jì)，這樣才有利于鹵水的合理采出，減小生產(chǎn)成本。

注 釋

① 新疆地礦局第二水文工程地質(zhì)大隊(duì).2006.新疆若羌縣羅北凹地鉀鹽礦詳查報(bào)告[R].