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

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

（1.中國地質大學（北京）水資源與環境學院，北京 100083；2.中國地質科學院水文地質環境地質研究所，河北 石家莊 050061）

水文地質參數對于區域地下水資源評價、數值模擬及預報、開發利用與保護及科學管理具有重要意義[1]。當實際含水介質的水文地質參數不易獲取時，通過已知水文地質結構、參數、源匯條件及觀測資料反演未知參數通常是一種十分有效的方法[1-3]。當前水文地質參數的反演方法主要包括解析法和數值法。隨著計算機技術的發展，地下水數值模擬技術逐漸成熟，數值法逐漸成為解決地下水相關問題的重要方法[4-5]，用于水文地質參數求取、地下水資源評價、預報等，為地下水開發利用及科學管理提供有效的工具和重要依據。

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

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

1 水文地質試驗

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

研究區承壓含水層與上部潛水—微承壓含水層之間的弱透水層（Q22）由多層淤泥層與砂層互層組成，發育總厚度約為40 m（圖1）。為準確反演淤泥質弱透水層的水文地質參數，評價淤泥層滲透性能，同時為使抽水試驗中，上下含水層較快產生響應，本文選擇了中更新統下段（Q12）承壓含水層之上最近的，同時也是與承壓含水層水力聯系最緊密的淤泥層（埋深74.2～75.8 m）作為研究對象，設計實施了抽水試驗和溶質運移試驗。

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

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

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

2 地下水數值模型

2.1 概念模型

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

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

試驗場距自然邊界較遠，且試驗時間為冬季，無人為開采影響，試驗點附近水位觀測井顯示區域地下水位較穩定，天然水力坡度約3‰，基于此將模型四周水平邊界條件概化為定水頭邊界。另外由于模擬期內無顯著降水，蒸發對研究地層的影響可以忽略，故不考慮除試驗中抽水井以外的源匯項條件。

模型涉及的水文地質參數主要包括滲透系數、儲水率、彌散系數[17]等。研究區水文地質條件簡單，水平方向未作參數分區，每層水文地質參數采用相同數值。

2.2 數學模型

根據上述水文地質條件概化，將研究區地下水系統概化為三維非穩定地下水流系統[18]，地下水流的控制方程為：

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

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

ε—源匯項/d-1；

Ss—儲水率/m-1。

溶質運移試驗中選擇的萘磺酸鈉溶質在地下水中不易發生化學反應，不易揮發和吸附[20]，因此模擬過程僅考慮溶質在地下水中的對流、彌散以及源匯項混合作用，溶質運移模型控制方程為：

式中：Ω—滲流區域；

C—溶質濃度/（mg·L-1）；

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

Dij—水動力彌散系數張量/（m2·d-1）；

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

2.3 數值模型

本文采用MODFLOW-USG 程序包建立地下水流數值模型。該程序包是美國地質調查局（USGS）2013年發布的新的MODFLOW 版本[21]，USG 是非結構網格Unstructured Grid 的縮寫（簡稱Ugrid）。在網格剖分方面非結構化網格比傳統結構化網格具有更強的靈活性[22]。同時基于控制體積有限差分（Control Volume Finite-difference，CVFD）的網格設計保證了網格幾何形狀靈活性，也提高了數值計算的穩定性。利用非結構化網格精細刻畫和靈活加密的特點，即可以實現對水文地質試驗點附近網格的精細刻畫，同時也因其在局部加密過程中產生的無效網格數量較少，可以大大降低模型計算負荷。

為實現試驗場地的精細刻畫，在水平方向上將約4 km2模擬區剖分為64 m×64 m 的矩形網格，在抽水井和觀測井附近進行非結構化局部加密，最細網格剖分精度達到1 m×1 m（圖2）。縱向上將模型剖分為3 層，1 006.8～1 008.5 m（標高，下同）為上部含水層，厚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 非結構網格剖分圖Fig.2 Diagram showing unstructured grid subdivision

本次數值模型模擬驗證期共45 d，劃分為91 個應力期。模擬期前2 d 水位變化明顯，按每小時一個應力期分為48 個應力期。2 d 后抽水井與觀測井中水位均基本穩定，按每天一個應力期分為43 個應力期。

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

模型中上下含水層水文地質參數依據經驗值給定[24]，弱透水層水文地質參數通過觀測數據與模型運行結果擬合進行反演。

在水流模型基礎上采用MT3DMS 模塊建立溶質運移模型。在上部含水層中將溶質輸入模型，并在下部含水層中設置溶質濃度觀測井。溶質運移模型中各層的流動特征參數依據經驗值給定[25-26]，縱向彌散度取10 m。

2.4 水文地質參數反演

基于地下水位與溶質濃度多源數據對弱透水層水文地質參數進行反演。總體反演思路可分為參數識別和參數驗證2 大步驟（圖3）：以水流模型識別過程作為反演參數的識別過程；溶質運移模型的識別過程作為反演參數的驗證過程。首先對水流模型參數進行識別，通過水位擬合反演出一組滿足水流模型的弱透水層水文地質參數；然后將這些參數代入溶質運移模型進行參數驗證，通過溶質運移模型的識別驗證反演出的水文地質參數。若參數驗證失敗（溶質濃度曲線無法擬合）便調整水流模型中含水層參數重新反演，反復迭代后得到同時滿足水流模型和溶質運移模型的弱透水層參數。基于水流模型和溶質運移模型的共同識別，可反演出可信度較高的弱透水層水文地質參數。

圖3 參數反演流程圖Fig.3 Flow chart showing parameter inversion

3 結果分析

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

通過地下水流模型與溶質運移模型的共同識別對弱透水層水文地質參數進行反演，最終在完成參數識別和參數驗證后，模型中的地下水位和溶質濃度與實際觀測值擬合效果見圖4。

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

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

雖然模型中未考慮到弱透水層水頭消散的滯后性，但抽水試驗前期水位下降滯后并不會影響到穩定水位及水位穩定后進行的溶質運移試驗。忽略弱透水層水位響應滯后性，抽水井與觀測井的后期水位變化曲線均與實際觀測數據有較好的擬合效果，溶質運移模型運行結果中溶質濃度峰值出現時間也與實際觀測基本一致，只是濃度峰值大小存在一定差異，這可能與水樣采集時刻及測量誤差有關，也可能因為溶質運移過程中發生的對流彌散外的其他作用。

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

表1 主要水文地質參數表Table 1 Main hydrogeological parameters

4 結論

（1）弱透水層滯后釋水會導致越流系統水位響應滯后。抽水試驗中，抽水井開始抽水19 h 后觀測井中才觀測到水位下降。

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

（3）模型中未對弱透水層滯后釋水過程進行詳細刻畫，導致觀測孔中計算水位與實測水位存在一定誤差。溶質運移模型中參數取值誤差也可能會對模型精度產生一定影響。針對弱透水層中滯后釋水、溶質運移等復雜水文地質現象，今后還需進行深入研究。