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

基于水平井技術的農田排水數值模擬研究

2017-03-21 05:41:25王曉燕劉元晴劉振英
中國農村水利水電 2017年9期
關鍵詞:模型

何 錦,王曉燕,劉元晴,劉振英

(1.吉林大學環境與資源學院,長春 130026;2.中國地質調查局水文地質環境地質調查中心,河北 保定 071051;3.河北地質大學水資源與環境學院,石家莊 050031)

農業工程排水是改造鹽堿地的有效措施之一,常見的工程排水方式主要有明渠、豎井、暗管等。其原理均是將地下水位控制在極限蒸發深度以下,防止地下水中的鹽分受地表蒸發的影響而上升至地表,導致表土積鹽。近些年,暗管排水技術因其施工簡便、占地少、便于后期維護,使用較為廣泛[1,2],但其屬于被動排水技術,存在著排水量較小、排水效率不高、控制面積有限等缺點。隨著非定向開挖技術的逐步成熟,水平井技術逐漸受到關注,其應用領域由原來的石油領域逐步擴展到地下水資源開發、地下水污染防治等方面[3-5]。

水平井是指井的濾水管呈水平放置的集水建筑物。由于它大大提高了濾水管與含水層的接觸面積,其采排效率要比普通豎井高得多[6]。目前,國內利用水平井建立農田排水工程尚不多見,對于其排水效率和工程參數也鮮有研究。本文選擇華北濱海滄州典型試驗地塊,利用水平井野外試驗基地觀測數據,采用Hydrus與Modflow耦合模擬技術,對水平井的控制范圍及在典型排灌方案下的排水效果進行模擬分析,以定量評估利用水平井的排水效率,為今后采用水平井排水改良土壤提供技術支撐。

1 材料與方法

1.1 試驗區概況

水平井排水試驗區位于河北省滄州市北部滄縣境內,地貌特征以濱海平原為主。海拔高程5~7 m,自然坡降0.008%。當地屬于暖溫大陸性季風型氣候,多年平均降雨量為540 mm,多年平均蒸發量為1 371 mm (標準蒸發器E601)。該區表層土壤以粉沙質黏土、黏土為主;鹽漬化程度為輕度,2015年地下水位埋深1.0~1.5 m;水質類型為Cl-Na型,TDS為3~5 g/L;試驗區種植結構為夏玉米單作,基本上無灌溉。

1.2 水平井排水系統及抽水試驗觀測

試驗區場地寬280 m,長200 m,面積約為5.33 hm2。對試驗區水平井附近JC-01鉆孔巖性適當歸并概化為5層,即0~0.3 cm為壤土,0.3~1.9 m為黏性粉土,1.9~5.0 m以粉土為主,夾有薄層粉砂,5~8 m主要為黏土,8~30 m主要為較厚層粉沙,夾有少量粉土。各層顆粒分析數據見表1。場地施工水平井3眼,水平間距110 m,施工深度為9 m,水平跨距200 m,濾水管內徑為210 mm,長度100 m(見圖1),同時施工深度15 m的觀測井3眼,其中濾水管為9~15 m,主要用于監測不同距離地下水位變化。

場地內架設小型氣象觀測站(Watchdog 2900 ET,USA),獲得降雨量、氣溫、風速等氣象資料。土壤水分監測采用TDR技術(WinTrase,USA),并用烘干稱量法進行標定,監測深度從地表向下為20、40、60、80、100和120 cm,監測頻率為每4 h一次。地下水水位采用荷蘭Eijkelkamp公司的Mini-Diver監測,采樣頻率為4 h一次,數據自動記錄。水平井抽水采用虹吸法[7],排水量采用超聲波流量計(FLEXIM F601)自動記錄。

表1 試驗區典型土壤剖面顆粒分析Tab.1 Soil particle analysis of the typical profile in the test area

圖1 試驗區水平井布置(單位:m)Fig.1 Sketch map of horizontal wells arrangement in the test area

1.3 研究方法

水平井所排水量不僅包括非飽和帶垂向入滲量,還包括由水平井所在含水層四周徑流而來的水量。因此,這就需要模型具有精細刻畫飽和帶和非飽和帶水流的功能。本次研究選用美國國家鹽改中心開發 的Hydrus -1D軟件來建立飽和帶-非飽和帶水分運移剖面數值模型,模擬分析典型地層結構下淺部地下水的入滲-蒸發特征,建立潛水水位埋深與地下水入滲補給量的關系;其次以潛水面為耦合界面,采用Hydrus -1D軟件計算的潛水面水分交換量替代Modflow模型的補給源匯項,將非飽和帶剖面模型與地下水流模型耦合,建立試驗區排水效果預測模型,對試驗區不同灌、排組合條件下水平井工程排水效果進行模擬分析。

