抗生素已在全球范围内得到广泛生产和使用^[1-2]，主要用于治疗、预防人类和动物疾病以及促进畜禽生长^[3-7]，我国抗生素的使用量也非常大，2013年抗生素的使用量达到16.2×10⁴ t^[8]。磺胺类抗生素作为常用的六大抗生素之一，已在诸多环境介质中被检出^[9-11]。与其他类型抗生素相比，磺胺类抗生素的吸附系数较低，具有良好的水溶性，因此具有更高的迁移性^[4]。Srinivasan等^[12]对新西兰的奶牛养殖区域牧草土壤进行采样，以研究磺胺类药物的迁移转化，发现磺胺类药物容易淋失出土壤剖面，运移至地下水，或经排水和地表径流进入地表水。目前，对于磺胺嘧啶在环境中的行为研究多集中在静态吸附，而对于磺胺嘧啶在土壤中的迁移，尤其是在原状土柱中的迁移特征及其模拟研究较少。

原状土柱基本保留土壤的原始结构，因而其中的溶质运移更能反映野外实际情况。已有研究表明，原状土柱中的优先流是异丙隆在土壤中迁移的主要机制^[13]。Zsolnay等^[14]研究发现阿特拉津在原状土柱中的迁移速度要远大于扰动土柱。国内外学者运用Hydrus-1D模型深入研究了不同溶质在土壤中的运移过程^[15-17]。任理和毛萌^[18]对阿特拉津在稳定流场饱和砂质土壤中的运移过程用化学非平衡的两点模型和两区模型进行模拟，发现化学非平衡对阿特拉津运移的影响较物理非平衡高。Dousset等^[19]对异丙隆在原状土柱和填装土柱中的迁移特征进行分析，并分别用化学平衡和非平衡的双渗透模型(Dual-permeability)模拟了其在两种土柱中的穿透曲线(Breakthrough Curve, BTC)，发现较高的优先流和较低的平衡吸附点位使得异丙隆在原状土柱中的迁移更快，化学非平衡的双渗透模型对其穿透曲线模拟效果最好。

为探讨磺胺类抗生素在土壤中的运移特征，本文以磺胺嘧啶(Sulfadiazine, SDZ)为研究对象，通过其在原状土柱中的出流试验，探讨其在土壤中的迁移过程，用Hydrus-1D软件对试验结果进行模拟，并进一步运用模拟得到的参数，对不同水流流速下磺胺嘧啶的迁移行为进行预测，以期为磺胺嘧啶对地下水的污染风险评价提供理论依据。

1 材料与方法 1.1 供试材料

磺胺嘧啶，纯度99%，购自北京百灵威科技有限公司；pK_a值为1.57/6.50^{[3, 20]}。

2016年3月、6月和8月三次在青岛市即墨市蓝村镇的同一块农田的0~15 cm和15~30 cm处取原状土土柱，分别设为原状土柱Ⅰ和Ⅱ，并同时取土样进行基本理化性质分析，见表 1。

1.2 示踪试验

以KBr为示踪剂，通过土柱混合置换试验，启动蠕动泵，先用0.01 molL^-1的CaCl₂溶液，自下而上饱和土柱，待土柱饱和，出流液稳定后，自上而下输入1个孔隙体积(Pore Volume，PV)的0.05 mol L^-1的KBr溶液，然后用0.01 mol L^-1的CaCl₂冲洗，用自动部分收集器收集出流液，直至无法检测到Br^-为止。控制流速尽量保持一致或相近，每组做两个重复。

孔隙体积PV=v×t/l，式中，v为平均孔隙水流速(cm min^-1)，t为时间(min)，l为土柱长度(cm)。

1.3 原状土运移试验

去除表面土壤，用内径5 cm、高15 cm的有机玻璃(下端为尖锐的切面)，平稳缓慢地压入土层，开挖剖面取出土柱，该过程对土壤的扰动性较小，可以忽略不计。以此方法分别在0~15 cm、15~30 cm深处取原状土柱，装入试验用有机玻璃柱，土柱两端分别填装石英砂，并放置与土柱内径相当的滤纸，以使溶液均匀渗入，并在土柱两端接口处，分别加入适量脱脂棉，以防止土壤颗粒堵塞。启动蠕动泵，先用0.01 mol L^-1的CaCl₂溶液自下而上饱和土柱，待土柱饱和并出流稳定后，再自上而下输入3个PV、pH为6的30 mg L^-1的SDZ，其后用0.01 mol L^-1的CaCl₂溶液冲洗，直至无SDZ检出。出流液稀释后，用紫外分光光度计(UV9100C PC，北京莱博泰科仪器有限公司)测定溶液中磺胺嘧啶的浓度^[21]。

