蘇南沿江地塊2，4-DCP遷移擴散的時空變化特征

尹 勇，薛 琦，張 華，蘇 星，吳佳楠，張文藝

（1. 常州大學(xué) 環(huán)境與安全工程學(xué)院，江蘇 常州 213164；2. 江蘇龍環(huán)環(huán)境科技有限公司，江蘇 常州 213031）

我國長江中下游沿江地區(qū)化工企業(yè)眾多，這些企業(yè)在非正常工況下（如火災(zāi)爆炸、罐體泄漏、管線滲漏等）極易發(fā)生污染物泄漏，對周圍環(huán)境造成破壞，這些污染物一旦擴散進入長江，則會對長江水體造成污染。2014年《全國土壤污染狀況調(diào)查公報》顯示：長江三角洲部分區(qū)域土壤污染問題較為突出［1］。近年來，隨著化工園區(qū)污染預(yù)防、控制和管理措施的逐步完善，發(fā)生環(huán)境污染事故的風(fēng)險大幅降低，但不可否認(rèn)沿江化工園區(qū)內(nèi)的污水處理站、危廢庫、儲罐區(qū)等仍然是潛在污染源，存在污染物泄漏風(fēng)險，必須予以高度重視。

地下水?dāng)?shù)值模擬是應(yīng)用模擬軟件建立數(shù)學(xué)模型來預(yù)測地下水及其中溶質(zhì)的運移特征和未來發(fā)展趨勢的一種方法，是當(dāng)前預(yù)測、掌控沿江地區(qū)地下水污染問題的主要手段之一［2-3］。目前常用的地下水?dāng)?shù)值模擬軟件有GMS、Visual MODFLOW、FEFLOW和TOUGH2等。其中，GMS軟件具有良好的使用界面、強大的前處理和后處理功能以及優(yōu)良的三維可視效果，應(yīng)用最為廣泛［4-6］。鄭玉虎等［7］用GMS軟件模擬了地下水中COD和氨氮的遷移特征，并在模擬過程中考慮了土壤吸附作用對溶質(zhì)運移的影響，結(jié)果表明，研究區(qū)域土壤對COD和氨氮的吸附效果顯著。高志鵬等［8］將HYDRUS-2D軟件和GMS軟件聯(lián)用，模擬了氮素在河流-非飽和帶-飽和含水層中的遷移特征，結(jié)果表明，氮素在地下水中的濃度逐漸升高，最終保持穩(wěn)定。KHAYYUN等［9］利用GMS軟件建立了MODFLOW和MT3DMS模型，研究了3種情景下工業(yè)廢水對伊拉克底格里斯河左右兩側(cè)地下水水質(zhì)的影響，結(jié)果表明底格里斯河左側(cè)地下水水質(zhì)受工業(yè)廢水影響較大。RAHNAMA等［10］等調(diào)查了伊朗平原地區(qū)地下水的砷污染狀況，并運用GMS軟件模擬預(yù)測砷濃度的變化，發(fā)現(xiàn)黏土層厚度與砷濃度之間存在地理相關(guān)性。

本工作以蘇南沿江某農(nóng)藥生產(chǎn)企業(yè)所在區(qū)域為研究對象，以該企業(yè)突發(fā)高濃度2，4-二氯苯酚（2，4-DCP）污水泄漏事故為情景，采用概念模型法，運用GMS軟件建立地下水流場模型，以實測水位數(shù)據(jù)擬合模擬水位，通過模型識別驗證獲得研究區(qū)塊的地下水流場模型，在此基礎(chǔ)上建立三維地下水溶質(zhì)運移模型，推演并分析地下水中2，4-DCP的遷移擴散和時空分布特征，判斷其對周邊環(huán)境及長江水體的影響，為沿江化工地塊土壤及地下水污染管控提供技術(shù)支持。

1 研究區(qū)概況

某農(nóng)藥生產(chǎn)企業(yè)位于蘇南沿江某工業(yè)園區(qū)（約10 km2），地貌單元屬長江三角洲沖積平原。地面坡度小于0.5%，地形平坦，西北稍高，東南略低，平均地形高程約4.50 m，最高 5.80 m，部分地區(qū)僅2～3 m。地質(zhì)構(gòu)造屬于揚子古陸東端的下?lián)P子白褶帶。第四系地層厚度160～200 m。