2 數值方程

2.1 非飽和帶水運動基本方程

根據野外試驗條件,將水分入滲過程概化采用一維垂向運動。

(1)水分運移方程。即:

(1)

式中:θ是體積含水率,cm3/cm3;t是時間,d;C為容水度,cm-1;K是導水度,cm/d;h是土壤壓力水頭,cm;z是位置水頭(向上為正),cm;Sa是作物根系吸水率,cm3/(cm3·d)。

(2)土壤水力特征方程。土壤水力特征方程采用van Genuchten公式來描述:

(2)

K(h)=KsSle[1-(1-S1/me)m]2

式中:θ為體積含水率,cm3/cm3;θr為殘留含水率,cm3/cm3;θs為飽和含水率,cm3/cm3;h為負壓,cm;α,n,m為控制土壤水分特征曲線形狀的參數;Ks為飽和導水率,cm/d;Se為相對飽和度。

(3)水分運移邊界條件和初始條件。初始條件:h(z,t)=h0(t),t=t0,z≥0;上邊界:h(z,0)=h0(t),t≥0,z=0;下邊界:h(z,0)=0,t>0,z=l。

2.2 飽和帶水運動基本方程

飽和帶水運動基本方程為:

(4)

式中:K為滲透系數,m/d;H為地下水位,m;w為源匯項強度,m/d;D為模擬區;μ為單位重力給水度;t為模擬時間;H0為初始地下水位,m;n為模擬區邊界外法線方向;τ為零流量邊界。

3 非飽和帶數值模擬與參數識別

3.1 初始條件及邊界條件

Hydrus模型剖面為地表以下整個非飽和帶范圍,按照實際地層劃分土壤剖面。為擬合土壤物理模型參數,非飽和帶初始含水量根據實際觀測數據賦值,上邊界為大氣邊界,降雨量由氣象站實測提供,蒸散量由氣象數據經過Penman-Monteith公式計算[8],不考慮植物蒸騰作用,下邊界為實測地下水位。

擬合數據時段為2015年5月1日-2015年12月31日,共計285 d。

3.2 模型參數

采用Hydrus-1D軟件中提供的Rosetta模型給定各層土壤水力參數初始值,然后通過試驗區不同層位實測數據進行參數擬合,利用模擬效率系數ENS和相對誤差RE來確定最優參數,最終土壤含水量擬合效果及選定參數見圖2和表2。

圖2 不同深度土壤實測體積含水量與模擬值對比Fig.2 Comparison chart of measured water and simulated values in soil at different depths

土壤類型θr/(cm3·cm-3)θs/(cm3·cm-3)α/cm-1nKs/(cm·d-1)NESRE/%壤土0.0780.4300.0361.56024.960.826.8粉黏0.0610.3910.0241.36410.000.863.2粉土0.0430.3630.0151.46321.320.724.3黏土0.0680.3800.0081.0902.000.751.4粉沙0.0420.3840.0391.592100.000.892.2

3.3 模擬情景

采用標定后的土壤動力學參數,考慮不同的邊界條件,其中上邊界按照(典型降雨、典型降雨+灌溉)2種情形設定:

(1)天然氣象條件。按豐水年P=578.7 mm(2013年日降雨量賦值)計算。

(2)灌溉條件。除天然降水量外,每年在3、4、5、10月份的最后一天進行灌溉,每次灌溉量為50 mm,年總灌溉量為200 mm。

作物種植模式為夏玉米單作,吸水模型采用Feddes模型[9],參數采用模型默認值,作物潛在蒸散量為文獻[10]所提供公式計算,其葉面積系數參考文獻[11]給出,下邊界為人工設定定水頭邊界,設定規則為水位埋深3 m以上間隔為0.5 m,水位埋深3 m以下間隔為1.0 m。設定目的是用于模擬不同埋深情況下試驗區淺層地下水入滲、蒸發特征。模擬時段為360 d(每個月30 d),采用變時間步長剖分方式。設定初始步長為0.001 d,最小步長為0.000 1 d,最大步長為0.2 d。

4 Hydrus-1D模擬結果分析

4.1 天然氣象環境下地下水入滲-蒸發特征與地下水位關系

