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

RZWQM2模型模擬地膜覆蓋時間對棉田氮素遷移特征和氮肥利用效率的影響

2022-02-23 08:17:58褚曉升丁奠元湯秋香
植物營養與肥料學報 2022年1期
關鍵詞:模型

夏 文,林 濤,褚曉升,丁奠元,顏 安,爾 晨,湯秋香,*

(1 新疆農業大學農學院/棉花教育部工程研究中心,新疆烏魯木齊 830052;2 新疆農業科學院經濟作物研究所,新疆烏魯木齊 830091;3 西北農林科技大學中國旱區節水農業研究院,陜西楊凌 712100;4 揚州大學水利科學與工程學院,江蘇揚州 225009;5 新疆草地修復與環境信息重點實驗室,新疆烏魯木齊 830052)

地膜覆蓋具有提高氮肥利用效率和作物產量等作用[1-2]。常年覆膜不僅加劇農田殘膜污染,還顯著增加土壤根區硝態氮的積累[3],使棉花養分吸收效率下降[4]。使用降解地膜可在一定程度上減輕污染[5],但在實際應用中存在降解過早導致保水保肥性能變弱等缺陷[6]。在保證棉花生長對土壤水熱需求的同時,適時揭膜可增加土壤透氣性,有利于有機質分解,提高土壤速效養分的含量,保證作物生長中后期對土壤養分的吸收利用[7-8]。RZWQM2 (Root Zone Water Quality Model-version 2)模型因操作簡便、模擬結果準確和代表性廣等優點,是模擬和評估農田氮素遷移轉化的核心方法[9]。李艷[10]對RZWQM2模型中土壤養分模塊進行了校正和驗證,結果表明土壤硝態氮模擬值剖面分布趨勢和實測值相同,冬小麥季土壤硝態氮含量在0—60 cm土層變化較大,而夏玉米季因大量降水致使土壤硝態氮淋洗出根層。也有研究表明,利用RZWQM2模型模擬預測不同施氮水平在不同降水年型下1 m土壤剖面硝態氮淋失量,發現大降水事件很大程度上增加了累積在農田土壤中的硝態氮的淋溶損失[11-12]。Ma等[13]通過RZWQM2評估施氮量和施氮時間對土壤氮平衡的影響,結果表明模型能夠充分模擬不同施氮量所對應的氮素損失,相關系數高于0.83,且對施氮量的響應比對施氮時間大。薛長亮等[14]利用RZWQM2模擬預測農田硝態氮流失和氨揮發,結果顯示玉米季氮素通過農田硝態氮流失和氨揮發的平均損失量分別占冬小麥-夏玉米輪作總施氮量的13.0%和17.2%。張芊等[15]利用RZWQM2模型提出減緩硝態氮淋溶的玉米施肥方案,認為輪作周期冬小麥施氮量為233 kg/hm2、玉米施氮量為212 kg/hm2時可以有效控制硝態氮的淋溶損失,且輪作氮肥偏生產力可達到23.2 kg/kg。前人的研究主要集中于利用RZWQM2模型探索不同生態條件下滿足作物生長的水肥運籌方案[16-17]。RZWQM2模型現有的覆膜模塊,未能從機理上體現對不同地膜覆蓋時間下棉田氮素平衡及植株氮養分吸收利用變化規律的模擬,限制了RZWQM2在地膜覆蓋潛在價值研究中的應用。因此,本研究依據積溫學說,借助RZWQM2模型更好地發揮地膜增產效益促進氮素高效利用,本研究利用RZWQM2模型定量評估地膜覆蓋時間對棉田的氮素遷移影響機制,進一步評估預測不同地膜覆蓋時間下土壤剖面氮素平衡及不同氣候條件下作物對氮素的吸收利用,綜合考慮地膜覆蓋時間與氮素吸收利用的關系,從而更好地運用RZWQM2模型模擬不同環境下土壤氮素運移,為棉田氮素合理利用、降解地膜生產標準和確定最佳地膜回收時期提供理論依據。

