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

截流緩釋生態(tài)工程對流域徑流影響

2023-07-20 09:26:06錢昆侖張漢辰張維江
中國農村水利水電 2023年7期
關鍵詞:生態(tài)工程模型

錢昆侖,張漢辰,張維江

(1. 寧夏大學土木與水利工程學院,寧夏 銀川 750021; 2. 寧夏大學地理科學與規(guī)劃學院,寧夏 銀川 750021)

1 研究背景

淤地壩是經多年治理實踐證明適合黃土高原的重要的水土保持措施[1]。在攔泥保土、減少入黃泥沙、防洪減災、淤地造田、鞏固退耕還林、保障生態(tài)安全、促進糧食生產和水資源合理利用及經濟社會穩(wěn)定發(fā)展等方面發(fā)揮了重要作用[2]。但是,目前淤地壩存在幾個問題:①黃土高原現(xiàn)存淤地壩大多建于20世紀60、70 年代,由于技術原因,大多淤地壩建設標準不高,配套設施不完善,存在極大的安全隱患[3,4]。②根據(jù)數(shù)量統(tǒng)計,黃土高原地區(qū)淤地壩已淤滿41 008 座,占總數(shù)58 776 座的69.77%[5],其中相當一部分淤地壩已經削弱甚至喪失繼續(xù)攔泥和滯洪能力。③已經建成的淤地壩缺乏必要的維護和保養(yǎng),致使存在垮塌致災的風險[6]。因此,面對這些問題,有必要探索淤地壩新的發(fā)展方向,繼續(xù)實現(xiàn)涵養(yǎng)水源、保持水土、恢復生態(tài)、維護流域生命共同體的目標,從源頭上為維護黃河生命安全量身定制理論和技術措施。

駱駝林流域位于寧夏南部山區(qū),該區(qū)域水土流失嚴重,生態(tài)環(huán)境脆弱,水力侵蝕強度高,屬于半干旱半濕潤地區(qū)。不僅如此,當?shù)厮Y源時間分布極不均衡,當?shù)刂餮雌诮邓空既杲邓康?9%~66%。流域內水資源工程和水保工程布局不盡合理,導致流域內降雨徑流一般以洪水形式出現(xiàn),河道基流量少,部分河道呈現(xiàn)季節(jié)性干涸形態(tài),洪水資源無法利用。為了改善這種狀況,在當?shù)匦藿艘幌盗薪亓骶忈屔鷳B(tài)壩,生態(tài)壩結合生物措施以及工程措施,植物和壩體可以在較寬的時空尺度上協(xié)同演化,一來可以減小溝壑侵蝕,減輕對下游河道沖刷,并可以起到削峰滯洪的作用,二來生物措施可以進一步加固壩體,可使工程免于人工維護。同時,河水被生態(tài)壩攔截后通過壩體滲透作用,形成穩(wěn)定的滲流,保證河道常流水。

CASC2D-SED 模型是一個結合地理信息系統(tǒng)、遙感以及計算機技術的優(yōu)點,基于物理基礎的分布式水文模型。李致家等[7]將模型產流模塊應用于無資料地區(qū)徑流模擬,發(fā)現(xiàn)在無資料地區(qū)模型具有適用性。張漢辰等[8]將模型應用于半干旱半濕潤地區(qū),發(fā)現(xiàn)模型在該地區(qū)具有良好的模擬結果。晁麗君等[9]通過模型應用發(fā)現(xiàn),模型在超滲產流流域與混合產流流域均有不錯的模擬精度。上述均為模型在一些地區(qū)的適用性研究,本文將截流緩釋生態(tài)工程建成后流域水文條件變化與模型徑流計算模塊耦合,耦合模型能更精確的表達出工程對流域徑流過程的影響。

本文以駱駝林流域為研究區(qū)域研究截流緩釋生態(tài)工程對流域水文的影響,通過工程未建時流域降雨-徑流資料構建流域CASC2D-SED 模型,在此模型基礎上,利用截流緩釋生態(tài)工程建設后土地利用參數(shù)、土壤參數(shù)以及水流動力參數(shù)的變化將工程與模型耦合,得到工程CASC2D-SED 模型,通過兩種模型對工程前后流域徑流的模擬分析,流域徑流系數(shù)和日均徑流量的計算分析,揭示截流緩釋生態(tài)工程對流域徑流的影響。最終可以驗證改良后的水土保持措施對小流域徑流的影響,為淤地壩新的發(fā)展提供借鑒。

2 數(shù)據(jù)來源及研究方法

