0 引言

环境介质中有机污染物浓度的时间变化已得到广泛研究, 其中出现较多的一个现象是部分有机污染物浓度的时间变化滞后于其排放历史.了解这种现象对研究有机污染物的迁移转化规律, 分析长期环境监测数据, 评估环境管理政策的效果, 以及预测环境介质中有机污染物的浓度有重要意义.目前, 已经在冰雪^[1]、冰芯^[2]和沉积物^[3-7]中发现这种滞后现象.

除了直接对环境样品取样测定外, 研究人员也从环境模型角度对其进行了研究, 但迄今为止国内报道较少. Czub和McLachlan^[8]使用食物链模型ACC-HUMAN讨论了生物体对环境浓度变化的响应, 指出不易降解的物质在人体中峰值浓度的滞后时间更长. Lim和Lastoskie^[9]使用湖泊和生物蓄积模型研究发现, 六溴联苯(hexabromobiphenyl, hexaBB)在环境介质和生物体中的时间变化均滞后于其排放历史, 营养级最高的鲑鱼体内hexaBB浓度峰值比排放峰值滞后4年. Quinn等^[10]使用CoZMoMAN模型研究了多氯联苯(polychlorinated biphenyls, PCBs)的代际差异, 发现婴儿出生前暴露于PCBs中的浓度峰值, 滞后于其排放约10年. Gouin和Wania^[11]使用全球分布模型Globo-POP的研究结果中, 三氟乙酸等主要分布在水中且不易降解的物质, 在北极地区的浓度变化相对于其排放滞后最明显. Choi和Wania^[12]使用区域污染物归趋模型CoZMoPoP 2, 系统研究了半衰期对有机物浓度时间变化滞后的影响, 提出了环境污染"可恢复"(reversibility)的概念.这个概念也被Persson等^[13]和MacLeod^[14]等用来从环境污染角度, 完善Rockström等^[15]提出的行星边界(Planetary Boundaries)理论.

虽然环境介质中一些有机污染物浓度的时间变化滞后现象已被研究人员所了解, 但其影响因素在国内外并未得到充分研究.为此, 选择131种有机物, 通过构建EQC (EQuilibrium Criterion)模型中标准环境的四级模型, 计算得出4种排放场景下环境介质中有机物浓度的时间变化, 在此基础上初步探讨滞后时间的影响因素.

1 材料与方法 1.1 EQC模型简介

选择EQC(EQuilibrium Criterion, 平衡标准)模型进行研究. EQC模型是一个通用模型, 于1996年提出^[16]后, 已被广泛用来进行有机化学品筛选和风险评估工作^[17-24].目前已有多款基于EQC模型的计算程序和软件可供使用, 例如加拿大Trent大学环境模型中心发布的Level Ⅲ 2.80软件^[25]和美国环境保护总署发布的EPI Suite软件^[26]等. Hughes等^[27]在2012年对EQC模型软件进行了更新, 使其得到了更科学更灵活的应用, 但环境参数并未更改.目前已有的EQC模型计算程序和软件中包括了一级(稳态、平衡、非流动系统)、二级(稳态、平衡、流动系统)和三级(稳态、非平衡、流动系统)模型.

在化学品筛选和风险评估工作中, 使用一种通用且典型的环境系统, 来对所有化学品进行无差别研究, 比使用实际环境系统更为方便^[28].这里选择EQC三级模型中的标准环境进行研究, 其面积为100 000 km$^{2}$, 相当于中国浙江省或江苏省的面积.现有的EQC模型计算程序和软件均只能计算出稳态的结果, 不能用其来描述系统的动力学特征.为了研究环境介质中有机物浓度的时间变化, 必须引入与时间有关的变量, 因此本研究需在EQC三级模型和模型中标准环境的基础上, 重新构建动力学模型, 即四级(非稳态、非平衡、流动系统)模型.

1.2 模型的构建

EQC标准环境包括水、大气、沉积物和土壤4个相.每个相中的微分质量平衡方程为: d$M$/d$t$=输入速率-输出速率.其中$M$为有机物质量, 单位为mol; $t$为时间, 单位为a.用逸度表示, $M$=ZVf, 其中$Z$为逸度容量, 单位为mol$\cdot $Pa$^{-1}\cdot $m$^{-3}$; $V$为相体积, 单位为m$^{3}$; $f$为逸度, 单位为Pa.根据每个相中有机物的输入和输出过程编写微分质量平衡方程, 共得到以下4个微分方程^[26]:

水(下标1):

