999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

兩種吸附模型對阿特拉津在壤質砂土中的模擬效果分析

2020-02-22 03:43:18李曉宇任仲宇李芳春魏明海
農業環境科學學報 2020年1期
關鍵詞:模型

李曉宇,任仲宇,李芳春,井 琦*,魏明海

(1.北京工業大學建筑工程學院,北京 100124;2.水質科學與水環境恢復工程北京市重點實驗室,北京 100124;3.同方知網(北京)技術有限公司,北京 100192;4.生態環境部環境規劃院,北京 100012)

阿特拉津(Atrazine)是目前世界上廣泛使用的除草劑之一,也是土壤和地下水中最常檢測到的農藥之一[1],已在中國的松花江流域、遼河流域和長江流域中檢測到[2-4]。阿特拉津被認定為有潛在的致癌作用并可影響人體的內分泌系統,美國、歐共體和日本等均把它列入內分泌干擾物名單[5-6]。阿特拉津及其代謝產物可在水和土壤中存在數十年[7],其殘留會對飲用水源和生態環境構成威脅,或通過食物鏈傳遞和生物富集作用對人體健康造成危害。因此,研究阿特拉津在多孔介質中的遷移與變化規律具有十分重要的現實意義。

國內外學者對影響阿特拉津遷移的主要過程(吸附與解吸)及作用機理的相關研究已有大量報道[8-9]。Kovaios等[10]研究了土壤對阿特拉津的吸附作用,發現阿特拉津的吸附是可逆的,并觀察了吸附和解吸過程。針對土壤中阿特拉津的易遷移性,許多學者進行了室內土柱試驗和田間試驗,并在此基礎上通過模型模擬阿特拉津的遷移,從而評價其對地下水的污染風險[4]。Mao等[11]、Zheng等[12]利用兩點(兩區)非平衡吸附模型模擬阿特拉津在土壤中的遷移,并取得了較好的模擬效果。Prata等[13]采用三點非平衡吸附模型成功地擬合了淋溶條件下阿特拉津在巴西一種氧化土中的穿透曲線。毛萌等[14]建立了包含吸附項和降解項的對流彌散方程,模擬室內滴灌施藥條件下阿特拉津在土壤中的遷移。Van Genuchten等[15]發展了兩點模型的控制方程,在原方程基礎上增加了降解項。Gamerdinger等[16]應用帶有降解項的兩點模型模擬了室內定流速條件下阿特拉津在土柱中的遷移。

綜上所述,已有研究對模型的改進主要表現在改進模型源匯項(例如考慮吸附、降解項等)或將兩點非平衡吸附模型改進為三點非平衡吸附模型。單點非平衡吸附模型是兩點非平衡吸附模型的特例[15],Selim等[17]曾利用單點非平衡吸附模型模擬分析2,4-D除草劑在飽和及非飽和分層土壤中的遷移,但未討論該模型的模擬效果。吸附-解吸模型可用來模擬病毒、細菌和膠體等在多孔介質中的遷移[18],Yu等[19]還應用吸附-解吸模型模擬了納米零價鐵在飽和砂柱中的遷移,杜青青等[20]應用吸附-解吸模型模擬了氨氮在污染場地的遷移,但有關吸附-解吸模型模擬分析阿特拉津等農藥在多孔介質中的遷移報道甚少。筆者在利用單點非平衡吸附模型模擬阿特拉津遷移過程中,發現在淋洗階段誤差較大,進而考慮采用吸附-解吸模型對阿特拉津遷移進行模擬,通過改進吸附、解吸的源匯項表達,以期消減此誤差。

本文通過構建室內土柱阿特拉津動態遷移試驗,獲取模擬所需參數,利用HYDRUS-1D軟件,分別建立了基于單點非平衡吸附模型和吸附-解吸模型的運移模型,采用整體調參和分段調參兩種模擬方式,模擬阿特拉津在飽和壤質砂土中的遷移與分布規律,深入分析了兩種模型模擬結果差異背后的物理機理,并對源匯項的分項逐項進行計算,定量分析了兩種吸附模型對模擬結果的影響過程和影響程度,可為深刻認識阿特拉津的遷移特征及其環境風險評價提供科學依據。

1 材料與方法

1.1 室內土柱試驗

1.1.1 材料與儀器