2.1 研究區(qū)概況

駱駝林流域作為全國18個集中連片的特殊貧困地區(qū)之一,處于涇河水系的頡河流域,位于寧夏涇源縣北28 km 處的六盤山鎮(zhèn),東經106°14'24″,北緯35°36'31″,處于六盤山段東麓,為六盤山土石山林區(qū),屬于黃土丘陵溝壑區(qū)的第二副區(qū)。流域的地形較為復雜,以黃土梁峁侵蝕地貌為主,溝壑縱橫,水土流失嚴重,多年侵蝕模數(shù)3 500 t∕km2,生態(tài)環(huán)境脆弱,立地條件較差,地勢西北高,東南低,海拔高度1 921~2 122 m。該流域所在六盤山地區(qū)多年平均降水量為504.0~537.4 mm,從五年平均滑動曲線來看,整個研究區(qū)年降水量呈波動式變化,約12~15 a 為一個大的波動周期。

2.2 截流緩釋生態(tài)工程

截流緩釋工程是對小型淤地壩的改良。其結合生物措施以及工程措施,工程部分是以土工織物作為充填袋,將級配砂石作為填充壘砌而成具有過濾透水作用的工程體;生物部分是選取適宜的草本植物以及木本植物構建植被帶,待植物長成后壩頂?shù)哪颈局参锖蜕舷掠蔚牟荼局参飳M成巨大的生物體[10]。生物措施和工程措施能夠在較寬的時空尺度協(xié)同演化,改變周邊水環(huán)境以及水力條件。由于結構簡單,其施工方便,因此建造速度快,具有一定的可推廣性。

2.2.1 截流緩釋生態(tài)工程布置情況

流域內共新建有39 座截流緩釋生態(tài)壩,壩高0.8~1.5 m,均于2018 年年初枯水季節(jié)全部建造完成,其中上海子(1、4 號水系,為流域干流)建有24 座,下海子(3、5 號水系,為流域支流)建有11座,刺槐溝(2號水系)建有4座。流域總面積5.7 km2,流域出口設置流量觀測站,流量站控制面積3.3 km2。具體布置情況如圖1。

圖1 工程布置情況Fig.1 Engineering layout

2.2.2 截流緩釋工程斷面特征概化

將截流緩釋工程假定為截流層、物理緩釋層和生物緩釋層3 個部分,其中截流層由河道斷面底部土或砂等顆粒較小的材料組成,緩釋層的下部是在河道斷面中上部分由砂或石等顆粒較大的材料組成的透水層,緩釋層的上部是在河道斷面頂部及河漫灘部分由喬木或灌木組成的植被層。其特征如圖2。

圖2 工程結構斷面特征Fig.2 Characteristics of Engineering Structural Sections

2.3 數(shù)據(jù)來源

研究中建立CASC2D-SED 水文模型對截流緩釋生態(tài)工程建設前后水沙情況進行模擬,模型輸入需要包括流域內降雨資料、流域DEM 數(shù)據(jù)、流域內土地利用類型、土壤類型、河道水系特征、河道網絡數(shù)據(jù)等資料。其中降雨資料以及流量資料采用2017-2019年的觀測資料。

(1)降雨資料。研究區(qū)域面積小,降雨可視為均勻分布,流域設置一臺型號為Vantage Pro2 的便攜式氣象站記錄降雨過程。

(2)DEM 數(shù)據(jù)。使用GPS 對駱駝林流域進行野外地形數(shù)據(jù)采集,獲得了1∶2 000 的流域地形圖,實測流域高程1 866~2 142 m,平均高程2 012 m,平均坡度10.6°。將其轉換為20 m精度的流域DEM圖,共生成14 267 個網格。

(3)土地利用和土壤類型。通過地理信息系統(tǒng)、遙感并結合現(xiàn)場調研獲得,其中土壤類型為黑壚土(湘黃土),主要土地利用類型為荒地、草地、林地。

(4)河道網絡數(shù)據(jù)、流域水系圖。通過ArcGIS 處理得到流域的DEM,并進行填洼、匯流等操作,生成水系圖、河道網絡圖;之后數(shù)據(jù)采用實地勘測進行對比修正。

(5)流量觀測。流域匯流出口處使用H-Flume 地表徑流觀測系統(tǒng)對徑流進行實時監(jiān)測。

(6)土壤缺水度。烘干法測定土壤含水率,土壤缺水度為土壤有效孔隙度減土壤含水率。

2.4 研究方法