天然條件下地下水蒸發、入滲以及綜合入滲補給強度隨埋深的變化關系曲線見圖3。由圖3可以看出,在地下水埋深較淺時(h<1.5 m)地下水蒸發強烈,隨埋深增大地下水蒸發量急劇減少,埋深達2.5 m后蒸發量減為0;淺埋地下水降水入滲隨埋深增大有所減少,埋深達1 m后趨于穩定;綜合補給強度表現為地下水埋深小于1.5 m時以蒸發為主,1.5~3.0 m是由蒸發為主向降水入滲為主過渡,大于3 m時以降水入滲為主,蒸發量減小到0。

地下水綜合補給量不僅隨埋深變化且隨時間呈動態變化(見圖4)。地下水位埋深較淺時(h<1.5 m)綜合補給強度隨季節變化強烈,1-5月份和10-12月份主要以蒸發為主,而6-9月份受到降雨增多的影響,綜合補給強度以入滲為主,且在埋深在1.5 m左右補給量達到隨著埋深的逐步增大,綜合補給量迅速增加。當埋深較大時(h>6 m),各季節段綜合補給強度均趨于穩定,綜合補給量大小排序為R6-9月>R10-12月>R1-5月。

圖3 天然條件下地下水補給強度~埋深關系曲線Fig.3 Relation graph of groundwater recharge intensity and depth under natural conditions

圖4 天然條件下不同季節段綜合補給量~埋深關系曲線Fig.4 Relation graph of comprehensive recharge and depth in different seasons under natural conditions

4.2 典型灌溉環境下地下水入滲、蒸發特征與地下水埋深的關系

在典型灌溉環境下,地下水蒸發、入滲以及綜合入滲補給量隨埋深的變化關系見圖5。由圖5可以看出,在地下水位埋深較淺時(h<1.5 m),典型灌溉環境下綜合補給強度隨埋深明顯增大,當地下水位增大到2.0 m以下時,年綜合補給強度也趨于穩定,其值由原來天然條件下的100 mm/a增加到180 mm/a。

通過分析典型灌溉條件下不同季節段補給強度與水位埋深關系(見圖6)可以看出,當地下水埋深較淺時(h<1.5 m),由于在3-5月份和10月份增加了灌溉量,1-5月和10-12月的綜合補給量與天然條件下相比增加了50%左右,6-9月份綜合補給強度與天然條件下相比相差不大。隨著水位埋深的增大,雨季的綜合補給量逐漸減小,而其他季節補給量增速放緩。當埋深較大時(h>6 m)各季節段綜合補給強度均趨于穩定 ,綜合補給量大小排序為R1-5月>R10-12月>R6-9月,各季節段補給量均較天然條件下補給量有所增大。

圖5 灌溉條件下地下水綜合補給強度~埋深關系曲線Fig.5 Relation graph of groundwater recharge intensity and depth under irrigation conditions

圖6 灌溉條件不同季節段綜合補給強度~埋深關系曲線Fig.6 Relation graph of comprehensive recharge and depth in different seasons under irrigation conditions

5 水平井排水效果模擬

5.1 Modflow-Hydrus耦合模型

Modflow-Hydrus耦合模型的建立是為了更加精細處理非飽和帶與飽和帶水分交換量的計算問題。傳統意義上Modflow在計算潛水補給量時采用降水入滲系數法[12]、在計算潛水蒸發量時采用經驗公式法[13],它們只能粗略地刻畫地下水的補給、蒸發規律。由前文模擬得知,試驗區地下水綜合補給量與水位埋深呈現非線性規律,在地下水埋藏較淺時(h<2.5 m),隨著埋深的增加蒸發量迅速下降,超過2.5 m后地下水蒸發的變化量減小,蒸發量逐漸變為0,而補給量在到達一定深度后為常量。因此與實際補給、蒸發情況相比,采用簡單線性關系處理非飽和-飽和帶水分交換誤差較大。因此,將不同季節階段綜合補給強度隨埋深變化規律擬合成綜合補給強度~埋深分段函數(見表3),利用迭代方式生成Modflow的面狀補給源數據文件,作為飽和水流模型的上部邊界條件,從而代替原有補給和蒸發模塊的處理方法,可以更加準確地模擬淺層地下水流場的變化[14]。

5.2 初始條件及邊界條件

