基于多源數(shù)據(jù)的弱透水層水文地質(zhì)參數(shù)反演研究
——以呼和浩特盆地某淤泥層為例

石鴻蕾，郝奇琛，邵景力，崔亞莉，張秋蘭

（1.中國(guó)地質(zhì)大學(xué)（北京）水資源與環(huán)境學(xué)院，北京 100083；2.中國(guó)地質(zhì)科學(xué)院水文地質(zhì)環(huán)境地質(zhì)研究所，河北 石家莊 050061）

水文地質(zhì)參數(shù)對(duì)于區(qū)域地下水資源評(píng)價(jià)、數(shù)值模擬及預(yù)報(bào)、開(kāi)發(fā)利用與保護(hù)及科學(xué)管理具有重要意義[1]。當(dāng)實(shí)際含水介質(zhì)的水文地質(zhì)參數(shù)不易獲取時(shí)，通過(guò)已知水文地質(zhì)結(jié)構(gòu)、參數(shù)、源匯條件及觀測(cè)資料反演未知參數(shù)通常是一種十分有效的方法[1-3]。當(dāng)前水文地質(zhì)參數(shù)的反演方法主要包括解析法和數(shù)值法。隨著計(jì)算機(jī)技術(shù)的發(fā)展，地下水?dāng)?shù)值模擬技術(shù)逐漸成熟，數(shù)值法逐漸成為解決地下水相關(guān)問(wèn)題的重要方法[4-5]，用于水文地質(zhì)參數(shù)求取、地下水資源評(píng)價(jià)、預(yù)報(bào)等，為地下水開(kāi)發(fā)利用及科學(xué)管理提供有效的工具和重要依據(jù)。

目前已有不少學(xué)者基于數(shù)值法進(jìn)行了水文地質(zhì)參數(shù)的反演研究[6-9]。毛喜云[10]結(jié)合群井抽水試驗(yàn)結(jié)果，建立地下水流數(shù)值模型對(duì)哈爾濱地區(qū)典型河漫灘二元結(jié)構(gòu)地層的水文地質(zhì)參數(shù)進(jìn)行反演。李貴仁[11]建立地下水流數(shù)值模型，通過(guò)水位擬合對(duì)安徽某鐵礦區(qū)巖溶裂隙含水層的水文地質(zhì)參數(shù)進(jìn)行反演。SzabóNP[12]建立了一種遺傳算法輔助反演方法反演非飽和帶地層的水文地質(zhì)參數(shù)。Moharir 等[13]利用MODFLOW軟件對(duì)玄武巖含水層參數(shù)進(jìn)行抽水試驗(yàn)反演建模。以往研究大都通過(guò)建立地下水流數(shù)值模擬模型，通過(guò)水位擬合反演參數(shù)。由于反演問(wèn)題是根據(jù)結(jié)果間接估計(jì)的，數(shù)值模型中可能存在異參同效的問(wèn)題，因此基于單一數(shù)據(jù)源的反演方法具有較強(qiáng)的不確定性。另外，以往研究大都是反演含水層的水文地質(zhì)參數(shù)，針對(duì)弱透水層的參數(shù)反演研究很少。但對(duì)于多層含水系統(tǒng)，弱透水層的水文地質(zhì)參數(shù)是計(jì)算層間越流量、預(yù)測(cè)和評(píng)價(jià)地面沉降的關(guān)鍵參數(shù)，對(duì)于指導(dǎo)地下水合理開(kāi)發(fā)利用具有重要意義[14]。

針對(duì)以上問(wèn)題，本文在野外抽水試驗(yàn)和溶質(zhì)運(yùn)移試驗(yàn)的基礎(chǔ)上，利用MODFLOW-USG 建立地下水流-溶質(zhì)運(yùn)移耦合模型，基于地下水位和溶質(zhì)濃度多源數(shù)據(jù)對(duì)弱透水層水文地質(zhì)參數(shù)進(jìn)行反演。相較于單一源數(shù)據(jù)反演，這種方法在一定程度上降低了模型異參同效性導(dǎo)致的不確定性[15]。

