寧 辛，林文婷，歐忠輝

(福建師范大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，福建 福州 350117)

隨著工業(yè)化和城市化進(jìn)程加快，土壤作為污染物的最終受體，污染排放已超過(guò)土壤環(huán)境容量[1].在發(fā)生“鎘米”、“毒地”等污染事件后，土壤污染才進(jìn)入公眾視野[2-6].Cd是一種廣泛的、不可生物降解的有毒重金屬，長(zhǎng)期穩(wěn)定地存在于環(huán)境中，并通過(guò)農(nóng)田生態(tài)系統(tǒng)和食物鏈累積對(duì)人類健康造成有害影響，例如對(duì)腎臟和骨骼造成嚴(yán)重?fù)p害，加劇了人類患癌的風(fēng)險(xiǎn)[7-12].目前可以通過(guò)生物炭、PH和有機(jī)水溶肥等方法來(lái)降低植物對(duì)鎘的吸收[13-15]，從而減少Cd在作物體內(nèi)的富集[16-18].

在土壤溶液中，除了自由離子外，Cd還與可溶性的有機(jī)配體形成絡(luò)合物[19-20].Barber模型可以用來(lái)描述養(yǎng)分離子(氮、磷、鉀) 的吸收機(jī)制[21]，但該模型只考慮了一種物質(zhì)(水溶態(tài))，不能正確描述Cd絡(luò)合物在根表面的解離和緩沖植物根系吸收Cd2+[22-23].對(duì)于Cd的絡(luò)合物的解離，Van[24]建立了自由離子運(yùn)動(dòng)模型，說(shuō)明絡(luò)合物解離速率取決于根表面的通量和根圍區(qū)域自由離子的擴(kuò)散.2002 年，Van 等[25]又建立了游離鎘離子(m)、 配體(l) 和絡(luò)合物(ml) 的守恒模型來(lái)模擬了m和ml在不同動(dòng)力學(xué)條件下的濃度分布，從而推導(dǎo)出絡(luò)合物處于完全惰性、不完全穩(wěn)定和部分不穩(wěn)定時(shí)根表面鎘離子的吸收通量，然后與Kouteck-Koryta(KK) 方法導(dǎo)出的通量進(jìn)行比較，從而確定了KK近似解的有效性范圍.Sterckeman等[26]將Barber模型應(yīng)用到作物對(duì)Cd的吸收.在Barber的基礎(chǔ)上，Lin等[27]建立了植物吸收鎘的機(jī)理模型，從而模擬出Cd在溶液中的擴(kuò)散和對(duì)流遷移、吸附和絡(luò)合動(dòng)力學(xué)以及根系吸收，但未給出模型的解析解.在Lin等[27]建立的模型的基礎(chǔ)上，Schneider等[28]考慮了完全惰性、完全不穩(wěn)定和部分不穩(wěn)定情況下模型的解析解且通過(guò)偏微分方程數(shù)值解驗(yàn)證解析解的有效性，但是他只給了根表面的Cd2+的濃度和通量，未給出整個(gè)根圍中的濃度分布.本文在Schneider等[28]工作的基礎(chǔ)上，計(jì)算出鎘的絡(luò)合物處于3種狀態(tài)下整個(gè)根圍鎘離子的濃度分布的全局解析解，再使用根圍鎘離子濃度的解析表達(dá)式來(lái)替代模型進(jìn)行計(jì)算和預(yù)測(cè)，并通過(guò)根長(zhǎng)密度估計(jì)整個(gè)根系對(duì)鎘的吸收，從而確定控制根吸收鎘的有效性因素.本文的工作將豐富有毒重金屬的污染生態(tài)理論，對(duì)土壤環(huán)境的可持續(xù)性提供理論依據(jù).

1 植物對(duì)鎘的吸收模型

