反映我國空間分異特性的多介質環境逸度模型的構建及十溴二苯醚的歸趨模擬

鄢世陽，王中鈺，陳景文, *，李雪花，于洋，林軍

1. 工業生態與環境工程教育部重點實驗室，大連理工大學環境學院，大連 116024 2. 生態環境部固體廢物與化學品管理技術中心，北京 100029

當前，在美國化學文摘社(Chemical Abstracts Service, CAS)注冊的化學品種類超過1.6億種[1]，且以每日約1.5萬種的速度增加。全球市場使用的化學品已達35萬種[2]。化學品在促進社會發展、改善人類生活質量的同時，也為人體健康及生態健康帶來威脅。化學品釋放后經由遷移、轉化等行為，廣泛分布于多種環境介質中。揭示不同環境介質中化學品的分布和歸趨，是進行化學品風險評價和管理的前提。

自然環境是一個多介質/界面、時空連續的整體，僅采用環境監測手段，難以反映污染物在環境中分布、水平和歸趨[3]。多介質環境模型能夠模擬和預測化學品在環境中的遷移、轉化、分布和歸趨，是進行化學品風險評價和管理的必要工具[4]。基于逸度的多介質環境模型最早由Mackay等研發，具有結構清晰、便于計算等特點[5]，在化學品的環境歸趨預測中被廣泛應用[6]。多介質環境逸度模型基于質量平衡原理，以逸度為平衡判據，建立化學品在不同介質內的質量守恒關系式[7]。逸度模型分為Ⅰ到Ⅳ級，其中，Ⅲ級逸度模型考慮穩態非平衡條件，模擬區域環境條件恒定情況下化學品的歸趨，常被用于化學品環境暴露與持久性的預測[8-12]。

傳統的多介質環境逸度模型不具備對空間異質性的表征能力[13]。Wania等開發了結合地理信息系統(GIS)和大氣傳輸模型的空間分異多介質逸度模型GloboPOP[14]，但GloboPOP的空間分辨率(15°×15°)仍較低。Zhu等[15-16]開發了反映我國環境屬性的具有空間分異特性的多介質環境模型SESAMe，并成功用于苯并[a]芘、三氯生等化學品的環境歸趨模擬。基于GIS的多介質環境逸度模型系統，往往需要特定地理區域的環境屬性參數和所模擬的目標化學物質詳盡的環境行為參數[17]。我國的國土面積大，空間異質性強，發展具有高分辨空間異質性的多介質環境逸度模型具有重要意義。

本研究發展了覆蓋我國的國土空間、50 km×50 km分辨率的Ⅲ級多介質環境逸度模型系統，并以一種典型化學品十溴二苯醚(BDE-209)為例，介紹了模型的模擬效果。BDE-209被廣泛用作阻燃劑，添加于電子電器、紡織品、化工、家具、建材和塑料制品中[18]。2017年我國消費BDE-209共1.57萬t，占全國阻燃劑市場份額10.8%[19]。BDE-209在我國不同省份的多種環境介質中被檢出[20]。由于BDE-209具有長距離遷移性、難降解性和生物蓄積性等特點[21]，且可轉化為毒性更大的五溴二苯醚和八溴二苯醚[22]，已于2017年被列入“關于持久性有機污染物的斯德哥爾摩公約”的管控/限制清單[23]。因此，模擬BDE-209的環境歸趨行為，具有理論與現實意義。

1 研究方法(Methods)

1.1 研究區域

研究區域為我國全境，為67°24'E～137°15'E，3°28'N～50°34'N的區域。使用50 km×50 km分辨率將研究區域分為東西方向100個、南北方向114個等距柵格的柵格系統(圖1)，使用ArcGIS軟件統計全國五級河流及湖泊水庫分布、水土面積占比、多年平均降水量和溫度等環境屬性參數[24]，進而對不同的柵格進行環境屬性參數的賦值。所模擬區域的土地面積為1.01×1013m2，水面積為3.02×1011m2。

圖1 研究區域圖Fig. 1 Map of the study area

1.2 BDE-209的環境排放量

