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

基于DNDC模型的南方典型稻田氮素去向及其優化管理

2022-04-14 02:38:02周王子劉原草
節水灌溉 2022年4期
關鍵詞:水稻模型

周王子,洪 軍,劉原草,董 斌

(1.武漢大學水資源與水電工程科學國家重點實驗室,武漢 430072;2.葛洲壩中固科技股份有限公司,武漢 430000)

0 引 言

水稻是我國重要的糧食作物,超過65%的人口主食水稻[1]。氮素作為水稻的主要營養元素之一,因其對水稻產量的促進作用而受到廣泛關注。稻田的氮素管理是水稻“兩高一優”生產研究的重點。不合理的氮素管理不僅無法促進水稻氮素吸收和提高產量,而且還造成氮素資源的浪費[2]。未被利用的氮素以地表徑流、深層滲漏與氣態的形式從農田中損失,導致面源污染與大氣環境惡化[3]。特別是對于氮素的反硝化氣態損失,稻田長期淹水條件為該項損失所需的厭氧環境創造了條件,其所產生的N2O是重要的溫室氣體,加劇了全球氣候變暖[4]。

稻田氮素去向的動態分析是實現氮素優化管理的基礎[5]。由于監測工作量大、試驗周期長、環境因素復雜等原因,全面涉及氮素各去向的田間試驗研究往往難以開展。近年來,生物地球化學模型已成為研究農田養分元素遷移動態的重要工具[6], 其中由Li 等[7,8]開發的DNDC (Denitrification-Decomposition)模型因其對氮素生物化學過程的良好描述、清晰的子模塊劃分以及友好的操作界面,從而在世界各地得到廣泛應用[9-11]。DNDC 模型最初用于農業生態系統中碳素與氮素氣態損失的模擬[12];隨著模型中水文與營養模塊的日益完善,該模型能夠有效反映作物對養分的吸收量[13]與養分淋失量[14],從而發展為綜合的養分平衡模型。

DNDC 模型作為有效的養分平衡分析工具,其應用目前多局限于農田中養分單項損失的模擬與驗證層面,缺乏對養分各去向特征的綜合把握。因此,本研究以2017年湖北省孝昌縣的田間試驗為基礎,在成功驗證DNDC 模型可靠性的基礎上,定量地分析了典型水稻氮素管理模式下氮素去向的動態特征,同時對不同的氮素管理模式進行情景模擬與評價。研究結果可為我國南方稻田的優化管理提供指導,以最終實現氮素高效利用與低氮污染。

1 材料與方法

1.1 研究區概況

研究區位于湖北省孝感市孝昌縣高崗村,地理位置為北緯32°23′,東經114°07′,屬亞熱帶季風氣候。全年平均氣溫為15.2°C,降雨量達1 130 mm,多年平均無霜期為246 d。土壤以水稻土與潮土為主,水稻土以潴育型水稻土為主。水稻為一季中稻,于四月上旬育秧,五月中上旬移栽,九月中旬收割。

1.2 大田觀測試驗

選擇研究區內三塊水稻田(編號為S1、S2 與S3)作為試驗觀測田(圖1),其主要施肥管理與土壤本底值如表1所示。三塊田的施肥管理存在一定區別,但其均為當地典型的施肥策略,對其采取只觀測不干預的原則。三塊田的施氮量分別為90 kg/hm2、120 kg/hm2與150 kg/hm2(以S1 到S3 的順序敘述,下同),按4∶3∶3 的結構分為基肥、分蘗肥與拔節肥施入稻田。磷肥與鉀肥均一次基施,施肥量分別為54 kg/hm2(折算為P2O5當量)與90 kg/hm2(折算為K2O當量)。在水稻種植前取得土壤樣品,根據標準方法[15]測定各項土壤本底值。

表1 試驗田的施肥管理與土壤本底值Tab.1 Fertilization management and soil properties in experimental fields

圖1 研究區內試驗田的地理位置Fig.1 Location of experimental fields in the study area