試驗材料:阿特拉津(純度98%,國家標準物質中心)、氯化鉀(分析純,北京化工廠)、甲醇(色譜純、國藥集團化學試劑有限公司)。供試土壤采自北京市永定河河灘,采樣深度40 cm,土壤取回后自然風干,過2 mm篩備用。用簡易比重計法進行了該土壤的粒徑分析,并按國際制將土壤質地定名為壤質砂土。分析土壤的砂粒含量為86.653%,粉粒含量為11.795%,黏粒含量為1.553%,有機質含量為4.230%,干容重為1.656 g·cm-3,有效孔隙度為0.130,pH值為6.300。

試驗儀器:電子天平(JA2003,上海恒平科學儀器有限公司)、電導率儀(sensION5,哈希水質分析儀器有限公司)、高效液相色譜儀(Agilent1100,安捷倫科技有限公司)等。

模擬裝置如圖1所示,由進水裝置、土柱主體、出水裝置三部分組成。進水裝置由馬氏瓶、進水管和升降支架組成,馬氏瓶用于維持土柱定水頭淹水狀態,升降支架用于承載馬氏瓶,可調節高度。土柱主體由有機玻璃柱和充填介質組成,有機玻璃柱工作長度36 cm,內徑6 cm,柱頂設置溢流孔,柱中及柱底分別設置取樣孔。出水裝置由出水管和集液箱組成。

圖1 柱試驗裝置示意圖Figure 1 Schematic diagramof column experimental device

1.1.2 土柱遷移試驗

在土柱底部鋪設5 cm厚粒徑為0.5~2 cm的卵礫石,然后將風干過篩后的土壤按10 cm層厚分層填入土柱,裝填過程中對土壤層間打毛并在表層土壤上部鋪設3 cm厚粒徑為0.5~2 cm的卵礫石。土柱填裝完成后對土柱進行排氣飽和并調整進水出水的水位差,得到穩定滲流速度,在維持水位差恒定條件下將濃度為22μg·mL-1的阿特拉津滲濾液從頂端注入柱內,注入時長200 min,總注入水量1650 cm3,約12.5個土壤有效孔隙體積(12.5 PV),注入結束后瞬時用去離子水淋洗土柱,淋洗時間150 min,總淋洗水量1238 cm3(9.4 PV)。定時從取樣孔取樣,并測定滲流液中阿特拉津濃度。每次每個觀測孔取樣2個,測定結果取樣品的平均值。

1.1.3 測定方法

應用安捷倫高效液相色譜儀測定阿特拉津濃度。試驗采用C18色譜柱,設置流動相(甲醇∶水=70∶30)的流速為1 mL·min-1,在225 nm波長條件下用UVVIS檢測器測定。通過比對峰面積/濃度標線,計算得到不同時刻阿特拉津濃度。

1.2 模型的建立

HYDRUS軟件可用于模擬分析非飽和、飽和多孔介質中的水和溶質運移[18]。本文利用HYDRUS-1D軟件建立土壤水流和溶質運移數值模型,模擬阿特拉津溶液入滲進入飽和土柱后的遷移過程以及隨后淋洗條件下阿特拉津在飽和土柱中的遷移過程,設定時間單位為min,質量單位為μg,長度單位為cm。根據土柱試驗的運行情況,將模擬時長定為350 min。

1.2.1 水流模型的基本方程

一維平衡水流運動采用Richards方程來描述[18,21]:

式中:h為土壤壓力水頭,cm;θ為土壤體積含水率;t為時間,min;S為源匯項,min-1;α為水流方向與垂直方向夾角(即:α=0°代表垂向流動,90°代表水平流動,0°<α<90°代表傾斜流);x為空間位置,cm;K為非飽和水力傳導度函數,cm·min-1,如式(2):

式中:Kr為相對水力傳導度;Ks為飽和水力傳導度,cm·min-1。

本研究水流模型概化為均質各項同性飽和一維垂向穩定流,因此方程(1)可簡化為下式:

式中各參數意義同上。

1.2.2 水流模型的初始條件及邊界條件

模型初始條件為整體飽和狀態,即初始含水率為飽和含水率。

模型上邊界為供水邊界,下邊界為排水邊界,且均為定水頭邊界,水頭差為5.17 cm,通過插值得到土柱內的初始流場。

1.2.3 溶質遷移模型的基本方程

通過改進的對流彌散方程來模擬阿特拉津在飽和多孔介質中的遷移,該方程[22]用一個源匯項表示阿特拉津與固體基質的反應(吸附和解吸過程)。

式中:n為孔隙度;C為溶質濃度,μg·mL-1;S為多孔介質中溶質固相濃度,μg·g-1;D為彌散系數,cm2·min-1;v為達西速度,cm·min-1;x為空間位置,cm;ρ為多孔介質干容重,g·cm-3;t為時間,min。