BDE-209在產品的生產、加工、使用、廢棄和回收的生命周期過程中，會釋放到大氣、水和土壤等環境介質中[25]。據統計，2017年我國BDE-209的總生產量為2.42萬t[19]。有研究表明，BDE-209的排放量與國民生產總值(GDP)呈正相關[26]，本研究按照每個柵格的GDP比值將BDE-209的總排放量分配至各個柵格單元中。每個柵格的GDP，使用ArcGIS軟件Spatial Analyst Tools統計的2010年全國第二次經濟普查1 km×1 km分辨率柵格人口分布標準化的GDP分布數據集[24]。

參照“歐洲自愿排放控制行動計劃(VECAP)”[27]提供的BDE-209全生命周期排放系數(表1)，應用排放因子法計算BDE-209在全生命周期的不同階段環境排放量，得到柵格單元內BDE-209向大氣、水和土壤的排放量，其計算公式如下：

表1 BDE-209生命周期排放系數(EF)[27]Table 1 Life cycle emission coefficient (EF) of BDE-209[27]

Mi,j=M×EFi,j(M=2.42×107kg)

(1)

Mi,total=∑Mi,j

(2)

Ei,total=Mi,total/8760

(3)

式中：Mi,j為在生命周期場景j下排入環境介質i的BDE-209量(kg)；EFi,j為排放因子(無量綱)；M為BDE-209的生產量(kg)；Mi,total為排入環境介質i的BDE-209總量(kg)，Ei,total為BDE-209的釋放速率(kg·h-1)，8 760為1年的小時數。

1.3 模型框架

本研究在非平衡、穩態的條件下，計算大氣、水、土壤和沉積物中BDE-209的環境濃度。基于Ⅲ級逸度模型對每個柵格建模[8]，其質量平衡方程如下：

大氣：E1+GA1cB1+f2D21+f3D31=f1(D12+D13+DR1+DA1)=f1DT1

(4)

水體：E2+GA2cB2+f1D12+f3D32+f4D42=f2(D21+D24+DR2+DA2)=f2DT2

(5)

土壤：E3+f1D13=f3(D31+D32+DR3)=f3DT3

(6)

沉積物：E4+f2D24=f4(D42+DR4+DA4)=f4DT4

為了求出tmin，本文采用蒙特卡羅算法，生成k個隨機整數xk，表示第k個時間段乘坐電梯的人數，然后根據電梯不用時所停放的層數、每一層人數的不同、電梯行駛每一層所消耗的時間、電梯最近接待原則等影響因素，考慮每個人等待電梯消耗的時間總和。

(7)

式中：Ei(i=1, 2, 3, 4)為BDE-209的排放速率(mol·h-1)；GAi為BDE-209隨介質平流輸入速率(m3·h-1)；cBi為平流輸入的BDE-209濃度(mol·m-3)；fi為介質i中BDE-209的逸度(Pa)；DRi為BDE-209在介質i的反應速率系數(mol·Pa-1·h-1)；DAi為從介質i平流輸出速率系數(mol·Pa-1·h-1)；Dij為從介質i向介質j的遷移速率系數(mol·Pa-1·h-1)；DTi為所有從介質i流失的D值總和(mol·Pa-1·h-1)。

借鑒Scheringer等[28]對柵格內化學品質量傳輸的設計，考慮柵格間的大氣、水體中化學品的傳輸過程。如圖2所示，每個柵格單元的周圍，有8個相鄰的柵格。各柵格單元通過氣流與水流連接，忽略大氣的氣流、洋流和地形因素。考慮周圍區域的排放和環境過程對每個柵格的影響，每個柵格中的大氣與水可與鄰近柵格進行質量交換。單個柵格內，大氣與水介質內1/8質量的化學品隨介質平流進入相鄰的柵格，同時獲取周圍柵格傳遞的平流介質內的化學品質量(圖2)，且周圍柵格獲得該單元傳輸的化學品質量是相同的。

圖2 各柵格內物質交換示意圖Fig. 2 Schematic diagram of substance inter-exchange among the grids

1.4 模型參數化及求解