1.4 运移模型

均质土柱中的水流和溶质运移常被视为垂向一维的，稳定流条件下，只考虑吸附作用，溶质在一维均匀介质中的运移过程可用对流弥散方程(Convection Dispersion Equation，CDE)来描述^[22]。方程可表示如下：

$ \frac{{\partial c}}{{\partial t}} = D\frac{{{\partial ^2}c}}{{\partial {x^2}}} - v\frac{{\partial c}}{{\partial x}} - \rho \frac{{\partial s}}{{\partial t}} $

(1)

式中，c为液相中的溶质浓度(mg L^-1)，t为时间(min)，v为平均孔隙流速(cm min^-1)，D=λv，D为弥散系数(cm² min^-1)，λ为弥散度(cm)，ρ为容重(g cm^-3)，x为距离(cm)，s为单位质量土壤吸附的溶质浓度(mg g^-1)。

对于非反应性溶质，如Br^-，∂s/∂t为0，因而可以通过土壤中Br^-、Cl^-等的穿透曲线，利用确定性平衡模型估算参数v和D。此时方程如下：

$ \frac{{\partial c}}{{\partial t}} = D\frac{{{\partial ^2}c}}{{\partial {x^2}}} - v\frac{{\partial c}}{{\partial x}} $

(2)

对于线性吸附，阻滞因子R_d的表达式为：

$ {R_{\rm{d}}} = 1 + \frac{{\rho {K_{\rm{d}}}}}{\theta } $

(3)

式中，K_d为线性吸附分配系数(g L^-1)，θ为体积含水量(cm³ cm^-3)。

非平衡单点模型(One-site Sorption Model，OSM)假定吸附过程为动力学吸附，控制方程为式(4)~式(6)。

非平衡两点吸附模型(Two-site Sorption Model，TSM)是用来描述化学非平衡的基本模型，该模型假定吸附点位可以分为两部分^[23-24]：类型Ⅰ假定吸附是瞬时的，用吸附平衡等温方程来描述；类型Ⅱ假定吸附是受速率限制(依赖于时间)的，并遵从一阶动力学方程^{[23, 25]}，控制方程为式(7)、式(8)。

非平衡两区模型(Two-region Model，TRM)是在CDE基础上建立起来的，主要考虑了不动水体在溶质运移过程中所起的作用，更真实的反映了土壤孔隙状况，被广泛应用于土壤中饱和及非饱和溶质运移规律研究^[26]。TRM模型假定液相被分为两个区域，即“可动(Dynamic)”区域和“不可动(Stagnant)”区域。溶质在土壤中的运移受土壤结构(如团聚体间孔隙和微孔隙等)限制，对流-弥散运移被限定在可动区域；可动与不可动区域间溶质的交换(吸附速率)受溶质扩散到不可动区域交换点的限制，用一级动力学方程来描述^[25-27]。在此模型中，认为所有吸附点的吸附总是处于平衡的(即吸附是瞬时的)，且吸附性溶质和非吸附性溶质都受物理非平衡的影响^[26]。控制方程为式(9)、式(10)^{[22-23, 26, 28-29]}：

$ \frac{{\partial \theta c}}{{\partial t}} + \rho \frac{{\partial s}}{{\partial t}} = \frac{\partial }{{\partial t}}\left( {\theta D\frac{{\partial c}}{{\partial x}}} \right) - \frac{{\partial qc}}{{\partial x}} $

(4)

$ \rho \frac{{\partial s}}{{\partial t}}=\alpha \rho \left[ {{s_{\rm{e}}} - s} \right] $

(5)

$ {s_{\rm{e}}} = {K_{\rm{d}}}c $

(6)

$ \left( {1 + \frac{{f\rho {K_{\rm{d}}}}}{\theta }} \right)\frac{{\partial c}}{{\partial t}} + \frac{\rho }{\theta }\frac{{\partial {s_2}}}{{\partial {x^2}}} = D\frac{{{\partial ^2}c}}{{\partial {x^2}}} - v\frac{{\partial c}}{{\partial x}} $

(7)

$ \frac{{\partial {s_2}}}{{\partial t}} = \alpha \left[ {\left( {1 - f} \right){K_{\rm{d}}}c - {s_2}} \right] $