蘇南沿江屬北亞熱帶濕潤性季風(fēng)氣候區(qū)，四季分明、雨量充沛、日照充足、無霜期較長，年平均氣溫約16.3 ℃，歷年平均降雨量1 096 mm、最大降雨量1 888 mm、最小降雨量592 mm。

根據(jù)地質(zhì)資料、現(xiàn)場調(diào)查和鉆孔數(shù)據(jù)，確定該區(qū)域的地層范圍屬第四系全新統(tǒng)及上更新統(tǒng)長江下游三角洲沖積層，巖性主要為填土、淤泥質(zhì)粉質(zhì)黏土、粉質(zhì)黏土、粉土和粉砂。水質(zhì)檢測結(jié)果顯示該區(qū)域地下水化學(xué)類型主要為HCO3-Ca·Na型水，礦化度0.5～1.2 g/L，pH約6.2～11.4。地下水類型按潛藏條件分為孔隙潛水和承壓水。其中，潛水主要埋藏于填土、淤泥質(zhì)粉質(zhì)黏土和粉質(zhì)黏土中，其水位埋深一般為0.3～2.0 m；承壓水主要埋藏于粉土和粉砂中，其埋深約10 m。蘇南長江段5～10月為豐水期，11月之后進入枯水期，全年水位標(biāo)高0.98～6.00 m，地下水水位會受其影響產(chǎn)生變化。潛水層水位變化主要受降水影響，年最大值大部分出現(xiàn)在7月，變化幅度0.55～1.28 m，年最小值常出現(xiàn)在枯水期，變化幅度0.06～0.30 m；承壓水層年最高水位出現(xiàn)在7月，變化幅度0.61～1.64 m，年最低水位出現(xiàn)在1～2月，變化幅度0.03～0.40 m。

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

2.1 水文地質(zhì)概念模型

研究區(qū)處于一個相對獨立的水文地質(zhì)單元，北部以長江為定水頭邊界；東西兩側(cè)河流與研究區(qū)地下水存在水力聯(lián)系，故可概化為與河流平均狀態(tài)相關(guān)的定水頭邊界；南部為自然的隔水邊界。垂向上，研究區(qū)潛水層上部邊界為水量交換邊界，接受大氣降雨入滲補給、河流入滲補給等，模型底部以不透水泥巖為隔水邊界，具體模擬范圍及邊界條件見圖1。根據(jù)研究區(qū)水文地質(zhì)條件和地層巖性特征，污染物進入地下主要污染潛水層，從環(huán)境最不利角度出發(fā)，模擬層位設(shè)為潛水層和承壓水層。

圖1 研究區(qū)范圍及邊界條件示意圖

本研究區(qū)內(nèi)地下水水流運動遵循質(zhì)量守恒原理，符合達西定律和滲流連續(xù)性方程。在空間維數(shù)上，以水平流為主，最終流入長江；垂向上，含水層之間存在水力流動，符合空間三維流特征。研究區(qū)水流系統(tǒng)的補徑排特征、水位和流速隨時間的變化產(chǎn)生差異，符合非穩(wěn)定流特征。

2.2 研究區(qū)地下水?dāng)?shù)學(xué)模型

依據(jù)區(qū)域水文地質(zhì)概念模型，本研究區(qū)的地下水流可概化為非均質(zhì)、各向異性、空間三維結(jié)構(gòu)、非穩(wěn)定的地下水流系統(tǒng)，其對應(yīng)的數(shù)學(xué)模型見式（1）。

式中：t為時間，d；Ω為模型模擬區(qū)；H為含水層水位，m；Kx、Ky、Kz分別為x、y、z方向上的滲透系數(shù)，m/d；μs為貯水率，m-1；W為含水層的源匯項，m3/d；H0為初始水位，m；Γ1為滲流區(qū)域的一類邊界；Γ2為滲流區(qū)域的二類邊界；n→為邊界Γ2的外法線方向；K為三維空間上的滲透系數(shù)張量，m/d；q為二類邊界上已知流量函數(shù)，流入為正、流出為負(fù)、隔水邊界為0。

2.3 初始條件

將模擬區(qū)域內(nèi)13個地下水監(jiān)測孔的水位（監(jiān)測初始時間為2021年1月）作為模擬預(yù)測的初始水位和模型識別驗證的標(biāo)準(zhǔn)。

2.4 源匯項

區(qū)域地下水補給來源主要為垂向補給，來自大氣降水入滲。研究區(qū)年均降雨量為1 096 mm，主要巖性為粉質(zhì)黏土，因此，降水入滲補給系數(shù)取值0.18［11］。