根據不同的吸附機理,源匯項方程式有不同的表達形式。本文分別采用單點非平衡吸附模型和吸附-解吸模型對比模擬分析阿特拉津在土柱中的遷移過程。

式中:α是一級反應物質交換率系數(表示單位時間液相與固相間溶質質量轉化大小),min-1;Kd是Henry線性吸附系數,mL·g-1。

式中:ka為吸附系數,min-1;kb為解吸系數,min-1;ψ 為溶質吸附函數;Smax是固相最大吸附量,μg·g-1;θ為土壤體積含水率,在飽和狀態下,砂土中θ可近似認為等于孔隙度。

1.2.4 溶質遷移模型的初始條件及邊界條件

模型初始條件為初始時刻土柱各位置溶質濃度為0。

模型上邊界為濃度通量邊界,下邊界設定0梯度邊界條件。佩克萊數Pe表示對流與擴散在溶質遷移中所占比例[24]。Parlange等[25]認為,當Pe>4時,近似認為多孔介質流體中溶質濃度與某時段內按出流量計算的平均溶質濃度相等,此時邊界條件為0濃度梯度邊界條件。Pe計算公式如下:

式中:L為柱長,cm;u為實際平均速度,u=v/n,cm·min-1;D為彌散系數,cm2·min-1。通過計算本試驗佩克萊數Pe=49.860>4,符合Parlange等[25]研究的判定條件,即擴散作用相對于對流作用在溶質遷移中所占的比例較小,此時出流溶質的濃度梯度近似為0,故下邊界設定為零梯度邊界。

1.3 模型參數的獲取

本模擬中需要的參數有土壤水力特征參數和溶質特征參數,土壤特征參數主要指孔隙度、飽和滲透系數等;溶質特征參數主要指彌散系數、吸附系數等。

1.3.1 土壤水力特征參數

土壤水力特性參數的準確性對于反映實際土壤水分運動過程具有重要意義[26]。本研究采用常水頭法測定土壤介質飽和滲透系數Ks,容積排水法測定土壤介質孔隙度。測量結果土壤飽和滲透系數為0.340 cm·min-1,孔隙度為0.349。

1.3.2 彌散系數

彌散系數是表征可溶性物質通過滲透介質時彌散現象強弱的指標,可通過實驗利用投加示蹤劑測定其穿透曲線的方式獲取。本研究在室內土柱模擬裝置中,以Cl-作為示蹤劑,在飽和一維穩定流條件下淋濾KCl溶液,并在土柱出口端定時測定Cl-濃度,利用公式(8)[27]計算彌散系數。計算得到彌散系數為0.600 cm2·min-1。

彌散系數計算表達式:

式中:t0.16、t0.5、t0.84為柱中一個固定點x處相對濃度,分別為0.16、0.5、0.84時對應的時間,min。

Cl-濃度測試采用電導率法。配制不同濃度的KCl溶液,并測定其對應的電導率值,繪制電導率與KCl濃度的關系曲線,試驗時用電導率儀測定樣品電導率,利用電導率與KCl濃度關系將電導率換算為濃度。

1.3.3 等溫吸附系數

吸附系數的確定由室內吸附試驗獲得,本研究中阿特拉津在土壤中的吸附行為符合Henry等溫線性吸附,線性吸附系數Kd為0.372 mL·g-1;將土壤在濃度為22μg·mL-1的阿特拉津溶液中恒溫振蕩200 min,得到吸附量為6.450μg·g-1,以此作為本試驗條件下的固相最大吸附量Smax。

2 結果與討論

2.1 模型參數率定

在水力試驗基礎上確定水流模型參數,建立土柱一維水流穩定流模型,調整邊界水頭以改變土柱內的流場,將觀測孔的實測測壓管水頭值和計算水頭值進行對比,反復進行三次,取誤差較小的情形所對應的參數為模型水流參數。模型中涉及的土壤水力特征參數如飽和滲透系數為0.340 cm·min-1,孔隙度為0.349。

水質模型中源匯項涉及的參數有彌散系數D、一級反應物質交換率系數α、線性吸附常數kd、吸附系數ka、解吸系數kb、固相最大吸附量Smax。彌散系數D、線性吸附常數kd、固相最大吸附量Smax通過實驗得到,并在模型中適當調整,最終得到這些參數值:D為0.600 cm2·min-1,kd為0.372 mL·g-1,Smax為6.450 μg·g-1;針對以上兩種吸附模型中不易通過試驗直接測得的參數如一級反應物質交換率系數α、吸附系數ka、解吸系數kb進行調參,得到綜合反映整個試驗過程的“平均的”α、ka和kb,分別為 0.200、0.179 min-1和 0.126×10-5min-1。