于2017年5月31日(分蘗期)至9月14日(收獲期)對涉及稻田氮素去向的數據進行觀測。其中,氮素隨地表水的流入(灌溉)與流出(徑流和滲漏)負荷通過進出水流量與水質總氮濃度計算而得。進出水流量由旋槳式流速儀測定,水樣總氮濃度采用堿性過硫酸鉀紫外分光光度法(HJ 363-2012)分析,監測與取樣頻率均為4 h 一次(灌溉或排水時段內)。

氮素滲漏損失通過田面水滲漏速率及其總氮濃度計算而得。滲漏速率由測滲筒測定[圖2(a)]。測滲筒為直徑16 cm,長度為70 cm 的圓柱形鋼管,將其垂直打入土表以下40 cm 深度,以保證其完全穿透犁底層。筒內水層深度的變化即為滲漏量,監測頻率為每天兩次。滲漏水采用自制的裝置采集[圖2(b)]。該裝置為L 形的聚乙烯管,其豎管高度70 cm,橫管長度50 cm,兩者內徑均為11 cm。橫管上部均勻布設直徑2 cm 的濾水小孔,并用紗布纏繞。將橫管埋入試驗田面以下40 cm,埋入后對裝置與土表的縫隙做防滲處理(鋪設防滲膜),以防田面水直接進入橫管內。采用隔膜泵從豎管將橫管內的滲漏水抽取上來,分析其總氮濃度。滲漏水取樣頻率為1 d 一次(施肥后7 d內)或6 d一次(非施肥影響期)。

圖2 測滲筒與滲漏水收集裝置示意圖Fig.2 Schematic diagrams of the permeameter and the leaching collection pipe

氮素的氨揮發損失采用通氣法[16]原位監測,該方法基于簡易的氨氣吸收裝置(圖3)。裝置由15 cm的聚氯乙烯圓管與兩塊等徑的海綿組成,每塊海綿充分吸收15 mL磷酸甘油,隨后固定在圓管內。低處的海綿用于吸收稻田揮發氨氣,頂部的海綿用于防止大氣中的氨被下部海綿吸收。將該裝置放置在試驗田表面,并以每天一次(施肥后7 d 內)或每3 d 一次(非施肥影響期)的頻率更換海綿。將被更換的低處海綿浸入200 mL 的氯化鉀溶液并振蕩1 h 時,采用納什試劑比色法(HJ 535-2009)測定浸提液的銨態氮濃度。氨揮發的速率可按照參考文獻[16]給出的公式計算而得。

圖3 通氣法測定氨揮發裝置Fig.3 Device for measuring ammonia volatilization by venting method

水稻收獲后對稻谷與稻草的吸氮量進行測定。在每個試驗田內隨機選擇三塊1 m2的區域,且在每塊區域內隨機割取三蔸水稻植株,將其帶回實驗室內進行考種與脫粒,隨后分別對籽粒與稻草烘干并稱重。采用過氧化氫消解-凱氏定氮法(NY/T 2017-2011)測定植物樣全氮含量。