區(qū)域地下水的主要排泄方式為蒸發(fā)，地下水蒸發(fā)量與其埋深有關(guān)。研究區(qū)年水面蒸發(fā)量為1 515.9 mm，地下水蒸發(fā)量比水面蒸發(fā)量小得多，研究區(qū)地下水位埋深約為2 m，主要巖性為粉質(zhì)黏土，因此，蒸發(fā)系數(shù)取值0.06［11］。

2.5 水文地質(zhì)參數(shù)選取

地下水模擬參數(shù)的取值見表1。參照《環(huán)境影響評價技術(shù)導(dǎo)則 地下水環(huán)境》（HJ 610—2016）［12］附錄表B.1和B.2，潛水層水平滲透系數(shù)取0.1 m/d，給水度取0.02；由于蘇南沿江為黏性土質(zhì)，降雨入滲補給系數(shù)取0.000 52 m/d［12］。弱透水層與承壓水層的滲透系數(shù)采用抽水試驗（采用一個抽水井帶兩個觀測井的非穩(wěn)定流抽水試驗方法）來確定，最終獲得弱透水層、承壓水層的水平滲透系數(shù)均為1.529 8 m/d。垂直滲透系數(shù)與水平滲透系數(shù)的比值按經(jīng)驗值取0.1。

表1 地下水模擬參數(shù)的取值

2.6 流場模型的建立與識別校正

在GMS軟件中將模型在XY方向上以80×80進行剖分，垂直方向上分為3層，每個單元格長寬均為50 m，共計12 800個單元格，有效活動單元格為5 728個。將表1的參數(shù)輸入GMS軟件，建立研究區(qū)的地下水流場模型。

模型的識別和驗證在地下水?dāng)?shù)值模擬過程中十分重要，是對所建立模型的一個校核過程，目的在于確保對該研究區(qū)的地下水?dāng)?shù)值模擬是正確的，所建立的模型是精確可靠的。主要原理是通過不斷修正降雨入滲補給系數(shù)、滲透系數(shù)等水文地質(zhì)參數(shù)，重復(fù)計算，直到模擬結(jié)果與水位觀測點的實測值之間的誤差在可接受范圍內(nèi)。GMS軟件提供了一套輔助校核工具，可以幫助用戶驗證實測水位與模擬水位的吻合性。

在建立地下水流場模型后，根據(jù)所監(jiān)測的水位資料對模型進行識別校正，最終得到校正后的地下水流場（見圖2）。由圖2可見：水流整體往東北方向流動，匯入長江，沿水流方向水位逐漸遞減，但減幅不大。由于研究區(qū)的中部小塊區(qū)域缺少部分實測點的水位數(shù)據(jù)，導(dǎo)致模擬值與實測值具有一定差別，但總體上看，校準(zhǔn)后的地下水流模型能夠較好地反映該區(qū)域地下水流場的特征。校核后的水文地質(zhì)參數(shù)見表2。

表2 校核后的水文地質(zhì)參數(shù)

圖2 校正后的地下水流場

選擇研究區(qū)域內(nèi)13個監(jiān)測井的實測水位與模擬水位進行擬合，設(shè)置實測值與計算值的極差為0.5 m，置信區(qū)間為95%，結(jié)果見圖3。由圖3可見：算點都均勻分布在45°線兩側(cè)，R2為0.973 6，表明擬合程度較高，誤差較小。

圖3 水位實測值與計算值的擬合曲線

圖4為水位擬合誤差情況圖，根據(jù)誤差棒顏色的不同，可以判斷誤差的大?。壕G色表示實測值與計算值的誤差在置信區(qū)間范圍內(nèi)；橙色表示誤差超出置信區(qū)間，但不超過置信區(qū)間的200%；紅色表示誤差超出置信區(qū)間200%。由圖4可見，只有一個點的誤差棒指示為橙色，其余均為綠色，沒有出現(xiàn)紅色。模擬結(jié)果中水位絕對誤差小于0.5 m的觀測點占已知水位觀測點的92%，遠高于《地下水資源管理模型工作要求》（GB/T 14497—1993）［13］的規(guī)定：預(yù)測水位值絕對誤差在0.5 m內(nèi)的觀測點須占所有已知水位觀測點的70%以上。這表明所建立的模型精度符合要求。

圖4 水位擬合誤差情況圖