(8)

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\left( {{\theta _{\rm{m}}} + f\rho {K_{\rm{d}}}} \right)\frac{{\partial {c_{\rm{m}}}}}{{\partial t}} + \left[ {{\theta _{{\rm{im}}}} + \left( {1 - f} \right)\rho {K_{\rm{d}}}} \right]\frac{{\partial {c_{{\rm{im}}}}}}{{\partial t}} = \\ {\theta _{\rm{m}}}{D_{\rm{m}}}\frac{{{\partial ^2}{c_{\rm{m}}}}}{{\partial {x^2}}} - q\frac{{\partial {c_{\rm{m}}}}}{{\partial x}} \end{array} $

(9)

$ \left[ {{\theta _{{\rm{im}}}} + \left( {1 - f} \right)\rho {K_{\rm{d}}}} \right]\frac{{\partial {c_{{\rm{im}}}}}}{{\partial t}} = {\alpha _{\rm{t}}}\left( {{c_{\rm{m}}} - {c_{{\rm{im}}}}} \right) $

(10)

式中，c为溶质浓度(mg cm^-3)，D为水动力弥散系数(cm² min^-1)，q为水流通量(cm min^-1)，s_e为吸附平衡时单位质量土壤吸附的溶质浓度(mg g^-1)，s为动力学吸附的溶质浓度(mg g^-1)，α为描述动力学吸附过程的一阶速率常数(min^-1)，s₂指类型Ⅱ吸附位点上的吸附浓度(mg g^-1)，f是在平衡时发生瞬时吸附的交换位点所占的分数，下标m和im分别指可动和不可动的区域，θ_m+θ_im=θ，q=θ_mV_m是容积水通量密度(Volumetric water flux density，cm min^-1)，α_t是描述在可动和不可动区域间溶质交换速率的一阶质量传递系数(min^-1)。

2 结果与讨论 2.1 Br^-的穿透曲线

由图 1可见，Br^-在原状土柱中的穿透曲线并不对称，且出流较早，在第一个流出液样品中，即15min(0.24 PV)内检测到Br^-，说明在原状土柱Ⅰ和Ⅱ中存在大孔隙优先流。由于连续大孔隙的存在，使得Br^-可以在原状土柱中快速迁移出流，原状土柱Ⅰ和Ⅱ的Br^-穿透曲线基本一致，说明此两层土壤结构性质相近，垂向空间异质性很弱。

2.2 磺胺嘧啶在原状土中的运移特征

对六组原状土土样进行土柱出流试验，其试验结果如图 2所示。第一组、第二组和第三组试验分别为青岛市蓝村镇3月份、6月份和8月份土样。

虽然三组土样均取于同一块农田，但因为土壤性质的空间变异性，SDZ在各土柱中迁移状况不同，因此所得的SDZ的穿透曲线也有所差异，对于第一组和第二组试验，原状土柱Ⅱ的穿透曲线相较于原状土柱Ⅰ发生左移，两组穿透曲线的峰值为0.8~0.9，不同深度穿透曲线峰值差别不大。根据表 1，随着采样深度的增加，土壤pH有所升高，阳离子交换量和土壤有机质含量降低。pH升高使得更多的SDZ呈现阴离子状态，不利于SDZ在土壤中的吸附，此外土壤有机质的降低也使得SDZ与土壤的作用力及交换能力减弱，因而SDZ在原状土柱Ⅱ中迁移更快。残留在土壤中的SDZ，随着深度的增加，更容易发生迁移，进而增加地下水污染的风险。第三组的穿透曲线峰值为1左右，且出流较快，很快便达到峰值，对土柱进行反冲洗时，SDZ浓度下降很快，推测造成这种现象的原因，一是取样前已对土地进行深耕，使得所取原状土柱含有大孔隙，产生优先流。另一个原因是土壤的pH变化所致，第一组原状土柱Ⅱ和第三组出流液的pH均较大(未列图)，其出流速度也较快。

2.3 Br^-穿透曲线参数确定

对于原状土柱，由于Br^-为非反应性物质，首先用LEA模型对其在原状土柱中的迁移进行模拟，获得参数θ_s，v和λ；由于原状土柱保持土壤的原有结构，存在不动水区域，因此用两区模型，在估算θ_s、v和λ的同时，估算不可动水所占分数θ_im，拟合的相关参数见表 2。由图 3可以看出，两个模型对Br^-在原状土柱中的拟合效果均较好，拟合曲线基本重合，拟合的决定系数R²可达0.997以上，均方根误差RMSE的值为0.021，其中不同模型拟合的不同深度的弥散度λ差别不大，对于平均孔隙流速v，原状土柱Ⅰ的LEA模型和原状土柱Ⅱ的TRM模型拟合值与实测值0.0625 cm min^-1较为接近，拟合得到的原状土柱Ⅰ和Ⅱ的不动水区域分别为0.152和0.105 cm³ cm^-3。