此外,還對研究區內試驗田的田面水位、農田管理(耕作、施肥、灌溉方式與時間)以及氣象數據進行監測,以完善DNDC 模型所需的輸入數據。其中,田面水位由HOBO U20L-04 自記水位計自動監測,監測頻率為每5 min 一次;逐小時氣象數據(降水強度、氣溫、氣壓、風速、日照時間)由中國氣象數據網(http://data.cma.cn/data/weatherBk.html)獲取。

1.3 DNDC模型及其設置、驗證與情景模擬

DNDC 模型是模擬農業生態系統中碳、氮循環的生物地球化學過程模型。該模型由6 個模塊組成:土壤氣候、植物生長、有機質分解、硝化、反硝化與發酵。前3個子模塊根據氣候、土壤、植被與人類活動的生態驅動因素計算土壤環境變量(土壤溫度、水分、pH、氧化還原電位和基質濃度),后3個子模塊使用計算所得的環境變量來模擬與氣體相關的生化反應(氮氣、氮氧化物以及甲烷的氣體釋放)。具體而言,土壤氣候子模塊通過整合氣候信息和土壤熱力與水力特性,得以計算土壤溫度、濕度以及氧化還原電位。該模塊還包含一個串聯的箱式模型,用于模擬土壤剖面內的水分入滲。植物生長子模塊根據積溫、氮需求量以及水分脅迫情況,得以計算作物生物量累積及其在谷物、秸稈與根系中的分配情況。不同作物的生長曲線與生物量分配系數被儲存在作物庫文件中,可被隨時調用。該作物庫文件還包含作物固氮的經驗系數,可實現對作物固氮的粗略估計。有機質分解子模塊通過整合氣候、土壤、植物與耕作等信息,得以計算土壤有機質向銨態氮、硝態氮與溶解性有機氮的轉化率。硝化作用模塊能夠模擬硝化菌活性以及銨態氮被硝化的過程,而反硝化作用模塊能夠模擬硝態氮被還原成氮氧化物與氮氣的過程。這兩個模塊同時運行,且兩者通過“厭氧球”這一概念實現耦合。發酵模塊量化了厭氧條件下甲烷的產生、氧化和運移。每個子模塊具體實現詳見Li等[7,8]的研究。

采用DNDC(9.5 版)對三塊試驗田的氮素去向(水稻吸氮量、滲漏損失、氨揮發、反硝化)進行模擬。模擬所需的基本數據包括氣象、土壤性狀、作物性狀與田間管理數據,這些數據大多已在田間觀測試驗中獲取。模型中所需的大氣CO2濃度被設置為390 mg/L,氮沉降指數為2.16 mg/L。作物生長參數(最高生物量、生物量比例)采用模型默認值而未被率定。由于本研究對作物吸氮量的關注主要是總量大小而非動態變化,故對總量進行驗證已經能夠滿足要求。模型中有機肥施入(Manure Amendment)、地膜覆蓋(Film mulch)以及放牧(Grazing)等與作物相關的模塊均未被設置,因為試驗田均未實施此類管理措施。

模型驗證通過對比氮素的滲漏損失動態、氨揮發動態與水稻植株吸氮量的實測值與模擬值來實現。其中,水稻吸氮量的模擬效果通過其與實測值的相對誤差反映,氮素損失動態的模擬效果由決定系數(R2)與納什系數(NSE)來反映,兩者計算公式為:

式中:Xobs(i)與Xsim(i)分別為第i個觀測值與模擬值;與分別為所有觀測值與模擬值的平均值;n為觀測數據的個數。R2與NSE的范圍分別為0~1 與-∞~1,兩者數值越接近1,表明模擬值對實測值的匹配程度越好,其中NSE=0 通常被作為判斷模型是否有效的閾值。

雖然氮素的徑流損失在田間試驗也得到了觀測,但該流失動態不用于模型的驗證。這是因為DNDC模型對徑流的模擬基于各土層含水量和田間持水量之間的水勢梯度,其與稻田徑流產生的實際機理存在較大偏差(稻田徑流的產生一般是由于田面水位高于田埂或排水口);DNDC 標準模型也鮮有用于稻田氮素徑流損失的模擬。后文將就模型對稻田徑流損失考慮不足所可能造成的不利影響進行討論。

DNDC 模型擁有明確物理意義的參數,參數不存在因模型輸入數據的不同而發生變化[17],因此模型對不同情景具有強大的預測能力。為了提供具體的氮素優化管理措施,需對不同施肥量與施氮方式下的氮素各去向負荷進行情景模擬。其中,施氮量在保持施肥結構不變的情況選從50 kg/hm2每次增加25 kg/hm2至200 kg/hm2;施氮方式則考慮尿素表施(試驗田現有方式)、尿素深施(20 cm 施肥深度)、廄肥混施與秸稈還田(均占基施化肥的一半)。獲得較高稻谷吸氮量(反映稻谷產量與品質)與較低氮素流失的氮素管理即為優化管理模式。

2 結果與分析

2.1 模型驗證

DNDC 模型對氮素滲漏損失的模擬效果如圖4所示。由圖4可知,模型能夠較為準確地模擬各試驗田的氮素滲漏損失動態。各試驗田內的滲漏損失模擬值與實測值擬合程度較為理想,R2與NSE均處于較高的水平(R2=0.860 2、0.695 5、0.649 3;NSE=0.786 2、0.621 4、0.513 5),其線性相關性均達到顯著水平(P<0.001),且高于判斷模型有效的NSE=0 這一閾值。盡管S1田塊在施分蘗肥時(移栽后37 d),其峰值滲漏損失的模擬值要低于實測值(模擬值為0.130 4 kg/(hm2·d),實測值為0.175 4 kg/(hm2·d),模擬低于實測25.66%),但考慮到試驗測量誤差等因素,這一峰值模擬差距并不會顯著影響氮素去向總體特征分析。

圖4 試驗田氮素滲漏損失模擬值與實測值的比較Fig.4 Comparison between simulated and measured values of nitrogen leaching losses in experimental fields

模型對氨揮發的模擬效果如圖5所示。由圖5可知,模型能夠較為準確地模擬氨揮發動態。模擬與實測的擬合度總體達理想水平(R2=0.766 9、0.725 5、0.783 9;NSE=0.590 1、0.673 5、0.755 9)。盡管各試驗田在分蘗肥使用時(S1 田塊為移栽后37 d,S2田塊為移栽后19 d,S3田塊為移栽后39 d)的模擬值與實測值均存在一定出入,但不會顯著影響對氮素總體特征的分析。

圖5 試驗田氨揮發模擬值與實測值的比較Fig.5 Comparison between simulated and measured values of ammonia volatilization in experimental fields

模型對水稻植株吸氮量的模擬效果如表2所示。各試驗田稻谷與稻草吸氮量的模擬值均較為接近實測值,實測值與模擬值的相對誤差絕對值大多低于10%;僅S1 田塊的稻草吸氮量模擬值高于實測值11.91%,但其相對誤差也僅略高于10%,誤差處在可接受的范圍內。雖然未能對植株吸氮的動態變化進行驗證,但考慮到模型內水稻生長參數具有很強的代表性,故認為模型對水稻植株吸氮量的模擬能夠用于氮素總體去向特征的分析。

表2 稻谷與稻草吸氮量的模擬值與實測值的比較Tab.2 Comparison between simulated and measured values of nitrogen uptakes by rice and straw

由于缺乏實測數據,模型未能對氮素的反硝化損失進行驗證。但是,模型中反硝化模擬與滲漏損失、氨揮發等氮素損失項模擬相互耦合,故上述驗證已能確保模型對包括反硝化在內的氮素各去向模擬的可靠性。

2.2 氮素去向特征

DNDC 模擬的試驗田氮素損失動態如圖6所示。由圖6可知,施肥是造成稻田氮素滲漏損失的重要誘因,各試驗田的氮素滲漏損失均在施肥后迅速升高,隨后5 d 內逐漸降至低水平。三塊試驗田在施肥后5 d 內的氮素滲漏損失占總滲漏損失的比例分別高達93.68%、81.72%與84.13%,可見該時間段是造成氮素滲漏的主要時期。

圖6 試驗田氮素滲漏損失、氨揮發與反硝化動態Fig.6 Dynamics of nitrogen leaching losses,ammonia volatilization and denitrification in experimental fields

氮素的氣態損失(氨揮發與反硝化)雖然在施肥后出現一定程度的上升趨勢,但其均對基肥施用的響應較弱。相比之下,兩者對分蘗肥的響應非常劇烈,特別是氨揮發在施分蘗肥后上升至全生育期內的最高水平;此外,氮素反硝化損失在水稻落干黃熟時期(分別為移栽后第133 d、118 d 與128 d)出現急劇上升的現象,甚至在S2 與S3 田塊中,該時期的峰值達到全生育期最高水平。由此可知施肥并不是造成氮素氣態損失高峰的唯一因素,氣候與水分均是影響氮素氣態損失的潛在因素。

表3給出了由DNDC 模擬的試驗田內氮素各去向負荷(除徑流損失)。徑流損失采用實測值,雖然觀測時段未能涵蓋返青期,但據調查這一時段并未發生暴雨與人為排水,故觀測期內的徑流損失基本涵蓋了整個生育期的情況。由表3可知,稻谷在S1、S2 與S3 田塊的吸氮量分別為49.51 kg/hm2、37.50 kg/hm2與46.98 kg/hm2;由于3個田塊的施氮量隨編號排序依次遞增,可知提升施氮量并不能保證更大的稻谷吸氮量。

試驗田的氮素氣態損失量(氨揮發與反硝化損失之和分別為10.35 kg/hm2、28.28 kg/hm2與34.69 kg/hm2)均高于隨水體流失量(徑流與滲漏損失之和分別為1.99 kg/hm2、4.34 kg/hm2與6.07 kg/hm2);前者占總損失量(12.35 kg/hm2、32.62 kg/hm2與40.76 kg/hm2)的比例分別高達83.86%、86.69% 與85.11%,而后者僅占16.14%、13.31%與14.89%。此外,施氮量較大的田塊造成的各項氮素損失亦較大,特別是對于氣態損失,其受施氮的促進效應遠高于施氮對氮素隨水流失的促進效應。

稻季后的土壤氮盈余量可根據氮素總輸入量與總輸出量計算而得,其結果也在表3給出。其中,氮素的總輸入除包含施氮以及田間觀測試驗所獲取的氮沉降與灌溉輸入外,還包含DNDC模型所采用作物相關系數法計算的生物固氮量,其值(10.09 kg/hm2、8.78 kg/hm2與10.12 kg/hm2)處于前人研究的范圍[18],具有可靠性。由表3可知,S1田塊的土壤氮含量在稻季后減少;S2 與S3 田塊則存在土壤氮盈余現象,且對于施氮量更大的S3田塊,其土壤氮盈余量也更大。

表3 試驗田氮素各去向途徑在水稻生育期內的總量 kg/hm2Tab.3 The amounts of different nitrogen outputs from experimental fields during a rice season

2.3 情景模擬

采用DNDC模型對不同氮素管理模式下的氮素各去向負荷(除徑流損失)進行模擬。其中,針對施氮量逐漸增加的情景模擬結果如圖7所示。由圖7可知,在施氮量低于100 kg/hm2時,施氮量增加促進稻谷對氮素的吸收;而當施氮量高于100 kg/hm2時,稻谷吸氮量基本保持穩定。可見,100 kg/hm2為研究區施氮閾值,超過這一閾值會導致氮素利用的下降。氮素氣態損失始終隨施氮量增加而增加,特別是氨揮發隨施氮量增加的斜率也在逐漸上升,這促使高水平施氮下的氨揮發損失達到極高水平,甚至可超過稻谷吸氮量。盡管施氮在一定程度上加劇氮素滲漏損失,但其值始終低于10 kg/hm2,維持在較低的水平。

圖7 試驗田不同施氮量下氮素各去向負荷Fig.7 Nitrogen output loads under different nitrogen application rates in experimental fields

針對不同施氮方式的情景模擬結果如表4所示。由表4可知,施氮方式對稻谷吸氮的影響僅在S1 田塊中有所體現,且秸稈還田方式下的吸氮量最高;S2 與S3 田塊的稻谷吸氮在四種施氮方式下均保持一致。施氮方式對氮素滲漏損失的影響較小,其值在不同施氮方式下的差異分別僅為0.49 kg/hm2、0.29 kg/hm2與2.28 kg/hm2。然而,施氮方式對氮素氣態損失的影響非常顯著。尿素深施相比于表施方式能夠略微減少氨揮發與反硝化的氮素損失負荷;農家肥基施雖然相比于尿素表施能夠減少氨揮發,但是卻加劇了反硝化的損失量;秸稈還田能夠同時抑制氨揮發與反硝化,其對氣態損失的控制最為顯著。

表4 試驗田不同施氮方式下氮素各去向負荷kg/hm2Tab.4 Nitrogen output loads under different nitrogen application methods in experimental fields

3 討 論

以往對DNDC模型的應用大多聚焦于養分的單項模擬,如夏文建等[19]僅考慮作物吸氮與氮素氣態損失,趙崢[6]僅考慮徑流與滲漏損失。本研究則全面考慮了氮素各去向的模擬,并且用原位觀測數據對模擬成功驗證。用于驗證氮素滲漏損失的實測數據分別多達27、30 與28 個,氨揮發實測數據分別達41、44 與44 個,兩項數據均滿足驗證數據量的要求。盡管模型未對氮素的徑流損失進行模擬,但實測數據表明其損失量處于極低的水平(表3),故其并不影響對氮素各去向總體特征分析。

氣態損失是氮素的主要損失途徑,這一發現與以往的研究結論相符[20]。氣態損失不僅造成氮素資源的浪費,還導致因氨揮發引起的大氣酸化[3]以及因反硝化產生氮氧化物引起的氣候變暖[4],故對氮素氣態損失的控制具有重要意義。分蘗肥導致氨揮發的顯著增高,故該施肥時期是控制氨揮發的關鍵時期。另外,在水稻臨近收獲的黃熟期,稻田自然落干,有利于N2O 的產生[21],因此三塊試驗田黃熟期的反硝化損失都出現了顯著的增大(圖6)。

稻田中氮素隨水體的流失較少,出現這一現象的原因在于稻田的水管理與土壤條件。具體而言,稻田中常修筑田埂以適應水稻生產對水的高需求。除了在水稻中期為曬田而進行的排水外,其他主動排水與降雨導致的徑流不多;此外,水稻耕作使得稻田土壤產生密實的犁底層,其顯著抑制了土壤水分的深層滲漏。因此,稻田排水不應當是農田面源污染的主要來源,在某種程度上其可被視為凈化水質的濕地[22,23];而被重點關注的農田面源污染可能主要源于旱田而非稻田[24]。

根據情景模擬可提出本研究區內稻田的氮素優化管理模式,即控制施氮量為100 kg/hm2,并結合秸稈還田方式。氮肥施用量在高施氮水平下的進一步提升并不會促進稻谷吸氮量增加,相反其存在加劇氮素氣態損失的風險。100 kg/hm2的施氮量是使得稻谷吸氮最大的施氮水平。研究區所習慣采用的90 kg/hm2施氮量為施氮不足,采用優化施氮水平可使稻谷吸氮量提升12.67%;研究區內120 kg/hm2與150 kg/hm2的施氮量為施氮過量,經優化后可在保證稻谷吸氮量不變的情況下分別減少33.74%與65.57%的氮素氣態損失。雖然DNDC(9.5版)難以直接模擬稻谷的生物量(其僅能模擬植稻谷對碳素與氮素吸收),但是稻谷吸氮量一般與產量呈顯著的正相關關系[25],因此100 kg/hm2可被認為是使得產量最高的最優施氮水平。秸稈還田能夠顯著降低施氮量,且同時抑制氨揮發與反硝化損失。當還田氮素占基施氮素的一半時,其與未還田處理相比分別削減17.64%、34.28%與29.71%的氮素氣態損失。

4 結 論

(1)經驗證后的DNDC模型能夠準確模擬中國南方典型稻田的氮素各去向負荷,包括水稻植株吸氮、氮素滲漏、氨揮發與反硝化,可用于稻田氮素去向結構分析。

(2)氮素的氣態損失是稻田內氮素損失的主要途徑。其中,對氨揮發與反硝化的控制要特別關注分蘗期施肥與水稻黃熟落干對該氣態損失的顯著促進作用。

(3)基于DNDC模型的情景模擬可提出具體的氮素優化管理模式。研究區內的氮素優化管理為控制施氮量100 kg/hm2與秸稈還田相結合,該優化管理模式可保證稻谷吸氮量與產量最高,且顯著抑制氮素氣態損失。

猜你喜歡
水稻模型
一半模型
什么是海水稻
有了這種合成酶 水稻可以耐鹽了
今日農業(2021年21期)2021-11-26 05:07:00
水稻種植60天就能收獲啦
軍事文摘(2021年22期)2021-11-26 00:43:51
油菜可以像水稻一樣實現機插
今日農業(2021年14期)2021-10-14 08:35:40
重要模型『一線三等角』
一季水稻
文苑(2020年6期)2020-06-22 08:41:52
重尾非線性自回歸模型自加權M-估計的漸近分布
水稻花
文苑(2019年22期)2019-12-07 05:29:00
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产精品污污在线观看网站| 欧美97欧美综合色伦图| 国产一区二区精品福利| 99久久99这里只有免费的精品| 亚洲综合精品第一页| 亚洲第一区在线| 国产性爱网站| 国产小视频在线高清播放| 免费又爽又刺激高潮网址| 亚洲一级毛片在线观播放| 久久青草精品一区二区三区| 日韩A级毛片一区二区三区| 亚洲综合片| 亚洲日韩AV无码一区二区三区人| 亚洲精品视频在线观看视频| 日本欧美中文字幕精品亚洲| 欧美乱妇高清无乱码免费| 久久久久青草大香线综合精品| 国产激情第一页| 亚洲高清中文字幕在线看不卡| 亚洲精品欧美重口| 中文字幕伦视频| 国产91在线|日本| 精品欧美一区二区三区久久久| 自拍中文字幕| 免费jjzz在在线播放国产| 国产在线观看91精品亚瑟| 国产在线视频欧美亚综合| 国产成人区在线观看视频| 国产在线观看91精品| 91年精品国产福利线观看久久 | 国产玖玖视频| 色婷婷在线播放| 免费国产在线精品一区| 亚洲欧洲日产无码AV| 精品少妇人妻无码久久| 五月天久久婷婷| 亚洲色图综合在线| 久久99精品久久久久纯品| 婷婷成人综合| 国产尤物jk自慰制服喷水| 欧美一级高清片欧美国产欧美| 伊人蕉久影院| 无码专区在线观看| 制服丝袜一区二区三区在线| 制服无码网站| 国产精品无码作爱| 亚洲成a人在线播放www| 青草视频久久| 54pao国产成人免费视频| 午夜三级在线| 国产在线精品99一区不卡| 亚洲国产精品不卡在线| 亚洲av无码久久无遮挡| 91综合色区亚洲熟妇p| 精品国产网站| 一本大道香蕉高清久久| 一本久道久久综合多人| 国产高清免费午夜在线视频| 久久中文字幕2021精品| 午夜国产大片免费观看| 欧美精品一区在线看| jizz亚洲高清在线观看| 26uuu国产精品视频| 亚洲A∨无码精品午夜在线观看| 波多野结衣第一页| 国产xx在线观看| 人人91人人澡人人妻人人爽 | 欧美一区二区精品久久久| 性欧美在线| 男女男免费视频网站国产| 国产午夜在线观看视频| 国产在线自在拍91精品黑人| 国产视频a| 久久久成年黄色视频| 国产在线观看一区二区三区| 一级毛片免费高清视频| 国产草草影院18成年视频| 久久婷婷国产综合尤物精品| 亚洲成人一区二区| 亚洲 成人国产| h视频在线观看网站|