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

金屬礦區重金屬遷移對水體影響的數值模擬

2015-02-22 07:01:50軒曉博李一平王仕彬
水資源保護 2015年2期
關鍵詞:水質模型

軒曉博 逄 勇,李一平,王仕彬,王 雪

(1.河海大學環境學院,江蘇 南京 210098;2.黃河勘測規劃設計有限公司,河南 鄭州 450003;3.河海大學淺水湖泊綜合治理與資源開發教育部重點實驗室,江蘇 南京 210098;4.福建省閩東南地質大隊,福建 泉州 362021)

金屬礦區重金屬遷移對水體影響的數值模擬

軒曉博1,2逄 勇1,3,李一平1,3,王仕彬4,王 雪1

(1.河海大學環境學院,江蘇 南京 210098;2.黃河勘測規劃設計有限公司,河南 鄭州 450003;3.河海大學淺水湖泊綜合治理與資源開發教育部重點實驗室,江蘇 南京 210098;4.福建省閩東南地質大隊,福建 泉州 362021)

為了解礦區土壤中重金屬對周圍地表水體的影響,以閩江支流浐溪河彭村水庫為例,基于MIKE21模型平臺構建二維重金屬預測模型,根據重金屬在水體和底泥中的轉移擴散、懸浮沉積及吸附解吸原理,采用土壤淋溶及浸泡試驗所獲取的庫區土壤重金屬含量等數據,模擬不同入庫流量下庫區Zn、Cd、Pb等3種重金屬污染分布規律。結果表明:小流量條件下,Zn的濃度為庫區中央最小,兩岸次之,靠近礦區區域最大;大流量下,庫區中央大于兩岸,下游大于上游。Pb的濃度為庫區中央大于兩岸,下游大于上游。Cd的濃度為庫區中央大于兩岸,下游大于上游。水庫蓄水后,重金屬濃度基本滿足GB 3838—2002《地表水環境質量標準》Ⅱ類水質標準。

金屬礦區;重金屬遷移; MIKE 21; 水質; 數值模擬

重金屬在土壤中的遷移轉化受金屬的化學特性、土壤的物理特性、生物特性和環境條件等因素影響,其遷移轉化過程分為物理遷移、化學遷移、物理化學遷移和生物遷移,其遷移轉化形式復雜多樣,是多種形式的錯綜結合[1-3]。

通常狀況下,土壤中的重金屬基本處于相對穩定的狀態,且難以遷出土體[4],當土壤物理特性以及環境條件發生變化時,重金屬就會發生遷移轉化。對于水庫,當其蓄水后,由于水體的浸泡作用,土壤中的重金屬除部分發生水解生成難溶物外,剩余重金屬多會析出,并且能以離子態在底泥以及水體中轉移擴散,從而有可能影響水體質量。多年來因重金屬污染物進入河流、水庫、湖泊等導致水環境污染釀成了不少悲劇,因此它是水庫(湖泊)、河流水環境評價的重要內容之一。

隨著計算機水平的快速發展,采用先進的計算軟件平臺來構建符合實際的重金屬預測模型,將成為一種新的預測手段。常見的水質預測模型有美國陸軍工程兵團水道實驗站的模型SELECT、CE-QUAL-R1和CE-QUAL-W2[4]等模型,美國國家環保署(EPA)開發的WASP[5]、EFDC[6]模型,美國農業部(USDA)開發的SWAT模型[7]等,以及丹麥水力學研究所(DHI)研制開發的MIKE軟件模型,這些模型在水庫水質模擬預測過程中都發揮了極其重要的作用。