模型的輸入參數包括環境屬性參數和BDE-209環境行為參數。環境屬性參數包含介質的面積(Ai)、水體深度(h2)、土壤/沉積物的有機碳含量以及環境介質的遷移速率等，其數值取自文獻[29-34]，如表2所示。本研究使用的部分環境屬性參數為Mackay等提供的推薦值，這部分參數的數據量龐大，獲取難度高，且以往研究顯示其對模型影響小[29]，故使用推薦值。截至2014年12月31日，全國城鎮土地總面積為890萬hm2[35]，約占我國國土面積的0.93%，本研究將城鎮地區歸于土壤介質，未單獨考慮。本研究的柵格單元模型應用的是經典Ⅲ級多介質環境逸度模型[7]，將森林區域亦歸并于土壤介質。

表2 本研究的部分環境屬性參數Table 2 Environmental attribute parameters used in the current study

各個柵格的介質流速、介質面積等，基于我國河道分布圖[24]使用ArcGIS進行統計。各個柵格的介質深度(包含大氣高度(h1)、沉積物厚度(h3)等)、土壤/沉積物的有機碳含量和傳質系數等，使用我國的均值及相關推薦值[33-34]。BDE-209的環境行為參數來自PubChem數據庫[36]與文獻[37-38]，如表3所示。

表3 BDE-209的物理化學性質及環境行為參數Table 3 Physical and chemical properties and environmental behavior parameters of BDE-209

使用C#語言編程求解模型，計算得到各地區不同環境介質的BDE-209濃度，將其與文獻報道的實測濃度對比，以驗證模型準確性。

1.5 靈敏度分析和不確定性分析

使用Morris法計算模型參數的靈敏度[39]，模型參數的靈敏度系數CS計算公式如下：

CS,i=(Y1.1,i-Y1.0,i)/(0.1×Y1.0,i)

(8)

CS,total=∑|CS,i|

(9)

式中：CS,i為模型參數在介質i的靈敏度系數，Y1.0,i和Y1.1,i分別為參數(如hi,Ai,KOW)取值1.1倍和1.0倍時介質i的濃度。CS,total為靈敏度系數CS,i絕對值加和。以CS,total>0.5為標準，篩選出靈敏度顯著的參數。

使用蒙特卡洛方法[40]分析模型結果不確定性，假定參數服從正態分布，隨機取值cS,total>0.5的參數，應用Crystal Ball軟件運行10 000次，計算各柵格不同環境介質的BDE-209濃度的變異系數及四分位差[40]。其計算公式如下：

CV=σ/μ

(10)

Q=Q3-Q1

(11)

式中：CV為變異系數，σ為標準差，μ為平均值；Q為四分位差(上四分位數與下四分位數的差)，Q3為上四分位數(即位于75%)，Q1為下四分位數(即位于25%)。

2 結果與討論(Results and discussion)

2.1 模型驗證

經計算獲得的BDE-209排放總量的空間分布如圖3所示。將BDE-209排放速率代入模型，得到BDE-209在我國各介質的濃度值分布(圖4)。將計算值與文獻中實測數據對比，考慮實際監測條件的復雜性及不確定性，實測值與絕大多數模型計算結果偏差在一個數量級內視為合理。比較各介質內環境預測濃度與文獻數據(表4)，可知模擬值與實測值吻合度較高，誤差在一個數量級的范圍內，說明模型結果可靠。

表4 中國不同地區BDE-209預測濃度值與實測濃度值的比較Table 4 Comparison of predicted BDE-209 concentrations with those measured values in different regions in China

圖3 BDE-209排放總量的空間分布Fig. 3 Spatial distribution of total emission of BDE-209

值得指出的是，以往的監測研究中，一些采樣點距離污染源比較近，或者距離城市等具有高BDE-209排放系數的區域比較近，導致相應的監測濃度較高，濃度的分布范圍廣。例如，Shi等[41]報告廣東地區沉積物中BDE-209的濃度范圍為32.7～2 105 ng·g-1。