1 水文地質(zhì)試驗(yàn)

呼和浩特盆地第四系含水層系統(tǒng)劃分為山前單一結(jié)構(gòu)含水層系統(tǒng)和平原區(qū)雙層結(jié)構(gòu)含水系統(tǒng)，雙層結(jié)構(gòu)含水系統(tǒng)主要分布于盆地中部。本文水文地質(zhì)試驗(yàn)點(diǎn)位于呼和浩特盆地中東部的大黑河沖湖積平原，為雙層結(jié)構(gòu)的第四系含水系統(tǒng)[16]，地表高程1 081 m，地勢(shì)較為平坦。以中更新統(tǒng)上段（Q22）淤泥質(zhì)層為分層標(biāo)志，地下含水系統(tǒng)劃分為淺層潛水—微承壓含水層、深層承壓含水層[17]。上部為上更新統(tǒng)至全新統(tǒng)（Q3-4）潛水—微承壓含水層，厚度10～40 m，含水層分布廣、埋藏淺、顆粒粗、水量豐富；下部為中更新統(tǒng)下段和下更新統(tǒng)（Q12— Q1）承壓含水層，平均厚度大于50 m，下更新統(tǒng)（Q1）承壓含水層因分布范圍小，埋藏較深，水量小。中更新統(tǒng)下段（Q12）承壓含水層分布廣而穩(wěn)定，厚度大，含水豐富。大黑河平原區(qū)承壓含水層與上部潛水—微承壓含水層間的淤泥質(zhì)層厚度自東向西由12 m 增至超過(guò)100 m，上下兩含水層可通過(guò)淤泥質(zhì)弱透水層發(fā)生水量交換[18]。以往研究表明，承壓水長(zhǎng)期超采造成水位持續(xù)下降，進(jìn)而潛水對(duì)承壓含水層越流增大，是呼和浩特盆地山前潛水含水層疏干的主要原因[18-19]。

研究區(qū)承壓含水層與上部潛水—微承壓含水層之間的弱透水層（Q22）由多層淤泥層與砂層互層組成，發(fā)育總厚度約為40 m（圖1）。為準(zhǔn)確反演淤泥質(zhì)弱透水層的水文地質(zhì)參數(shù)，評(píng)價(jià)淤泥層滲透性能，同時(shí)為使抽水試驗(yàn)中，上下含水層較快產(chǎn)生響應(yīng)，本文選擇了中更新統(tǒng)下段（Q12）承壓含水層之上最近的，同時(shí)也是與承壓含水層水力聯(lián)系最緊密的淤泥層（埋深74.2～75.8 m）作為研究對(duì)象，設(shè)計(jì)實(shí)施了抽水試驗(yàn)和溶質(zhì)運(yùn)移試驗(yàn)。

在試驗(yàn)場(chǎng)施工2 眼相距2 m 的水文地質(zhì)鉆孔（圖1），其中抽水孔深82.5 m，在76～82 m 處安裝濾水管，抽水試驗(yàn)過(guò)程中對(duì)下部含水層（III 層）進(jìn)行抽水并觀測(cè)其水位；觀測(cè)孔深74.2 m，在72.5～74.2 m 處安裝濾水管，觀測(cè)該層（I 層）地下水位。兩含水層之間夾1.6 m厚的淤泥層（II 層）即為本次研究的弱透水層。