以上各類模型軟件都有各自的側重點和局限性,并且在模擬水流流場以及水質分布方面存在一定的差異。如SELECT模型,只能預測已知密度水庫的垂向水質分布,且對一些參數進行了物質化處理;CE-QUAL-R1模型與SELECT模型一樣,也是確定水質參數垂向位置隨時間變化的值,它將水庫概化為垂直方向一系列的水平層,且每一次的熱能與雜質都均勻分布;另外該模型對氮、磷等營養元素的遷移轉化過程能較好地模擬,但對重金屬等應用較少;SWAT模型主要針對農業面源污染進行模擬,等等。近年來經過國內外許多工程實際運用,認為MIKE模型在預測水庫水質方面具備明顯優勢,且能模擬接近實際的水流條件,尤其對難溶解的污染物如重金屬的遷移擴散等能進行較好地模擬。MIKE模型分為MIKE11、MIKE21以及MIKE Basin等,MIKE21是應用較為廣泛的平面二維數學模型,該模型可用于模擬河流、湖泊、河口、海灣的水流、波浪、泥沙及環境場。MIKE21模型包括水動力學模塊和水質模塊兩大部分,可定義多種類型的水邊界條件,如流量、水位或流速等,能進行、干、濕節點和干、濕單元的設置,能較方便地進行灘地水流的模擬。其中水動力學模塊是MIKE 21軟件最核心的基礎模塊,可以模擬由于各種作用力的作用而產生的水位及水流變化,它包括了廣泛的水力現象,可用于任何忽略分層的二維自由表面流的模擬,該模塊為環境模擬提供了水動力學的計算基礎。

目前該軟件在國內的應用發展很快,并在一些大型工程中被廣泛應用。如渾河流域采用了MIKE11模型進行水質預測,較好地分析了不同發展模式的合理性[8];上海世博園白蓮涇區域水量水質預測模型依據MIKE11進行構建,并模擬了“6·29大雨”白蓮涇河區域的河道水量水質情況[9];劉偉等[10]利用MIKE11模型計算了松花江干流松江大橋至通氣河口段的NH3-N納污能力,對不同水文條件下的流量、流速進行實時模擬,避免了因排污口概化造成的計算偏差;此外,MIKE模型還應用于長江口綜合治理工程、太湖富營養化治理工程、香港新機場工程建設等[11]。

通過查詢相關課題研究成果及設計報告,并檢索了相關文獻,目前在國內對MIKE軟件的使用,多數集中在天然河流、湖泊以及人工水庫的模擬計算,但對于土壤中重金屬遷移對水庫的影響研究案例較少。本文以閩江支流浐溪河彭村水庫為例,基于MIKE21模型平臺構建二維重金屬預測模型,采用土壤淋溶及浸泡試驗所獲取的庫區土壤重金屬含量等數據,模擬不同入庫流量下庫區重金屬污染分布規律,參照GB 3838—2002《地表水環境質量標準》對水庫水質進行影響評價。

1 數據獲取與分析

1.1 數據來源

彭村水庫地處福建省德化縣,位于閩江下游的支流——大樟溪上游段,水庫正常蓄水位為642.0 m,總庫容7 871萬m3,具有多年調節性能。水庫壩址控制流域面積144.50 km2,占浐溪全流域的14.7%。該水庫主要功能以供水、防洪為主,兼顧發電等,具有較好的調節性能,可對流域來水“蓄豐補枯”,對下游起到一定的補償作用。彭村水庫庫區上游段河道右岸分布有未開發的鉛鋅礦,礦體分布高程725 m以上,高于正常蓄水位約83 m(圖1)。

圖1 彭村水庫工程位置及與礦區位置關系

本次評價所采取的源數據以及相關參數根據現場試驗獲取,試驗內容包括淋溶試驗、浸泡試驗及試坑滲水試驗,其中浸泡試驗包括水泡試驗和酸泡試驗。現場采樣根據不同的高程按等高距10 m結合100 m×20 m網度進行取樣控制,自下游向上游垂直于主河道在庫區左岸布置16條測線,測線間距100 m,共采取50個樣品,樣品測試因子為Pb、Zn、Cd。

本次試驗,采用不同的pH值分別進行試驗,試驗結果隨著pH的不同而有所不同。

土壤淋溶試驗。測定Zn、Cd、Pb的質量濃度分別為: 2.61~16.70 mg/m3、0.01~0.91 mg/m3、0.49~3.97 mg/m3。