本研究中，多介質環境逸度模型給出的濃度，均是各個柵格空間內的平均水平。此外，以往監測的大氣中BDE-209濃度，普遍高于本研究BDE-209在大氣中預測濃度，應該也是上述原因所導致。另外，在監測大氣中BDE-209濃度時，在低于方法檢測限的情況下，相應的濃度值往往被記錄為“未檢出”，本研究未選擇預測濃度與“未檢出”數據進行比較。

2.2 BDE-209的多介質歸趨與分布

模型預測BDE-209的濃度范圍：大氣中1.95×10-11～7.67×10 μg·m-3，水中1.89×10-9～9.19×104μg·m-3，土壤中1.97×10-12～7.8 μg·kg-1，沉積物中1.04×10-8～5.05×105μg·kg-1。上海、山東和廣東等東部地區大氣與土壤中BDE-209濃度較高，陜西、山西中部地區水與沉積物內BDE-209濃度較高。原因是東部地區BDE-209排放量較大，導致該地區BDE-209的環境介質濃度高；中西部地區在BDE-209排放水平較高的情況下，水資源量較少，其水介質與沉積物內BDE-209濃度高。

統計各柵格數據，得到BDE-209在我國的整體遷移、轉化狀況(圖5)，其在大氣、水、土壤和沉積物中的平均濃度分別為2.02×10-6μg·m-3、6.64×10-6μg·m-3、1.93 μg·kg-1和3.65×10 μg·kg-1。在穩態條件下，土壤中BDE-209的含量約占總量的79.31%，是BDE-209的主要匯；沉積物中BDE-209的含量占總量的19.34%，為BDE-209濃度最高的環境介質；BDE-209在水中的含量占總量的0.85%；BDE-209在大氣中的含量占總量的0.003%。BDE-209的相間遷移主要是由水體向沉積物遷移，為32.19 kg·h-1，占相間總遷移量的60.39%。由此可見，土壤和沉積物是BDE-209主要的匯。這是由于BDE-209的有機碳吸附系數(KOC)較高，易吸附分配于土壤與沉積物中。

圖5 BDE-209在我國環境的分配、遷移和轉化情況Fig. 5 Distribution, migration and transformation of BDE-209 in China

2.3 靈敏度分析

靈敏度分析結果如圖6所示。大氣高度(h1)與大氣停留時間(TA)的靈敏度系數較大，說明h1與TA對BDE-209在大氣中的濃度影響較大。水面積(AW)、水停留時間(TW)和水深(h2)等水環境屬性參數對水體中BDE-209濃度影響較大，表明水環境屬性參數(AW,TW,h2)是影響水相濃度的重要因素，AW越大、TW越長，越易增強BDE-209在大氣-水間的交換作用，使其從大氣向水體遷移，使其在水中的濃度增大。

圖6 參數靈敏度系數比較Fig. 6 Comparison of parameter sensitivity coefficients

BDE-209的KOC、土壤中降解半減期(t1/2(soil))和土壤面積(AE)這3個因素，對其在土壤中濃度影響較大。BDE-209的KOC較大，導致BDE-209易富集于土壤的有機質中。土壤中BDE-209的t1/2(soil)影響其在土壤中的持久性，BDE-209在土壤中具有較長的t1/2(soil)，是造成其在土壤中的累積的重要因素之一。土壤AE越大，越易增強BDE-209在大氣-土壤間的交換作用，使其在土壤濃度增大。

沉積物面積(AS)、沉積物介質停留時間(TS)和沉積物厚度(h4)等環境屬性參數與正辛醇-水分配系數(KOW)、KOC等BDE-209環境行為參數這5個因素，對其在沉積物中的濃度有較大影響，BDE-209具有較高的KOW(logKOW>4)[3]值，具有很強的疏水性，易積累于沉積物上，使其在沉積物濃度增大。沉積物的環境屬性參數(AS,TS,h4)影響BDE-209在水和沉積物間的交換作用，可導致其在沉積物中的濃度變化。

2.4 不確定性分析