1 材料與方法

1.1 試驗區概況

2017—2018年在新疆阿克蘇阿瓦提棉花綜合試驗基地 (40°06′N,80°44′E)進行田間試驗。該區屬于溫帶大陸性干旱氣候,海拔1025 m,日照時數2679 h,≥10℃年積溫3987.7℃,無霜期211天,年均降水量46.4 mm,蒸發量2900 mm,農業生產完全依賴于灌溉。2017—2018年棉花生育期內氣溫和太陽輻射量動態變化見圖1。試驗區土壤質地為砂壤土,播前0—80 cm土層硝態氮含量如表1所示。

圖1 2017和2018年棉花生育期內氣溫和太陽輻射量Fig. 1 Temperature and solar radiation during the growth period of cotton in 2017 and 2018

表1 試驗地播前0—80 cm土層硝態氮含量(mg/kg)Table 1 NO3--N content in 0-80 cm soil profile before sowing

1.2 田間試驗設計

試驗共設置6個地膜覆蓋時間處理:40、55、70、85、100 天和全生育期覆蓋對照,分別表示為D40、D55、D70、D85、D100和CK。單因子隨機區組設計,每處理重復3次。采用無色普通PE地膜(新疆獨山子天利高新技術股份有限公司生產),一膜兩管6行的種植模式,株距11 cm。2017年4月8日播種,4月30日出苗,9月21日收獲,全生育期167 天;2018年4月15日播種,5月10日出苗,9月25日收獲,全生育期164 天。采用“一水一肥”的施肥方式,施肥量(折合純氮)為300 kg/hm2,基追比為2∶8。其他管理方式同大田。

1.3 測定項目及方法

1.3.1 氣象參數 RZWQM2模型所需氣象參數由包含083E空氣濕度傳感器、0766強制通風罩、034B風速風向傳感器、TE525雨量筒、LI190SB光合有效輻射儀、LI200X總輻射儀的高精度自動氣象站測定(Campbell Scientific Inc.,USA)。該氣象站采用10HZ頻率采集原始數據,并提供每30 min的計算均值。

1.3.2 土壤溫度 覆膜后將自動地溫計(Micro Lite 5032PU,Israel)埋入各處理寬行膜下10 cm處,其測量范圍-30℃~80℃,分辨率0.1℃。采集時間為覆膜后至棉花收獲期,采樣間隔為1 h。

1.3.3 土壤貯水量 土壤貯水量可以系統地評價不同地膜覆蓋時間下的土壤水分利用情況。80 cm深土層貯水量計算公式[18]:

式中:W為土壤貯水量(mm);Wi為第i層土壤體積含水量(cm3/cm3);Hi為第i層土層厚度(cm);n為80 cm深土層含水量測定層數。

1.3.4 土壤硝態氮 苗期至吐絮期用土鉆于每處理的第2幅膜的窄行、寬行和裸行3個位點,取0—80 cm的土樣,每20 cm為一層。采用CaCl2浸提法獲得提取液,紫外分光光度計法測定硝態氮的含量[19]。

1.3.5 棉株氮養分吸收量 每小區隨機選取有代表性的3株棉花取棉株地上部分,烘干后粉碎,過0.5 mm篩,經H2SO4-H2O2消解定容后,采用奈氏比色法測定棉株的氮含量[20]。

棉株氮養分吸收量=實收單位面積株數×成熟單株干質量×成熟單株含氮量(kg/hm2)。

1.4 模型校正與驗證

RZWQM2 模型輸入的初始數據包括氣象數據(每日最高溫度、最低溫度、平均相對濕度、風速、降水和凈輻射量)、試驗田塊的物理性質(土壤質地、容重、飽和含水率和萎焉系數等)、土壤化學性質(硝態氮含量、銨態氮含量、土壤有機質等)和試驗中田間管理措施(播種時間、播種量、收獲時間、灌溉處理、施肥處理和耕作處理等)。