綜上，從地下水流場和觀測孔水位對比兩方面看，水位觀測孔擬合結(jié)果良好，地下水流場宏觀擬合效果良好，模型識別的水文地質(zhì)參數(shù)與水文地質(zhì)條件基本符合，能夠反映地下水流動特征，具有較高的精度。因此，可基于該地下水流模型模擬2，4-DCP的遷移轉(zhuǎn)化特征。

3 2，4-DCP的遷移轉(zhuǎn)化

3.1 2，4-DCP運移模型的控制方程

模擬2，4-DCP的遷移時只考慮對流、彌散和吸附作用，描述三維地下水流系統(tǒng)中某種污染物遷移運動的控制方程見式（2）。

式中：R為滯后因子；θ為地下介質(zhì)的孔隙度；ρ為溶質(zhì)的質(zhì)量濃度，mg/m3；xi，xj分別為沿X軸和Y軸的距離，m；Dij為水動力彌散系數(shù)張量，m2/d；vi為滲流或線性孔隙水流速度，m/d，與單位流量或達西流量qi有關(guān)，vi=qi/θ；qs為單位體積含水層源（正值）和匯（負(fù)值）的體積流量，1/d；ρs為源匯流中溶質(zhì)的質(zhì)量濃度，mg/m3；λs為吸附相的反應(yīng)速率常數(shù)，1/d；λ為溶解相的反應(yīng)速率常數(shù)，1/d；ρb為多孔介質(zhì)密度，mg/m3。

3.2 污染情景及參數(shù)設(shè)置

本次模擬情景設(shè)置為該農(nóng)藥廠地塊內(nèi)污水收集管破損，導(dǎo)致工業(yè)廢水通過土壤包氣帶進入地下水造成污染。假設(shè)裂口面積約為0.01 m2，污水滲漏量約為0.7 m3/d，污染類型為連續(xù)入滲。經(jīng)調(diào)查取樣分析，該工業(yè)廢水中污染物以苯系物為主，其中2，4-DCP質(zhì)量濃度較高，為107 mg/L。因此，選擇2，4-DCP作為典型污染因子進行污染物遷移規(guī)律研究。

模擬期為7 300 d（20 a），模擬輸出時間段為100 d，1 000 d，1 825 d（5 a），3 650 d（10 a），5 475 d（15 a），7 300 d（20 a）。本次模擬不考慮化學(xué)反應(yīng)、生物降解以及外部源匯項（溶質(zhì)通過源進入或者通過匯離開研究區(qū)域），即忽略擴散過程中2，4-DCP的生物降解和化學(xué)分解等因素。

模擬過程中相關(guān)參數(shù)的設(shè)置見表3。其中：2，4-DCP的干預(yù)值參照荷蘭干預(yù)值標(biāo)準(zhǔn)［14］；吸附常數(shù)經(jīng)實驗室等溫吸附實驗測得；縱向彌散度通過實驗室土柱實驗和現(xiàn)場試驗獲得，橫向彌散度通常比縱向彌散度低一個數(shù)量級，垂向彌散度比橫向彌散度低一個數(shù)量級［15-16］，橫縱彌散度比值取0.2，垂縱彌散度比值取0.05；孔隙度采用沿江地區(qū)的經(jīng)驗值［17］。

表3 模擬過程中相關(guān)參數(shù)的設(shè)置

3.3 2，4-DCP遷移擴散過程中的時空變化特征

初始化三維地下水溶質(zhì)運移模型，將相關(guān)參數(shù)導(dǎo)入已建立的地下水流場模型，運行后得到2，4-DCP在不同模擬時間條件下的時空變化特征，結(jié)果見圖5和表4。如圖5可見：2，4-DCP沿地下水水流朝東北方向呈橢圓形擴散，污染面積逐漸擴大（圖5a）；當(dāng)模擬時間為5 475 d（15 a）時，2，4-DCP在水平方向上已經(jīng)擴散至長江江岸（圖5e）；當(dāng)模擬時間為7 300 d（20 a）時，2，4-DCP在水平方向上已經(jīng)擴散進入長江（圖5f）。

表4 模擬過程中2，4-DCP的遷移情況

圖5 不同模擬時間條件下2，4-DCP在第一層含水層中的時空分布特征（埋深0.2～1.5 m）