不確定性分析結果如圖7所示，BDE-209在大氣、水、土壤和沉積物濃度的Q分別為7.60×10-7μg·m-3、1.98×10-6μg·m-3、1.08×10-7μg·kg-1和1.39×10-5μg·kg-1；CV分別為0.16、0.19、0.21和0.24。不確定性分析結果表明，沉積物的Q最大，且其CV結果最大，說明沉積物中BDE-209濃度離散程度較大，即BDE-209的沉積物濃度存在較大不確定性，因此，需要進一步規范沉積物參數以提高計算結果準確性。

圖7 BDE-209在各介質內的預測濃度分布Fig. 7 Probability distribution of simulated concentration of BDE-209 in each media

2.5 模型對比

Abbasi等[60]使用BETR-Global模型預測了全球BDE-209的濃度，其所預測的我國空氣中BDE-209的平均濃度為18.5 pg·m-3。本研究模型預測我國空氣中BDE-209的平均濃度為2.02 pg·m-3，與Abbasi等[60]報告的濃度值在一個數量級。BETR-Global模型的空間分辨率是15°×15°，而本研究中模型的分辨率是50 km×50 km (相當于0.5°×0.5°)。BETR-Global模型的一些單柵格計算濃度包含其他國家的排放情況，只考慮不同地區大氣與海水間的污染物交換，未考慮內陸水體對污染物傳輸的影響。這些主要因素，導致了2個模型預測值的差異性。

O’Driscoll等[61]使用基于Ⅲ級逸度模型的EQC模型計算了BDE-209在中國臺灣地區的歸趨，其結論是BDE-209主要富集于土壤與沉積物中(占總量的95.3%)，與本研究的結果一致。由于EQC模型不具有空間分異特異性，只能計算整體區域的污染物環境平均濃度，無法展示不同區域的污染物分布。

本研究對一些地表區域進行了簡化，將城市與森林區域均歸并于土壤介質。近年來，我國城市化快速發展，城市地表硬化面積快速增大，城市的建筑物表面及地面和土壤相比，對污染物具有不同的吸附能力，進一步的研究有必要考慮城市建筑物表面和地面對污染物賦存與歸趨的影響。此外，森林系統富含各種天然有機質，有可能富集各種疏水性有機污染物，進一步的研究，也需要考慮我國森林覆蓋率變化和森林結構變化對污染物賦存與歸趨的影響。基于GDP分配各網格BDE-209排放量時，忽略人口數量的影響；應用VECAP提供的全生命周期系數計算BDE-209的排放量時，其使用場景與我國存在一定差異，這些會導致BDE-209排放量計算與實際情況有偏差。我國境內河流的流向復雜，在50 km×50 km分辨率下逐個確定網格內河流的流向具有難度。本研究將大氣與水體的流向進行概化，設定柵格間的大氣與水體的交換是無方向，這也是以往多介質環境模型研究中的通行做法[15-16]。然而，不考慮流動介質的流向，會導致模擬結果的誤差。盡管如此，與以往的多介質環境模型相比，本研究所構建的模型具有空間分異特性，嵌入了我國地理環境的屬性數據，具有精細的內陸水文數據，可更詳細地描述化學品在我國不同區域的環境歸趨，有助于化學品的分區域監管。

本研究基于Ⅲ級逸度模型開發了適用于我國區域環境的空間分辨率為50 km×50 km的多介質環境逸度模型，利用該模型預測了BDE-209在我國環境介質的分布與濃度，預測結果與實測值吻合。上海、山東和廣東等東部地區大氣與土壤內BDE-209濃度較高，陜西、山西等中西部地區水與沉積物內BDE-209濃度較高。土壤和沉積物是BDE-209主要的匯。環境屬性參數(水面積、水停留時間和水深)及BDE-209的環境行為參數(正辛醇/水分配系數、有機碳吸附系數以及土壤中降解半減期等)，對模型結果影響較大。BDE-209在沉積物中的濃度存在較大不確定性。所構建的具有空間分辨特性且反映我國環境屬性特征的Ⅲ級多介質環境逸度模型，適用于我國環境中化學品的分布與賦存預測，有助于化學品的分區域風險預測與管理。

◆