2.4.1 CASC2D-SED模型簡介

CASC2D 模型結構最初源于Julien 教授對二維地面徑流計算方法的發(fā)展,采用APL 語言編寫計算程序[7],后來加入Green-Ampt 下滲計算、一維顯示擴散波河道計算、二位土壤侵蝕算法和地下徑流模擬等,模型得以更加完善、系統(tǒng),模型名稱也變更為CASC2D-SED。

模型水流演算部分包括[11]:①降雨計算,采用距離平方倒數(shù)法描述降雨分布。②截流計算,降雨等扣除截流等初始損失后開始入滲。③產流計算,采用Green-Ampt產流模式計算柵格瞬時下滲率和累積下滲量。④坡面匯流計算,采用擴散波的二維顯式有限差分格式計算。⑤河道匯流計算,采用一維顯式有限差分擴散波方法計算。

CASC2D-SED 模型通過土地利用參數(shù)(如:降雨截留系數(shù))反映降雨的截留過程,通過土壤參數(shù)(如:飽和毛管水頭、飽和水力傳導度、初始土壤缺水量等)反映降雨的產流過程,通過水流運動參數(shù)(如:曼寧系數(shù)、糙率等)反映匯流過程。

2.4.2 工程對流域徑流影響分析方法

流域匯流出口設有徑流觀測系統(tǒng),可以測定流域出口斷面流量,進行逐月的降雨徑流分析、徑流系數(shù)計算、徑流量以及徑流系數(shù)變化趨勢分析,闡明截流緩釋生態(tài)工程對流域徑流的影響。

根據(jù)17 年工程未建時的流域降雨徑流實測資料建立駱駝林流域CASC2D-SED 模型(以下簡稱流域模型),可通過此模型推求流域未受工程影響時的降雨徑流過程。根據(jù)18 年后實測資料,將建成的工程與模型耦合,耦合過程中是在已建成的流域CASC2D-SED 模型基礎上,保持其他模型輸入不變,將變化的土地利用參數(shù)、土壤參數(shù)、水流運動參數(shù)以及DEM 數(shù)據(jù)代替原有數(shù)據(jù),作為新的模型輸入。土地利用參數(shù)、DEM 根據(jù)實測資料獲?。挥捎诹饔蛲寥李愋筒蛔?,土壤參數(shù)中變化要素主要為缺水度,采用烘干法測定土壤含水量,后根據(jù)模型中土壤含水量與缺水度的關系計算模型輸入缺水度;采用水準儀測量工程建設前后溝道坡降變化值;水流運動參數(shù)結合參數(shù)率定得到結果。耦合過程中模型的截流、產流、匯流過程發(fā)生變化,最終得到與工程耦合的CASC2D-SED 模型(以下簡稱工程模型)可以模擬流域受截流緩釋生態(tài)工程影響下的徑流過程。

研究通過截流緩釋生態(tài)工程建成前后流域徑流量、徑流系數(shù)對比分析,結合CASC2D-SED 模型對工程建成前后流域徑流模擬對比分析,得到截流緩釋生態(tài)工程對流域徑流的影響結果。

2.4.3 精度評價標準

根據(jù)現(xiàn)行《水文情報預報規(guī)范》研究采用確定性系數(shù)、洪峰相對誤差來評定徑流模擬結果[12]。

(1)確定性系數(shù)Dc。模擬流量與實測流量之間的吻合程度用確定性系數(shù)作為指標,該系數(shù)越接近1,說明模擬與實測流量之間擬合程度越好,計算公式為:

式中:n為模擬時間段數(shù);Qc為模擬流量,m3∕s;Q0為實測流量,m3∕s;為實測流量平均值,m3∕s。

(2)相對誤差Df,%。模擬誤差除以實測值為相對誤差,以百分數(shù)表示。文中用來衡量模擬洪峰與實測洪峰之間的偏差程度,計算公式為:

式中:Q0為實測洪峰流量,m3∕s;Qc為模擬洪峰流量,m3∕s。

3 結果與分析

3.1 參數(shù)率定

選用駱駝林流域2017 年水文觀測數(shù)據(jù),建立流域CASC2D-SED 模型。使用人工試錯法進行參數(shù)率定,用確定性系數(shù),相對誤差衡量模擬結果的精度,在參數(shù)率定中,發(fā)現(xiàn)截流深度、飽和水力傳導度、飽和毛管水頭、初始土壤缺水量均為高敏感參數(shù),糙率為極高敏感度參數(shù)。土地利用類型參數(shù)和土壤類型參數(shù)率定結果見表1、2。