土壤水泡試驗。采用蒸餾水進行試驗,浸取結果顯示:Zn、Cd、Pb的質量濃度分別為30.44~212.30 mg/m3、0.40~1.21 mg/m3、1.18~21.62 mg/m3;Zn、Cd達到GB 3838—2002《地表水環境質量標準》Ⅰ~Ⅱ類標準,Pb符合Ⅰ類標準占40%,其余符合Ⅱ~Ⅲ類標準。

土壤酸泡試驗。采用稀硫酸和硝酸進行試驗,酸泡試驗pH值選取5.2、5.5、5.8、6.0,之所以選取這4個值,是根據彭村水庫下游龍門灘水庫2012年1月3日以來的逐日監測報表,于當年7月25日出現最低pH值5.2確定的。浸取結果顯示:Zn、Cd、Pb的質量濃度分別為 31.51~153.80 mg/m3、0.18~6.71 mg/m3、0.07~2.42 mg/m3;均小于GB 5085.3—2007《危險廢物鑒別標準-浸出毒性鑒別》的標準中所列的濃度限值,該土壤不具有浸出毒性特征的危險廢物。

試坑滲水試驗。試驗得出該區域土壤滲透系數為0.099 4~0.305 7 m/d。

1.2 數據可靠性分析

本次預測源數據通過現場布點、進行一系列的浸溶試驗得到,布點時考慮了水庫水位上漲與消落的過程,取不同高程點進行垂向布置,并根據礦區與水庫的位置關系,依據“近密遠疏”的原則,進行了縱向與橫向的布置,試驗的每一個過程嚴格按照GB/T 14158—93《區域水文地質工程地質環境地質綜合勘查規范》(1∶50 000)試驗規范要求進行,因此,所取得的數據是可靠的。

2 模型構建

所構建的模型包括水動力學模塊和水質模塊兩個部分,水動力學模塊主要模擬水流條件,水質模塊主要模擬污染物的遷移轉化過程。

2.1 水動力學模塊

MIKE 21 FM水動力部分的輸入數據分為以下幾個部分:①計算域和相關時間參數,包括網格地形及時間設置;②校準要素,包括底床阻力、渦黏系數和風摩擦阻力系數;③初始條件,如水面高程;④邊界條件,分為開邊界條件和閉邊界條件,這里主要為開邊界條件,包括流量、水位等;⑤其他驅動力,包括風速、風向、源匯項和波浪輻射應力等。

該模塊初始條件根據工程設計的有關技術參數設定,驅動力則由氣候統計資料以及現場調查所得。本次研究選擇預測范圍內3個邊界的流量以及水位時間序列作為開邊界條件。由于一般水庫豐水時段水質較差,因此計算時采用典型年份(P=10%的1974—1975年的水文年份)豐水期的日流量及水位資料,模型參數率定采用水質與水文同步監測的數據。

該模塊網格概化見圖2,地形內插結果參見圖3,模型基本設置如下:網格數量為18 981個,最小網格尺度為15 m, 時間步長為300 s, 渦黏系數為0.4,曼寧數為32,風摩擦阻力系數為0.001 3。

圖2 彭村水庫水動力學模塊網格概化

2.2 水質模型

2.2.1 模型基本方程

流域水質模型主要模擬以下污染因子:溶解態重金屬、吸附態重金屬、底泥孔隙水中溶解態重金屬、底泥中吸附態重金屬、懸浮顆粒物及沉積物量。重金屬在水中和底泥中轉移擴散、懸浮沉積及吸附解吸的原理如圖4所示,據此原理推算出重金屬模型見式(1)。

(1)

圖3 彭村水庫水動力學模塊地形內插結果

圖4 重金屬在水體和底泥中的遷移轉化示意圖

其中Aadss=kwKdSHMXSS

Ddess=kwXHM