本研究利用2017—2018年全生育期覆蓋處理土壤貯水量、硝態氮含量及植株氮養分吸收量的實測數據依次校正RZWQM2模型的水分模塊和養分模塊。用5個覆蓋時期處理的數據對模型進行驗證。分別用統計檢驗中的均方根誤差(root mean square error,RMSE)和平均相對誤差(mean relative error,MRE)評估模型的可靠性。RMSE和MRE的值越小表明模擬值與實測值的差異越小,當RMSE達到最小值為優,MRE趨近于0模擬效果為優,其計算如式 (2)~(3)[21]。

式中:N為模擬值或實測值個數;Pi和Qi分別為第i個模擬值和實測值。

1.5 模型應用

借助校正與驗證后的RZWQM2模型,將地膜覆蓋時間55~105 天細化到每3 天為一間隔,模擬計算氮素吸收量。用此模型計算每個地膜處理施肥和未施肥下的棉花氮素吸收量,按照下列公式,計算不同地膜覆蓋時間處理下的氮肥利用效率。

氮肥農學效率(kg/kg)=(施氮區籽棉產量-不施氮區籽棉產量)/施氮量;氮生理利用率(kg/kg)=(施氮區籽棉產量-不施氮區籽棉產量)/(施氮區吸收氮量-不施氮區吸收氮量);氮肥回收率(%)=施氮區吸收氮量-不施氮區吸收氮量)/施氮量×100[22]。

1.6 數據分析

采用DPS Version 7.05 軟件和Excel 2016軟件進行數據的整理與分析,采用最小顯著性差異(LSD)法進行顯著性測驗。

2 結果與分析

2.1 模型建立與校正

2.1.1 水分模塊的建立與校正 利用2017和2018年全生育期覆蓋處理下,土壤貯水量的實測數據調整模型的土壤水分參數。校正后的土壤貯水量的模擬值與實測值的變化趨勢一致(圖2),土壤貯水量模擬值與實測值的MRE值的變化范圍為5.67%~7.42%,RMSE值變化范圍為1.41~1.51 mm,達到允許范圍。

圖2 棉花不同生長天數土壤貯水量實測值與模擬值Fig. 2 Measured and simulated soil water storage quantity at different growth days of cotton

2.1.2 養分模塊的建立與校正 利用棉花主要生育時期實際測定的0—80 cm土層土壤硝態氮含量調整模型參數,直至模擬值與實測值之間的MRE值范圍在0.12%~0.37%,RMSE值變化范圍為2.72~9.47 mg/kg,0—40 cm土壤硝態氮含量實測值與模擬值的離散程度大于40—80 cm相應指標,模擬值大于實測值(圖3)。

圖3 采用全生育期覆蓋對照處理對土壤硝態氮含量實測值進行的模擬值校正Fig. 3 Calibration of simulated values using measured soil nitrate nitrogen content under film mulching across the growing period of cotton (CK)

圖4顯示,2017和2018年全生育期地膜覆蓋條件下,棉株氮養分吸收量的模擬值與實測值之間MRE的變化范圍為4.82%~6.52%,RMSE值范圍為7.69~12.13 kg/hm2,均在合理誤差范圍內。

2.2 模型驗證

表2顯示,2017和2018年各處理土壤貯水量的模擬值與實測值的MRE分別為6.61%~20.08%和6.74%~16.17%,RMSE值在1.14~2.87和1.38~2.76 mm,總體來看,以D100處理的模擬效果最好。

表2 不同地膜覆蓋時間下土壤貯水量實測值與模擬值(mm)的擬合度Table 2 Agreement of measured and simulated soil water storage as affected by film mulching days

圖4 2017年和2018年地膜覆蓋出苗后棉株不同生長天數氮素吸收量的模擬值與實測值Fig. 4 The simulated and measured N uptake amount of cotton at different growing days under film mulching in 2017 and 2018