表1 土地利用類型參數(shù)表Fig.1 Engineering layout

表2 土壤類型參數(shù)表Fig.2 Characteristics of Engineering Structural Sections

流域匯流出口處進行徑流觀測,設有H-Flume 地表徑流觀測系統(tǒng)。測量時通過時間間隔超聲波液位傳感器捕捉水槽水位,水位與水槽定義截面積乘積可以計算出出口處流量。流量根據(jù)系統(tǒng)設置五分鐘測量一次。根據(jù)地表徑流觀測系統(tǒng)對流域徑流的觀測結果,擇了其中的5 場進行模擬。其中前3 場對模型參數(shù)進行率定,后2場對模型進行驗證,產流模擬相關值見表,模擬值及實測值對比圖如圖3。

表3為率定期和驗證期相對誤差以及確定性系數(shù)值。由表3可知,洪峰相對誤差絕對值在20%以內,說明模擬洪峰接近于實測洪峰;確定性系數(shù)均在0.7 以上,說明模擬流量與實測流量之間吻合度高??傮w說明CASC2D模型在研究區(qū)具有適用性。

表3 參數(shù)率定及驗證Tab.3 Comparison of simulation effects in Luotuolin watershed

3.2 結果與分析

3.2.1 逐月降雨-徑流過程分析

對駱駝林流域2017-2019 年間具有明顯降雨徑流過程的12個月經行降雨-徑流過程分析,對流域凈降雨量變化趨勢、徑流量變化趨勢進行分析,對單位降雨量、單位徑流量、徑流系數(shù)進行計算。計算結果如表4。

表4 駱駝林流域降雨-徑流特征表Tab.4 Rainfall-Runoff characteristic table in Luotuolin Watershed

對流域內2017-2019 年12 個月的日均降雨量和日均流量進行了計算,流域平均降雨量為0.116 mm∕d,從2017 年的0.144 mm∕d 下降為2019 年的0.100 mm∕d。日均徑流量有明顯增加的趨勢,2017 年徑流相對較小僅為33.71 m3∕h,2018 和2019 年徑流量相對較大,分別為47.29 m3∕h 和28.28 m3∕h。2017、2018、2019 年平均徑流系數(shù)分別為11.83%、13.55%、15.26%,具有逐年增大的趨勢,如圖4。對于月均徑流系數(shù),通過徑流系數(shù)與時間的線性擬合可發(fā)現(xiàn),當去掉偏差較大的兩個點后R2幾乎達到0.8,相關系數(shù)為0.892,如圖5,說明徑流系數(shù)與工程建設時間長短的相關程度高,月均徑流呈現(xiàn)顯著的增加趨勢。

圖4 日均降雨徑流量趨勢圖Fig.4 Trend chart of daily average rainfall runoff

圖5 去除偏差較大點后徑流系數(shù)圖Fig.5 Runoff coefficient diagram after removing points with large deviation

3.2.2 典型洪水降雨-徑流過程分析

對2017 和2018、2019 年各選取一場典型降雨徑流過程,采用流域模型和工程模型分別進行模擬,對比分析截流緩釋生態(tài)工程對洪水的調控過程。首先對2018 和2019 年的降雨-徑流過程進行工程CASC2D-SED 模型模擬,后與流域CASC2D-SED模型對相同年份模擬結果進行對比;2017年的洪水過程進一步與構建的工程CASC2D-SED 模型對同一年洪水模擬結果進行對比,其結果如表5。

表5 工程對典型洪水過程影響結果Tab.5 Impact results of closure and slow-release ecological system on typical flood process

201707050000 次洪水實測最大流量為0.023 m3∕s,流域CASC2D-SED 模型模擬最大流量為0.025 m3∕s,通過工程CASC2D-SED 模型模擬該次洪水過程受截流緩釋生態(tài)工程影響后的徑流過程,模擬最大流量為0.015 m3∕s。截流緩釋生態(tài)工程影響下最大流量減少了34.78%。實測值與模擬值如圖6。

201808100000 次洪水實測最大流量為0.012 m3∕s,工程CASC2D-SED 模型模擬最大流量為0.012 m3∕s,通過流域CASC2D-SED 模型模擬該次洪水過程不受截流緩釋生態(tài)工程影響前的徑流過程,模擬最大流量為0.019 m3∕s。截流緩釋生態(tài)工程影響下最大流量減少了36.84%。實測值與模擬值如圖7。