實際情況下,隨注入、淋洗階段時間的增長,模型中的一級反應物質交換率系數、吸附和解吸系數會發生變化。為了探究實際情況下各反應系數的變化過程及對模擬效果的影響,本研究以觀測孔1和觀測孔2為調參檢驗點,結合阿特拉津穿透曲線形態,將模擬按照穿透曲線的“上升段”“峰值段”“下降段”分段進行調參模擬,得到各階段參數值如表1所示。

2.2 模擬結果

圖2、圖3分別為取樣孔1和取樣孔2阿特拉津濃度實測值與整體調參(未分段調參)擬合的單點非平衡吸附模型、吸附-解吸模型模擬計算結果對比圖。圖4、圖5分別為取樣孔1和取樣孔2阿特拉津濃度實測值與分段調參擬合的單點非平衡吸附模型、吸附-解吸模型模擬計算結果對比圖。

單點非平衡吸附模型在阿特拉津注入過程中擬合情況較好,各觀測點濃度模擬值和實測值較為一致,但在淋洗過程中,模擬結果相比于試驗結果有較為明顯的誤差,不能較好地重現整個過程(見圖2和圖3);吸附-解吸模型在整個過程中模擬值與實測值趨勢更為一致,且誤差較小,可以較好地描述阿特拉津在土壤中的遷移過程。通過觀察圖4、圖5并與圖2、圖3對比可知,分段調參擬合的單點非平衡吸附模型的擬合度在注入階段提升不明顯,在淋洗階段有顯著提升;分段調參擬合的吸附-解吸模型的擬合度在注入階段和淋洗階段均提升很小。

2.3 模型結果檢驗

選用Nash-Suttcliffe模擬效果系數(NSC)和均方根誤差(RMSE)評價模型計算結果,計算公式分別為:

表1 各階段一級反應物質交換率系數、吸附和解吸系數擬合值Table 1 Fitting values of exchange rate coefficient of first order reactive substance,adsorption and desorption coefficients at each stage

圖2 觀測孔1整體調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 2 Integral parameter fitting of observation hole 1:The simulation result compartion chart of one-site non-equilibriumadsorption model and adsorption-desorption model

圖3 觀測孔2整體調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 3 Integral parameter fitting of observation hole 2:The simulation result compartion chart of one-site non-equilibriumadsorption model and adsorption-desorption model

圖4 觀測孔1分段調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 4 Segmental parameter fitting of observation hole 1:The simulation result compartion chart of one-site non-equilibriumadsorption model and adsorption-desorption model

式中:Xobs為實測值,Xmodel為模型計算值,Xˉobs為實測值的算術平均值。當模型計算值和實際監測值相等時,NSC=1,模擬效果最好;通常NSC在0~1之間,NSC越大,說明計算值與觀測值匹配程度越好,一般當NSC>0.75時,表示模擬結果可以接受。

式中:N表示實測值個數,Xobs,i和Xmodel,i為試驗觀測值與模型計算值。RMSE越小,說明計算值與觀測值匹配程度越好。

分整體調參擬合和分段調參擬合兩種情形,按照阿特拉津注入階段、淋洗階段和注入-淋洗全階段,計算單點非平衡吸附模型和吸附-解吸模型的NSC和RMSE值,如表2所示。

由表2可知,整體調參擬合情形下,單點非平衡吸附模型模擬注入階段(吸附階段)NSCs與RMSEs分別為0.909~0.922和2.752~3.167,吻合較好;而對于淋洗階段,NSCs為負,說明擬合度很低;單點非平衡吸附模型模擬阿特拉津注入-淋洗的全過程的NSCs與RMSEs分別為0.138~0.593和6.180~8.395。分段調參擬合后,注入階段的NSCs提高到0.919~0.941,RMSEs減小到2.405~2.986,淋洗階段的NSCs提高到0.927~0.940,RMSEs減小到2.309~2.036,注入-淋洗全階段的NSCs提高到 0.921~0.941,RMSEs減小到 2.204~2.717。說明用分段擬合的單相非平衡吸附模型模擬阿特拉津注入、淋洗和注入-淋洗的全過程均是適宜的。