式中:XHM為水體中吸附態重金屬質量濃度,g/m3;Asdss為吸附過程中溶解態重金屬濃度,g/(m3·d);Ddess為解吸過程中溶解態重金屬濃度,g/(m3·d);Ssev為沉積過程吸附態重金屬濃度,g/(m3·d);Rresv為起懸過程吸附態重金屬濃度,g/(m3·d);t為計算時段,s;kw為水體中的解吸率,d-1;kd為重金屬在顆粒態和水之間的分配系數;SHM為水體中溶解態重金屬質量濃度,g/m3;XSS為水體中懸浮顆粒質量濃度,g/m3;Vvsm為SS的沉降速度,m/d;Z為計算層厚度,m;Rresrat為SS的再懸浮率,g/(m2·d);XHMS為底泥里吸附態重金屬濃度,g/m2;XSED為沉降通量,g/m2。

由于吸附態重金屬占總量的95%以上,所以一般不考慮水體中溶解態重金屬,則上式有時也可簡化為

2.2.2 模型邊界條件

模型邊界條件包括水動力邊界條件以及水質邊界條件。根據MIKE 21 模型設置要求,水質模型邊界主要為開邊界,與水動力學模塊開邊界一樣,此次研究中共有3個水質開邊界,包括1個上游開邊界和2個下游開邊界。從偏保守角度考慮,本次模擬時的邊界濃度取最不利情況下的數值。風場邊界的確定根據德化縣30年長系列氣象觀測資料,考慮概化為空間均勻的風場,利用風速和風向因子,數據格式為dfs0文件格式。

2.2.3 模型率定和驗證

對浐溪流域水體監測斷面進行了污染物計算值與實測值驗證,結果見圖5。采用2012年2月浐溪流域的實測水質資料進行模型參數率定。

圖5 浐溪河監測斷面污染物質量濃度計算值與實測值驗證

水動力模型參數率定過程中,針對模型的糙率進行調整,通過對水位和流量的率定確定模型糙率值為0.03。模型驗證過程中,水位和流量誤差均在25%以內,表明參數選取較為合理。水質模型參數主要利用ECO Lab中的水質模塊和重金屬模塊進行率定驗證,通過反復的試算與比較,最終確定參數值。模型中部分參數的變化對整個模型的影響相對較小,故采用默認值。最終確定的主要參數見表1。

3 結果及分析

3.1 模型預測結果

采用已建模型對蓄水后礦區附近水庫水質進行模擬計算,分大、小流量不同條件進行計算,大流量采取典型年份豐水期的平均值,小流量采取典型年份枯水期的平均值。通過計算得出不同條件下的污染分布圖,結果見圖6。

由圖6可知不同流量條件下不同重金屬污染分布規律如下:

a. Zn的質量濃度分布規律。在小流量情況下,庫區內靠近礦區一岸的狹長區域污染濃度較大,大致介于5.281 2×10-3~5.281 6×10-3mg/L,在水庫中央以及下游濃度最小,基本介于5.276 8×10-3~5.277 6×10-3mg/L,其他靠近河岸邊的區域濃度居中;在大流量情況下,整體的濃度要較小流量條件下的大,庫區內的分布規律與小流量條件下的恰好相反,在庫區中央地帶濃度較大,基本處于5.8×10-3~5.83×10-3mg/L,靠近兩岸處較小,為5.78×10-3~5.81×10-3mg/L。不管是小流量條件下,還是大流量條件下,庫區蓄水后,Zn的質量濃度整體滿足GB 3838—2002《地表水環境質量標準》Ⅱ類標準要求。

表1 水質模型參數取值范圍

圖6 彭村水庫中重金屬在不同流量條件下的污染分布

b. Pb的質量濃度分布規律。在小流量條件下,庫區中央的濃度較大,基本在9.782 8×10-3~9.783 3×10-3mg/L之間,越靠近兩岸越小,基本介于9.781 2×10-3~9.779 8×10-3mg/L之間;在大流量條件下,庫區整體分布規律與小流量基本一致,中央區域大,兩岸較小,與小流量條件下分布不同的是,靠近兩岸的濃度更小,基本介于9.74×10-3~9.742×10-3mg/L之間。可以看出,兩種流量條件下,下游濃度大于上游濃度,整體都滿足相關水質標準要求。

