• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    基于數字柵格的河網小流域水質數學模型構建與應用

    2021-07-09 05:18:58尹海龍林夷媛趙東華石澤敏
    同濟大學學報(自然科學版) 2021年6期
    關鍵詞:武進點源柵格

    尹海龍,林夷媛,趙東華,石澤敏

    (1.同濟大學環(huán)境科學與工程學院,上海200092;2.中交上海航道勘察設計研究院有限公司,上海200120)

    我國南方河網地區(qū)河流縱橫,降雨豐富。隨著污水收集率的不斷提高,在城鎮(zhèn)點源污染排放得到有效治理的情況下,村鎮(zhèn)小流域地區(qū)的點源、面源污染排放成為水環(huán)境治理關注的問題。為此,需要建立耦合點源、面源排放與河道水質的數學模型,評估小流域內各類污染排放造成的水質動態(tài)響應,為污染物減排方案的制定提供科學依據[1-4]。

    由于河網模擬范圍較大,大多數情況下只能采用數值方法進行模擬[5]。河道水動力水質模型可以全面描述水體中的水動力—水質過程,但是難以準確獲得面源污染的動態(tài)輸入。因此,非常有必要將降雨產匯流模型和河道水動力水質模型耦合,以此建立河流和流域間的水力聯(lián)系。Liu等[6]集成運用HSPF模型和EFDC模型模擬圣路易斯灣河口的水力水質情況,發(fā)現(xiàn)圣路易斯灣河口對密西西比河的水質具有重大影響;Li等[7]采用SWAT和EFDC耦合方法,研究了長潭水庫1952—2010年間58年降水資料的頻率響應與水質模擬分析;Thompson等[8]將水文模型MIKE SHE和水力學模型MIKE11進行耦合,每個時間段都進行數據的交換和迭代計算,通過這種模擬方式反映了坡面與河道洪泛區(qū)之間的相互作用關系;Narasimhan等[9]基于SWAT和WASP建立了水庫水質評價綜合管理模型,結果表明需要削減35%的TN和TP才能明顯減少水庫葉綠素a的濃度;朱瑤[10]等耦合運用SWAT與WASP模型,模擬苕溪流域面源污染負荷與水質響應情況。然而,物理性分布式水文模型還存在水文過程描述的過參數化問題,對空間信息的精度要求也很大程度上限制了模型的整體效果。

    為了克服水文過程過參數帶來的不確定問題并提高空間尺度上的模擬精度,本文采用柵格基礎的數字高程模型(Raster-based Digital Elevation Model,DEM)構建方法,自動提取水系和劃分子流域,以DEM柵格單元為基本模擬單元進行產流和產污計算,精細化描述降雨徑流導致的面源污染排放與受納水體間的空間響應關系,且應用方便。在自主研發(fā)的柵格化徑流污染模型的基礎上耦合開放式的河道水動力水質模型。開放式模型的優(yōu)勢在于有助于克服傳統(tǒng)商業(yè)軟件建模面臨的智能通訊和系統(tǒng)集成性能差的問題,可進一步耦合智能算法,為實現(xiàn)參數自動率定、污染物精準溯源等奠定基礎[11-12]。以常州市武進港小流域為例開展實證研究,對河網水質進行日尺度的動態(tài)模擬和精準預測。

    1 模型的基本原理

    1.1 水動力水質模型

    1.1.1 一維非恒定流模型

    對于河網地區(qū)而言,河道縱向長度遠大于其寬度和深度,可將河網概化為一維問題。根據質量守恒和動量守恒定律,可建立描述河網水動力過程的一維圣維南方程組,如下所示:

    式中:z為河道水位;x為河道縱向距離;t為時間;Q為斷面過流流量;A為橫斷面面積;Bt為河道橫斷面的非流動調蓄面積;ql為單位長度上的旁側入流量,包括點源產生的污水排放量和降雨徑流入河量,其中降雨徑流入河量由柵格化降雨徑流污染模型計算確定,如1.2節(jié)所述;g為重力加速度;R為河道橫斷面的水力半徑;n為河道粗糙系數。

    1.1.2 水質模型

    河網水質模型控制方程由如下公式描述:

    式中:φ為水質組分質量濃度;E為縱向離散系數;SE為點源和面源污染負荷排放匯入,其中面源污染負荷由柵格化面源模型計算確定,如1.2節(jié)所述;SI是描述水體中生物化學反應動力學過程,其作用機理架構如圖1所示。

    圖1 水質模擬組分的耦合作用關系描述Fig.1 Description of coupling effect of components in water quality simulation

    在圖1中,水質模型描述了水體溶解氧平衡為核心的模型組分之間的相互作用關系。溶解氧的消耗考慮了含碳有機物生化降解耗氧和氨氮硝化作用耗氧以及河道底泥有機質分解耗氧;溶解氧的補充考慮了大氣向水體中的復氧過程。其中氨氮的轉化又涉及到氮的循環(huán)過程,即有機氮、氨氮、硝酸鹽之間的轉化。因此,水質模擬組分主要包括碳化BOD(CBOD)、有機氮、氨氮、硝酸鹽和溶解氧及其之間的降解、轉化速率等。

    考慮到水動力和水質模型的開放性以及將來進一步建立智慧水務平臺的需要,水動力和水質模型的模擬計算由開放式的模擬工具HEC-RAS模型系統(tǒng)(Hydrologic Engineering Center’s River Analysis System)實現(xiàn)[13-14]。在此基礎上實現(xiàn)與自主開發(fā)的柵格化降雨徑流污染模型的集成,滿足河網地區(qū)水質動態(tài)模擬和水質動態(tài)達標分析的功能。

    1.2 柵格化降雨徑流污染模型

    柵格化降雨徑流污染模型以數字等高模型(Digital Elevation Model,DEM)為基礎,自動提取水系和劃分子流域,然后以子流域為基礎計算降雨徑流產流、產污量,并最終分配到對應的河道計算單元。模型的主要功能包括柵格降雨量插值計算提取、柵格徑流深度計算、柵格徑流污染物量計算和柵格水流匯流計算4個模塊,計算思路如圖2a所示。

    圖2 柵格化降雨徑流污染模型原理Fig.2 Principles of precipitation-runoff nonpoint pollution model based on DEM raster

    (1)降雨量柵格分布計算。由于降雨量站點的實測數據僅為站點所處位置的降雨量,采用反向距離加權法(IDW)進行插值獲得整個流域內每個柵格的降雨量[15]。

    (2)降雨徑流深度柵格分布計算。采用SCS曲線值法計算對應不同用地類型的柵格徑流深度。SCS法為降雨徑流關系法,能反映不同土壤類型、不同土地利用方式及前期土壤含水量對降雨徑流的影響,具有參數少、簡單易行等特點[16]。SCS曲線值法的方程為

    其中

    式中:S為潛在最大入滲量;P為一次降雨總量;Q為實際徑流深;CN為該日的曲線值;Ia為由于地表儲存、截留和下滲而導致的地表徑流的初始損失量,與S呈一定的正比關系,通??扇ˇ藶?.2。

    (3)降雨徑流污染物柵格分布計算。進一步,柵格像元上的降雨徑流攜帶污染物量采用如下公式計算:

    式中:L是柵格像元的污染物排放量;Cland為柵格像元的對應土地利用類型污染物徑流質量濃度。

    (4)柵格徑流和污染物匯流計算。采用8點重力流向算法(deterministic eight-node,D8),基于數字地形高差進行匯流流向判斷,自動提取水系,劃分匯水區(qū)域[17]。為避免匯流過程模擬的失真,確定水流方向前需要進行DEM數據無凹陷處理。D8算法的基本思想是:在3×3窗口中,計算中心網格與鄰域8個格網之間的距離權落差,具有最大權落差值的鄰域格網方向即定義為水流方向,并且規(guī)定一個柵格單元的水流方向用一個特征碼表示,如圖2b所示。距離權落差計算公式如下:

    式中:Δhi是中心柵格與鄰域柵格的高程差值,柵格間的距離d與方向有關,在對角線方向為2倍柵格間距,在其他方向為1倍柵格間距。

    根據流向計算結果分析每一柵格單元對應的上游匯水面積值。通過從上游到下游的匯水面積累加,建立河道計算單元與匯水面積的對應關系。結合柵格徑流量和污染物計算結果,將降雨徑流污染分配到對應的河道計算單元。

    自主開發(fā)的降雨徑流污染模型優(yōu)勢在于:①產流計算基于柵格式SCS模型,同時考慮了降雨和下墊面條件的空間不均勻性,具有空間尺度上計算精度高的特點;②產污計算基于輸出系數模型,避免了面源污染形成機理的復雜性,所需參數少、操作簡便。

    2 模型的建立與率定

    2.1 研究區(qū)域概況

    武進地處長三角地理中心,位于常州、無錫及宜興三市交界處。武進港是太湖梅梁灣的主要入湖河流之一,是連通京杭大運河和太湖的主要水上通道。武進港全長29 km,水深2~3 m,河寬25~30 m,流入太湖水量枯水年為2.72×108m3·a-1,豐水年為4.36×108m3·a-1。

    本文將模型從武進港干流擴展到整個流域范圍內的河網地區(qū),研究區(qū)域土地利用類型以農業(yè)用地、居住用地、工業(yè)用地為主。概化河道包括武進港及武進港兩岸采菱港、武南河、永安河、錫漂漕河、雅浦港5條主要河流,另外還有禮嘉大河、小留河、虎臣河、政平河、東陽岸河等28條支流。武進港及其流域概況如圖3所示。

    圖3 研究區(qū)域概況Fig.3 Description of the site studied

    2.2 模型邊界條件

    水動力模型上游和下游邊界均采用水位過程線,數據來自坊前、常州、百瀆口等水文站2018年日監(jiān)測數據。

    水質模型邊界條件包括點源和面源污染負荷的輸入,點源又可分為集中式點源和分散式點源。研究區(qū)域內集中式點源主要包括7座城鎮(zhèn)污水處理廠尾水排放。分散式點源主要包括未經處理直排的農村生活污水和畜禽養(yǎng)殖廢水,依據《太湖流域主要入湖河流水環(huán)境綜合整治規(guī)劃編制技術規(guī)范》[18],采用排污系數法對分散式點源的化學需氧量(COD)、氨氮(NH3-N)、總氮(TN)、總磷(TP)負荷量進行估算,計算公式如下:

    式中:W為分散式點源污染負荷量;Np為人口數;NL為畜禽養(yǎng)殖量;αp為農村生活污水排污系數,其中COD取27 g·人-1·d-1,NH3-N取4 g·人-1·d-1,TN取6 g·人-1·d-1,TP取0.2 g·人-1·d-1;αL為畜禽養(yǎng)殖排污系數,換算成生豬計算,COD取40.55 g·頭-1·d-1,NH3-N取1.31 g·頭-1·d-1,TN取3.05 g·頭-1·d-1,TP取0.38 g·頭-1·d-1;βP為農村生活污水入河系數,取0.7;βL為畜禽養(yǎng)殖入河系數,取0.6。

    面源邊界條件輸入采用降雨徑流污染模型的模擬結果,降雨徑流污染模型中污染排放系數取值如表1所示。需要注意的是,污染負荷統(tǒng)計參數COD需與水質模型輸入參數CBOD進行換算,換算系數依據經驗取值:取5日生化需氧量BOD5與碳的生化需氧量CBOD的比值為0.70,BOD5與COD的比值為0.50,可 推 算 得CBOD與COD的 比 值為0.71[19]。

    表1 研究區(qū)域不同土地利用類型徑流面源污染輸出系數Tab.1 Runoff nonpoint pollution export coefficients from different land use types in the catchments studied

    計算得研究區(qū)域排放總流量為1.6×108t·a-1,COD、NH3-N、TN、TP污染負荷排放總量分別為9 024.1 t·a-1、621.4 t·a-1、1 275.8 t·a-1、80.8 t·a-1,圖4給出了集中式點源、分散式點源和徑流面源污染的負荷排放占比。可見雖然分散式點源排放流量很小,但排放污染負荷量占比卻很高,COD、NH3-N、TN、TP占 比 分 別 為42.1%、46.3%、58.5%、45.2%。對于NH3-N,降雨徑流污染占比也較高。

    圖4 研究區(qū)域各污染類型組成比例Fig.4 Discharge and proportion of each pollution type in the catchments studied

    圖5進一步展示了武進港流域以高精度DEM 5m分辨率柵格為單位的氨氮產生量以及不同河段的污染物空間分布,可以看出,雨季(8月)的降雨徑流污染負荷入河量明顯高于旱季(1月)。

    2.3 模型率定驗證

    水質模型率定包括以下幾個步驟:①輸入上下游水動力邊界及旁側入流量,通過圣維南方程建立河道水動力模型;②輸入點源和面源污染負荷以及各水質參數以模擬河道中水質的時空變化;③比較河道斷面水質的實測值和模擬值,必要時調整水質參數。

    分別采用納什效率系數(Nash-Sutcliffe efficiency,NSE)和決定系數(R2)來評價水動力模型和水質模型質量的好壞,計算公式如下:

    式中:Qsimi、Qobsi分別為日均流量模擬值和實測值;為日均流量實測值的平均值;n代表總共匹配的天數。NSE取值為負無窮至1,NSE越接近1,表示模型質量越好,模型可信度高。通常認為NSE達到0.75以上代表模擬結果很理想,NSE值在0.75~0.36之間模擬結果也是可以接受的,當NSE低于0.36代表模擬結果不理想[20]。

    式中:fsimi、fobsi分別為月均水質濃度模擬值和月水質濃度實測值分別為月均水質濃度模擬值和實測值的平均值;m代表計算時間段內的月份數。R2取值為0至1,R2越接近1,表示模擬精度越好。通常認為R2大于0.6代表模擬結果理想,R2在0.6~0.5之間認為模擬結果總體上可以接受[21]。

    研究區(qū)域河網水動力模型采用黃埝橋水文站2018年1—12月的日監(jiān)測流量數據進行率定,如圖6所示??梢钥闯?,模擬值與實測值吻合情況較好,計算得NSE值為0.71,模擬誤差在允許范圍之內,率定后的河道粗糙系數取值為0.030~0.035。

    圖6 黃埝橋斷面流量率定結果Fig.6 Modeling calibration result for flow discharge at Huangnian Bridge station

    水質模型采用武進港月采樣監(jiān)測數據(2018年4—12月),率定溶解氧(DO)、氨氮(NH3-N)、高錳酸鹽(CODMn)和總磷(TP)4個水質指標,同樣水質監(jiān)測指標CODMn需與水質模型中的率定參數CBOD進行換算,取BOD5/CBOD=0.70、BOD5/CODMn=0.89,可推算得:CODMn/CBOD=0.79[22]。選取慈瀆橋(CDQ2)和武南河下游(WN7)2個典型斷面進行年尺度模擬。CDQ2斷面位于武進港上游,主要污染源為上游客水,未受支流匯入污染影響;WN7斷面位于武南河和武進港交匯處,水質受采菱港、武南河等匯入支流污染影響,總體而言水質差于CDQ2斷面。由率定結果(圖7)可以看出,模擬值與實測值吻合情況較好。依據式(9)計算得:①針對CDQ2斷面,DO、NH3-N、CODMn和TP的R2值分別為0.82、0.69、0.61和0.71;②針對WN7斷面,DO、NH3-N、CODMn和TP的R2值 分 別 為0.88、0.63、0.59和0.66,說明模擬誤差均在允許范圍之內,水質模型中主要參數率定結果如表2。率定結果顯示,降雨強度較大時,受降雨徑流匯入負荷的沖擊,河道斷面水質短時間內會出現(xiàn)急劇升高的現(xiàn)象,降雨結束后又恢復至常態(tài)。

    圖7 CDQ2、WN7斷面水質參數率定結果Fig.7 Calibration results of parameters for water quality modeling at CDQ2 and WN7 station

    表2 水質模型主要參數率定結果Tab.2 Calibrated results of main parameters for water quality modeling

    3 水質改善方案評估

    3.1 污染削減方案比選

    針對流域內分散點源占比高的特點,提出優(yōu)先削減分散實點源污染排放,分析污染削減效果與入湖考核斷面之間的水質響應關系。

    由2.2節(jié)討論可知,雖然分散式點源排放流量很小,但排放污染負荷量占比卻最高。因此優(yōu)化考慮農村分散式點源污染的削減方案,并評估污染削減前后的河道水質動態(tài)達標頻次。

    3.2 畜禽養(yǎng)殖污染控制

    通過全面清理整頓非法和不符合規(guī)范標準的養(yǎng)殖場,進行養(yǎng)殖專業(yè)化,可實現(xiàn)畜禽養(yǎng)殖接近“零排放”。模擬研究片區(qū)內畜禽養(yǎng)殖全部關?;蛲ㄟ^還田再利用實現(xiàn)畜禽污染零排河,即片區(qū)內畜禽養(yǎng)殖污染負荷削減率100%的情況下姚巷橋斷面水質的響應情況,模擬結果見圖8和表3。

    圖8 畜禽養(yǎng)殖負荷削減前后入湖斷面水質響應情況Fig.8 Predicted water quality at lake inflow station before and after livestock pollution reduction

    削減畜禽養(yǎng)殖負荷后水質得到了一定程度的提升,特別是CODMn改善顯著,年達標頻次(地表水III類標準,簡稱III類)從23.0%提升至64.9%,但是NH3-N仍無法達標。

    3.3 農村污水處理

    收集未經處理的農村生活污水,經處理后以一級A標準排放。在畜禽養(yǎng)殖“零排放”基礎上,結合流域內污染物空間分布特征,比選出對水質影響大的污染優(yōu)先控制村落群,優(yōu)先削減該村落群80%農村生活污水,以期進一步提高水質年達標頻次。

    共設計了3種污染優(yōu)先控制村落比選方案:①工況1。選取片區(qū)內農村生活污水排放量最高的區(qū)域7個村落為污染優(yōu)先控制村落,共削減COD 2 591.7t·a-1、NH3-N 120.4 t·a-1、TP 12.9 t·a-1;②工況2。在工況1的基礎上,增加排入武進港干流的10個 村 落,共 削 減COD 3 330.4t·a-1、NH3-N 178.9 t·a-1、TP 30.3 t·a-1;③工況3。在工況1的基礎上,增加排入上游來水京杭大運河的7個村落,共削 減COD 3 205.4 t·a-1、NH3-N 170.9 t·a-1、TP 29.1 t·a-1。分別模擬3種工況下姚巷橋斷面水質的響應情況,模擬結果見圖9和表3。

    圖9 不同農村污水削減工況下入湖斷面水質響應情況Fig.9 Predicted water quality at lake inflow station of different rural sewage reduction schemes

    表3 不同負荷削減工況下姚巷橋斷面水質達標情況分析Tab.3 Annual attainment rates of NH3-N,CODMn,and TP at Yaoxiang Bridge station of different load reduction schemes

    模擬結果可知,進一步削減污染優(yōu)先控制村落相應農村生活污水后,姚巷橋斷面水質得到明顯改善。針對工況1,CODMn和TP改善效果較明顯,CODMn年均值降至5.86 mg·L-1,年達天數(III類)從84 d提升至252d;TP的年均值降至0.22 mg·L-1,年達天數(III類)也從60 d提升至142 d,但是NH3-N仍無法達到地表水III類標準。工況2、3在工況1基礎上分別削減排入武進港和京杭大運河村落群的入河負荷,針對CODMn和TP,工況2、3模擬結果相差不大;但是針對NH3-N,工況3的模擬結果明顯優(yōu)于工況2,NH3-N年均值降低至1.06 mg·L-1,年達標頻次(III類)達66.3%,年達標天數(III類)達242 d。綜上可得工況3為相對最優(yōu)工況。

    綜上,改善客水水質、加強農村生活污水處理設施建設以及分散畜禽養(yǎng)殖污染治理可有效改善武進港水質;其中改善上游來水京杭大運河對提高姚巷橋斷面年達標率效果尤為顯著。相對最優(yōu)工況下仍存在不達標天數主要是因為雨天受降雨徑流匯入負荷的沖擊,河道斷面水質短時間內會出現(xiàn)急劇升高的現(xiàn)象。因此從水質動態(tài)達標的角度考慮,若要實現(xiàn)考核斷面全面、穩(wěn)定達標,需進一步加強對農田徑流和地表徑流污染的控制,如推進面源污染的河道生態(tài)攔截等工程措施,從而確保在雨后一段時間內能盡快恢復水質。

    4 結論

    建立了適用于河網小流域水環(huán)境系統(tǒng)模擬的水文—水動力—水質耦合數學模型,并在常州市武進港流域進行了實證研究。

    (1)建立的耦合模型中,降雨徑流污染模型采用自主開發(fā)的模型系統(tǒng),以DEM柵格單元為基本模擬單位,可以對降雨徑流導致的面源污染時空分布和匯流去向進行精準模擬。在此基礎上,河道水動力水質模型采用開源模型HEC-RAS,構建了流域污染物排放與河道水質的響應關系,實現(xiàn)河網地區(qū)日尺度的水質模擬和動態(tài)預報。模型系統(tǒng)同時具有開放性,能夠為小流域智慧監(jiān)管平臺提供技術支撐。

    (2)針對武進港小流域分散式點源污染負荷高、部分匯入支流水質較差、長期難以達標的問題,將水動力水質模型從武進港干流擴展到整個流域范圍內的河網地區(qū),從而能夠對武進港流域內的污染總量控制進行定量評估。結合流域內入河污染源組分和空間分布特點,模擬不同負荷削減工況下入湖考核斷面姚巷橋的水質響應情況。從動態(tài)達標角度分析考慮,對不同水質方案進行評估。最優(yōu)工況下,針對水質指標CODMn、NH3-N和TP,考核斷面姚巷橋的年達標頻次(III類)分別從23.0%、0、16.4%升高至71.8%、66.3%、75.9%。

    作者貢獻說明:

    尹海龍:數學模型建立、論文撰寫。

    林夷媛:數學模型建立、論文撰寫。

    趙東華:技術和材料支持。

    石澤敏:技術和材料支持。

    猜你喜歡
    武進點源柵格
    常州武進:真情服務開創(chuàng)僑務工作新局面
    華人時刊(2023年7期)2023-05-17 09:05:26
    Listing Tables:a Strategy of Problem Solving in Maths from Image Thinking to Logical Reasoning
    立足武進 僑連海外
    華人時刊(2021年17期)2021-12-02 03:26:10
    基于鄰域柵格篩選的點云邊緣點提取方法*
    關于脈沖積累對雙點源干擾影響研究
    靜止軌道閃電探測性能實驗室驗證技術研究
    基于標準化點源敏感性的鏡面視寧度評價
    武進淹城遺址功能新考
    東方考古(2016年0期)2016-07-31 17:45:44
    不同剖面形狀的柵格壁對柵格翼氣動特性的影響
    基于CVT排布的非周期柵格密度加權陣設計
    雷達學報(2014年4期)2014-04-23 07:43:13
    国产黄a三级三级三级人| 黄色毛片三级朝国网站| 可以在线观看毛片的网站| 亚洲狠狠婷婷综合久久图片| 久久香蕉国产精品| 淫秽高清视频在线观看| 亚洲人成77777在线视频| 亚洲国产欧美人成| 国产一区二区激情短视频| 五月玫瑰六月丁香| 波多野结衣高清作品| 美女高潮喷水抽搐中文字幕| 亚洲av电影在线进入| 88av欧美| 熟妇人妻久久中文字幕3abv| 又黄又粗又硬又大视频| 国产伦一二天堂av在线观看| 亚洲专区国产一区二区| 变态另类成人亚洲欧美熟女| 亚洲中文字幕日韩| 一进一出抽搐gif免费好疼| 男人舔女人的私密视频| av欧美777| 99精品欧美一区二区三区四区| 日日干狠狠操夜夜爽| 日本一二三区视频观看| 高清毛片免费观看视频网站| 亚洲五月婷婷丁香| 久久久国产成人免费| 午夜精品久久久久久毛片777| 国产亚洲av嫩草精品影院| 99在线视频只有这里精品首页| 热99re8久久精品国产| 免费在线观看成人毛片| 看免费av毛片| 国产爱豆传媒在线观看 | 亚洲国产欧美一区二区综合| www.熟女人妻精品国产| 在线观看www视频免费| 他把我摸到了高潮在线观看| 国产亚洲av嫩草精品影院| 成人18禁高潮啪啪吃奶动态图| 亚洲人成网站在线播放欧美日韩| 日本精品一区二区三区蜜桃| 成人av在线播放网站| 一本精品99久久精品77| 国产亚洲av高清不卡| 99久久精品热视频| 久久性视频一级片| 真人做人爱边吃奶动态| 亚洲美女黄片视频| 亚洲七黄色美女视频| 中出人妻视频一区二区| 亚洲一区中文字幕在线| 亚洲成人国产一区在线观看| 成人手机av| 成人18禁高潮啪啪吃奶动态图| 全区人妻精品视频| 黄色视频不卡| 成年免费大片在线观看| 老司机在亚洲福利影院| 99热只有精品国产| 亚洲人与动物交配视频| 成人精品一区二区免费| 无遮挡黄片免费观看| 久久人妻福利社区极品人妻图片| 熟女少妇亚洲综合色aaa.| 不卡av一区二区三区| av国产免费在线观看| 亚洲七黄色美女视频| 久久草成人影院| 免费在线观看视频国产中文字幕亚洲| 国产精品香港三级国产av潘金莲| 男人的好看免费观看在线视频 | 久久久久久国产a免费观看| 99在线视频只有这里精品首页| 久久久精品大字幕| 一进一出抽搐gif免费好疼| 亚洲精品粉嫩美女一区| 成人av在线播放网站| 高潮久久久久久久久久久不卡| 夜夜看夜夜爽夜夜摸| 深夜精品福利| 热99re8久久精品国产| 国产熟女xx| 日韩成人在线观看一区二区三区| 99国产精品一区二区蜜桃av| 老司机午夜十八禁免费视频| 国产精品野战在线观看| 日韩欧美精品v在线| 亚洲七黄色美女视频| 熟女少妇亚洲综合色aaa.| 国产精品1区2区在线观看.| 久久99热这里只有精品18| 丝袜人妻中文字幕| 久久精品国产99精品国产亚洲性色| 欧美丝袜亚洲另类 | 此物有八面人人有两片| 搡老岳熟女国产| 天堂动漫精品| 国产不卡一卡二| 国产99久久九九免费精品| 夜夜躁狠狠躁天天躁| 成人一区二区视频在线观看| 日韩精品中文字幕看吧| 日韩免费av在线播放| 久久99热这里只有精品18| 制服丝袜大香蕉在线| www国产在线视频色| 超碰成人久久| 国产精品,欧美在线| 一本一本综合久久| 搡老熟女国产l中国老女人| 成人三级黄色视频| 黄片大片在线免费观看| 欧美高清成人免费视频www| 亚洲乱码一区二区免费版| 成人午夜高清在线视频| 久久久久久久久免费视频了| 伊人久久大香线蕉亚洲五| 午夜福利免费观看在线| 人妻久久中文字幕网| 变态另类成人亚洲欧美熟女| 淫秽高清视频在线观看| 国产精品亚洲av一区麻豆| 亚洲第一电影网av| 无限看片的www在线观看| 国产又色又爽无遮挡免费看| 国产av又大| 黑人欧美特级aaaaaa片| 日韩免费av在线播放| av在线天堂中文字幕| 特大巨黑吊av在线直播| 亚洲中文字幕一区二区三区有码在线看 | 久久久国产欧美日韩av| 十八禁人妻一区二区| 欧美最黄视频在线播放免费| 午夜福利欧美成人| 亚洲人成电影免费在线| 欧美日韩国产亚洲二区| 国产精品,欧美在线| 麻豆成人午夜福利视频| 成年版毛片免费区| 女生性感内裤真人,穿戴方法视频| 在线看三级毛片| 岛国在线观看网站| 亚洲成人精品中文字幕电影| 日韩精品青青久久久久久| 窝窝影院91人妻| 亚洲国产欧美网| 一二三四在线观看免费中文在| 国产精品影院久久| 黄片小视频在线播放| 少妇裸体淫交视频免费看高清 | 国产97色在线日韩免费| 久久久久久久久免费视频了| 白带黄色成豆腐渣| 亚洲 欧美 日韩 在线 免费| 巨乳人妻的诱惑在线观看| 五月玫瑰六月丁香| 2021天堂中文幕一二区在线观| 在线观看www视频免费| 久久久国产欧美日韩av| 久久精品国产清高在天天线| 亚洲精品国产精品久久久不卡| 国模一区二区三区四区视频 | 国产区一区二久久| 啦啦啦观看免费观看视频高清| 黑人巨大精品欧美一区二区mp4| 嫩草影视91久久| 日韩欧美三级三区| 亚洲性夜色夜夜综合| 国产高清有码在线观看视频 | 中文字幕av在线有码专区| av天堂在线播放| 老鸭窝网址在线观看| 制服人妻中文乱码| 精品电影一区二区在线| 午夜亚洲福利在线播放| 中文字幕久久专区| 国产亚洲精品第一综合不卡| 欧美成狂野欧美在线观看| 欧美高清成人免费视频www| 男人舔女人下体高潮全视频| 麻豆国产av国片精品| 久久精品亚洲精品国产色婷小说| 变态另类丝袜制服| 欧美日韩福利视频一区二区| 男女那种视频在线观看| www日本在线高清视频| 色哟哟哟哟哟哟| 久久婷婷人人爽人人干人人爱| 老鸭窝网址在线观看| 啦啦啦观看免费观看视频高清| 亚洲欧美日韩无卡精品| 成人手机av| 又黄又粗又硬又大视频| 久久久久久久午夜电影| 夜夜夜夜夜久久久久| 亚洲中文日韩欧美视频| 我的老师免费观看完整版| 丰满的人妻完整版| 老司机午夜十八禁免费视频| 日本黄大片高清| 制服丝袜大香蕉在线| 国产在线精品亚洲第一网站| 国产不卡一卡二| 99精品久久久久人妻精品| 日韩国内少妇激情av| 91麻豆精品激情在线观看国产| 黄色片一级片一级黄色片| 白带黄色成豆腐渣| av有码第一页| 99热这里只有精品一区 | 国产野战对白在线观看| tocl精华| 婷婷精品国产亚洲av| 琪琪午夜伦伦电影理论片6080| 老司机福利观看| 国产精品综合久久久久久久免费| 国模一区二区三区四区视频 | 日韩成人在线观看一区二区三区| 国产精品美女特级片免费视频播放器 | av免费在线观看网站| 亚洲精品国产精品久久久不卡| aaaaa片日本免费| 午夜福利成人在线免费观看| 精品国产乱码久久久久久男人| 国产成人影院久久av| 亚洲一区高清亚洲精品| 亚洲一区二区三区不卡视频| 欧美+亚洲+日韩+国产| 国产精品爽爽va在线观看网站| 欧美成人免费av一区二区三区| 亚洲一区二区三区色噜噜| 亚洲欧美一区二区三区黑人| 成人一区二区视频在线观看| 日本免费一区二区三区高清不卡| 国产不卡一卡二| 成人国产一区最新在线观看| 欧美黑人欧美精品刺激| www.www免费av| 久久伊人香网站| 91在线观看av| 亚洲 国产 在线| 国产欧美日韩精品亚洲av| 日韩高清综合在线| 正在播放国产对白刺激| 曰老女人黄片| 丁香欧美五月| 99国产综合亚洲精品| 成人av在线播放网站| 欧美大码av| 国产精品亚洲一级av第二区| 天天躁狠狠躁夜夜躁狠狠躁| 国产三级黄色录像| 午夜老司机福利片| 国产视频一区二区在线看| 精品国产乱码久久久久久男人| 人妻丰满熟妇av一区二区三区| 日韩 欧美 亚洲 中文字幕| 真人一进一出gif抽搐免费| 亚洲一区中文字幕在线| 中文字幕熟女人妻在线| 精品国产乱子伦一区二区三区| 美女免费视频网站| 免费在线观看成人毛片| 天天躁夜夜躁狠狠躁躁| 最近在线观看免费完整版| 欧美精品亚洲一区二区| 制服丝袜大香蕉在线| 黄色毛片三级朝国网站| 哪里可以看免费的av片| 欧美大码av| 中文字幕人成人乱码亚洲影| www国产在线视频色| 久久久久久九九精品二区国产 | 久久久久国产精品人妻aⅴ院| 国产精品日韩av在线免费观看| 观看免费一级毛片| 少妇被粗大的猛进出69影院| 老司机福利观看| 午夜免费成人在线视频| 热99re8久久精品国产| 国产一级毛片七仙女欲春2| 免费搜索国产男女视频| 日韩国内少妇激情av| a级毛片a级免费在线| 两个人看的免费小视频| 三级国产精品欧美在线观看 | 美女扒开内裤让男人捅视频| 99精品在免费线老司机午夜| 少妇人妻一区二区三区视频| 香蕉丝袜av| 日本在线视频免费播放| 国产精品野战在线观看| 欧美+亚洲+日韩+国产| 99国产综合亚洲精品| e午夜精品久久久久久久| 精品不卡国产一区二区三区| 丝袜美腿诱惑在线| 国产人伦9x9x在线观看| 欧美又色又爽又黄视频| 国产亚洲精品久久久久5区| 亚洲avbb在线观看| 国产精品99久久99久久久不卡| 国产av一区二区精品久久| 国产视频内射| 最好的美女福利视频网| 久久久精品国产亚洲av高清涩受| 老熟妇仑乱视频hdxx| 禁无遮挡网站| 精品国产亚洲在线| 欧美成狂野欧美在线观看| 女同久久另类99精品国产91| 一边摸一边抽搐一进一小说| 手机成人av网站| 91av网站免费观看| 国产免费男女视频| 亚洲欧美一区二区三区黑人| 国产精品乱码一区二三区的特点| 久久人妻福利社区极品人妻图片| 亚洲专区字幕在线| 91成年电影在线观看| 窝窝影院91人妻| www.999成人在线观看| 久久久久精品国产欧美久久久| 日韩精品免费视频一区二区三区| 亚洲国产中文字幕在线视频| 国产主播在线观看一区二区| 99热这里只有精品一区 | 男女下面进入的视频免费午夜| 女人被狂操c到高潮| 国产av在哪里看| 无人区码免费观看不卡| 少妇粗大呻吟视频| 午夜精品久久久久久毛片777| 亚洲专区国产一区二区| 国产三级黄色录像| 国产蜜桃级精品一区二区三区| 欧美成人性av电影在线观看| 午夜久久久久精精品| 久久婷婷成人综合色麻豆| 久久香蕉国产精品| 日本撒尿小便嘘嘘汇集6| 超碰成人久久| 香蕉av资源在线| 免费在线观看黄色视频的| 一本大道久久a久久精品| 在线观看一区二区三区| 男插女下体视频免费在线播放| 日本三级黄在线观看| av免费在线观看网站| 深夜精品福利| 别揉我奶头~嗯~啊~动态视频| 黄色a级毛片大全视频| 色哟哟哟哟哟哟| 舔av片在线| 在线播放国产精品三级| 久久精品人妻少妇| 国产成人精品无人区| 黄色a级毛片大全视频| 欧美av亚洲av综合av国产av| 青草久久国产| 香蕉久久夜色| 亚洲性夜色夜夜综合| 亚洲国产日韩欧美精品在线观看 | 午夜福利18| 亚洲免费av在线视频| 婷婷精品国产亚洲av在线| 精品一区二区三区四区五区乱码| 欧美日本亚洲视频在线播放| 少妇的丰满在线观看| 成人三级黄色视频| 国产精品久久电影中文字幕| 国产欧美日韩一区二区三| 国产片内射在线| 91麻豆精品激情在线观看国产| 国产高清videossex| 久久婷婷人人爽人人干人人爱| 成人18禁在线播放| 午夜免费观看网址| 国内精品久久久久精免费| 久久久久久久久久黄片| 国产精品电影一区二区三区| a级毛片a级免费在线| 在线观看舔阴道视频| 叶爱在线成人免费视频播放| 国产午夜精品论理片| 一个人免费在线观看电影 | 国产在线观看jvid| 国产精品久久电影中文字幕| 宅男免费午夜| 久久久精品大字幕| 色综合站精品国产| 一本综合久久免费| 免费搜索国产男女视频| 欧美成人午夜精品| 99精品久久久久人妻精品| 18禁观看日本| 97超级碰碰碰精品色视频在线观看| 午夜福利18| 特大巨黑吊av在线直播| 他把我摸到了高潮在线观看| 高潮久久久久久久久久久不卡| 久久久水蜜桃国产精品网| 日日夜夜操网爽| 久久久久久大精品| 两性午夜刺激爽爽歪歪视频在线观看 | 50天的宝宝边吃奶边哭怎么回事| 国产精品日韩av在线免费观看| 在线国产一区二区在线| 在线观看美女被高潮喷水网站 | 国产精品精品国产色婷婷| 免费在线观看视频国产中文字幕亚洲| 久久亚洲真实| 精品一区二区三区视频在线观看免费| 别揉我奶头~嗯~啊~动态视频| 嫩草影视91久久| 最新在线观看一区二区三区| 亚洲片人在线观看| 久久午夜综合久久蜜桃| 亚洲精品国产精品久久久不卡| 中文字幕av在线有码专区| 免费看美女性在线毛片视频| 777久久人妻少妇嫩草av网站| 国产成人影院久久av| 三级国产精品欧美在线观看 | 亚洲成av人片在线播放无| 在线观看美女被高潮喷水网站 | 91av网站免费观看| 国产精品98久久久久久宅男小说| 欧美成狂野欧美在线观看| 亚洲专区中文字幕在线| www国产在线视频色| 国产主播在线观看一区二区| 淫妇啪啪啪对白视频| 黄色视频不卡| 熟妇人妻久久中文字幕3abv| 1024视频免费在线观看| 老熟妇仑乱视频hdxx| 桃色一区二区三区在线观看| 亚洲欧美精品综合一区二区三区| 禁无遮挡网站| 成人国产一区最新在线观看| 99久久综合精品五月天人人| 宅男免费午夜| 一级毛片高清免费大全| 亚洲av第一区精品v没综合| 老司机在亚洲福利影院| 1024手机看黄色片| 2021天堂中文幕一二区在线观| 欧美日韩亚洲国产一区二区在线观看| 免费电影在线观看免费观看| 亚洲精品av麻豆狂野| 香蕉丝袜av| 婷婷精品国产亚洲av| 国产av又大| 国产av在哪里看| 亚洲精品国产一区二区精华液| 在线免费观看的www视频| 日日夜夜操网爽| 国内毛片毛片毛片毛片毛片| 啪啪无遮挡十八禁网站| 久久婷婷成人综合色麻豆| aaaaa片日本免费| 国产不卡一卡二| xxx96com| 欧美性猛交╳xxx乱大交人| 亚洲av五月六月丁香网| 国产欧美日韩一区二区精品| 精品久久久久久,| 久久久久精品国产欧美久久久| 天堂√8在线中文| 国产午夜精品论理片| 一级片免费观看大全| 亚洲精品中文字幕在线视频| 欧美3d第一页| 国产私拍福利视频在线观看| 国产成人精品久久二区二区免费| 日韩av在线大香蕉| 久久久久久免费高清国产稀缺| 亚洲欧美激情综合另类| 少妇被粗大的猛进出69影院| 真人做人爱边吃奶动态| cao死你这个sao货| 不卡一级毛片| 成人国语在线视频| 久久天堂一区二区三区四区| 中文字幕高清在线视频| 国产精品久久久久久人妻精品电影| 全区人妻精品视频| 99久久精品国产亚洲精品| 国产高清激情床上av| 国产99久久九九免费精品| 在线观看美女被高潮喷水网站 | 免费人成视频x8x8入口观看| 麻豆久久精品国产亚洲av| 欧美3d第一页| 欧美黄色片欧美黄色片| 国产伦在线观看视频一区| 又大又爽又粗| 亚洲 国产 在线| 亚洲avbb在线观看| 国产精品九九99| 又黄又爽又免费观看的视频| 久久久国产成人免费| 一级a爱片免费观看的视频| 日本一本二区三区精品| 日韩欧美免费精品| 一区二区三区激情视频| 亚洲色图 男人天堂 中文字幕| 成人av一区二区三区在线看| 麻豆成人午夜福利视频| 好男人电影高清在线观看| 97人妻精品一区二区三区麻豆| 亚洲人成伊人成综合网2020| 午夜福利18| 99在线视频只有这里精品首页| 色综合婷婷激情| 久久香蕉国产精品| 看片在线看免费视频| 国产精品自产拍在线观看55亚洲| 国产亚洲精品综合一区在线观看 | 欧美久久黑人一区二区| 国产欧美日韩精品亚洲av| 一级毛片高清免费大全| 亚洲九九香蕉| 国语自产精品视频在线第100页| 欧美成人性av电影在线观看| 国产黄片美女视频| cao死你这个sao货| 精品不卡国产一区二区三区| 别揉我奶头~嗯~啊~动态视频| av中文乱码字幕在线| 国产视频内射| 国产又色又爽无遮挡免费看| 首页视频小说图片口味搜索| 国产精品香港三级国产av潘金莲| 国内久久婷婷六月综合欲色啪| 两个人免费观看高清视频| 成人18禁高潮啪啪吃奶动态图| 哪里可以看免费的av片| 亚洲国产欧美网| 日本a在线网址| 老司机在亚洲福利影院| 午夜免费激情av| 老熟妇仑乱视频hdxx| 女警被强在线播放| 一区二区三区高清视频在线| 免费在线观看亚洲国产| 国产精品久久久av美女十八| 一本精品99久久精品77| 亚洲中文av在线| 一本综合久久免费| 亚洲国产欧美一区二区综合| 久久久久久免费高清国产稀缺| 18禁黄网站禁片午夜丰满| 国产精品久久久av美女十八| 国产高清videossex| АⅤ资源中文在线天堂| 久久久久久久久中文| 香蕉丝袜av| 欧美激情久久久久久爽电影| 国产亚洲精品久久久久5区| 久久久久国产精品人妻aⅴ院| 啪啪无遮挡十八禁网站| 亚洲人成网站在线播放欧美日韩| 久久九九热精品免费| www.自偷自拍.com| 欧美3d第一页| 久久精品夜夜夜夜夜久久蜜豆 | 国产视频一区二区在线看| 午夜精品在线福利| 久久精品人妻少妇| 久久草成人影院| 国产av一区二区精品久久| 亚洲在线自拍视频| 90打野战视频偷拍视频| 两性午夜刺激爽爽歪歪视频在线观看 | 91麻豆av在线| 色综合亚洲欧美另类图片| 高潮久久久久久久久久久不卡| 91成年电影在线观看| 亚洲最大成人中文| a在线观看视频网站| 日本免费一区二区三区高清不卡| 欧美丝袜亚洲另类 | 五月伊人婷婷丁香| 国产激情欧美一区二区| 久久精品成人免费网站| 国产野战对白在线观看| 亚洲 欧美一区二区三区| 日本免费a在线| 在线a可以看的网站| 国产欧美日韩一区二区精品| 欧美日本亚洲视频在线播放| 亚洲国产中文字幕在线视频| 亚洲男人天堂网一区| 国产精品电影一区二区三区| 听说在线观看完整版免费高清| 欧美性长视频在线观看| 久久天躁狠狠躁夜夜2o2o| 麻豆久久精品国产亚洲av| xxxwww97欧美| 国产高清激情床上av| 中文字幕久久专区| 久久中文字幕一级| 亚洲专区国产一区二区| 久久香蕉激情| 蜜桃久久精品国产亚洲av| 国产av一区二区精品久久| 成人一区二区视频在线观看| 欧美日韩乱码在线| 熟女电影av网|