試驗區地下水流場采用Processing Modflow軟件進行模擬,模型有效計算范圍為東西向2 650 m,南北向2 500 m的矩形區域,以50 m矩形網格對計算區進行均勻離散,水平井附近采用10 m格距進行加密剖分。經試算該范圍四周邊界在3水平井大降深(如5 m)抽水時平均降深小于1 cm,可視為無窮遠邊界,四周邊界均定為隔水邊界。上部源匯項(即潛水的綜合補給量)采用前文總結的不同時間的分段函數賦值(見表3),初始水位條件根據2015年初地下水位給定,滲透系數K和重力給水度μd根據已有的抽水試驗數據賦值。模擬時間段為1 800 d(每年360 d,每月30 d)。水平井采用drainage模塊,根據實際抽水量按照單元格個數平均分配排水量。

表3 研究區天然氣象環境下地下水位埋深與潛水綜合補給量函數關系 mm

5.3 水平井控制范圍模擬分析

利用Modflow-Hydrus耦合模型分別對試驗區天然氣象條件與典型灌溉條件進行了模擬,每種環境又分為連續排水和間歇性排水,以分析不同年排水量和工程管理方式對單井控制范圍的影響。其中連續排水為模擬期內每天以定流量進行排水,間隙排水為每年只有3、5、10月份排水,共計90 d排水,其他時間段均不排水。

在天然氣象環境下,對單條水平井排水強度為0.39 m3/(d·m)情況下進行連續抽水模擬。模擬結果顯示(見圖7),在連續排水條件下,3~4個月以后排水井兩側200 m以內范圍,潛水埋深值由現狀平均值1.25 m,增加到2.0m以上,潛水蒸發量大大減少。由圖8可以看出,排水井兩側300 m處平均水位較初始平均水位降落達0.75 m,即平均埋深達到2.0 m以上,地下水環境可得到明顯改善,推測單井控制兩側范圍不小于300 m;且排水井兩側400 m處平均水位較初始水位降落亦能達到0.45 m;若水平井平行排布,根據疊加原理可知,疊加后的降深不小于0.8 m,單井控制兩側范圍可達800 m。

圖7 距水平井不同距離潛水位變化過程曲線Fig.7 Curve of groundwater level change process in different distance from horizontal wells

圖8 模擬穩定后距水平井不同距離潛水位動態曲線Fig.8 Dynamic curve of phreatic water level in different distance from horizontal wells when simulation stability

連續排水條件下,水平井排水效果明顯、控制范圍大,但可操作性差。通過模擬可以看出(見圖9),若采用間歇性排水方式,單水平井控制兩側400 m范圍內潛水位埋深由初始的1.25 m降落至2.33 m,中心點降落至3.56 m;第2年進一步下降至2.44 m和3.71 m,水位降落十分明顯,疏干效果較好。

表4為水平排水井分別以排水強度0.39、0.3 m3/(d·m)連續排水和0.75 m3/(d·m)間歇排水情況下的單井控制范圍(以水平井兩側平均水位降至2.0 m為準)數據。結果表明,天然環境下年排水量不同,單排水井的控制范圍不同,單井控制范圍與年排水量成正相關關系,且與連續性排水相比,間歇性排水控制范圍明顯減小。

表4 天然環境下不同年排水量單井控制范圍 m

圖9 試驗區單井間歇排水條件下不同排水期潛水位埋深曲線[排水強度為1.0 m3/(d·m)]Fig.9 Curve of phreatic water level in different drainage period of single well intermittent drainage condition in test area [Drainage strength 1.0 m3/(d·m)]

表5為典型灌溉環境下,水平排水井分別以排水強度0.54、0.45 m3/(d·m)連續排水和0.9 m3/(d·m)間歇排水情況下的單井控制范圍數據。結果表明,由于增加了田間灌溉量,導致潛水補給量增加,單井排水控制范圍減小,而間歇性排水控制范圍減小程度更大。

表5 典型灌溉環境下不同年排水量單井控制范圍 m

5.4 試驗區排水對包氣帶水分通量的影響

通過對試驗區水平井進行間歇性排水后的降雨入滲量和地下水蒸發量進行分區統計(見圖10),結果可以看出:在現狀天然環境下(h=1.25 m),試驗區年均降水入滲量和蒸發量基本保持平衡,在未進行排水情況下,包氣帶中水分通量平均是靜止的。當排水工程運行后,對降雨入滲量和蒸發量都有明顯影響。按照地下水年均蒸發量與入滲量比值大小來劃分其影響范圍。