植物根系只吸收溶液中游離的鎘離子.在溶液中，m和l的絡(luò)合反應(yīng)決定了m和ml之間的相互轉(zhuǎn)換.在絡(luò)合反應(yīng)中，溶液中絡(luò)合物的絡(luò)合速率ka和解離速率kd決定了自由離子和絡(luò)合物的相互轉(zhuǎn)換速率.m在溶液中以Dm的擴(kuò)散速率向根表面遷移，而絡(luò)合物的擴(kuò)散速率比自由離子慢得多，Dm遠(yuǎn)遠(yuǎn)大于Dml.在反應(yīng)層(μ) 中，當(dāng)根表面的m的濃度降低時(shí)，ml會(huì)進(jìn)行解離并以自由離子的形式被根吸收，從而絡(luò)合物解離的速率也取決于絡(luò)合的動(dòng)力學(xué)和根表面自由離子的濃度.在絡(luò)合反應(yīng)下，土壤溶液中含有m、l和ml.在土壤中，鎘固態(tài) ms的解吸也可以緩沖自由離子的生物利用而減少在擴(kuò)散層(δ) 中自由離子的擴(kuò)散(如圖1).

圖1 Cd的絡(luò)合物在土壤溶液中的解吸、吸附、對(duì)流、擴(kuò)散和Cd2+的吸收過(guò)程[28]Fig.1 The desorption, adsorption, convection, diffusion of Cd complex in soil solution and absorption of Cd2+[28]

Schneider 等[28]根據(jù)鎘離子以對(duì)流和擴(kuò)散的遷移方式，建立了游離鎘離子、配體和絡(luò)合物在土壤中的反應(yīng)擴(kuò)散的模型:

(1)

其中ci(i=m，l，ml) 是土壤溶液中游離鎘離子、配體和絡(luò)合物的濃度，bi是土壤緩沖力，t是時(shí)間，r是距離根中心的距離，f是液態(tài)下土壤的阻抗，θ是土壤含水量，Di是擴(kuò)散系數(shù)，r0是根半徑，v0是根表面水通量，ka為溶液中絡(luò)合物的絡(luò)合速率常數(shù)，kd為溶液中絡(luò)合物的分離速率常數(shù).

通過(guò)米氏動(dòng)力學(xué)來(lái)實(shí)現(xiàn)植物在根系表面吸收鎘離子

其中Jm是植物在根系表面吸收鎘離子的通量，Imax是金屬離子被根系吸收的最大速率，Km是米氏常數(shù).

根表面與初始濃度之間的區(qū)域(即δ(t)) 稱為擴(kuò)散層，Schneider等[28]把右邊界設(shè)定在濃度等于初始濃度的位置，該位置隨時(shí)間逐漸遠(yuǎn)離根表面

Schneider等[28]根據(jù)絡(luò)合物在土壤溶液中處于完全惰性、完全不穩(wěn)定和部分不穩(wěn)定分別建立子模型，并得到了穩(wěn)態(tài)的通解、根表面的濃度和通量解析表達(dá)式.由于根表面的濃度表達(dá)式復(fù)雜，Schneider等[28]沒有進(jìn)一步利用左右邊界條件確定通解的待定系數(shù)和根圍中濃度的表達(dá)式.

當(dāng)溶液中絡(luò)合物處于完全惰性時(shí)，由于絡(luò)合物的解離速率很慢，暫時(shí)不考慮配體和絡(luò)合物對(duì)根系吸收鎘離子的影響，建立了絡(luò)合物處于完全惰性時(shí)的子模型:

(2)

m的遷移通量等于根表面吸收m的通量

(3)

完全惰性下的右邊界為[28]

(4)

其中bm,t是土壤對(duì)鎘離子和絡(luò)合物的緩沖力，

此時(shí)，m和ml的遷移通量等于根表面吸收m的通量，即

(5)

(6)

完全不穩(wěn)定下的右邊界為[28]