c. Cd的質量濃度分布規律。在小流量條件下,庫區中央濃度較大,基本介于8.118 5×10-4~8.119×10-4mg/L,兩岸較小,且靠近礦區一側的較對岸濃度稍大,基本介于8.115 5×10-4~8.116 2×10-4mg/L,且中央區域為下游大于上游,兩岸區域為下游和上游小于礦區附近;在大流量條件下,整體分布規律與小流量條件下基本一致,中央區域濃度較小流量條件下小,除靠近礦區一側區域,兩岸區域濃度比小流量條件下稍大,介于8.106×10-4~8.11×10-4mg/L。可以看出,在彭村水庫蓄水后,水庫的Cd質量濃度均可以滿足水質標準要求。

綜上可知,在不同流量條件下,礦區重金屬污染范圍均較小,按GB 3838—2002《地表水環境質量標準》進行評價,均滿足Ⅱ類標準要求。

3.2 預測結果合理性分析

本模型是在MIKE21平臺基礎上,充分考慮污染因子的釋放及轉移規律而建立起來。首先按照2012年2月份水質實測資料進行模型率定和驗證,并參照太湖流域水質預測模型中有關對重金屬預測的參數取值來確定本次預測的參數值。同時在模型計算條件設置中充分考慮了氧化還原、水體紊動強度等方面的因素,使得模擬計算的過程更趨近于實際,因此預測結果基本合理。

4 結 語

筆者基于MIKE21模型,結合該工程區的特點以及本水庫特性,構建了重金屬二維預測模型,進一步梳理了重金屬在水體、土壤、底泥過程中的遷移轉化過程,對未開發礦區水庫蓄水所產生的重金屬污染影響狀況進行了模擬計算。計算結果表明:水庫蓄水后,Zn濃度分布規律為:小流量條件下,庫區中央最小,兩岸次之,靠近礦區區域最大,大流量下,中央大于兩岸,下游大于上游;Pb濃度分布規律為:中央大于兩岸,下游大于上游;Cd濃度分布規律為:中央大于兩岸,下游大于上游。整體上重金屬濃度均滿足GB 3838—2002《地表水環境質量標準》Ⅱ類標準要求。由此可知,該未開發的金屬礦,在水庫蓄水后不會對水庫水質產生影響。通過對污染分布規律的研究,為水庫設計調度提供優化建議,并為類似水庫水質預測提供參考。

[1] 孫鐵珩.污染生態學[M].北京: 科學出版社, 2002: 18-24.

[2] 鮑桐, 廉梅花,孫麗娜,等.重金屬污染土壤植物修復研究進展[J].生態環境,2008,l7(2): 858-865.(BAO Tong,LIAN Meihua, SUN Lina,et al.Research progress on the phytoremediation of soil contaminated by heavy metals[J].Ecology and Environment,2008,17(2) 2008:858-865.(in Chinese))

[3] 李麗,王富華,王旭,等.韶關土壤重金屬污染狀況[J].農業環境與發展, 2010, 27(1): 24-25.(LI Li,WANG Fuhua,WANG Xu,et al.Heavy metal contamination of soils in Shaoguan[J].Agricultural Environment and Development,2010,27(1):24-25(in Chinese))

[4] 王濟,王世杰.土壤中重金屬環境污染元素的來源及作物效應[J].貴州師范大學學報:自然科學版,2005,23(2):111-113.(WANG Ji,WANG ShiJie.The sources and crops effect of heavy metal elements of contamination in soil[J].Journal of Guizhou Normal University:Natural Sciences,2005,23(.2):111-113.(in Chinese))