$ \begin{align} Z_1V_1\dfrac{{\rm d}f_1}{{\rm d}t}=E_1+f_{2}D_{21}+f_{3}D_{31}+f_{4}D_{41}-f_1(D_{12}+D_{13}+D_{1{\rm R}}+D_{\rm 1out}), \end{align} $

(1)

大气(下标2):

$ \begin{align} Z_2V_2\dfrac{{\rm d}f_2}{{\rm d}t}=E_2+f_{1}D_{12}+f_{4}D_{42}-f_{2}(D_{21}+ D_{24}+D_{2{\rm R}}+D_{\rm 2out}), \end{align} $

(2)

沉积物(下标3):

$ \begin{align} Z_3V_3\dfrac{{\rm d}f_3}{{\rm d}t}=E_3+f_{1}D_{13}-f_3(D_{31}+D_{3{\rm R}}), \end{align} $

(3)

土壤(下标4):

$ \begin{align} Z_4V_4\dfrac{{\rm d}f_4}{{\rm d}t}=E_4+f_{2}D_{24}-f_4(D_{41}+D_{42}+D_{4{\rm R}}). \end{align} $

(4)

方程中各参数的定义见表 1. Mackay在其著作^[29]中详细讨论了各参数的计算方法, 这里不再列出.

1.3 排放速率

一般来说, 一种有机物排放进入环境中的量会呈现出先升高后降低的时间分布规律^[12].为了模拟这种排放场景, 选择如下的排放速率函数进行计算^[9]:

$ \begin{align} E_n=C\exp\Big(\dfrac{t_{\max}-t}{A}\Big)\Big/\Big(1+\exp\Big(\dfrac{t_{\max}-t}{A}\Big)\Big)^2. \end{align} $

(5)

其中$E_{n} (n=1$, 2, 3, 4)为排放速率, 单位为mol$\cdot $a$^{-1}$; $t_{\max}$为最大排放速率对应的时间, 取$t_{\max}$=10 a; $t$为时间, 单位为a; $A$为常数, 取$t<t_{\max}$时$A=0.4$, $t>t_{\max}$时$A=0.85$; $C$也为常数, 单位为mol$\cdot $a$^{-1}$, 与最大排放速率相关, 由于本研究只关注时间变化, $C$的取值对此无影响.

1.4 排放场景

有机物有可能直接排放进入的环境介质有水、大气和土壤, 一般不直接进入沉积物中, 因此设定4种不同的排放场景进行研究, 具体见表 2.其中场景1-3分别为有机物仅通过水、大气和土壤排放, 用来模拟单一的排放场景; 场景4为有机物同时向水、大气和土壤中排放, 用来模拟复合排放场景.所有场景中向沉积物的排放速率均设为0.

1.5 模型的计算

模型的输入参数总结见表 3. 1.2节中的方程均属于一阶线性常微分方程, 使用MATLAB 2012b软件中的Simulink工具进行建模和计算^[30].计算时, 使用ode15s方法, 最大步长设为0.1 h, 其余计算参数均取默认值.计算得出环境介质中有机物浓度时间变化曲线后, 再编写程序读取有机物峰值浓度对应的时间, 减去排放峰值浓度对应的时间, 即为环境介质中有机物峰值浓度的滞后时间.

1.6 有机物的选择

选择131种有机物, 包括8种多环芳烃、8种PCBs、5种二噁英和26种农药等常见有机污染物.这些有机物的理化性质和半衰期大部分取自文献[29], 其余根据EPI Suite软件^[26]计算其理化性质, 根据相关文献来设定其半衰期^[31-33].所有有机物的理化性质均取25 ℃下的值.

1.7 数据统计

使用Minitab 17软件进行正态分布检验(使用Ryan-Joiner检验方法)、相关分析和回归分析, $P<$0.05认为具有统计学意义.

2 结果与讨论 2.1 典型有机物的峰值浓度滞后时间