圖7 201808010000 次洪水降雨-徑流過程實測模擬對比Fig.7 Comparison of measured and simulated flood rainfall-runoff process (201808010000)

201908020100 次洪水實測最大流量為0.008 m3∕s,工程CASC2D-SED 模型模擬最大流量為0.008 m3∕s,通過流域CASC2D-SED 模型模擬該次洪水過程不受截流緩釋生態(tài)工程影響前的徑流過程,模擬最大流量為0.011 m3∕s。截流緩釋生態(tài)工程影響下最大流量減少了27.27%。實測值與模擬值如圖8。

圖8 201908020100 次洪水降雨-徑流過程實測模擬對比Fig.8 Comparison of measured and simulated flood rainfall-runoff process (201908020100)

4 結 論

在構建的流域CASC2D-SED 模型基礎上,利用截流緩釋生態(tài)工程建設后土地利用參數(shù)、土壤參數(shù)以及水流動力參數(shù)的變化將工程與模型耦合,得到工程CASC2D-SED 模型,通過兩種模型對工程前后流域徑流的模擬分析,流域徑流系數(shù)和日均徑流量的計算分析,得到以下結論。

(1)在3 年的觀測期內,在截流緩釋生態(tài)工程的影響下,流域日均徑流量顯著增加;去除偏差較大值后,徑流系數(shù)與工程建設時間長短呈顯著相關關系,相關系數(shù)達到0.892,從2017年到2019 年,年平均徑流系數(shù)從11.83%增加到15.26%。截流緩釋生態(tài)工程建設后流域水資源時間分布更加合理,溝道基流量顯著增加。

(2)通過截流緩釋生態(tài)工程建設前后洪峰流量模擬對比分析,結合流量實測資料,發(fā)現(xiàn)受截流緩釋生態(tài)工程影響的流域洪峰流量相較于未受截流緩釋生態(tài)工程影響的流域洪峰流量有較為顯著的削減,且能夠減少洪峰流量的30%左右。

5 展 望

流域下墊面是復雜的,工程建設對流域下墊面影響也是復雜的,在工程與模型耦合過程中所考慮的土地利用參數(shù)、土壤參數(shù)以及水動力參數(shù)等不足以概括流域下墊面的全部變化,因此模擬精度會受到一定的影響。從峰現(xiàn)時間來看,工程建設前后峰現(xiàn)時間并無顯著差別,若是減小模型輸出步長,峰現(xiàn)時間之差可能會更加顯著。由于截流緩釋生態(tài)工程建設完成時間較短,所得到的實測資料較少,因此數(shù)據(jù)的代表性上可能會有欠缺。

國內外眾多研究結果表明,淤地壩可以明顯削減洪峰流量,減少洪水總量[13,14]。而淤地壩削減洪峰主要通過壩前水庫以及擴大河道過水斷面降低洪水流速,削減洪水[15]。截流緩釋生態(tài)工程同樣可以顯著的減小洪水總量、洪洪峰流量,其原因主要在于:在無降雨發(fā)生或降雨量較小時,截流緩釋工程的截流層和物理緩釋層起主導作用,水流通常不會超出河道,截流層的材料顆粒較小,水流僅可能通過滲透的方式通過過;物理緩釋層的材料顆粒較大,能夠讓水流通過,減緩水流流速。在降雨較大時,由于徑流量變大,河道水流在河道底部攔截的作用下會漫上河漫灘,這種情況下截流緩釋工程的生物緩釋層也發(fā)揮重要作用,匯流斷面由原本的河道斷面變成了由河道和河漫灘組成的混合斷面,同時由于植被的作用,斷面形態(tài)特征和水力條件均發(fā)生變化,匯流糙率明顯增加,在相同流量的條件下,截流緩釋工程的匯流斷面面積更大,必然會減緩水流流速,削減洪峰流量。 由于駱駝林流域面積較小,流域對洪水調蓄能力較弱,對降雨的響應較為敏感,產匯流速度較快,因此截流緩釋工程對該流域峰現(xiàn)時間的影響不夠明顯。

淤地壩建壩初期以攔蓄水資源為主,一些研究結果表明,淤地壩的存在顯著減少了流域徑流量[16,17],降低流域徑流系數(shù)[18,19]。但是截流緩釋工程存在一定的差異,主要因為截流緩釋生態(tài)工程不截斷水流,因此會降低水資源因為攔蓄所帶來損失,另一方面,由于工程作用,一部分地表水以入滲的方式轉換為壤中水和地下水,非雨期流量隨之增加,即增加了溝道基流量,影響水資源再分配。地下水比例增加也可以減少蒸發(fā)損失,因此在截流緩釋生態(tài)工程作用下流域徑流量和徑流系數(shù)均有增加趨勢。因此在截流緩釋生態(tài)工程的作用下,流域由原來的季節(jié)性斷流,改變?yōu)槿缃竦某A魉?/p>