同樣通過觀察吸附-解吸模型的孔1和孔2的NSC和RMSE值可知,無論整體調參模擬和分段調參模擬,吸附-解吸模型模擬阿特拉津注入、淋洗和注入-淋洗全過程均是較為適宜的。分段擬合的吸附-解吸模型模擬阿特拉津注入、注入-淋洗全過程的模擬效果系數有微弱提高。

2.4 兩種模型機理對比分析

圖5 觀測孔2分段調參擬合:單點非平衡吸附模型與吸附-解吸模型模擬結果對比圖Figure 5 Segmental parameter fitting of observation hole 2:The simulation result compartion chart of one-site non-equilibrium adsorption model and adsorption-desorption model

單點非平衡吸附模型認為固相中的吸附量S與一級反應物質交換率系數α和線性吸附常數與溶質液相濃度的乘積KdC有關。S和KdC隨著吸附過程的發生不斷發生變化。根據單點非平衡吸附模型,可將阿特拉津在土柱中的遷移、吸附過程描述如下:初始時刻柱中各位置阿特拉津固相及液相濃度均為0,阿特拉津隨水流遷移進入土柱,土柱中液相濃度升高,同時發生吸附作用,固相濃度隨之增加,達到吸附平衡后停止吸附。隨后,注入去離子水,吸附在固相中的阿特拉津逐漸解吸到水中隨水流出。由圖2、圖3可知,阿特拉津連續注入階段(吸附階段)模擬值與實測值偏差較小,但在去離子水注入階段(淋洗階段)則偏差很大。這是因為整個注入-淋洗階段的一級反應物質交換率系數α為同一常數,淋洗階段由于采用了與注入階段相同的α,導致淋洗階段實測值與模擬值差異較大。分段擬合后,單點非平衡吸附模型由于在注入階段和淋洗階段采用了與對應階段相匹配的α,因此注入階段和淋洗階段的擬合度均較高。

表2 各階段模型模擬效果系數(NSC)和均方根誤差(RMSE)值Table 2 Simulation effect coefficient(NSC)and root-mean-squared error(RMSE)values of models at each stage

為了進一步說明兩種吸附模型的結果差異,本文對兩種吸附模型在不同階段的源匯項值進行了準確計算,分別列示于表3與表4。

表3列出了觀測孔1單點非平衡吸附模型源匯項計算值。由表3可知:隨著阿特拉津溶液的持續注入,模型計算得到的液相濃度和固相濃度逐漸增大,均為正值且逐漸減小,到阿特拉津注為0.003 μg·cm-3·min-1,此時固相濃度仍有較小的正向增長;在淋洗階段,隨著去離子水的持續注入,模型計算得到的液相濃度和固相濃變為負值(負值表示解吸過程)且絕對值與注入阿特拉津時的源匯項值相近,即吸附到固相中的阿特拉津逐漸解吸到溶液中,解吸速率與注入阿特拉津溶液時的吸附速率大致相當,因此采用單點非平衡吸附模型在淋洗階段模擬值遠超出實測值。

吸附-解吸模型認為吸附和解吸速率不同,將模型分為吸附項和解吸項分別描述,用ka和kb代表吸附和解吸的速率。在吸附項中,吸附量與溶質液相濃度C和最大固相濃度Smax有關,當吸附量達到Smax時停止吸附。根據吸附-解吸模型,可將阿特拉津在土柱中遷移、吸附過程描述如下:初始時刻柱中各位置阿特拉津液相和固相濃度均為0,隨著阿特拉津注入土柱,土柱中液相濃度升高,同時發生吸附作用,固相濃度隨之增加,隨著土柱中各位置吸附量達到最大值后,吸附不再發生,整個過程中阿特拉津也會發生少量的解吸;對于淋洗階段,同樣存在著吸附和解吸作用,解吸系數kb在數值上很小,說明吸附到壤質砂土中的阿特拉津解吸速率很慢,這也是吸附-解吸模型計算值遠低于單點非平衡模型計算值的原因。