研究涉及的有机物有131种, 限于篇幅仅列举典型有机物的计算结果.为了选取典型有机物, 需要考查模型输入参数(见表 3)和模型结构.除各相中的排放速率函数为假设的函数外, 模型输入参数可分为两类, 第一类参数描述有机物的固有性质, 在任何情况下均为常数, 与环境条件无关, 这类参数包括蒸汽压、水溶解度、lg$K_{OW}$、熔点和分子量; 第二类参数为4种环境介质中的半衰期, 其并非常数且与环境条件有关.在第一类参数中, 蒸汽压、水溶解度和lg$K_{OW}$在模型计算中参与程度较高, 均参与了各介质中逸度容量的计算; 而熔点和分子量参与计算的程度较低, 其中熔点只是用来判断模型中气溶胶逸度容量的计算是否使用过冷蒸汽压, 分子量只是用来对水溶解度数据进行单位换算, 因此这两个参数在选择典型有机物时不予考虑.使用模型中参与计算程度较高的7种参数对有机物进行分类研究.为方便分类, 将蒸汽压、水溶解度、lg$K_{OW}$和持久性(表示半衰期的长短)分为"高"和"低"两个水平, 其中"高"分别表示值超过10 Pa、100 mg/L、3和150 d (即3 600 h)的水相中半衰期.蒸汽压高于10 Pa是常用的判断挥发性有机物的标准^[34], 水溶解度高于100 mg/L的物质则常被认为是易溶于水^[35], lg$K_{OW}$高于3的物质常被认为具有一定的生物蓄积性^[36], 水相中半衰期超过150 d的物质常被认为不可降解^[37].按照这种方法, 有机物可被分为16类.其中的11类可在所选择的131种物质中确定对应的有机物, 剩余5类不能确定对应的持久性为"高"的有机物, 因此对文献[29]中的半衰期进行调整, 详见表 4.在不同的环境条件下, 有机物的半衰期有可能改变, 因此这种调整是有意义的. 表 4中16类典型有机物在环境介质中峰值浓度的滞后时间见表 5.

从表 4和表 5总结可见:

(1) 大多数情况下, 环境介质内有机物浓度的时间变化出现滞后, 表明这种滞后现象在环境系统中普遍存在.

(2) 不同介质中的滞后时间不同.在大多数情况下, 土壤相和沉积物相中的滞后时间相对于水相和大气相中的滞后时间更长, 这是因为EQC模型中对水相和大气相加入了平流过程, 使其传质过程较土壤相和沉积物相更快.

(3) 不同类型有机物的滞后时间不同.有机物性质会影响其滞后时间, 不同性质的有机物参与的主要传质过程不同, 因此其滞后时间也会有差异.其中lg$K_{OW}$和持久性对滞后时间的影响, 比水溶解度和蒸汽压更大. lg$K_{OW}$和持久性都为"高"的有机物(分类1、5、7和13)有最长的滞后时间.目前报道的环境介质中浓度变化出现明显滞后的有机污染物均为这种物质^[1-7].

2.2 滞后现象的影响因素

为了研究滞后时间的影响因素, 对蒸汽压、水溶解度、lg$K_{OW}$、持久性(用水相中半衰期表示)和滞后时间进行相关分析.正态分布检验结果表明, 大部分数据不服从正态分布, 并且由于数据中存在0值, 不能进行Box-Cox变换^[38], 因此使用Spearman相关系数, 结果见表 6.从中可见: ①蒸汽压和水溶解度与滞后时间均无相关性. ② lg$K_{OW}$只与沉积物相的滞后时间呈极显著正相关, 而与其他相的滞后时间无相关性. ③场景1沉积物相中, 持久性与滞后时间呈显著正相关; 其他情况下持久性与滞后时间均呈极显著正相关.总结以上3点, 水相、大气相和土壤相中滞后时间的影响因素主要为持久性, 而沉积物相中滞后时间的影响因素主要为lg$K_{OW}$和持久性.这是因为大气相中的传质过程较快, 而水相和土壤相均与大气相接触, 相界面的传质过程也较快, 所以在这3个相中, 传质速率较少受介质间的传质过程(用蒸汽压、水溶解度和lg$K_{OW}$来计算)制约, 而介质内的传质过程(用持久性来计算)成了整个传质过程的瓶颈; 沉积物相仅与水相接触, 其介质间和介质内的传质过程都较缓慢, 所以用来计算水-沉积物相界面传质过程的参数lg$K_{OW}$, 以及计算沉积物相内传质过程的参数持久性都会对影响沉积物相中的滞后时间.

2.3 滞后现象的预测 2.3.1 滞后时间的预测公式

在131组有机物的数据基础上, 增加表 4和表 5中的五组半衰期调整后的数据, 共使用136组数据, 来对滞后时间和模型输入参数进行逐步线性回归.由于蒸汽压、水溶解度和半衰期的取值范围较大, 对其先取对数后再进行回归分析.每个相的预测变量仅选择与相内或相间传质过程相关的参数.所得结果见表 7.

