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

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

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

0 引言

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

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

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

2 模型概化及建立

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

研究區(qū)地下水流整體上以水平運(yùn)動為主、垂向運(yùn)動為輔，地下水系統(tǒng)符合質(zhì)量守恒定律和能量守恒定律；在常溫常壓下地下水運(yùn)動符合達(dá)西定律；考慮含水層之間的流量交換，地下水運(yùn)動可以概化為空間三維流；地下水系統(tǒng)的垂向運(yùn)動是由層間水頭差異引起的；地下水系統(tǒng)的輸入輸出隨時間、空間變化，故地下水為非穩(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—自由面以下含水層儲水率(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軟件對上述數(shù)學(xué)方程進(jìn)行求解，地下水模擬系統(tǒng)(Groundwater Modeling System)，簡稱GMS，是由美國Brigham young University的環(huán)境模型研究實(shí)驗(yàn)室在綜合Modflow、Modpath等已有地下水模型基礎(chǔ)上研發(fā)而成的，是一個具有綜合性、用于地下水模擬的圖形界面軟件。基于GMS軟件的Solid模塊建立的水文地質(zhì)結(jié)構(gòu)模型(圖2、圖3)，為分析研究區(qū)的沉積特征提供了方法和思路，為建立地下水?dāng)?shù)值模擬模型奠定了基礎(chǔ)。

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

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

在建模工作中，首先根據(jù)水文地質(zhì)條件和野外工作獲得的水文地質(zhì)參數(shù)，按參數(shù)分區(qū)給定參數(shù)初值，通過水位擬合進(jìn)行參數(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)有各均衡要素資料和水位動態(tài)資料情況，確定2019年11月—2020年11月這一個完整的水文年為模擬期，以月為單位共劃分為12個應(yīng)力期。

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

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

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

圖9 模擬期末潛水含水層流場擬合圖

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

圖10 模擬期末承壓含水層流場擬合圖

圖11 潛水觀測孔水位擬合曲線圖

圖12 承壓水觀測孔水位擬合曲線圖

3 鹵水開采預(yù)測

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

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

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

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

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

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

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

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

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

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

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

經(jīng)3年開采后，研究區(qū)承壓水總體降深0.00～4.00 m左右，其中個別開采區(qū)中心的最大降深可達(dá)到10.00 m以上。降深面積分布較為廣泛，開采區(qū)降深4.00 m以上區(qū)域面積達(dá)194.11 km2。而開采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ì)參數(shù)以及研究區(qū)地下水系統(tǒng)均衡分析來看，研究區(qū)地下水流數(shù)值模擬模型反映了羅北凹地地下水流動規(guī)律和特征，符合研究區(qū)實(shí)際的水文地質(zhì)條件，反映了研究區(qū)內(nèi)地下水的運(yùn)動特點(diǎn)和動態(tài)變化趨勢，基本達(dá)到模型精度要求，模型識別和均衡計(jì)算結(jié)果與實(shí)際水文地質(zhì)條件相符。

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

注 釋

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