由表3可知,D40、D55、D70、D85和D100處理的MRE值變化范圍分別為0.17%~0.44%、0.11%~0.41%、0.10%~0.33%、0.06%~0.31%和0.07%~0.57%,RMSE值分別為4.81~16.38、3.74~13.45、2.73~8.64、2.05~8.58 和 3.75~17.28 mg/kg。從生育時期來看,花鈴期模擬結果較吐絮期差,這可能是由于花鈴期棉花處于生殖生長階段,所需的氮素營養較多,不同地膜覆蓋時間下棉花氮素吸收存在差異,而模型未能對這種差異做出很好地響應。

表3 不同地膜覆蓋時間下棉花不同生育期土壤硝態氮模擬值與實測值的擬合程度Table 3 The agreement between simulated and measured soil nitrate nitrogen contents at different growth stages of cotton as affected by film mulching days

棉株氮吸收量的模擬值與實測值吻合度均較高(表4),地膜覆蓋時間對棉花氮素積累量影響顯著,不同地膜覆蓋時間下棉花氮素積累量隨覆膜時間的增加而增加,在出苗后100天達到最大值。驗證結果表明棉株氮吸收量的MRE值為1.95%~8.39%,RMSE 值在 4.39~24.31 kg/hm2。

表4 不同地膜覆蓋時間下棉花氮素吸收量模擬值與實測值(kg/hm2)Table 4 The simulatedand measured Nabsorptionof cotton as affected by film mulching days

綜上所述,RZWQM2模型可以較準確地模擬不同地膜覆蓋時間下棉田土壤貯水量、土壤硝態氮及棉株氮養分吸收量,模擬結果總體趨勢在可信范圍內。

2.3 基于模擬值的棉田氮素去向分析

利用校正后的RZWQM2模型,從地膜覆蓋55天開始,每3天為一間隔,模擬計算了棉田0—80 cm土壤剖面氮素平衡(表5)。表5表明,2年氮素礦化量和固定量隨地膜覆蓋時間的增加而增加,其峰值分別為91~93和97~99天,但礦化量的增加量顯著高于固定量。隨著地膜覆蓋時間的延長,氮素總損失量(氨揮發+反硝化+滲漏)隨地膜覆蓋時間的增加而減少,其中土壤氨揮發和反硝化氮量無明顯差異,滲漏量隨地膜覆蓋時間的增加而減少,因此,滲漏是氮素損失的主要途徑,其在總氮損失中的占比最大。 2017—2018年氮素滲漏最小值分別在地膜覆蓋94~96和97~99天,較地膜覆蓋55~57天分別減少26.81%和44.65%。從減少氮素滲漏的角度看,地膜覆蓋時間以94~96和97~99天為宜。

表5 地膜覆蓋55~105天棉田0—80 cm土壤剖面氮素平衡Table 5 Nitrogen balance of soil profile in 0-80 cm cotton field after 55-105 days of film mulching

2年棉花氮素吸收量峰值分別出現在94~96和97~99天,分別為383.81和392.51 kg/hm2。綜上所述,地膜覆蓋時間過長,不利于棉株在不同生育時期氮養分的吸收利用,地膜覆蓋時間過短,由于氮素凈礦化較少而氮素損失較多造成供植株吸收利用的氮素較少。

2.4 基于模擬值計算的棉花氮素利用率

利用校正后的RZWQM2模型,從地膜覆蓋55天開始,每3天為一間隔,模擬計算植株氮養分吸收量和籽棉產量,計算農學效率、生理利用率和回收率(表6)。農學效率、生理利用率和回收率呈先上升后下降趨勢。