2.4 磺胺嘧啶在原状土柱中的穿透曲线及参数确定

根据表 2用LEA模型拟合得到的原状土柱相关参数，固定θ_s、v和λ，分别用OSM和TSM对磺胺嘧啶在原状土柱中的迁移进行模拟；固定θ_im、θ_s、v和λ，用TRM对磺胺嘧啶在原状土柱中的迁移进行模拟，对不同模型的模拟结果进行对比分析，拟合结果见图 4和表 3。

由表 3可见，三个模型均能较好地拟合SDZ在原状土柱中的迁移，除第二组的原状土柱Ⅰ之外，决定系数R²均大于0.974，均方根误差RMSE均在0.052以内，拟合的BTCs基本重合，TRM拟合结果最好，OSM的拟合结果较TSM好。TSM拟合得到的f值除第一组的原状土柱Ⅱ外均很小，基本上可以忽略不计，说明SDZ在原状土柱中基本上不存在瞬时平衡吸附。推测可能是因为大孔隙的存在，增强了SDZ的流动性，从而减少了土柱对SDZ的瞬时吸附。由于SDZ在土壤中降解缓慢且试验周期不长，模型不考虑SDZ的降解，反演模拟得到的在可动区域吸附位点所占分数f均在0.154以内，说明不可动区域也是SDZ在土壤中吸附的重要部分。第一组和第二组拟合得到的吸附系数K_d值较第三组大，且用TRM反演模拟得到的K_d要大于OSM和TSM，这与SDZ在原状土柱中的BTCs相一致。吸附平衡经验系数β的变化趋势与之相反，α的数值均较小，且无明显变化规律。可见，TRM对SDZ在原状土柱中的模拟效果最好。

2.5 磺胺嘧啶在土壤中运移的预测

由表 3可知，TRM能更好地模拟磺胺嘧啶在原状土柱中的迁移过程，第二组试验结果较第一组和第三组更接近实际情况，因此，利用第二组原状土柱反演模拟得到的相关参数，对不同流速下SDZ在土壤剖面中的运移进行模拟预测，以分析其对浅部地下水可能产生的污染风险。设SDZ的输入浓度为30 mg L^-1，输入持续时间为10 PV，流速分别为0.017、0.030和0.100 cm min^-1；土壤剖面深度为200 cm，分为两层，其中第一层为0~15 cm，相关参数与第二组的原状土柱Ⅰ相同，第二层为15~200 cm，相关参数与原状土柱Ⅱ相同；分别在土壤剖面的15、30、50、100、150、200 cm处设置观察点，预测结果如图 5，出流特征见表 4。

由表 4可见，在同一深度土层中，随着流速的增加，SDZ出流时间、到达峰值时间及出流持续时间均逐渐变小，SDZ的出流浓度逐渐增大，说明低流速时，SDZ与土壤接触时间长，增加了吸附量。从图 5可以看出，同一流速时，土壤剖面中SDZ的浓度峰值随着深度增加而减小，且出流时间逐渐增大。低流速时，SDZ的穿透曲线具有较好的对称性，当流速增加到0.100 cm min^-1时，穿透曲线不对称，并具有不同程度的“拖尾”现象，说明高流速时，SDZ在土壤中的解吸过程较为缓慢。在200 cm处，水流流速为0.100 cm min^-1时，SDZ在3.90 PV时即可出流，可见SDZ在土壤中的迁移非常快，当水流流速较高时，可快速穿过土层进入地下水，对地下水存在较高的污染风险。

3 结论

磺胺嘧啶在原状土柱Ⅱ中迁移速度要大于原状土柱Ⅰ，这一方面是受优先流的影响，另一方面是受土壤理化性质，如pH、CEC和有机质含量等的影响。运用Hydrus-1D中的非平衡单点模型(OSM)、非平衡两点模型(TSM)和非平衡两区模型(TRM)对SDZ穿透曲线进行反演模拟发现，TRM可更好地模拟SDZ在原状土柱中的迁移过程，拟合的R² > 0.94，RMSE < 0.046。在此基础上，用TRM对不同水流流速时SDZ在200 cm厚的土壤剖面中的迁移过程进行了预测。结果显示，同一深度处，随着流速从0.017 cm min^-1增加到0.030 cm min^-1再到0.100 cm min^-1，磺胺嘧啶的出流时间越来越短，而出流的浓度峰值越来越大，当流速为0.100 cm min^-1时，SDZ可快速穿过整个土壤剖面。可见，高流速时，SDZ对地下水污染存在较大风险。