由表4可見：在0～100 d內(nèi)，2，4-DCP的水平遷移距離為121.15 m，擴散速率最快，平均擴散速率為1.211 5 m/d；100～1 000 d時，2，4-DCP的擴散速率有所降低，水平遷移距離為64.40 m，平均擴散速率為0.071 6 m/d；1 000～1 825 d時，水平遷移距離為53.19 m，平均擴散速率為0.064 5 m/d；1 825～3 650 d時，水平遷移距離為94.41 m，平均擴散速率為0.051 7 m/d；3 650～5 475 d時，水平遷移距離為76.17 m，平均擴散速率為0.049 3 m/d；5 475～7 300 d時，2，4-DCP的擴散速率最低，水平遷移距離為57.35 m，平均擴散速率降低至0.037 1 m/d。綜上，隨著模擬時間的延長，2，4-DCP的遷移擴散速率在不斷降低。如果不采取任何污染控制措施，到7 300 d（20 a）后，2，4-DCP在水平方向上的最大遷移距離為466.67 m；污染面積將達到110 982 m2，而且，隨著時間的推移，污染面積還會持續(xù)增大。由表4還可看出：在污染前期（0～3 650 d），隨著模擬時間的延長，污染羽中心2，4-DCP的質(zhì)量濃度逐漸增大；當(dāng)模擬時間為7 300 d時，污染羽中心的2，4-DCP質(zhì)量濃度達到最大，為41.453 mg/L，超過荷蘭干預(yù)值標(biāo)準(zhǔn)（0.03 mg/L）1 375倍。而且，隨著時間的推移，污染羽中心2，4-DCP的濃度還會持續(xù)增大。

圖6為模擬時間7 300 d時垂直方向上地下水中2，4-DCP的分布情況。由圖6可見：在垂直方向上，2，4-DCP由潛水層先后穿透弱透水層和承壓水層，擴散至模型底部，最大污染深度可達40 m；隨著地下水深度增加，污染物擴散速度明顯減慢，由初始的-5 m處的41.453 mg/L擴散至-25 m處的4.607 mg/L大約經(jīng)歷7 300 d（約20 a），表明2，4-DCP的垂直擴散速率相對較低，這與江水豐水期的頂托作用有關(guān)。

圖6 模擬時間為7 300 d（20 a）時垂直方向上地下水中2，4-DCP的分布

以上分析表明：土壤吸附對2，4-DCP的滯留作用有限，2，4-DCP很容易穿過粉質(zhì)黏土層，對地下水造成污染；2，4-DCP的凈化主要依靠含水層的稀釋作用，這意味著該區(qū)域地下水一旦受到污染，只依靠地下水的循環(huán)是難以實現(xiàn)自凈的。由此可見，該廠工業(yè)廢水的泄漏對周邊地區(qū)的土壤地下水環(huán)境安全以及長江生態(tài)存在一定威脅，應(yīng)及時采取有效措施防止2，4-DCP的擴散。

4 結(jié)論

由于沿江地塊場地土壤及水文地質(zhì)條件的復(fù)雜性和不確定性，很難建立起科學(xué)、準(zhǔn)確反映污染物彌散、吸附和降解等環(huán)境行為的數(shù)學(xué)模型?；诖?，本工作采用一種間接模擬方式，即先利用GMS軟件構(gòu)建能反映研究區(qū)域的地下水流場模型，然后利用地下水觀測井實測數(shù)據(jù)檢驗、校對地下水流場模型，從而建立起能準(zhǔn)確反映本研究區(qū)域地下水流場的模型。在此基礎(chǔ)上，利用建立的流場模型構(gòu)建溶質(zhì)運移模型，揭示污染物擴散規(guī)律，即推演地下水中2，4-DCP在遷移擴散過程中的時空分布特征，判斷其是否會對周邊環(huán)境及江水造成污染影響。

a）模擬結(jié)果表明：2，4-DCP將會隨著地下水流往長江方向遷移，100 d后，2，4-DCP的水平遷移距離為121.15 m，1 000 d后為185.55 m，5 a后為238.74 m，10 a后為333.15 m，15 a后為409.32 m（屆時2，4-DCP將擴散至長江江岸），20 a后為466.67 m；20 a后污染羽中心2，4-DCP的質(zhì)量濃度將達到41.453 mg/L，最大污染深度為40 m，污染面積為110 982 m2，2，4-DCP將擴散至長江水體?？梢姡搹S工業(yè)廢水的泄漏對周邊地區(qū)的土壤地下水環(huán)境安全以及長江生態(tài)存在一定威脅，應(yīng)及時采取有效措施防止2，4-DCP的擴散。

b）本文所建立的沿江地區(qū)污染物擴散與預(yù)測方法，對于沿江化工地塊土壤及地下水污染管控有一定的參考價值。