圖10 天然氣象環境下間歇性排水時地下水埋深與蒸發量和潛水入滲量比值關系[排水強度1.0 m3/(d·m)]Fig.10 The relation graph of groundwater level and the ratio of evaporation and infiltration in intermittent drainage condition under natural meteorological environment [Drainage strength 1.0 m3/(d·m)]

(1)試驗區。排水工程運行后,在區內年均水位埋深在2.5 m以下,地下水蒸發量僅占入滲量的10%,水份通量以絕對優勢向下運行,使淺部含水層水質有很強的淡化趨勢。

(2)強烈影響區。年均水位埋深在1.75~2.5 m,蒸發量僅占入滲量30%,面積約為試驗區的4倍。該區地下水鹽環境明顯改善,水份通量向下運行占優勢,使淺部含水層水質有明顯淡化趨勢。

(3)明顯影響區。年均水位埋深在1.50~1.75 m,蒸發量約占入滲量50%,面積約為水平井試驗區的3倍。該區部分非排水季節以蒸發地下水為主,雨季和灌溉季以入滲補給地下水為主。以水文年的宏觀時間尺度看,年均水份通量總體保持向下運行,可一定程度上抑制土壤鹽漬化,使淺部含水層水質有一定的淡化趨勢。

按水分通量劃分影響分區,試驗區已施工的3眼水平井排水工程運行后,可使8倍試驗區面積(試驗區1倍、強烈影響區4倍、明顯影響區3倍) 地下水鹽運移環境得到改善,在8倍試驗區面積影響區內,水分通量總體向下運行,總蒸發量與入滲量比值不大于50%。

5.5 水平井排水的經濟效益分析

以降低示范區面積(約6 hm2)內農田地下水位至地下水臨界深度(地下埋深為2.5 m)所需的工程量為標準,確定不同排水方式所采用的井型數量以及配套工程的相關費用,并考慮管理成本綜合來計算2者的經濟成本。其中水平井施工工程參數依據前文所定,暗管工程設計參數參考前人研究成果[15],排水管埋深為2.5 m,材質為PVC打孔波紋管,間距為25 m,總鋪管長度為2 800 m。

經對比分析(見表6),試驗區水平井排水系統投資預算為3.54 萬元/hm2,傳統的暗管排水系統投資預算為1.85 萬元/hm2。

表6 暗管排水及水平井排水投資效益對比Tab.6 The Investment Benefit between subsurface drain with horizontal well drain

此結果略高于彭成山[16]等人開展暗管改堿工程投資0.59~1.46 萬元/hm2。通過比較可以看出水平井由于施工及管材費用較高,投資比暗管排水高出50%~100%,但是水平井排水具有施工簡便、水位調節能力強、后期維護方便等優點,是一種比較有推廣前景的新排水技術。

6 結 語

(1)Hydrus模型可以較好模擬試驗區土壤水分變化過程。利用率定后的模型研究了地下水位埋深與地下水綜合補給量的關系,確定了天然條件下試驗區地下水蒸發深度在2.5 m左右。

(2)以潛水面為耦合界面,建立了Hydrus-Modflow耦合模型,并基于確定的地下水蒸發深度預測了不同排灌模式下試驗區的地下水位變化情況。結果顯示,在天然氣象環境連續排水條件下,單一水平井控制范圍可達800 m,間歇排水條件下可達200 m;在典型灌溉環境下連續排水條件下可達500 m,間歇排水條件下可達100 m。單井控制范圍與年排水量成正相關關系。水平井平行排布時,水平排水井的控制范圍明顯增大。

(3)通過模擬顯示,試驗區內3條水平井在間歇抽水情況下,可使相當于試驗區8倍的土地面積中包氣帶水分通量具有下移趨勢,潛水含水層補給量增大,淺部水土環境逐步得到改善。

[1] 遲道才,程世國,張玉龍,等. 國內外暗管排水的發展現狀與動態[J].沈陽農業大學學報,2003,34(3):312-316.

[2] 陳 陽,張展羽,馮根祥,等.濱海鹽堿地暗管排水除鹽效果試驗研究[J].灌溉排水學報,2014,33(3):38-41.

[3] Zhang H. Analytical study of capture time to a horizontal well[J].Journal of hydrology,1999,217(1-2):46-54.

[4] Anggle D G. A horizontal well recovery system to capture LNAPL and affected groundwater[J].Ground Water,1994,32(5):847-848.

[5] Conger R M. A groundwater pumping application for remediation of a chlorinated hydrocarbon plume with horizontal well technology[J].Ground Water Management,1993,15:47-60.