綜合上述研究成果,截流緩釋生態(tài)工程一方面可以削減洪峰,另一方面可以增加溝道基流量,這對于類似于寧南山區(qū)這種缺水地區(qū)的水資源利用有著重要的意義。工程結合工程措施和生物措施,得到人工干預下的降雨-徑流響應機制,將雨洪資源區(qū)域性攔蓄,使溝道由原來的季節(jié)性斷流變?yōu)殚L流水,增加水資源的利用效率,挖掘水資源潛力。同時截流緩釋生態(tài)壩結構簡單,易于施工,且筑壩材料易于就地取材。因此,該工程可應用于降雨集中,降雨強度大的水資源匱乏區(qū)域,由于成本低廉、施工簡單、材料易得,因此該工程具有可推廣性。本文所涉及的僅僅為淤地壩改良的一種思路,期待有更優(yōu)秀的工程型式可以挖掘出缺水地區(qū)的水資源潛力,實現(xiàn)人水和諧。

猜你喜歡
生態(tài)工程模型
一半模型
“生態(tài)養(yǎng)生”娛晚年
重要模型『一線三等角』
住進呆萌生態(tài)房
學生天地(2020年36期)2020-06-09 03:12:30
重尾非線性自回歸模型自加權M-估計的漸近分布
生態(tài)之旅
3D打印中的模型分割與打包
子午工程
太空探索(2016年6期)2016-07-10 12:09:06
工程
生態(tài)
領導文萃(2015年4期)2015-02-28 09:19:05
主站蜘蛛池模板: 美女高潮全身流白浆福利区| 亚洲男人天堂网址| 一本大道AV人久久综合| 久草国产在线观看| 激情综合网址| 中文无码毛片又爽又刺激| 99资源在线| 欧美激情视频二区| 69综合网| 亚洲Av激情网五月天| 九色免费视频| 国产一区二区免费播放| 成色7777精品在线| 国产成人亚洲精品蜜芽影院| 国产欧美视频在线观看| 波多野结衣久久高清免费| 亚洲午夜久久久精品电影院| 亚洲国产精品成人久久综合影院| 久久天天躁狠狠躁夜夜2020一| 国产成人精品在线| 国国产a国产片免费麻豆| 欧美日韩导航| 91网址在线播放| 无码啪啪精品天堂浪潮av| 国产视频a| 99精品伊人久久久大香线蕉| 成人综合在线观看| 3D动漫精品啪啪一区二区下载| 欧美69视频在线| 又污又黄又无遮挡网站| 久久91精品牛牛| 日日拍夜夜操| 久久国产精品麻豆系列| 露脸一二三区国语对白| 久久频这里精品99香蕉久网址| 久久久久88色偷偷| 国产91丝袜在线播放动漫 | 国产成人a毛片在线| 992Tv视频国产精品| 中国国产高清免费AV片| 麻豆国产精品一二三在线观看| 三上悠亚在线精品二区| 亚洲国产天堂在线观看| 欧美福利在线播放| 成人年鲁鲁在线观看视频| 欧美福利在线播放| 欧美午夜网| 免费在线视频a| 国产三级国产精品国产普男人| 精品国产99久久| 亚洲VA中文字幕| 波多野结衣国产精品| 国产免费黄| 亚洲天堂高清| 亚瑟天堂久久一区二区影院| 尤物成AV人片在线观看| 高清无码不卡视频| 国产在线观看91精品| 一本色道久久88| 人妻精品全国免费视频| 国产人成乱码视频免费观看| 男女精品视频| 色婷婷在线播放| 91在线播放免费不卡无毒| 日韩欧美91| 精品无码日韩国产不卡av| 黄色片中文字幕| 欧美成一级| 黄色网在线免费观看| 真人免费一级毛片一区二区| 国产91小视频在线观看 | 少妇高潮惨叫久久久久久| 日本影院一区| 四虎国产精品永久一区| 六月婷婷激情综合| 爆操波多野结衣| 免费A级毛片无码无遮挡| 无码日韩视频| 国产成人久久综合一区| 亚洲免费成人网| 波多野结衣无码视频在线观看| 亚洲中文字幕手机在线第一页|