[5] 孫穎,陳肇和,范曉娜,等.河流及水庫水質模型與通用軟件綜述[J].水資源保護,2001(2):9-10.(SUN Ying,CHEN Zhaohe,FAN Xiaona,et al.A review of river and reservoir water quality models and general software[J].Water Resources Protection, 2001(2):9-10(in Chinese))

[6] 劉夏明,李俊清,豆小敏,等.EFDC 模型在河口水環境模擬中的應用及進展[J].環境科學與技術,2011,34(6G):136-140.(LIU Xiaming,LI Junqing,DOU Xiaomin,et al.The application and advance of environmental fluid dynamics code (EFDC) in estuarine water environment[J].Environmental Science & Technology,2011,34(6G):136-140. (in Chinese))

[7] 賴格英,吳敦銀,鐘業喜,等.SWAT模型的開發及應用進展[J].河海大學學報:自然科學版,2012, 40(3):243-247.(LAI Geying,WU Dunyin,ZHONG Yexi,et al.Progress in development and applications of SWAT model[J].Journal of Hohai University:Natural Sciences,2012,40(3):243-247.(in Chinese))

[8] 常旭,王黎,李芬,等.MIKE11模型在渾河流域水質預測中的應用[J].水電能源科學,2013,31(6):60-61.(CHANG Xu,WANG Li,LI Fen,et al,Application of MIKE11 model in water quality prediction of Hunhe River Basin[J].Water Resources and Power,2013,31(6):60-61.(in Chinese))

[9] 黃琳煜,聶秋月,周全,等.基于MIKE11的白蓮涇區域水量水質模型研究[J].水電能源科學,2011,29(8):21-24.(HUANG Linyu,NIE Qiuyue,ZHOU Quan,et al.Study of water quantity and water quality model of Bailianjing Region based on MIKE11[J].Water Resources and Power,2011,29(8):21-24.(in Chinese))

[10] 劉偉,劉洪超,徐海巖.基于MIKE11模型計算河流水功能區納污能力方法[J].東北水利水電,2009(8):69-70.(LIU Wei,LIU Hongchao,XU Haiyan.Calculation method of water environment capacity for water function area based on MIKE 11 model[J].Water Resources & Hydropower of Northeast,2009(8):69-70.(in Chinese))

[11] 王領元.應用mike對河流一、二維的數值模擬[D].大連:大連理工大學,2007:33-40.

Numerical simulation of influence of heavy metal migration on water in metallic mining areas

XUAN Xiaobo1,2, PANG Yong1,3, LI Yiping1,3, WANG Shibin4, WANG Xue1

(1.CollegeofEnvironment,HohaiUniversity,Nanjing210098,China; 2.YellowRiverEngineeringConsultingCo.Ltd.,Zhengzhou450003,China; 3.KeyLaboratoryofIntegratedRegulationandResourcesDevelopmentonShallowLakes,MinistryofEducation,HohaiUniversity,Nanjing210098,China;4.SoutheastFujianGeologicalBrigadeofFujianProvince,Quanzhou362021,China)

In order to make it clear the effects of heavy metals in the soil of mining area to its surface water around, taking the Pengcun Reservoir on the Chanxi River which is a branch of Minjiang River as an example, a two-dimension heavy metal prediction model was set up based on the MIKE21 model. According to the principles of migration-diffusion, suspension-deposition and adsorption-desorption of heavy metals in water and sediment, and by using the concentration data of heavy metals in the reservoir area obtained from the experiment of soil leaching and immersion, the contamination distribution of heavy metals such as Zn, Cd and Pb under different reservoir inflow was simulated. The results show that the distribution order from small to big of Zn concentration under a low flow condition is the central zone of reservoir,river banks,the places close to mining area, while under a high flow condition, the order from big to small is the central zone of reservoir, river banks and downstream, upstream; the distribution order from big to small of Pb concentration is the central zone of reservoir, river banks and downstream, upstream; the distribution order from big to small of Cd concentration is the central zone of reservoir, river banks and downstream, upstream; after reservoir impounding, the concentration of heavy metals in the reservoir can basically meet the water quality standard of Class II in GB 3838—2002 surface water environment quality standard.