2017年,農學效率和生理利用率隨地膜覆蓋時間增加而增加的快速增長期分別為地膜覆蓋67~93和58~96天,分別較地膜覆蓋55~57天提高126.19%和69.25%,隨后農學效率和生理利用率均出現緩慢下降,至地膜覆蓋103~105天分別下降6.15%和3.28%。回收率隨地膜覆蓋時間增加呈平緩增加的趨勢,達到峰值后回收率較地膜覆蓋55~57天提高22.23%,地膜覆蓋97天后回收率隨之下降,至地膜覆蓋103~105天回收率下降5.12%。2018年,農學效率、生理利用率和回收率隨地膜覆蓋時間增加而增加的快速增長期分別為地膜覆蓋67~96、79~96和70~99天,分別較地膜覆蓋55~57天提高66.67%、113.20%和18.83%,隨后呈緩慢下降趨勢,至地膜覆蓋103~105天分別下降9.71%、21.43%和0.88%。綜上所述,適當的地膜覆蓋時長有利于提高氮素的農學效率、生理利用率和回收率,但過長的地膜覆蓋時間反而降低棉花的農學效率和生理利用率。基于2017和2018年度數據模擬的最佳覆蓋時間為94~96天。

表6 基于RZWQM2模型模擬的不同地膜覆蓋時間棉花氮素利用率Table 6 Nitrogen utilization efficiency based on simulated N uptake with RZWQM2 model under different film mulching time

2.5 基于模型對適宜地膜覆蓋時間的評估

以2019和2020年氣象數據為基礎,利用校正后的RZWQM2模型模擬植株氮素吸收量、氮素農學效率、氮素生理利用率和氮素回收率隨地膜覆蓋時間的變化(圖5)。2019年植株氮素吸收量隨地膜覆蓋時間的變化均呈先上升后下降的趨勢,在地膜覆蓋100 天達到峰值,隨后降低。氮素農學效率和氮素生理利用率在覆膜40 天最低,在地膜覆蓋50~70天緩慢上升,至85~90天達到峰值,之后緩慢下降,2020年變化幅度高于2019年。氮素回收率在地膜覆蓋40天最低,在地膜覆蓋50~80天緩慢上升,至85~90天達到峰值,之后緩慢下降,2019年變化幅度高于2020年。由此可以看出,RZWQM2模型模擬地膜覆蓋時間與棉花氮素利用效率的關系充分考慮了氣象因素的影響,評估結果更加客觀可靠。

圖5 RZWQM2模型模擬的植株氮素吸收及氮肥效率隨棉花生長天數的變化Fig. 5 Simulated changes in cotton N uptake and N use efficiency across the growing days using RZWQM2 model

3 討論

3.1 RZWQM2模擬地膜覆蓋時間的適用性

本研究運用RZWQM2模型對不同地膜覆蓋時間下新疆棉田土壤氮素遷移及吸收利用進行了模擬和預測。通過依次對水分模塊和養分模塊參數進行校正,其2017和2018年對不同地膜覆蓋時間下全生育期土壤貯水量的模擬結果和實測結果之間的MRE和RMSE值分別在6.61%~20.08%、6.74%~16.17%和1.14~2.87、1.38~2.76 mm,變化誤差處于可接受的范圍內,表明了模擬結果與實測結果有較為理想的擬合度。模擬數據的變化也與棉花生育期對水分和養分的需求規律相吻合。如出苗后60天是棉花對水分需求較大的時期,此時氣溫升高,土壤水分蒸發量也較大,故80 cm土層內貯水量變化較大,該時期應為棉田頭水灌溉的時間[23]。模型對不同地膜覆蓋時間下土壤剖面硝態氮模擬效果有一定的偏差,MRE和RMSE值變化范圍分別為0.06%~0.57%和2.05~17.28 mg/kg,可能因為借助RZWQM2模型模擬氮素的運移和轉化對灌溉和降水等環境因素的響應較弱[24]。不同地膜覆蓋時間下棉株氮素吸收量模擬值與實測值吻合度較高,隨著地膜覆蓋時間的增加棉株氮素吸收量呈先增大后減小的趨勢,且D100處理在整個棉花生育過程中氮素吸收量最大,這與Wang等[25]研究結果一致,適宜的地膜覆蓋時間可保證作物生長中后期對土壤養分充分吸收利用,提高養分積累量,但地膜覆蓋時間過長也不利于棉花對氮素的吸收利用。模擬結果印證RZWQM2模型已經可以應用于新疆地區不同地膜覆蓋時間下氮素遷移轉化規律研究,為氮素高效利用提供依據。