研究中选择的131种有机物中, 分类1、3、5、7、9、11和15中的有机物数据较少, 预测公式主要根据剩余9类有机物的数据获得.由于样本量较少, 公式对上述7类有机物的滞后时间预测是否准确, 仍需进一步验证.

研究得出的预测方法仅基于EQC模型和本研究假设的排放函数, 而在实际环境系统中, 温度、流速等环境参数变化较大, 无论使用何种模型, 都很难对其中的滞后时间进行准确预测.但预测方法在环境管理和风险评估中仍可用来进行初步的使用, 如用来筛选可能对生态环境产生长期甚至不可恢复影响的物质.

2.3.2 沉积物中滞后现象的预测

在4种环境介质中有机污染物浓度的长期连续监测研究中, 水相、大气相和土壤相中已发表的数据很少^[39-41], 并且数据量和时间跨度暂不满足与模拟结果的对比讨论; 而沉积物相的样品可通过$^{210}$Pb等方法定年, 相关研究相对较多.这里选择沉积物相, 来对其中有机物的滞后现象做进一步分析.

沉积物相中的主要传质过程包括水相对沉积物相的输入、沉积物相向水相的输出和沉积物相内的转化输出.沉积物埋藏也是一种传质过程, 但由于其传质速率非常缓慢, 可忽略.在模型中, 描述这3种传质过程的$D$值分别为$D_{13}$、$D_{31}$和$D_{3{\rm R}}$.计算输入过程和输出过程$D$值的差, 见表 8.从中可见, 沉积物中滞后时间较长的分类中(分类1、5、7和13), ($D_{13}-D_{31}-D_{3{\rm R}})$的值均为正值.因此, 可通过($D_{13}-D_{31}-D_{3{\rm R}})$这个参数, 来判断沉积物相中是否发生明显滞后.物理意义上, ($D_{13}-D_{31}-D_{3{\rm R}})$的符号可用来指示沉积物相传质过程的方向, 符号为正表示沉积物相属于有机物的"汇", 有机物更易在沉积物相中积累, 且这个积累过程更长, 超过了排放量增长的过程, 因此有机物浓度变化产生滞后.

在已有的报道中, 沉积物相中出现滞后现象的物质有DDT和总PCBs^[3-7].计算其($D_{13}-D_{31}-D_{3{\rm R}})$值, 分别为5.1$\times $10$^{12}$和1.5$\times $10$^{12}$, 均为正值, 可见用($D_{13}-D_{31}-D_{3{\rm R}})$的符号来判断沉积物中是否发生明显滞后现象是可行的.另一方面, ($D_{13}-D_{31}-D_{3{\rm R}})$的符号有明确的物理意义, 使用范围较表 7中的公式更广泛, 不局限于EQC模型.

虽然通过($D_{13}-D_{31}-D_{3{\rm R}})$的符号可初步判断沉积物相中是否发生明显的滞后现象, 但使用这种方法对有机物直接排入的相进行研究时发现, 由于排放函数的引入, 并不能获得有意义的结论, 这也是本研究中使用的常微分方程组数值解法的劣势.若使用常微分方程组的解析解法, 可对滞后现象进行更深入研究^[42].但由于本研究中构建的模型包括的方程数量较多, 模型解析解的推导过程也较复杂, 对模型的解析解理论仍需开展进一步研究.

3 结论及展望

(1) 本研究通过构建EQC标准环境的四级多介质模型, 计算得出了131种有机物在4种排放场景下, 环境介质中有机物浓度的时间变化.

(2) 在模拟结果中, 大多数情况下环境介质内有机物浓度的时间变化, 相对于排放历史出现滞后, 表明这种滞后现象普遍存在.土壤相和沉积物相中的滞后时间相对水相和大气相更长, lg$K_{OW}$和持久性高的有机物浓度变化滞后最明显.

(3) 持久性为水相、大气相和土壤相中滞后时间的主要影响因素, lg$K_{OW}$和持久性为沉积物相中滞后时间的主要影响因素.对滞后时间和模型输入参数进行逐步线性回归, 得到了EQC模型中滞后时间预测公式.对于沉积物相, 可用($D_{13}-D_{31}-D_{3{\rm R}})$的符号来判断其中有机物浓度变化是否会产生明显的滞后.

(4) 模型的解析解法可对滞后现象进行更深入的研究, 但需对其求解方法进行完善.

(5) 目前, 在环境监测数据分析, 环境管理政策制定和效果评估, 有机污染物在环境介质中的浓度预测等工作中, 均未考虑环境介质中有机污染物浓度时间变化的滞后现象.未来需要对这种滞后现象展开更广泛更深入的研究.