metallic mining area;heavy metal migration;MIKE21;water quality;numerical simulation

10.3880/j.issn.1004-6933.2015.02.006

國家自然科學基金(51179053)

軒曉博(1979—),男,工程師,碩士研究生,主要從事水資源保護規劃及環境影響評價研究。E-mail:xbxuan@126.com

X524

A

1004-6933(2015)02-0030-06

2014-10-18 編輯:徐 娟)

猜你喜歡
水質模型
一半模型
水質抽檢豈容造假
環境(2023年5期)2023-06-30 01:20:01
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
一月冬棚養蝦常見水質渾濁,要如何解決?這9大原因及處理方法你要知曉
當代水產(2019年1期)2019-05-16 02:42:04
這條魚供不應求!蝦蟹養殖戶、垂釣者的最愛,不用投喂,還能凈化水質
當代水產(2019年3期)2019-05-14 05:42:48
圖像識別在水質檢測中的應用
電子制作(2018年14期)2018-08-21 01:38:16
3D打印中的模型分割與打包
濟下水庫徑流水質和垂向水質分析及評價
FLUKA幾何模型到CAD幾何模型轉換方法初步研究
主站蜘蛛池模板: 国产成人免费高清AⅤ| 国产伦片中文免费观看| 国产精品视频观看裸模| 亚瑟天堂久久一区二区影院| 99久久无色码中文字幕| 国产成人三级| 97超碰精品成人国产| 又黄又湿又爽的视频| 午夜福利网址| 亚洲人妖在线| 日韩国产无码一区| 在线看免费无码av天堂的| 成人字幕网视频在线观看| 久草视频福利在线观看| 精品剧情v国产在线观看| 国产永久无码观看在线| 国产爽爽视频| 亚洲欧美h| 中文字幕无码中文字幕有码在线| 国产精品欧美激情| 亚洲人成网18禁| 一区二区三区国产| 国产精品污视频| AV老司机AV天堂| 亚洲av无码片一区二区三区| 国产精品无码AⅤ在线观看播放| 日本人妻一区二区三区不卡影院| 亚洲成a人片在线观看88| 伊人久久久久久久| 国产精品55夜色66夜色| 国产成人综合亚洲网址| 国产特级毛片aaaaaaa高清| 性欧美在线| jizz在线观看| 无码AV动漫| 99久久精品免费看国产电影| a毛片在线| 成人国产精品一级毛片天堂 | 亚洲第一视频免费在线| 嫩草国产在线| 久久久精品无码一二三区| 欧洲成人在线观看| 免费a级毛片18以上观看精品| 91欧美亚洲国产五月天| 色九九视频| 午夜久久影院| 亚洲福利视频一区二区| 亚洲欧美日韩中文字幕一区二区三区| 婷婷六月在线| 毛片网站在线看| 99久久精品视香蕉蕉| 国产H片无码不卡在线视频| 99久久精品视香蕉蕉| 视频在线观看一区二区| 在线国产91| 一级成人a毛片免费播放| 全免费a级毛片免费看不卡| 国产精品自在自线免费观看| 亚洲精品第1页| 久久精品国产国语对白| 国产成人综合欧美精品久久| 精品国产成人三级在线观看| 青青青伊人色综合久久| 久久国产热| 国产一级视频久久| 亚洲三级电影在线播放| 国产一级无码不卡视频| 亚欧乱色视频网站大全| 91精品国产无线乱码在线| 激情五月婷婷综合网| 亚洲精品在线观看91| 色婷婷狠狠干| 波多野结衣的av一区二区三区| 色网站在线视频| 亚洲一区二区日韩欧美gif| 亚洲毛片网站| 老司机aⅴ在线精品导航| 特级欧美视频aaaaaa| 久久精品国产999大香线焦| 亚洲欧美日本国产综合在线| 蜜芽国产尤物av尤物在线看| 永久在线精品免费视频观看|