3.2 基于RZWQM2的棉田氮素變化過程對不同地膜覆蓋時間的響應

影響土壤氮素礦化的最主要環境因子是土壤溫度和水分[26]。地膜覆蓋下可增加土壤溫度和保持水分,改善土壤微生物呼吸,提高土壤氮素的礦化速率,使活性有機氮庫下降[8],進而對土壤氮素有效性產生影響[27]。也有研究表明,地膜覆蓋通過提高土壤微生物的活性,進而加快土壤氮素礦化進程[28]。筆者利用驗證后的模型將地膜覆蓋時間處理以3天為間隔進行細化,模擬研究發現土壤氮素礦化隨地膜覆蓋時間增加總體呈上升趨勢,在地膜覆蓋時間超過99天后有少量下降,其原因可能是地膜覆蓋時間未超過99天時,對土壤水分和溫度均產生了正效應,地膜覆蓋大于99天后,植株耗水量增加,土壤水分下降,而溫度仍在逐漸升高,在綜合影響下土壤氮素礦化速率下降。

氮素的去向主要分為氨揮發、反硝化、滲漏和作物吸收4個方面,其中土壤氮素滲漏的損失取決于土壤剖面上部氮的荷載量與土壤剖面內下行水流的強弱。Liu等[29]研究表明,地膜覆蓋可降低氮素的淋溶損失,但Hai等[30]通過對土壤無機氮平衡狀況的估算推測,地膜覆蓋可能增加表層土壤氮的淋溶損失。Berger等[31]指出,地膜覆蓋在一定程度上會影響反硝化作用,但是其氮氧化物排放對地膜覆蓋下土壤氮素損失的貢獻非常小。本研究結果表明,隨地膜覆蓋時間增加,土壤氮素滲漏情況總體呈下降趨勢,其原因可能是地膜覆蓋時間延長促進了作物對水分和氮素的有效利用,進而降低了氮素的滲漏;土壤反硝化在氮素損失總量上占比極低,這與Berger等[31]研究結論一致。適宜的地膜覆蓋時間明顯增加棉株對氮素的吸收,但該吸收量的來源是土壤氮還是肥料氮仍需要進一步研究。綜上,合理的地膜覆蓋時間可以促進氮肥向植物可吸收氮進行轉化,減少氮素損失。

3.3 RZWQM2模型對不同年份氮素吸收利用分析

利用2017和2018年田間試驗對模型進行校正和驗證,確定RZWQM2模型可較好的模擬不同地膜覆蓋時間下氮素吸收量及氮素利用效率,以2019—2020年的氣候條件模擬不同地膜覆蓋時間下氮素吸收利用情況。4年間氮素吸收量和氮素利用率隨地膜覆蓋時間的增加總體呈上升趨勢,地膜覆蓋時間超過100 天后略有下降,這與多數大田試驗報道相同[32-34]。但有學者研究表明,地膜覆蓋主要提高了作物對土壤本身養分的吸收,而不是來自當季施用的肥料[35],其原因可能是地膜覆蓋下由于水熱條件的改變,作物會優先吸收原本已經存在的比較穩定的有機氮。本試驗條件下,在利用RZWQM2模型模擬2019和2020年棉花氮素利用時發現氮素農學效率和氮生理利用率均為2020年變化幅度遠大于2019年,而氮肥回收率恰恰相反。比較2年氣象數據,2020年較2019年氣溫高、降水少,筆者認為在高溫缺水的條件下,地膜覆蓋時間的增加使保墑效果更明顯,有較高的水氮利用率,與不施氮條件下的產量差異較大。然而在該條件下棉花氮吸收量差額降低,造成氮肥回收率變化幅度較小。要精確模擬預測棉株氮素吸收利用對地膜覆蓋時間的響應情況,應將長時間序列的氣象數據與不同栽培措施結合起來,在今后的研究中通過田間試驗進一步探究。

