反映我国空间分异特性的多介质环境逸度模型的构建及十溴二苯醚的归趋模拟

鄢世阳，王中钰，陈景文, *，李雪花，于洋，林军

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在沉积物中的浓度存在较大不确定性。所构建的具有空间分辨特性且反映我国环境属性特征的Ⅲ级多介质环境逸度模型，适用于我国环境中化学品的分布与赋存预测，有助于化学品的分区域风险预测与管理。

◆