[6] Maurer W C. Recent advances in horizontal drilling [J]. J Canadian Pet Technol,1995,34:25-33.

[7] 付 雷,何 錦,安永會,等. 滄州地區虹吸水平井開采弱滲透微咸水技術研究[J].節水灌溉,2015,10(1):83-87.

[8] Allen R G, Pereira L S, Raes D, et al. Crop evapotranspiration guidelines for computing crop water requirement[R]∥ FAO Irrigation and Drainage Paper 56. Rome, Italy, 1998.

[9] Feddes R A, Kowalik P J, Zaradny H. Simulation of field water use and crop yield[M]. New York: John Wiley and Sons, 1978.

[10] Prasad R. A linear root water uptake model[J]. Journal of Hydrology, 1988,99(3-4):297-306.

[11] 基于Hydrus-1D 模型的玉米根系吸水影響因素分析[J].農業工程學報,2011,27(sup 2):66-72.

[12] 張蔚榛. 地下水與土壤水動力學[M]. 北京:中國水利水電出版社,1997:23-24.

[13] 唐海行,蘇逸深,張和平. 潛水蒸發的實驗研究及其經驗公式的改進[J]. 水利學報,1989,(10):37-44.

[14] 王曉燕.淺層水平井技術用于水利土壤改良排水效果數值模擬研究[D].石家莊:石家莊經濟學院,2013.

[15] 彭成山,楊玉珍,鄭存虎,等. 黃河三角洲暗管改堿工程技術實驗與研究[M] 鄭州:黃河水利出版社,2006:9-10.

[16] 張月珍,張展羽,張宙云,等. 濱海鹽堿地暗管工程設計參數研究[J]. 灌溉排水學報,2011,30(4):96-99.

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 激情国产精品一区| 99精品一区二区免费视频| 亚洲动漫h| 国产剧情一区二区| 91精选国产大片| 亚洲黄色视频在线观看一区| 99激情网| 日韩大片免费观看视频播放| 99激情网| 国产成人高清在线精品| 国产精品3p视频| 中国丰满人妻无码束缚啪啪| 久青草网站| 国产精品亚洲а∨天堂免下载| 第一区免费在线观看| 亚洲h视频在线| 欧美全免费aaaaaa特黄在线| 亚洲国产91人成在线| 久久6免费视频| 欧美三級片黃色三級片黃色1| 久热中文字幕在线观看| 色丁丁毛片在线观看| 91视频青青草| 国产视频久久久久| 亚洲无码精彩视频在线观看| 亚洲综合亚洲国产尤物| 国产不卡一级毛片视频| 国产精品爽爽va在线无码观看| 91精品久久久久久无码人妻| 国产女人在线观看| V一区无码内射国产| 国产精品亚洲专区一区| 日本一区二区三区精品国产| 国产精品开放后亚洲| 中文国产成人精品久久| 在线播放真实国产乱子伦| 日韩精品免费在线视频| 国产毛片基地| 99这里只有精品在线| 99久久国产综合精品女同| 国产精品视频猛进猛出| 男人的天堂久久精品激情| 免费a级毛片视频| 粉嫩国产白浆在线观看| 欧美日韩国产成人高清视频| 日本在线国产| 亚洲狠狠婷婷综合久久久久| 国产美女一级毛片| 在线五月婷婷| 国产福利不卡视频| 日韩视频福利| 亚洲成A人V欧美综合| 日韩欧美中文字幕一本| 手机看片1024久久精品你懂的| 久久www视频| 国产精品免费电影| 色综合成人| 日韩高清在线观看不卡一区二区| 日韩精品久久无码中文字幕色欲| 亚洲国产在一区二区三区| 国产成人精品2021欧美日韩| 成年人国产视频| 国产精品一区二区国产主播| 久久综合一个色综合网| 国产精品漂亮美女在线观看| 国产在线观看一区精品| 免费女人18毛片a级毛片视频| 九色在线观看视频| 亚洲美女一区| 国产一区成人| 久久亚洲欧美综合| 成人福利在线视频免费观看| 欧美午夜视频| 久久久噜噜噜久久中文字幕色伊伊| 香蕉在线视频网站| 国产凹凸一区在线观看视频| 伊人色天堂| 亚洲中文字幕国产av| 国产精品香蕉| 拍国产真实乱人偷精品| 国产女人18毛片水真多1| 国产欧美日韩资源在线观看|