曹文庚，楊會峰，高媛媛，南 天，王 哲，徐素娟

（1.中國地質科學院 水文地質環(huán)境地質研究所，河北 石家莊 050061；2.自然資源部 地下水科學與工程科普基地，河北 石家莊 050061；3.水利部南水北調規(guī)劃設計管理局，北京 100038；4.海河水利委員會水文局，天津 300170；5.河北省地礦局地質勘查技術中心，河北 石家莊 050081）

1 研究背景

南水北調中線受水區(qū)保定平原近40年來，由于長期、持續(xù)、大規(guī)模過度開采地下水，區(qū)域地下水水位持續(xù)下降，超采嚴重[1]，部分含水層被疏干或枯竭，引發(fā)了地面沉降、地裂縫、水質惡化等一系列生態(tài)與環(huán)境地質問題，嚴重制約當地經濟社會可持續(xù)發(fā)展和生態(tài)文明建設，危及生態(tài)安全和供水安全。南水北調東中線一期工程經過10 多年的建設分別于2013年、2014年順利通水。隨著工程的通水，2018年南水北調工程開始對保定市的拒馬河、白洋淀進行生態(tài)補水，地下水超采狀況有所緩解，保定平原受水區(qū)供水格局已發(fā)生改變，在緩解水資源緊張局面的同時，水位恢復不當也會引起一系列的環(huán)境地質問題[2]，特別是由于水位上升，在水巖交換作用下，會引發(fā)地下水化學組份的變化，局部地區(qū)地下水存在潛在污染的可能。因此本文基于校正的地下水流模型，利用美國地調局開發(fā)的PHREEQC 軟件，模擬2030年淺層地下水位上升后，水巖交換作用下保定平原地下水主要化學組分變化趨勢，分析地下水質量演變機制。

目前，基于PHREEQC 的水巖交換模擬主要采用水文地球化學系統(tǒng)的化學平衡模擬方法[3]，國內外學者利用PHREEQC 開發(fā)出多種模型和方法進行水巖交換反應模擬，Claus 利用PHREEQC 對不同深度的地下水化學組份進行了建模，定義了一種一維的反應輸運模型，該模型考慮了地下水位上升過程中部分氧化的含硫黃鐵礦沉積物與地下水中所發(fā)生的水文地球化學反應模擬[4]；王廣才團隊建立了耦合的Monte Carlo-WATEQ4F模型對水-巖反應平衡狀態(tài)進行模擬和分析[5]；Stigter 團隊建立正向地球化學模型方法，模擬了地表水與地下水混合灌溉對地下水環(huán)境的影響機制[6]；湯潔團隊[7]及趙江濤等[8]利用反向水文地球化學模擬方法解釋了地下水化學成分的形成、分布及演化機理。通過梳理總結，目前國內外在水巖交換方面較為成熟和流行的模型包括：（1）表面吸附模型（SCMs），主要基于水相及固相的化學特性，用數學的方法刻畫水巖交換作用中的離子吸附作用、競爭及協同效應，尤其擅長模擬礦物表面化學鍵，例如絡合反應[9-10]；（2）擴散雙層模型（DLL），是顯示的計算擴散層組成和非靜電層的表面絡合的模型，可以模擬水溶液中選擇性地吸附某種離子[11]。（3）電位分布-多點位絡合模型（CD-MUSIC），可以描述鐵氧化物及氫氧化物表面的吸附行為和離子形態(tài)分布[12-13]，還能預測環(huán)境條件改變下的離子形態(tài)的響應特征[14]，同時可以捕捉主要離子的非線性電荷效應[15]。因此，CD-MUSIC模型非常適合本次模擬的主題，故選擇CD-MUSIC模型模擬地下水位恢復后的地下水化學組分的演化趨勢。

2 研究區(qū)概況

2.1 自然地理概況南水北調受水區(qū)保定段主體為保定市平原區(qū)，北臨京廊、西依太行、東臨滄衡、南鄰石家莊，行政區(qū)劃上包括保定市區(qū)、滿城區(qū)、清苑區(qū)、徐水區(qū)、雄安新區(qū)、涿州市、定州市、高陽縣、蠡縣、高碑店、定興、望都和安國，平原區(qū)面積約1.1 萬km2（圖1）。該地區(qū)地勢平坦，海拔多在50 m 以下，自山前向白洋淀微傾斜。地表水系屬大清河流域，流經的主要河流有拒馬河、易水河、漕河、界河和唐河等。