4 結論

經過校正與驗證,根區水質模型RZWQM2模擬的土壤氮素動態和作物氮素吸收利用值與實測值的均方根誤差(RMSE)和平均相對誤差(MRE)值均在允許范圍之內。利用該模型模擬的棉花氮素吸收量和氮肥利用率隨地膜覆蓋時間變化的差異反映了氣象和土壤因素的影響,因而用該模型評估適宜新疆棉花生產的地膜覆蓋時間,是一種客觀和可靠的方法。綜合考慮,在新疆地區地膜覆蓋時間達94~99 d時可提高植株氮素的吸收利用能力,減少資源浪費,減輕殘膜污染。

猜你喜歡
模型
一半模型
一種去中心化的域名服務本地化模型
適用于BDS-3 PPP的隨機模型
提煉模型 突破難點
函數模型及應用
p150Glued在帕金森病模型中的表達及分布
函數模型及應用
重要模型『一線三等角』
重尾非線性自回歸模型自加權M-估計的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: 不卡无码网| 538精品在线观看| 91av成人日本不卡三区| 2021国产精品自拍| 99热这里只有免费国产精品| 日本黄色不卡视频| 久久黄色免费电影| 特级做a爰片毛片免费69| 69精品在线观看| 国产农村精品一级毛片视频| 午夜限制老子影院888| 亚洲欧美一区二区三区麻豆| 一级在线毛片| 亚洲精品桃花岛av在线| 国产乱子伦无码精品小说| 日本黄色a视频| 国产91丝袜在线观看| 国产精品综合色区在线观看| 国产精品天干天干在线观看| 中文无码精品a∨在线观看| 亚洲精品少妇熟女| 任我操在线视频| 污污网站在线观看| 激情無極限的亚洲一区免费| 熟女日韩精品2区| 国产黄在线观看| 国产传媒一区二区三区四区五区| 久久香蕉国产线| 日韩东京热无码人妻| 欧美日韩在线第一页| 午夜高清国产拍精品| 国产在线欧美| 亚洲中文字幕日产无码2021| 粉嫩国产白浆在线观看| 国内精品视频| 国产人人乐人人爱| 九九线精品视频在线观看| 免费中文字幕在在线不卡| 在线播放国产一区| 精品无码国产一区二区三区AV| 欧美精品成人一区二区在线观看| 国产青青操| 亚欧美国产综合| 午夜视频www| 国产精品无码影视久久久久久久| 亚洲 成人国产| 啪啪国产视频| 国产男女XX00免费观看| 久久夜色撩人精品国产| 亚洲人成网站观看在线观看| 日韩国产精品无码一区二区三区| 国产成人综合亚洲欧美在| 波多野结衣中文字幕久久| 91精品国产综合久久不国产大片| 欧美日本激情| 天堂成人在线视频| a天堂视频| 拍国产真实乱人偷精品| 91欧美在线| 一级毛片a女人刺激视频免费| 国产剧情伊人| 国产亚洲美日韩AV中文字幕无码成人| 人人91人人澡人人妻人人爽| 暴力调教一区二区三区| 国产精品页| 中文字幕在线欧美| 制服丝袜在线视频香蕉| 91视频首页| 毛片卡一卡二| 草草影院国产第一页| 欧美黄网站免费观看| 亚洲国产天堂久久综合| 亚洲伊人天堂| 免费观看欧美性一级| 久久精品娱乐亚洲领先| 亚洲免费毛片| 亚洲综合第一区| 成人一级免费视频| 欧美在线三级| 青青操国产视频| 日韩在线1| 久久精品人人做人人爽97|