試驗(yàn)前在抽水孔和觀測(cè)孔中分別放置荷蘭vanEssen公司的水位自動(dòng)記錄儀（Mini-Diver），該水位儀可自動(dòng)監(jiān)測(cè)并記錄地下水的水位及溫度。將水位記錄儀放置在靜水位以下約10 m，保證抽水過(guò)程中水位計(jì)始終位于水面以下。抽水前的初始水位埋深為：抽水井50.06 m、觀測(cè)井50.02 m。整個(gè)抽水試驗(yàn)采用100 m3/d的定流量抽水，持續(xù)抽水40 d。開(kāi)始抽水后設(shè)置抽水井中水位記錄器記錄間隔0.5 h，觀測(cè)孔中記錄間隔1 h。

待上下含水層水位穩(wěn)定后，測(cè)得上下含水層水頭差為3.10 m。選用萘磺酸鈉作為示蹤劑進(jìn)行溶質(zhì)運(yùn)移試驗(yàn)。在開(kāi)始抽水后8 d，一次性向觀測(cè)孔中投入59 kg萘磺酸鈉粉末，溶質(zhì)投入后每天從抽水井取一次水樣，利用天津港東F-280 型熒光分光光度計(jì)對(duì)采集水樣進(jìn)行測(cè)試。累計(jì)測(cè)試了40 d 抽水井中水樣溶質(zhì)濃度。

2 地下水?dāng)?shù)值模型

2.1 概念模型

為確保抽水井抽水不影響邊界水位，進(jìn)而影響模型運(yùn)行結(jié)果，本次模擬以抽水試驗(yàn)和示蹤試驗(yàn)場(chǎng)地為中心，確定長(zhǎng)寬均為2 km 的矩形區(qū)域?yàn)槟M區(qū)。本次反演的是埋深74.2～75.8 m、厚度1.6 m 淤泥層的水文地質(zhì)參數(shù)。根據(jù)地層剖面（圖1）選擇埋深72.5～82.5 m的含水巖組作為模擬對(duì)象。垂向上將模型概化為3 層，上部中細(xì)砂層，透水性和富水性良好，概化為上部承壓含水層；中部為淤泥層，透水性差，概化為弱透水層；下部亦為中細(xì)砂層，透水性和富水性良好，概化為下部承壓含水層。

圖1 地層剖面示意圖Fig.1 Schematic stratigraphic section

試驗(yàn)場(chǎng)距自然邊界較遠(yuǎn)，且試驗(yàn)時(shí)間為冬季，無(wú)人為開(kāi)采影響，試驗(yàn)點(diǎn)附近水位觀測(cè)井顯示區(qū)域地下水位較穩(wěn)定，天然水力坡度約3‰，基于此將模型四周水平邊界條件概化為定水頭邊界。另外由于模擬期內(nèi)無(wú)顯著降水，蒸發(fā)對(duì)研究地層的影響可以忽略，故不考慮除試驗(yàn)中抽水井以外的源匯項(xiàng)條件。

模型涉及的水文地質(zhì)參數(shù)主要包括滲透系數(shù)、儲(chǔ)水率、彌散系數(shù)[17]等。研究區(qū)水文地質(zhì)條件簡(jiǎn)單，水平方向未作參數(shù)分區(qū)，每層水文地質(zhì)參數(shù)采用相同數(shù)值。

2.2 數(shù)學(xué)模型

根據(jù)上述水文地質(zhì)條件概化，將研究區(qū)地下水系統(tǒng)概化為三維非穩(wěn)定地下水流系統(tǒng)[18]，地下水流的控制方程為：

式中：h—地下水位/m；

K、Kz—含水介質(zhì)的水平和垂向滲透系數(shù)/（m·d-1）；

ε—源匯項(xiàng)/d-1；

Ss—儲(chǔ)水率/m-1。

溶質(zhì)運(yùn)移試驗(yàn)中選擇的萘磺酸鈉溶質(zhì)在地下水中不易發(fā)生化學(xué)反應(yīng)，不易揮發(fā)和吸附[20]，因此模擬過(guò)程僅考慮溶質(zhì)在地下水中的對(duì)流、彌散以及源匯項(xiàng)混合作用，溶質(zhì)運(yùn)移模型控制方程為：