研究區(qū)地貌單元主要包括太行山山前沖洪積傾斜平原、沖積湖積微傾斜平原及平原東部的沖湖積低平原。山前沖洪積傾斜平原分布于易縣、滿城、唐縣、望都、定州、曲陽和行唐一帶，由拒馬河、易水河、漕河、府河、唐河等河流沖洪積扇構成，地面高程80～30 m，地形坡降4‰～1‰，地勢大致自西南向東北傾斜。沖積湖積微傾斜平原分布于定興西南、容城以北、徐水、清苑、博野和安國一帶，位于山前沖洪積傾斜平原東部。上部為近代河流沖積層或扇前洼地相堆積物，下伏沖洪積層。地勢相對平坦，地面高程30～10 m，地形坡降2‰～1‰。沖湖積低平原分布于以白洋淀為中心的周邊地帶，主要為近代河流沖積相和湖沼相共同組成的平原，地勢低平、地面高程10～7 m，地形坡降一般0.7‰～0.2‰。

2.2 水資源狀況及環(huán)境地質問題據《保定市第二次水資源評價報告》顯示，保定平原多年平均地表水資源量為1.49億m3/a，地下水資源量為15.79億m3/a，淺層地下水可開采資源量為15.39億m3/a。依據2016年保定市水資源公報數據，2016年保定平原總供水量22.98億m3，其中地下水開采量達到20.34億m3，包括淺層地下水19.37億m3，深層地下水0.97億m3。保定水資源整體狀況匱乏，自1980年代開始大規(guī)模開采地下水，多年的地下水超采造成地下水水位持續(xù)下降，形成多個淺層地下水位降落漏斗，依據2016年6月淺層地下水位等值線（圖1（b）），形成了“保定徐水區(qū)-容城淺層地下水漏斗”、“雄縣地下水漏斗”以及“高陽-蠡縣-清苑淺層地下水漏斗”，最大水位埋深為49.02 m，漏斗總面積達到4144.4 km3。

2.3 南水北調保定平原區(qū)外調水量包括南水北調和引黃入冀工程，南水北調計劃分配給保定水量為5.0015億m3，引黃入冀工程計劃分配保定水量2.55億m3，白洋淀生態(tài)補水1.1億m3，2017年開始試通水。2018年9月水利部選擇河北省境內的滹沱河、滏陽河、南拒馬河3條典型河流的重點河段進行生態(tài)補水，2020年計劃將生態(tài)補水范圍擴大到14條河流。南水北調持續(xù)來水，保定平原全區(qū)皆為南水北調受水區(qū)，一定程度上替代保定平原全區(qū)現有的地下水開采，大幅度壓減了地下水開采量，從而引起區(qū)域地下水流場的變化，地下水位上升明顯。

圖1 保定平原水文地質鉆探取樣點及淺層地下水流場

3 研究方法

3.1 地下水流場預測獲取2030年淺層地下水位基于已校正的地下水流模型[16]，依據國務院批復的《海河流域綜合規(guī)劃（2012—2030）》，南水北調河北省2020年分配水量30.4億m3，保定市年均分配江水量5.05億m3；2030年南水北調河北省規(guī)劃分配水量分別為42.3億m3，同時考慮雄安新區(qū)建設，保定2030年南水北調分配水量預計達到8.0億m3。在預測模型中，模擬預測期為2016.6—2030.12，共175個應力期，以月為應力期。源匯項部分，白洋淀入滲、河流、降水等采用全部平水年的數據；邊界流量采用現狀模型的定流量處理。依據預測結果（圖1（c）（d）），2030年保定平原全區(qū)淺層地下水位呈現不同程度的上升，上升幅度最大的為原地下水漏斗區(qū)域的雄縣、高陽及蠡縣，水位最大上升幅度達到35 m，山前地區(qū)水位上升幅度相對較小，一般小于15 m（圖1（d））。