當(dāng)溶液中絡(luò)合物處于部分不穩(wěn)定時(shí)，需要考慮cl，并且cml和cm不再保持比例關(guān)系，因此需要用到式(1).由于式(1)是在柱坐標(biāo)下，且又是多個(gè)方程，導(dǎo)致這個(gè)子模型是不可解.采用近似方法進(jìn)行求解: 先將式(1)化為平面坐標(biāo)下的方程，求出部分不穩(wěn)定和完全惰性下穩(wěn)態(tài)解的比值，再通過(guò)比值乘以柱坐標(biāo)下完全惰性的穩(wěn)態(tài)解來(lái)近似表示柱坐標(biāo)下部分不穩(wěn)定的穩(wěn)態(tài)解.右邊界δpl(t)=(δin(t)-δfl(t))μ/δin(t)+δfl(t) 采用線性近似的方法得到，其中μ是反應(yīng)層厚度.近似方法的合理性需要數(shù)值模擬驗(yàn)證.

Schneider等[28]將通解代入式(3)和式(5)后，可以得到根表面的濃度，進(jìn)而可以計(jì)算Cd2+在根表面上的通量.但是根表面的濃度表達(dá)式已經(jīng)很復(fù)雜，再代入方程的通解后，表達(dá)式變得更加復(fù)雜，計(jì)算難度很大，并未給出整個(gè)根圍中的Cd2+濃度分布.本文的主要工作是借助Mathematica給出通解中待定常數(shù)的表達(dá)式和根表面濃度的完整表達(dá)式，分別求出絡(luò)合物處于完全惰性和完全不穩(wěn)定下模型的穩(wěn)態(tài)解; 采用柱坐標(biāo)下完全惰性的穩(wěn)態(tài)解加權(quán)從而近似表示絡(luò)合物處于部分不穩(wěn)定的穩(wěn)態(tài)解.

2 模型求解

(7)

其中A1=Imax-Kmv0，A2=Imax+Kmv0，A3=[r0/(r0+δin)]Pe，A4=[r/(r0+δin)]Pe.

將r=r0代入式(7)，得到根表面的鎘離子濃度

(8)

絡(luò)合物處于完全不穩(wěn)定下當(dāng)絡(luò)合物處于完全不穩(wěn)定下，在根表面的Cd2+濃度為

+C1C8

(9)

當(dāng)r=r0時(shí)，由式(9)得到根表面的鎘離子濃度

[-C1(C5+C7)(1+K*)+Imax(C5-C6)(δpl+λK*tanh(δpl/λ))]/

式(7) 和式(9)是絡(luò)合物分別處于完全惰性和完全不穩(wěn)定時(shí)的穩(wěn)態(tài)解，本文用數(shù)值解與穩(wěn)態(tài)解比較來(lái)驗(yàn)證解析方法的正確性; 式(9)是用柱坐標(biāo)下完全惰性的穩(wěn)態(tài)解加權(quán)后近似表示柱坐標(biāo)下部分不穩(wěn)定的穩(wěn)態(tài)解，并不是子模型(1)的穩(wěn)態(tài)解，精確穩(wěn)態(tài)解和近似穩(wěn)態(tài)解是否充分接近，也需要用數(shù)值解來(lái)驗(yàn)證.數(shù)值解還用來(lái)驗(yàn)證本文右邊界確定的合理性.

3 解析解的驗(yàn)證和數(shù)值解的計(jì)算

表1 模型參數(shù)和參數(shù)值[27]Tab.1 Model parameters and parameter[27]

圖2顯示3種情況鎘離子濃度的變化曲線.在圖2中，t0=1 d，t1=5 d和t2=10 d內(nèi)右邊界的平均移動(dòng)速度約為: 8.10×10-8，3.47×10-8，2.43×10-8cm·s-1，這說(shuō)明右邊界的移動(dòng)速度與絡(luò)合物的狀態(tài)無(wú)關(guān)，但隨著時(shí)間的增加而減慢.在某一時(shí)刻，鎘離子的濃度隨距離根表面的距離增加而增加.

圖2 絡(luò)合物處于完全惰性、完全不穩(wěn)定和部分不穩(wěn)定情況下Cd2+濃度解析解和數(shù)值解Fig.2 Comparison of analytical and numerical solutions of Cd2+ concentration when the complex is fully inert,fully labile,or partially labile