式中：Ω—滲流區(qū)域；

C—溶質(zhì)濃度/（mg·L-1）；

vi—各方向上地下水實(shí)際流速/（m·d-1），通過(guò)地下水流方程求得；

Dij—水動(dòng)力彌散系數(shù)張量/（m2·d-1）；

CS—固體顆粒吸附的溶質(zhì)濃度/（mg·L-1）。

2.3 數(shù)值模型

本文采用MODFLOW-USG 程序包建立地下水流數(shù)值模型。該程序包是美國(guó)地質(zhì)調(diào)查局（USGS）2013年發(fā)布的新的MODFLOW 版本[21]，USG 是非結(jié)構(gòu)網(wǎng)格Unstructured Grid 的縮寫(xiě)（簡(jiǎn)稱(chēng)Ugrid）。在網(wǎng)格剖分方面非結(jié)構(gòu)化網(wǎng)格比傳統(tǒng)結(jié)構(gòu)化網(wǎng)格具有更強(qiáng)的靈活性[22]。同時(shí)基于控制體積有限差分（Control Volume Finite-difference，CVFD）的網(wǎng)格設(shè)計(jì)保證了網(wǎng)格幾何形狀靈活性，也提高了數(shù)值計(jì)算的穩(wěn)定性。利用非結(jié)構(gòu)化網(wǎng)格精細(xì)刻畫(huà)和靈活加密的特點(diǎn)，即可以實(shí)現(xiàn)對(duì)水文地質(zhì)試驗(yàn)點(diǎn)附近網(wǎng)格的精細(xì)刻畫(huà)，同時(shí)也因其在局部加密過(guò)程中產(chǎn)生的無(wú)效網(wǎng)格數(shù)量較少，可以大大降低模型計(jì)算負(fù)荷。

為實(shí)現(xiàn)試驗(yàn)場(chǎng)地的精細(xì)刻畫(huà)，在水平方向上將約4 km2模擬區(qū)剖分為64 m×64 m 的矩形網(wǎng)格，在抽水井和觀測(cè)井附近進(jìn)行非結(jié)構(gòu)化局部加密，最細(xì)網(wǎng)格剖分精度達(dá)到1 m×1 m（圖2）?？v向上將模型剖分為3 層，1 006.8～1 008.5 m（標(biāo)高，下同）為上部含水層，厚1.7 m；1 005.2～1 006.8 m 為中部弱透水層，厚1.6 m；998.5 ～1 005.2 m 為下部含水層，厚6.7 m；模型垂向總厚度10 m。

圖2 非結(jié)構(gòu)網(wǎng)格剖分圖Fig.2 Diagram showing unstructured grid subdivision

本次數(shù)值模型模擬驗(yàn)證期共45 d，劃分為91 個(gè)應(yīng)力期。模擬期前2 d 水位變化明顯，按每小時(shí)一個(gè)應(yīng)力期分為48 個(gè)應(yīng)力期。2 d 后抽水井與觀測(cè)井中水位均基本穩(wěn)定，按每天一個(gè)應(yīng)力期分為43 個(gè)應(yīng)力期。

根據(jù)概念模型中對(duì)模型的概化，設(shè)定模型初始條件與邊界條件。根據(jù)初始水力坡度，將模型四周邊界設(shè)為定水頭邊界；根據(jù)實(shí)測(cè)水位數(shù)據(jù)設(shè)定模型上下含水層的初始水位，弱透水層初始水位取上下含水層水位平均值[23]。在抽水井中設(shè)定100 m3/d 定流量抽水。

模型中上下含水層水文地質(zhì)參數(shù)依據(jù)經(jīng)驗(yàn)值給定[24]，弱透水層水文地質(zhì)參數(shù)通過(guò)觀測(cè)數(shù)據(jù)與模型運(yùn)行結(jié)果擬合進(jìn)行反演。