3.2 地下水化學組分模擬點位及鉆探巖心采集深度選取在全區(qū)選取32個水文地質鉆探點，點位的分布基本控制全區(qū)所有地質地貌單元（圖1（a）），進行地下水樣品及鉆孔巖心樣品采集，巖心的采集深度范圍與2016 至2030年淺層地下水位上升高度范圍一致，見圖1（d），樣品采集間隔為每10 cm 采集巖心樣品50 mg，然后進行混合、過篩晾曬后送往河北省地質實驗測試中心進行土樣的易溶鹽測試，測試項目包括：Na、Ca、Mg、HCO3、Cl、SO4含量。由于每個地下水樣品的地下水環(huán)境不同，鉆孔易溶鹽含量也存在較大的差異，因此本次模擬針對每個地下水樣品點都建立獨立的水文地球化學模型，分別進行32組水巖相互作用的地球化學模擬。

3.3 水文地球化學模擬方法本次模型基于美國地調局開發(fā)的PHREEQC 軟件，采用的水文地球化學模擬方法是建立在平衡常數基礎上的平衡常數法。通常，水文地球化學模擬方法分為正向地球化學模擬和反向地球化學模擬，前者是依據假設反應及動力學參數來預測地下水化學組分及水-巖間的質量交換；后者是依據水、土樣品的化學組分數據來反向推測水巖交換反應的過程。

（1）本次研究已經掌握的數據包括：①2016年淺層地下水化學數據；②2016年水樣點附近巖心土樣的化學元素含量數據。

（2）模擬預測的目標為：①2030年淺層地下水化學組分；②2030年淺層地下水化學類型。

本次模擬預測工作采用正向水文地球化學模擬，然而，正向模擬需要提供室內試驗條件下的反應動力學方程以及相關模型參數。但是，由于本次工作缺乏室內試驗，無法通過試驗獲取這些參數，因此在正向模擬之前，還需要加入反向地球化學模擬工作。依據2016年已知的地下水化學組分和巖心樣品化學組分，通過反向模擬，獲取了水文地球化學模擬的相關參數。然后，將參數重新代入正向地球化學模擬過程中，利用CD-MUSIC模型提供的反應方程，采用水文地球化學混合作用模型、表面絡合作用模型、離子交換模型，進行正向地球化學模擬，預測得出2030年淺層地下水恢復條件下水化學組分（圖2）。

圖2 南水北調受水區(qū)保定段水位恢復地下水化學組分模擬技術路線

正向地球化學模擬主要基于以下4個方程式。

利用PHREEQC 進行地下水組分存在形式的計算，各種表面絡合組分可有以下方程確定：

式中：Cj為第j中絡合組分的活度；Kj為第j種絡合組分的熱力學平衡常數；αk為第k種組分的活度；Pk.j為第k種組分相對第j種絡合組分的化學計量數。上式即質量作用方程，對于地下水中的主要組分，若第k種主要組分的化學分析濃度為Ak，則質量守恒方程為：

從物理意義上講，ak和Cj必須滿足如下條件：

正向模擬反應的基本步驟為：起點水溶液的化學組分+反應相→終點水溶液的化學組分+生成相。可用下式表示：

式中：p為礦物相（反應相+生成相）數目；ap為第p種礦物相進入或離開每升溶液的摩爾數；bp，k為第p種礦物相中第k種元素的化學計量數；mT.k為溶液中第k種元素的總摩爾濃度。

將32個水位點的化學組分數據以及相關熱力學參數代入以上方程式進行迭代，分別建立32個水化學模擬模型，得出最終的評價結果。