圖2也可以說(shuō)明3種情況下的鎘離子濃度變化速率不同.在距離根表面0.055 cm處，當(dāng)絡(luò)合物處于完全惰性時(shí)，1-5 d的鎘離子濃度變化速率大約為4.48×10-13μmol cm-3·s-1，5-10 d的鎘離子濃度變化速率大約為1.04×10-13μmol cm-3·s-1; 當(dāng)絡(luò)合物處于完全不穩(wěn)定時(shí)，1-5 d的鎘離子濃度變化速率大約為3.56×10-12μmol·cm-3·s-1，5-10 d的鎘離子濃度變化速率大約為9.03×10-13μmol·cm-3·s-1; 當(dāng)絡(luò)合物處于部分不穩(wěn)定時(shí)，1-5 d的鎘離子濃度變化速率大約為4.21×10-13μmol·cm-3·s-1，5-10 d的鎘離子濃度變化速率大約為9.95×10-14μmol·cm-3·s-1，這說(shuō)明鎘離子濃度的變化速率隨時(shí)間的增加而下降且絡(luò)合物處于完全不穩(wěn)定時(shí)鎘離子濃度變化速率最大.在圖2中可以看到，當(dāng)鎘離子濃度達(dá)到穩(wěn)態(tài)時(shí)，根表面濃度已經(jīng)接近于0，而右邊界的濃度始終為初始濃度，穩(wěn)態(tài)解幾乎呈直線，解析解和數(shù)值解在根圍中間部分有偏離.

圖3 不同初始濃度下根表面Cd2+濃度隨時(shí)間的變化Fig.3 Changes of Cd2+ concentration on root surface with time under different initial concentrations

圖3顯示了在不同初始濃度下，根表面鎘離子濃度的變化.當(dāng)鎘離子初始濃度為3.94×10-7，1.58×10-6，3.15×10-6，6.3×10-6μmol·cm-3，經(jīng)過(guò)30 min的富集作用，溶液中的鎘離子濃度分別為3.16×10-8，1.28×10-7，2.59×10-7，5.33×10-7μmol·cm-3.說(shuō)明了根圍對(duì)鎘離子的凈化速率隨著鎘離子初始濃度的增加而增加，即鎘離子初始濃度與溶液中鎘離子的去除相關(guān).2000年，陳英旭等[29]采用蘿卜根系為材料對(duì)環(huán)境中重金屬鉛、鎘的富集修復(fù)作用進(jìn)行了研究，說(shuō)明了隨時(shí)間的推移，溶液中鎘濃度逐漸下降，即鎘的去除量升高，同時(shí)根系對(duì)溶液中鎘的去除也與初始溶液濃度有關(guān).

解析解和數(shù)值解的變化趨勢(shì)基本上是一致的，驗(yàn)證了解析方法的合理性.本文使用根圍鎘離子濃度的解析解來(lái)代替模型進(jìn)行計(jì)算和預(yù)測(cè)，描述了整個(gè)根系的鎘離子濃度.

4 結(jié)論

本文在文[28] 的工作基礎(chǔ)上，將絡(luò)合物的穩(wěn)定性分成3種情況，采用先求通解再求待定常數(shù)和近似方法，求解了具有自由邊界條件模型的穩(wěn)態(tài)解(即整個(gè)根圍區(qū)域鎘離子濃度的解析表達(dá)式).并通過(guò)數(shù)值解驗(yàn)證了解析解的正確性和求解穩(wěn)態(tài)解方法的合理性.

從解析結(jié)果和數(shù)值模擬發(fā)現(xiàn)右邊界初始濃度的降低時(shí)，Cd的絡(luò)合物會(huì)進(jìn)行解離并以自由離子的形式被根吸收，說(shuō)明了Cd2+初始濃度是控制根吸收鎘的有效性因素.進(jìn)而說(shuō)明了將初始濃度控制在一定范圍之內(nèi)，可有效降低植物對(duì)Cd的吸收和植物體內(nèi)的鎘含量.