在水流模型基礎(chǔ)上采用MT3DMS 模塊建立溶質(zhì)運(yùn)移模型。在上部含水層中將溶質(zhì)輸入模型，并在下部含水層中設(shè)置溶質(zhì)濃度觀測(cè)井。溶質(zhì)運(yùn)移模型中各層的流動(dòng)特征參數(shù)依據(jù)經(jīng)驗(yàn)值給定[25-26]，縱向彌散度取10 m。

2.4 水文地質(zhì)參數(shù)反演

基于地下水位與溶質(zhì)濃度多源數(shù)據(jù)對(duì)弱透水層水文地質(zhì)參數(shù)進(jìn)行反演?？傮w反演思路可分為參數(shù)識(shí)別和參數(shù)驗(yàn)證2 大步驟（圖3）：以水流模型識(shí)別過(guò)程作為反演參數(shù)的識(shí)別過(guò)程；溶質(zhì)運(yùn)移模型的識(shí)別過(guò)程作為反演參數(shù)的驗(yàn)證過(guò)程。首先對(duì)水流模型參數(shù)進(jìn)行識(shí)別，通過(guò)水位擬合反演出一組滿(mǎn)足水流模型的弱透水層水文地質(zhì)參數(shù)；然后將這些參數(shù)代入溶質(zhì)運(yùn)移模型進(jìn)行參數(shù)驗(yàn)證，通過(guò)溶質(zhì)運(yùn)移模型的識(shí)別驗(yàn)證反演出的水文地質(zhì)參數(shù)。若參數(shù)驗(yàn)證失?。ㄈ苜|(zhì)濃度曲線(xiàn)無(wú)法擬合）便調(diào)整水流模型中含水層參數(shù)重新反演，反復(fù)迭代后得到同時(shí)滿(mǎn)足水流模型和溶質(zhì)運(yùn)移模型的弱透水層參數(shù)。基于水流模型和溶質(zhì)運(yùn)移模型的共同識(shí)別，可反演出可信度較高的弱透水層水文地質(zhì)參數(shù)。

圖3 參數(shù)反演流程圖Fig.3 Flow chart showing parameter inversion

3 結(jié)果分析

水文地質(zhì)試驗(yàn)中觀測(cè)到：抽水井在開(kāi)始抽水后13 h 基本達(dá)到穩(wěn)定，穩(wěn)定降深3.76 m；觀測(cè)井中水位下降存在滯后，在開(kāi)始抽水后19 h 觀測(cè)孔水位開(kāi)始下降，36 h 觀測(cè)孔水位基本穩(wěn)定，穩(wěn)定降深0.66 m。由于抽水泵的影響，抽水井中后期水位觀測(cè)值出現(xiàn)小幅波動(dòng)，與總體水位變化相比波動(dòng)不大，不會(huì)對(duì)穩(wěn)定水位造成影響。溶質(zhì)運(yùn)移試驗(yàn)顯示，在開(kāi)始抽水24 d后，采集水樣中觀測(cè)到了萘磺酸鈉濃度峰值（3.38×10-8mol/L）。

通過(guò)地下水流模型與溶質(zhì)運(yùn)移模型的共同識(shí)別對(duì)弱透水層水文地質(zhì)參數(shù)進(jìn)行反演，最終在完成參數(shù)識(shí)別和參數(shù)驗(yàn)證后，模型中的地下水位和溶質(zhì)濃度與實(shí)際觀測(cè)值擬合效果見(jiàn)圖4。

圖4 含水層水位和抽水井中溶質(zhì)濃度擬合效果Fig.4 Fitting effect of groundwater level in the aquifer and solute concentration in the pumping well