3.4 模型參數選擇基于本項目已掌握的數據、模擬目標及方法思路，結合保定平原淺層含水層的礦物環(huán)境特征及地下水的水化學特征，確定模型的三種參數，具體內容及選擇理由為：①鐵錳氧化物及氫氧化物的吸附參數，依據本區(qū)的1∶5萬區(qū)域地質調查成果，此區(qū)域淺層地下水所賦存的晚更新世地層中，重礦物主要為針鐵礦、磁鐵礦、鈦鐵礦等，因此本次模型考慮針鐵礦中鐵錳氧化物的吸附性能，因此反應方程采用針鐵礦和各離子吸附解析的平衡方程（見表1）。表1中，lgK為收集已有的化學反應的熱力學常數；lgKopti為進行參數優(yōu)化后的熱力學常數；②表面吸附節(jié)點的參數，主要考慮表面吸附的物理吸附和微觀結構特征，本次依據CD-MUSIC模型定義了針鐵礦的兩種吸附點位類型，第一種是只有一個鐵原子的吸附點位，即Geo_uni，第二種為有三個鐵原子的吸附點位，即Geo_tri，基于固相的表觀形貌特征來反映吸附能力；③方解石、白云石、石膏等飽和指數。依據保定市1∶20萬水文地質圖及其說明書等資料，晚更新世是太行山保定段山前沖洪積扇發(fā)育最為強烈的時期，由于山區(qū)廣泛分布的為薊縣系霧迷山組和鐵嶺組的白云巖、白云質灰?guī)r等，受山區(qū)匯水區(qū)溶濾作用影響，平原區(qū)地下水陽離子主要為Ca、Mg 等離子，陰離子以HCO3及SO4為主，因此主要考慮方解石（CaCO3）及白云石（CaMg（CO3）2）以及石膏（CaSO4·2H2O）礦物的溶解及沉淀過程。

3.5 水文地球化學反應方程確定依據模型參數的選擇，本次模擬的反應方程選取具體見表1。

3.6 模擬結果檢驗本次通過PHREEQC 模擬預測獲取的2030年保定平原淺層地下水化學數據，需要對結果進行可靠性檢查，從而保證水文地球化學模擬的質量。

本次采用陰陽離子平衡檢查，天然地下水的化學成分中，陰陽離子是主要的組成部分，因此在水化學分析前進行陰陽離子的平衡檢查是非常必要的。將所有陰陽離子的單位由原來的“mg/L”換算為毫克當量濃度（meq/L），通過計算水中陰陽離子的相對誤差來判斷預測的2030年水化學數據的可靠性。

離子平衡的檢查公式為：

式中：E為相對誤差，%；Mc及Ma為陽離子及陰離子的毫克當量數，meq/L。

4 結果與分析

4.1 反向水文地球化學模擬結果保定平原從1970年代開始大規(guī)模開采地下水，地下水位不斷下降，本次選擇某一個地下水樣品點，例如以SH0538C 點為例，此點位于保定市博野縣內，2016年水位標高為-9.45 m，預測2030年水位標高恢復至20 m，水位恢復接近30 m。本文通過鉆探采集SH0538C 水樣點附近巖心，選擇水位恢復段巖心樣品后，進行過篩、混合，對混合樣品測試易溶鹽組分，獲取水位恢復段的巖土易溶鹽組分特征，由于易溶鹽是采用水溶法，即用4倍的純水進行溶解，并測試溶解出來主要離子含量。所以，易溶鹽組分特征實際上是巖土的水溶態(tài)。

表1 主要離子表面絡合反應

基于以上條件，可以得到以下假設：水位恢復段的巖土水溶態(tài)含量基本等于自20世紀70年代水位下降過程中，表面絡合到巖土顆粒上的吸附量?；诖思僭O，嘗試建立反向地球化學模型，以巖土易溶鹽含量作為模擬的終點，反向模擬強開采條件下水位下降過程中水巖交換反應，對研究方法中3個參數進行調整，包括：表面吸附點位參數、鐵氧化物和氫氧化物的吸附參數以及方解石、石膏的飽和指數。

假設驗證：通過調整參數，對Ca、Mg、Na、HCO3、SO4、Cl 共6種離子進行擬合，擬合度結果見表2。從表2可以看出，整體擬合度較好，Ca、Na、HCO3、SO4擬合度最優(yōu)（圖3及圖4），Cl 擬合度最差，其原因是Cl 的化學性質較為穩(wěn)定，在水巖交換過程中不發(fā)生任何交換反應，溶解度較高，因此造成Cl含量實測值與預測值擬合度較低，Cl 的擬合度基本可以滿足本次模擬的精度。

表2 水位恢復段巖土易溶鹽含量指標與反向水文地球化學模擬計算的表面吸附量擬合度評價指標