表4列出了觀測孔1吸附-解吸模型源匯項計算值。由表4可知:隨著阿特拉津溶液的持續注入,模型計算的液相濃度和固相濃度逐漸增大,分別達到注入濃度和固相最大吸附量后維持不變,源匯項隨著固相濃度的持續增大逐漸減小,達到最大固相吸附量后變成負值且維持恒定,此時的源匯項為-1.388×10-5μg·cm-3·min-1,可知固相濃度雖然在逐漸降低,但下降速率近似為0;在淋洗階段,隨著去離子水的持續注入,模型計算的液相濃度逐漸降低,源匯項為-1.388×10-5μg·cm-3·min-1,表明吸附在固相上的阿特拉津解吸速率極慢。通過觀察阿特拉津液相濃度實測穿透曲線發現沒有拖尾現象,也說明了其實際解吸速率很小。對比表3、表4的源匯項值計算結果,在解吸階段單點非平衡吸附模型的源匯項計算值比吸附-解吸模型的源匯項計算值高出4個數量級,這是引起兩種模型計算結果差異較大的原因。

表3 觀測孔1單點非平衡吸附模型源匯項計算值Table 3 Calculation values of source and sink terms of the onesite non-equilibriumadsorption model of the observation hole 1

表4 觀測孔1吸附-解吸模型源匯項計算值Table 4 Calculation values of source and sink terms of the adsorption-desorption model of the observation hole 1

3 結論

(1)整體調參下,單點非平衡吸附模型在連續注入阿特拉津的條件下,飽和壤質砂土中阿特拉津遷移模擬的NSCs系數為 0.909~0.922,RMSEs為 2.752~3.167,采用單點非平衡吸附模型是適合的;淋洗條件下阿特拉津遷移模擬的NSCs為負值,說明單點非平衡吸附模型不適用于淋洗條件下阿特拉津遷移模擬。

(2)整體調參下,吸附-解吸模型模擬連續注入、淋洗過程中阿特拉津在飽和壤質砂土中的遷移模擬的NSCs為 0.901~0.954,RMSEs為 2.037~3.289,表明吸附-解吸模型模擬連續注入、淋洗過程中阿特拉津遷移均是適合的。

(3)分段調參可提高單點非平衡吸附模型模擬阿特拉津連續注入、淋洗過程中在飽和壤質砂土中遷移的擬合度。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 九九香蕉视频| 色综合久久综合网| 亚洲区一区| 狠狠v日韩v欧美v| 亚洲天堂网在线观看视频| 国模私拍一区二区三区| 精品国产免费观看| 天堂va亚洲va欧美va国产| 国产91视频免费观看| 国产色偷丝袜婷婷无码麻豆制服| 性色一区| 污污网站在线观看| 小13箩利洗澡无码视频免费网站| 婷婷伊人久久| 欧美特黄一免在线观看| 在线国产欧美| 91麻豆久久久| 亚洲综合专区| 亚洲三级电影在线播放| 99精品在线看| 扒开粉嫩的小缝隙喷白浆视频| 免费AV在线播放观看18禁强制| 在线无码九区| 亚洲av片在线免费观看| 欧美国产日产一区二区| 91久久精品日日躁夜夜躁欧美| 99久久精品国产综合婷婷| 成人午夜久久| 国产美女一级毛片| 成年人久久黄色网站| 欧美日韩中文字幕二区三区| 波多野结衣AV无码久久一区| 国产丝袜啪啪| 美女毛片在线| 99re在线免费视频| 亚洲最新地址| 日本三区视频| 国产主播喷水| 国产免费福利网站| 大陆精大陆国产国语精品1024| 亚洲中文字幕无码mv| 99re经典视频在线| 亚洲欧美自拍一区| 欧美成在线视频| 人妻免费无码不卡视频| 噜噜噜久久| 久操线在视频在线观看| 久久久波多野结衣av一区二区| 91视频精品| 国产午夜精品一区二区三区软件| 992Tv视频国产精品| 青青操国产视频| 伊人久久精品无码麻豆精品| 国产三级成人| 中文字幕有乳无码| 熟女日韩精品2区| 9久久伊人精品综合| 国产女人水多毛片18| 亚洲成aⅴ人片在线影院八| 国产精品永久不卡免费视频| 欧美日韩成人在线观看| 香蕉久久永久视频| 国产成人免费高清AⅤ| 日韩不卡免费视频| 国产精品视频999| 扒开粉嫩的小缝隙喷白浆视频| 综合网久久| 国产精品理论片| 国产精品视频观看裸模| 一区二区三区四区在线| 亚洲综合色吧| 超清无码一区二区三区| 一区二区三区四区在线| www亚洲天堂| 亚洲精品欧美日本中文字幕 | 亚洲AⅤ波多系列中文字幕| 国产色爱av资源综合区| 91丨九色丨首页在线播放| 日本成人福利视频| 亚洲欧美日韩综合二区三区| 女人18毛片久久| 日韩在线视频网|