從水位及溶質(zhì)濃度擬合效果中可以看出，抽水井中的水位和溶質(zhì)濃度與實(shí)際觀測(cè)數(shù)據(jù)擬合較好，觀測(cè)孔中模擬水位變化明顯提前于實(shí)際水位變化，這是由于地下水流模型中未考慮到弱透水層中存在的滯后釋水現(xiàn)象[14]。模型中下部含水層水頭變化會(huì)直接導(dǎo)致弱透水層水頭下降，從而引發(fā)越流；而實(shí)際情況是下部含水層水頭下降后，會(huì)先導(dǎo)致弱透水層水頭由下自上逐漸消散，待水頭消散影響到弱透水層的上部邊界時(shí)才會(huì)引起上部含水層向下越流。實(shí)測(cè)水位波動(dòng)比較大是由于區(qū)域地下水位的變化導(dǎo)致的。

雖然模型中未考慮到弱透水層水頭消散的滯后性，但抽水試驗(yàn)前期水位下降滯后并不會(huì)影響到穩(wěn)定水位及水位穩(wěn)定后進(jìn)行的溶質(zhì)運(yùn)移試驗(yàn)。忽略弱透水層水位響應(yīng)滯后性，抽水井與觀測(cè)井的后期水位變化曲線(xiàn)均與實(shí)際觀測(cè)數(shù)據(jù)有較好的擬合效果，溶質(zhì)運(yùn)移模型運(yùn)行結(jié)果中溶質(zhì)濃度峰值出現(xiàn)時(shí)間也與實(shí)際觀測(cè)基本一致，只是濃度峰值大小存在一定差異，這可能與水樣采集時(shí)刻及測(cè)量誤差有關(guān)，也可能因?yàn)槿苜|(zhì)運(yùn)移過(guò)程中發(fā)生的對(duì)流彌散外的其他作用。

總體而言，模型運(yùn)行結(jié)果與觀測(cè)數(shù)據(jù)的擬合效果說(shuō)明模型能較好地反映實(shí)際水文地質(zhì)條件。模型識(shí)別后的主要參數(shù)見(jiàn)表1。整體來(lái)看，各參數(shù)取值滿(mǎn)足經(jīng)驗(yàn)范圍，抽水含水層水文地質(zhì)參數(shù)取值與解析法計(jì)算出的結(jié)果相近（假設(shè)不考慮越流條件下，根據(jù)承壓井Dupuit 公式計(jì)算出抽水含水層滲透系數(shù)K=4.17 m/d）。最終反演得出的淤泥層的水文地質(zhì)參數(shù)為：垂向滲透系數(shù)Kz=1.2×10-4m/d；儲(chǔ)水率Ss=1.0×10-5m-1。反演結(jié)果與該地區(qū)以往研究成果[17]具有較好的一致性，可信度較高。

表1 主要水文地質(zhì)參數(shù)表Table 1 Main hydrogeological parameters

4 結(jié)論

（1）弱透水層滯后釋水會(huì)導(dǎo)致越流系統(tǒng)水位響應(yīng)滯后。抽水試驗(yàn)中，抽水井開(kāi)始抽水19 h 后觀測(cè)井中才觀測(cè)到水位下降。

（2）基于地下水位與溶質(zhì)濃度多源數(shù)據(jù)對(duì)弱透水層水文地質(zhì)參數(shù)進(jìn)行反演，得出淤泥質(zhì)弱透水層水文地質(zhì)參數(shù)：垂向滲透系數(shù)Kz=1.2×10-4m/d；儲(chǔ)水率Ss=1.0×10-5m-1。

（3）模型中未對(duì)弱透水層滯后釋水過(guò)程進(jìn)行詳細(xì)刻畫(huà)，導(dǎo)致觀測(cè)孔中計(jì)算水位與實(shí)測(cè)水位存在一定誤差。溶質(zhì)運(yùn)移模型中參數(shù)取值誤差也可能會(huì)對(duì)模型精度產(chǎn)生一定影響。針對(duì)弱透水層中滯后釋水、溶質(zhì)運(yùn)移等復(fù)雜水文地質(zhì)現(xiàn)象，今后還需進(jìn)行深入研究。