圖3 易溶鹽指標含量實測值與反向水文地球化學模擬計算的表面吸附量對比

4.2 正向水文地球化學預測結果基于反向地球化學模擬得出的表面吸附節(jié)點參數、鐵錳氧化物吸附參數及方解石石膏的飽和指數，代入正向地球化學模型中，進行正向地球化學預測工作，基于方程（1）—（3），導入反向地球化學模型模擬得出表面絡合參數、鐵錳吸附參數以及礦物的飽和指數，在給定參數后就可計算出地下水中各種組分存在形式的濃度（活度）。建立水位恢復上升的一維流水-巖作用模型，基于PHREEQC 本身自帶表面吸附模型和離子交換模型，對2016年每個地下水樣品點進行正向模擬，分別輸入各點水位上升前水化學數據，分別確定各自表面吸附過程、離子交換過程及礦物的溶解沉淀過程。最終獲取保定平原全區(qū)2030年淺層地下水化學組分（表3）。

采用研究方法中3.6 的模擬檢驗公式，對預測出的2030年淺層地下水化學組分進行檢驗，32組水化學數據電荷平衡誤差皆小于10%，本次預測結果較為可靠。

5 討論

依據模擬預測結果，2030年南水北調受水區(qū)保定平原域淺層地下水位恢復明顯，區(qū)內地下水化學組分也隨之發(fā)生變化。為進一步分析南水北調供水對保定平原淺層地下水環(huán)境的影響，對保定平原2016年及2030年的地下水化學規(guī)律進行對比研究，分析水化學類型及影響地下水質量的主要指標的變化趨勢。

5.1 水化學類型變化趨勢依據2016年及2030年保定平原淺層地下水化學數據，分別制作了2016年淺層地下水化學類型分區(qū)圖以及2030年淺層地下水化學類型分區(qū)圖，對比二者可以看出：至2030年，南水北調供水條件下，由于水位大面積抬升，全區(qū)淺層地下水化學類型發(fā)生了較大的變化。

表3 2030年保定平原淺層地下水化學組分預測 （單位：mg/L）

總體上，淺層地下水逐漸向鹽化方向發(fā)展，HCO3-SO4型水代替HCO3型水成為保定平原主要的陰離子水化學類型，平原東部的地下水排泄區(qū)出現了Cl型水；Na-Ca型水成為平原中部主要的陽離子水化學類型（圖4）。

分帶特征變化明顯，由山前至平原中部，2016年淺層地下水化學類型由山前的Ca-Mg-HCO3型水，至平原中部逐漸演化為Na-Ca-HCO3-Cl型水；2030年由山前至平原中部，山前沖洪積平原的Ca-Mg-HCO3型水面積大幅下降，由2016年6749 km2下降至3351 km2，平原中部Na-Mg-HCO3-SO4型水面積由2016年的921 km2擴大至4109 km2（表4及圖4）。大部分地區(qū)，Na和Cl 是地下水化學組分中增長比例最大的離子，造成此問題的原因是由于水位恢復段巖心存在的大量Na+及Cl-，巖鹽的溶解是造程平原東部大部分地區(qū)Cl元素含量上升的主要原因，高蠡清淺層地下水漏斗是Cl含量增幅最大的區(qū)域，經過預測，至2030年，在平原東部地下水漏斗中心，行政區(qū)位置屬于高陽縣東北部，形成了Na-Cl型水，面積達到121.99 km2（表4）。

圖4 2016年及2030年水位恢復條件下保定平原淺層地下水水化學類型分區(qū)對比

表4 保定平原2016年及2030年淺層地下水水化學類型面積對比

5.2 TDS變化趨勢分析前人多項成果已經表明，地下水超采以及地下水的快速下降，會造成地下水TDS含量的增加。依據本次模型，水位回升條件控制下，2030年全區(qū)淺層地下水TDS含量變化表現出明顯的水文地質分帶性。依據2016年及2030年淺層地下水TDS含量制作的TDS含量變差圖（圖5）可以看出，自山前向平原東部，TDS含量由增加逐漸變?yōu)闇p少，保定平原北部永定河沖洪積平原、平原中部的易水河沖洪積平原、平原南部的唐河-沙河沖洪積扇，都表現為TDS 上升的趨勢，Cl及SO4含量上升是造成TDS 增大的主要原因；保定市區(qū)、安新、雄縣地區(qū)TDS含量表現為下降趨勢，其分布范圍與白洋淀附近的高水位區(qū)基本一致，TDS含量下降最大處分別位于保定市區(qū)及白洋淀周邊地區(qū)。

5.3 硬度變化趨勢分析硬度是研究地下水水質的重要指標之一，也是區(qū)域上影響地下水質量分級的主控因素。雖然目前已有大量的文獻對硬度與地下水位、地下水開采的關系進行了分析論證，通常認為水位下降、地下水超采是造成地下水硬度升高的主要原因。隨著南水北調供水、地下水開采減少，水位逐漸恢復，包氣帶厚度減小，部分包氣帶轉化為飽和帶，地下水環(huán)境由氧化環(huán)境逐漸變?yōu)檫€原環(huán)境，地下水入滲途徑縮短，大大減少了白云石等礦物溶解帶來的Ca2+及Mg2+含量。依據2016年至2030年淺層地下水硬度變差圖（圖5）可以看出，全區(qū)大部分地區(qū)，淺層地下水硬度下降明顯。至2030年，硬度下降幅度最大的區(qū)域主要為雄安新區(qū)及其周邊地區(qū)，硬度下降幅度超過200 mg/L，雄安新區(qū)及其附近水質有逐漸好轉的趨勢；而清苑南部以及保定市北部徐水容城一帶，淺層地下水硬度有小幅度上升。

圖5 2016—2030年保定平原淺層地下水TDS及硬度含量變幅

6 結論

針南水北調受水區(qū)保定平原開展地下水位回升條件下的水質演化預測研究，基于水巖相互作用過程中的表面吸附作用、鐵錳氧化物及氫氧化物的吸附作用、離子的競爭與協同吸附效應，以及不同礦物的溶解平衡，分析得出：

（1）總體上，南水北調受水區(qū)保定平原淺層地下水逐漸向鹽化方向發(fā)展，HCO3-SO4型水代替HCO3型水成為保定平原主要的陰離子水化學類型，平原東部高陽-蠡縣地區(qū)出現了Cl型水，Na-Ca型水成為保定平原中東部地區(qū)主要的陽離子水化學類型，巖鹽的溶解是造成平原東部大部分地區(qū)Cl元素含量上升的主要原因，高陽-蠡縣淺層地下水漏斗區(qū)是Cl含量增幅最大的區(qū)域。

（2）至2030年淺層地下水位整體回升后，保定平原的東西部地下水質量演化趨勢存在較大的差異。西部的山前平原淺層地下水TDS含量增大明顯，其中拒馬河沖洪積扇、易水河沖洪積扇以及唐河-沙河沖洪積扇增幅達到100 mg/L 以上，水位上升后地層中巖鹽及石膏的礦物溶解是造成TDS 上升的主要原因。中東部平原的保定市區(qū)、安新、雄縣、高陽等地區(qū)TDS含量表現為下降趨勢，其分布范圍與淺層地下水降落漏斗基本一致；至2030年，全區(qū)淺層地下水硬度下降明顯，硬度降幅最大的區(qū)域主要為地下水漏斗中心的雄縣、高陽及蠡縣地區(qū)，硬度下降幅度超過200 mg/L，陽離子交換作用是造成Ca、Mg含量下降的主要原因。南水北調供水及壓采后，水位恢復對地下水漏斗區(qū)域水質改善具有一定的積極作用。

（3）地表水的混合主要發(fā)生在山前的幾條補水河道中，這種混合作用對局部河段影響帶地下水質量的影響是比較大，但是由于補水期較短、補水量較小，且未見長期補水規(guī)劃，對區(qū)域地下水質量影響相對較小。因此，本次水文地球化學模擬工作主要考慮的是“南水北調來水后，本地的地下水壓采產生的水位抬升”，未能考慮地人工回補的地表水和地下水的混合作用，這是本次工作的缺憾，筆者下一步工作會開展相關研究，在模型中加入地表水混合作用。