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

    用SahysMod 模型研究不同灌排管理情景土壤水鹽動態(tài)

    2020-07-22 14:37:02黃亞捷卓志清黃元仿
    農(nóng)業(yè)工程學(xué)報 2020年11期
    關(guān)鍵詞:根區(qū)荒地排水溝

    黃亞捷,李 貞,卓志清,興 安,黃元仿※

    (1.中國農(nóng)業(yè)大學(xué)土地科學(xué)與技術(shù)學(xué)院/農(nóng)業(yè)農(nóng)村部華北耕地保育重點實驗室/自然資源部農(nóng)用地質(zhì)量與監(jiān)控重點實驗室,北京 100193;2.中國農(nóng)業(yè)科學(xué)院農(nóng)業(yè)資源與農(nóng)業(yè)區(qū)劃研究所,北京 100081)

    0 引 言

    銀北灌區(qū)在西北干旱與半干旱地區(qū)農(nóng)業(yè)生產(chǎn)中占有重要地位,它是寧夏土地整治和高標(biāo)準(zhǔn)灌溉綠洲農(nóng)田建設(shè)的重點區(qū)域。灌區(qū)耕荒地交錯分布、土壤鹽漬化嚴(yán)重[1]。在灌區(qū)4 413 km2耕地中,鹽漬化耕地21 480 km2,鹽堿荒地560 km2。銀北的石嘴山市、平羅縣、惠農(nóng)區(qū)等地和銀南的部分地區(qū)耕地土壤鹽漬化十分嚴(yán)重[2]。周德[3]認(rèn)為雖然銀北地區(qū)的土壤鹽漬化面積總體呈下降趨勢,但由于不合理的灌排方式等影響,土壤鹽漬化風(fēng)險仍然很大。

    已有研究[4]表明,在土地整治過程中簡單平整造地導(dǎo)致的灌排不配套與水鹽不耦和是造成區(qū)域土壤鹽漬化的重要原因。以水為中心,如在土地整治過程中以水定地,建立完整、配套的灌排系統(tǒng),是解決土壤鹽漬化的有效措施[4]。Konukcu 等[5]也指出健全灌排系統(tǒng)、控制地下水位是防止鹽堿化的重要措施。在灌水方面,基于合理的灌水量及不同灌水方式(如畦灌、溝灌)淋洗鹽分,可加快土壤脫鹽速度,有利于降低土壤鹽漬化[6]。對于水資源約束地區(qū),如銀北灌區(qū),利用咸水、微咸水灌溉是解決水資源短缺與降低土壤鹽漬化的有效措施。在排水方面,在土地整治過程中綜合考慮排水溝的設(shè)計規(guī)格、走向和與地勢的關(guān)系,排水溝和地下水位、土層結(jié)構(gòu)的關(guān)系,排水溝的排水、排鹽效果等,通過排水溝設(shè)置控制地下水埋深,最終可實現(xiàn)降低土壤鹽漬化的目的[7-8]。但是,對于如何通過土地整治加強灌區(qū)灌排管理,減輕或避免灌區(qū)土壤鹽漬化仍是灌區(qū)面臨的主要問題之一。

    通過模型模擬及預(yù)測,研究在土地整治過程中不同灌排管理下區(qū)域土壤水鹽運動規(guī)律和運行特征,對實現(xiàn)土壤鹽漬化的預(yù)測與防治具有重要的指導(dǎo)意義[5,9-11]。遺憾的是,區(qū)域內(nèi)土壤鹽漬化過程與水鹽運移規(guī)律較復(fù)雜,土壤水鹽具有明顯的空間變異性,以往對于灌排管理及鹽漬化治理的研究多集中于實驗室和田間小區(qū)試驗[5,12],而在土地整治過程中綜合考慮區(qū)域耕地及荒地土壤水鹽空間變異性的灌排管理研究還相對較少;同時,常用的模型(如HYDRUS 模型等)只能模擬情景比較簡單的水分溶質(zhì)運移,不能在同一構(gòu)建模型中給定不同類型的大氣邊界條件,這就導(dǎo)致不同景觀單元(如耕地與荒地)的土壤屬性輸入不能區(qū)分,也不能同時考慮不同的作物類型,根系不能跨區(qū)[13];除此之外,以往的研究較少關(guān)注在土地整治過程中不同灌排管理情景下耕地與荒地土壤水鹽在未來的長期動態(tài)變化。

    SahysMod(Spatial Agro-hydro-salinity Model)是荷蘭土地開墾和改良國際研究所(ILRI)的 Boonstra 和Oosterbaan 教授以水分和鹽分平衡原理為基礎(chǔ),集成SaltMod(Soil Salinity Model)和 SGMP(Standard Groundwater Model Program)形成的三維模型[14]。被逐步應(yīng)用于模擬與預(yù)測土壤水(鹽)分、地下水和排水的鹽分、地下水埋深、排水量等方面[15-17]。SahysMod 在模擬及預(yù)測未來土壤水鹽時充分考慮區(qū)域土壤屬性的空間變異性,以及作物種植制度的不同所引起的灌溉與排水的差異性[18]。因此,本研究基于SahysMod 模型同時考慮耕地及荒地土壤水鹽的空間變異性,采用平衡方程研究土地整治片區(qū)未來10 a 內(nèi)不同灌排管理情景下區(qū)域土壤水鹽動態(tài)變化,以期通過優(yōu)化灌排管理模式緩解灌區(qū)土壤鹽漬化。

    1 材料與方法

    1.1 研究區(qū)概況

    寧夏平羅縣西大灘(38°49′25"N,106°25′54"E)是銀北灌區(qū)典型的土壤鹽漬化分布區(qū)及土地整治片區(qū)(圖1),該區(qū)域?qū)儋R蘭山東麓洪積平原與黃河沖擊平原過渡地帶,面積為6.89 km2,其中耕地、荒地面積分別為4.90、1.99 km2。區(qū)域內(nèi)土壤膠體吸附大量的Na+,土壤堿化度在15%~60%之間,pH 值為8.0~10.4。區(qū)內(nèi)地勢平坦,荒地地勢低洼,高程值為1 077~1 122 m。

    圖1 研究區(qū)位置及采樣點空間位置分布圖 Fig.1 Location of study area and spatial distribution of soil sampling sites

    研究區(qū)屬于黃河中上游灌溉地區(qū),隸屬于第三支溝,田間農(nóng)渠相鄰或相間布置。區(qū)域主要依靠于黃河水、地下井抽水及排水再利用進(jìn)行灌溉,具體的灌溉制度見表 1,總灌水量約為670 mm。由于灌水所攜帶的一定數(shù)量的可溶鹽在灌區(qū)積累(灌溉水電導(dǎo)率約為1.05 dS/m),土壤含鹽量增高,部分荒地的土壤鹽漬化更為嚴(yán)重。研究區(qū)田間排水主要依靠農(nóng)溝排水,由表2 所示,斗排主要接收來自田間排水溝的排水;農(nóng)溝多為1.5 m 深,部分農(nóng)溝1.8 m 深。區(qū)域內(nèi)的排水標(biāo)準(zhǔn)一般為1.5 mm/d。區(qū)內(nèi)排水設(shè)施較為落后,溝道淤積、排水不暢等也加重了土壤鹽漬化。

    表1 研究區(qū)基本灌溉制度 Table 1 Basic irrigation regime in study area

    表2 排水系統(tǒng)的基本特征 Table 2 Characteristics of drainage system

    綜合考慮作物種植制度、灌排系統(tǒng)及高程等因素,本研究以240 m 為采樣間距,共布設(shè)70 個監(jiān)測樣點,具體布點如圖1 所示。土壤樣本采集在每次灌水前,在每個采樣點用土鉆分層取樣,每個樣品是由每個布設(shè)樣點的鄰近直徑10 m 范圍內(nèi)3 個土樣混合而成。采樣層次分別為0~10、>10~20、>20~30、>30~40、>40~50 cm,采樣時間為2015 年5 月、10 月,2016 年5 月、10 月及2017 年5 月?;?~10、>10~20、>20~30 cm 的土壤樣品混合形成0~30 cm 深度的土樣;基于30~40、>40~50 cm 深度的樣品混合形成30~50 cm 的土樣。參照[19-21]測定土壤含水率、容重、田間持水量、機械組成、電導(dǎo)率、土壤與地下水全鹽含量、地下水礦化度。

    1.2 研究方法

    1.2.1 SahysMod 模型介紹

    SahysMod 模型是通過多邊形網(wǎng)格(包括內(nèi)部及外部多邊形網(wǎng)格)對區(qū)域內(nèi)土壤水分、鹽分空間變異進(jìn)行劃分的3D 平衡模型[22]。它包括水平衡模塊、地下水模塊及鹽平衡模塊。根據(jù)當(dāng)?shù)氐臍夂驐l件、作物生長,SahysMod模型可分為1~4 個模擬季節(jié),并從垂直方向上分為4 層研究水鹽平衡,即地表、根區(qū)、過渡層和含水層。水平衡方程如下所示

    式(1)是地表水平衡方程。Pp為降雨量,Ig為灌水量,λ0為從根區(qū)進(jìn)入地表的水量,E0為地表水蒸發(fā)量,λi為從地表進(jìn)入根區(qū)的水量,S0為地表徑流量,ΔWs為儲存在地表的水量變化量。式(2)是根區(qū)土壤水平衡方程。Rr為進(jìn)入根區(qū)的毛細(xì)管上升水,Era為蒸騰耗水量,Lr為根區(qū)滲漏水量,ΔWf為根區(qū)在非飽和態(tài)儲存的水量。當(dāng)根區(qū)水量處于田間持水量與完全飽和之間時,ΔWr為此狀態(tài)下的儲水量。式(3)是過渡層水平衡方程。Lc為灌溉渠系滲漏水,Vr從含水層垂直進(jìn)入過渡層的水量,VL從過渡層進(jìn)入含水層的水量,Gd為通過排水溝及管道等排出的水量,△Wx為儲存在過渡層的水量變化量。式(4)是含水層水平衡方程。Gi為含水層的進(jìn)水量,G0為含水層的出水量,Gw為地下抽水量,ΔWq為儲存在含水層的水量變化量。式中的變量單位均為mm。

    鹽平衡方程是基于上述各個層次的水平衡方程及其所攜帶的鹽分含量建立。地下水流動是基于有限差分法確定。詳細(xì)的SahysMod 模型對區(qū)域土壤水鹽平衡方程及地下水流動的計算可參考以往研究[14,23-24]。

    基于SahysMod 模型確定研究區(qū)垂直方向上的各層(即地表、根區(qū)、過渡層和含水層)厚度?;赟ahysMod Working Group of ILRI[23-24]對地表厚度的定義,認(rèn)為區(qū)域內(nèi)地表層厚度為0。如圖1 所示,研究區(qū)內(nèi)主要作物類型為玉米,玉米根系平均長度為0.3 m[25]。SahysMod Working Group of ILRI[23-24]認(rèn)為根系平均長度即為根區(qū)厚度,因此確定根區(qū)厚度為0.3 m;過渡層厚度是基于最大地下水埋深確定。根據(jù)以往研究[26],研究區(qū)內(nèi)非生育期內(nèi)的耕地及荒地最大地下水埋深分別為 2.0 及1.8 m,耕地與荒地平均高度差為0.2 m。將最大地下水埋深減去根區(qū)厚度后,確定耕地與荒地的過渡層厚度分別為1.7 及1.5 m;含水層只考慮到潛水含水層,厚度為20 m[27-28]。

    1.2.2 SahysMod 模型參數(shù)設(shè)置

    考慮到土地平整、高程、作物類型以及灌排設(shè)施等差異對土壤鹽分的影響,在SahysMod 模型中共設(shè)置138個多邊形網(wǎng)格,其中內(nèi)部多邊形網(wǎng)格98 個(對應(yīng)網(wǎng)格編號1-98)、外部多邊形網(wǎng)格40 個(對應(yīng)網(wǎng)格編號99-138),每個網(wǎng)格的面積為270×270 m2(圖2)。SahysMod 模型認(rèn)為每個多邊形網(wǎng)格均是1 個均質(zhì)的單元,主要考慮不同網(wǎng)格之前的土壤水鹽運動。內(nèi)部多邊形網(wǎng)格為研究區(qū)域,每個內(nèi)部多邊形網(wǎng)格參數(shù)一致,如果在同1 個網(wǎng)格內(nèi)存在2 種或2 種以上的作物類型(如玉米與水稻),在計算網(wǎng)格內(nèi)的參數(shù)時(如灌水量)是基于不同作物的面積所占比例進(jìn)行確定。外部多邊形所處的位置為研究區(qū)邊界,由圖1 可知,研究區(qū)邊界為封閉灌渠及排水溝,根據(jù)以往研究認(rèn)為此條件下的研究區(qū)外邊界條件為定水頭邊界條件[29]。

    SahysMod 模型中比例尺設(shè)定為1∶20 000,預(yù)測周期為未來10 a。本研究將每年分為2 個模擬時期,即生育期(5—9 月,共5 個月)和非生育期(10 月—次年4 月,共7 個月)?;赑enman—Monteith 公式[30]、FAO 在灌溉與排水中提出的作物系數(shù)[31]及銀北灌區(qū)潛水蒸發(fā)量計算方法[32],區(qū)域過去10 a(2006—2016 年)中生育期平均降雨量為185 mm,平均潛水蒸發(fā)量為917 mm,非生育期的平均降雨量25 mm,平均潛水蒸發(fā)量為585 mm。由于研究區(qū)域地勢平坦、徑流極其微弱,可認(rèn)為地表徑流為0[33]。研究區(qū)地下水埋深及含鹽量數(shù)據(jù)由地下水觀測井常年觀測獲取。值得注意的是,采用SahysMod 模型計算的電導(dǎo)率為田間土壤飽和電導(dǎo)率,即EC 值[23]。本研究基于土壤飽和浸提液的電導(dǎo)率ECe 與田間土壤飽和條件下電導(dǎo)率EC 的關(guān)系(EC=2ECe)[23],將SahysMod 模型中輸出的EC 值均已換算為ECe 值。下文所提及的土壤電導(dǎo)率均指ECe。

    圖2 SahysMod 模型中多邊形網(wǎng)格設(shè)置 Fig.2 Nodal network dividing the experimental sites for SahysMod model

    SahysMod 主要輸入?yún)?shù)包括氣象、土壤、作物、地下水、灌溉以及排水等,主要輸出數(shù)據(jù)包括土壤鹽分、排水和地下水的礦化度、地下水埋深、排水量等。區(qū)域中各季度的氣象數(shù)據(jù)、土壤鹽分、作物類型、灌排水、地下水埋深以及礦化度等基礎(chǔ)參數(shù)通過實際監(jiān)測及有關(guān)文獻(xiàn)的取值范圍獲取,部分中間過程參數(shù)值采用模型默認(rèn)值,具體的數(shù)據(jù)來源及參數(shù)值見表3 和表4。

    表3 研究區(qū)主要數(shù)據(jù)來源 Table 3 Source of main data in study area

    考慮空間變異性的參數(shù)主要包括根區(qū)初始土壤鹽分、過渡層初始土壤鹽分、根區(qū)總孔隙度、根區(qū)有效孔隙度、根區(qū)土壤容重、根區(qū)田間持水量、根區(qū)淋洗效率,各參數(shù)值見圖3 所示。其中,在SahysMod 模型中,基于ECe 與土壤鹽分的強正相關(guān)性,采用ECe 反映土壤鹽漬化情況[23-24]。用2015—2016 年土壤鹽分?jǐn)?shù)據(jù)進(jìn)行模型參數(shù)率定,2017 年土壤鹽分?jǐn)?shù)據(jù)進(jìn)行驗證。

    1.2.3 SahysMod 模型參數(shù)敏感性分析

    通過敏感性分析識別不同輸入?yún)?shù)對SahysMod 模型中土壤水鹽的影響,在進(jìn)行敏感性分析的參數(shù)主要包括根區(qū)淋洗效率(Flr)、過渡層淋洗效率(Flx)、含水層 淋洗效率(Flq)及含水層導(dǎo)水率(Kaq)。其中,F(xiàn)lr 為根區(qū)滲漏水的鹽分質(zhì)量濃度與飽和土壤水的平均鹽分質(zhì)量濃度的比值,取值范圍為0~1[18,23]。Flx 為過渡層滲漏水的鹽分質(zhì)量濃度與飽和土壤水的平均鹽分質(zhì)量濃度的比值,取值范圍為0~1[18,23]。Flq 是指從含水層滲漏出的溶液鹽分濃度與含水層飽和時的平均鹽分濃度的比值。研究表明[18],F(xiàn)lq 的取值范圍在0.01~2 之間,值越大表明淋洗效率越高。另外,F(xiàn)lr 在率定時考慮了空間變異性;由于研究區(qū)面積較小,過渡層和含水層性質(zhì)較為均一,因而對Flx、Flq 及Kaq 的率定均未考慮變異性。參考以往研究[18,37],采用參數(shù)±20%及±50%來評定其敏感性。在進(jìn)行敏感性分析時認(rèn)為SahysMod 模型中其他參數(shù)基本不變。采樣均方根誤差(Root Mean Squared Error,RMSE)指標(biāo)進(jìn)行模型參數(shù)率定評價。

    表4 SahysMod 模型參數(shù)值 Table 4 Values of parameters for SahysMod model

    圖3 SahysMod 模型中輸入?yún)?shù)值 Fig.3 Input parameters for SahysMod model

    1.2.4 灌排管理模式情景設(shè)置

    通過綜合考慮土地整治前后區(qū)域土壤水鹽運動規(guī)律,最終設(shè)定4 種情景研究在土地整治過程中不同灌排模式對土壤水鹽長期動態(tài)變化的影響(表5)。研究區(qū)排水主要為田間農(nóng)溝排水,多數(shù)深為1.5 m,部分深1.8 m,考慮到1.8 m 的農(nóng)溝占比很小,因此SahysMod 模型中認(rèn)為現(xiàn)有模式下的排水溝深為1.5 m。FAO[31]認(rèn)為當(dāng)土壤飽和浸提液電導(dǎo)率為1.7 dS/m 時,玉米產(chǎn)生鹽分脅迫導(dǎo)致減產(chǎn)。由于研究區(qū)內(nèi)主要作物類型為玉米,因此認(rèn)為當(dāng)土壤電導(dǎo)率為1.7 dS/m 時,耕地根區(qū)土壤鹽分積累達(dá)到障礙水平。同時,在SahysMod 模型模擬不同灌排管理模式對土壤鹽分的影響時,僅僅考慮耕地網(wǎng)格在生育期未來10 a 內(nèi)的土壤鹽分變化。

    表5 不同灌排管理模式情景設(shè)置 Table 5 Setting of different irrigation and drainage management scenarios

    2 結(jié)果與分析

    2.1 SahysMod 模型率定及驗證

    取不同的根區(qū)淋洗效率Flr,模擬計算根區(qū)土壤電導(dǎo)率,將根區(qū)土壤電導(dǎo)率預(yù)測值與實測值進(jìn)行比較,吻合最好的根區(qū)淋洗效率即為實際的淋洗效率。如表6 所示,在現(xiàn)有灌排管理情景下,隨機選取荒地多邊形網(wǎng)格1、26、60 及耕地網(wǎng)格7、61、76、87、89 進(jìn)行比較,當(dāng)設(shè)定Flr的范圍為0.509~0.698,不同網(wǎng)格的RMSE 值都接近于0,說明土壤電導(dǎo)率的預(yù)測值與實測值基本一致,根區(qū)淋洗效率確實存在空間變異性,如耕地網(wǎng)格編號61 的Flr 值為0.659,荒地網(wǎng)格編號26 對應(yīng)的Flr 為0.612。另外,耕地Flr 一般要大于荒地,這和以往的研究基本一致。已有研究表明Flr 的值主要取決于土壤質(zhì)地、孔隙度及灌溉條件等,土壤較黏重、土壤有效孔隙度低及用堿性灌溉水灌溉會導(dǎo)致Flr 低[18,38]。Yao 等[18]也表明濱海農(nóng)田Flr為0.60~0.85,Jia 等[39]認(rèn)為寧夏的銀南地區(qū)的根區(qū)的淋洗效率為0.55。

    表6 基于SahysMod 模型的根區(qū)淋洗效率(Flr)率定 Table 6 Calibration of leaching efficiency of root zone (Flr) for SahysMod

    如表7 所示,選取不同F(xiàn)lq 值(Flq 取0.6~1.8)比較實際地下水埋深與預(yù)測地下水埋深。結(jié)果發(fā)現(xiàn),當(dāng)Flq 取1.2 時,均方根誤差RMSE 值最小(RMSE=0.020 m),地下水埋深的模擬值與實測值吻合最好?;谌《ú煌腇lx 值(Flx 取0.4~0.9)模擬地下水埋深。結(jié)果發(fā)現(xiàn),當(dāng)Flx 為0.8 時,RMSE 值最小,為0.020 m,地下水埋深的模擬值與實測值吻合最好?;谌《ú煌琄aq模擬計算地下水埋深,結(jié)果發(fā)現(xiàn),當(dāng)Kaq 取10 m/d,RMSE 值最?。≧MSE=0.020 m),地下水埋深的模擬值與實測值吻合程度最好,這和以往的銀北灌區(qū)的含水層導(dǎo)水率的值基本相近[40-41]

    以2015—2016 年土壤電導(dǎo)率數(shù)據(jù)進(jìn)行SahysMod 模型參數(shù)率定,以2017 年數(shù)據(jù)進(jìn)行驗證(表8),結(jié)果發(fā)現(xiàn)無論是耕地或者荒地土壤,模型預(yù)測平均值與實測平均值比較接近,RMSE 值較小,表明SahysMod 模型可以用于模擬及預(yù)測研究區(qū)內(nèi)不同灌排情景下的土壤鹽分、地下水埋深等動態(tài)變化規(guī)律。

    2.2 SahysMod 模型參數(shù)敏感性分析

    如圖4 所示,考慮到Flr 空間變異性,選取網(wǎng)格編號為26、60、76、87,研究Flr 對土壤電導(dǎo)率的敏感性。結(jié)果發(fā)現(xiàn)當(dāng)Flr 增加時(如Flr+20%及Flr+50%),根區(qū)土壤電導(dǎo)率降低;當(dāng)Flr 減小時(如Flr-20%及Flr-50%),根區(qū)土壤電導(dǎo)率增加,F(xiàn)lr 對根區(qū)土壤電導(dǎo)率比較敏感。另外,F(xiàn)lr值對地下水埋深無影響,這與以往的研究基本一致[18]。

    表7 過渡層淋洗效率、含水層淋洗效率及含水層水平導(dǎo)水率的率定 Table 7 Calibration of leaching efficiency of transition zone, aquifer zone and horizontal hydraulic conductivity

    如圖5 所示,不同的Flx、Flq 對地下水埋深具有較低的敏感性,Kaq 對地下水埋深具有較高的敏感性。當(dāng)Flx減小時(Flx-20%,即0.6;Flx-50%,即0.4),地下水埋深略增加;當(dāng)Flx 增加時(Flx+20%,即1.0),地下水埋深略減少。由于Flx 率定值為0.8,其范圍為0~1.0,故沒有考慮Flx 增加50%時地下水埋深的變化。當(dāng)Flq 增加時(Flq+20%,即1.44;Flq+50%,即1.80),地下水埋深略有減?。划?dāng)Flq 減小時(Flq-20%,即0.96;Flq-50%,即0.60),地下水埋深略有增加。這和以往研究略有差異,如Yao 等[18]認(rèn)為Flq 對地下水埋深沒有影響,這可能是由于本研究區(qū)面積(灌區(qū)尺度)相比較于Yao 的研究區(qū)域(農(nóng)田尺度)較大,研究區(qū)內(nèi)的地下水埋深具有一定的差異,并且由于耕地與荒地的土壤水鹽交換,使得Flq 對地下水埋深有微弱影響。最后,當(dāng)Kaq 增加時(Kaq+20%,即12 m/d;Kaq+50%,即15 m/d),地下水埋深增加,當(dāng)Kaq減小時(Kaq-20%,即8 m/d;Kaq-50%,即5 m/d),地下水埋深減小,這和以往研究相一致。Singh 等[17]指出在地下水埋深淺的鹽漬化區(qū)域,Kaq 的值比較小。Yao 等[18]認(rèn)為加大土壤導(dǎo)水率Kaq 有利于地下水埋深增加。

    表8 基于SahysMod 模型的土壤電導(dǎo)率預(yù)測值與實測值比較 Table 8 Comparison of measured and predicted soil electrical conductivity based on SahysMod model (dS·m-1)

    圖4 根區(qū)淋洗效率對土壤電導(dǎo)率敏感性分析 Fig.4 Sensitivity analysis of rootzone leaching efficiency on soil electrical conductivity

    圖5 過渡層淋洗效率、含水層淋洗效率及含水層水平導(dǎo)水率對地下水埋深的敏感性分析 Fig.5 Sensitivity analysis of leaching efficiency of transition zone, leaching efficiency of aquifer and horizontal hydraulic conductivity of aquifer on groundwater depth

    2.3 不同情景下土壤水鹽長期動態(tài)變化規(guī)律

    2.3.1 現(xiàn)有灌排管理模式

    基于SahysMod 預(yù)測在現(xiàn)有灌排管理情景下(即土地整治前)未來10 a 內(nèi)的耕地與荒地土壤電導(dǎo)率變化。如圖6 所示,隨機選取研究區(qū)內(nèi)的耕地與荒地網(wǎng)格(如荒地編號26、60,耕地編號76、87)。結(jié)果發(fā)現(xiàn),在預(yù)測初期(2017—2022 年)荒地土壤鹽分逐年升高,在預(yù)測后期(2023—2027 年)荒地土壤鹽分變化平緩且趨向于最大值,這可能是由于地下水埋深變淺(圖6),荒地區(qū)域在受到附近的灌溉耕地土壤水分水平滲透及地下水的垂直補給(荒地地下水電導(dǎo)率為5.6 dS/m),并且在強蒸發(fā)的作用下,水分?jǐn)y帶鹽分逐漸向土壤表層運移,水分被蒸發(fā)而鹽分最終積聚在荒地表層土壤之中,造成荒地土壤鹽分逐年升高。在預(yù)測后期荒地鹽分的積累可能逐步趨向于最大值,這是因為鹽分積累使得土壤結(jié)構(gòu)變化導(dǎo)致導(dǎo)水率降低,土面蒸發(fā)量減少,反過來影響了荒地土壤鹽分的進(jìn)一步積累。研究結(jié)果和以往的研究相一致,郭文聰?shù)萚42]認(rèn)為由于鈉離子的分散作用、降水的淋洗作用以及溶質(zhì)勢梯度作用下的鹽分?jǐn)U散和彌散作用對荒地鹽分積累的抑制作用,使得荒地鹽分積累存在最大值。

    耕地土壤鹽分在預(yù)測初期(2017—2022 年)變化緩慢,在預(yù)測后期(2023—2027 年)逐年增加。這是因為耕地的土壤鹽分初始值比較低(網(wǎng)格76 及87 初始土壤電導(dǎo)率分別為0.260 及0.243 dS/m),雖然地下水的垂直補給及灌溉帶來了大量鹽分(耕地地下水電導(dǎo)率為4.2 dS/m,灌溉水電導(dǎo)率為1.05 dS/m),但是地下水不斷得到灌溉的補充,灌溉的耕地與低洼的荒地形成地下水位差,受地下徑流的作用,高地勢區(qū)的耕地土壤鹽分最終通過灌溉淋洗作用不斷運移到荒地及周圍的田間排水溝,所以耕地土壤鹽分低且變化范圍不明顯。在預(yù)測后期,常年的灌溉導(dǎo)致地下水埋深變淺(由圖7 可知,生育期的地下水埋深由2017 年的1.65 m 到2027 年的1.22 m;非生育期的地下水埋深由2017 年的2.00 m 到2027 年的1.69 m),由于地下水的垂直補給及不斷的灌溉,給耕地帶來了大量鹽分。同時,在預(yù)測后期,灌溉耕地與荒地、排水溝之間的水力梯度變小,荒地在后期的積鹽量下降。盡管排水溝排鹽量仍在逐年增加(圖7),但是排水不僅要排走灌溉及地下水補給帶來的鹽分,而且也排走耕地土壤歷史上殘余鹽分,因此,在預(yù)測后期只依靠現(xiàn)有灌溉排水模式已不足以保證耕地土壤鹽分的排出效果。

    圖6 現(xiàn)有灌排管理下的未來10 a 土壤電導(dǎo)率變化 Fig.6 Change of existing irrigation and drainage on soil electrical conductivity in next ten years

    圖7 未來10 a 內(nèi)排水電導(dǎo)率及地下水埋深變化 Fig.7 Change of electrical conductivity of drained water and groundwater depth in next ten years

    2.3.2 不同灌水量

    如圖8 所示,隨機選取網(wǎng)格編號為76、87 的耕地,結(jié)果顯示在原有灌水量(670 mm)條件下,2024 年耕地根區(qū)土壤鹽分積累到障礙水平。與原有灌水量(670 mm)相比,當(dāng)減少灌水量時,耕地土壤鹽分逐年增加。灌水量從600 mm 減小到400 mm,土壤鹽分增加幅度提升,耕地根區(qū)土壤鹽分提前積累到障礙水平。當(dāng)增加灌水量,土壤鹽分每年的增加量減小,并且可以延遲耕地根區(qū)土壤鹽分積累到障礙水平的時間。灌水量從800 mm 增加到1 000 mm,耕地土壤鹽分在預(yù)測初期基本保持不變,在預(yù)測后期逐年緩慢上升,并且上升幅度逐年減小。

    2.3.3 不同灌溉水質(zhì)

    如圖9 所示,耕地根區(qū)土壤鹽分隨著灌溉水電導(dǎo)率的增大而增加,并且土壤鹽分提前積累到障礙水平。這可能是由于灌溉水帶入的鹽分大于排水溝排出的鹽分及荒地積累的鹽分,根區(qū)的鹽分迅速積累,造成耕地土壤鹽分上升。值得注意的是,在2025—2027 年耕地根區(qū)土壤鹽分雖然有所增加,但是增加幅度變小。這可能是由于土壤水鹽溶液已經(jīng)逐漸趨向于飽和,鹽分溶解量逐漸減小。

    圖8 灌溉水量對季節(jié)1 耕地土壤電導(dǎo)率的影響 Fig.8 Effect of irrigation amount on soil electrical conductivity of cultivated land in season 1

    圖9 灌溉水電導(dǎo)率對季節(jié)1 耕地土壤電導(dǎo)率的影響 Fig.9 Effect of electrical conductivity of irrigation water on soil electrical conductivity of cultivated land in season 1

    在不同灌溉水水質(zhì)下,以灌溉水的電導(dǎo)率為2.10 dS/m的耕地土壤鹽分增加量最大,在預(yù)測第10 年,網(wǎng)格編號為76 及87 的耕地土壤電導(dǎo)率分別達(dá)5.97 及6.13 dS/m,根據(jù)銀北灌區(qū)土地鹽漬化程度分級標(biāo)準(zhǔn)[43-44],此時的土壤已為鹽土,耕地已經(jīng)不適合作物種植。當(dāng)灌溉水電導(dǎo)率為0.60 dS/m,耕地土壤鹽分在未來10 a 內(nèi)基本保持不變,土壤電導(dǎo)率最大值約0.5 dS/m,遠(yuǎn)低于1.7 dS/m,基于FAO對玉米產(chǎn)生鹽分脅迫導(dǎo)致減產(chǎn)的程度[31],因此認(rèn)為在此情景下可以滿足作物的正常生長需求。

    2.3.4 不同排水溝深度

    如圖10 所示,當(dāng)排水溝深為1.5 m 時,耕地根區(qū)土壤鹽分增加最快,并且在2024 年耕地土壤鹽分積累到障礙水平。在預(yù)測后期,根區(qū)土壤鹽分增加幅度逐年降低,這可能是由于土壤水鹽溶液已經(jīng)逐漸趨向于飽和,鹽分溶解量逐漸減小,進(jìn)而導(dǎo)致耕地土壤鹽分的增加幅度逐年降低。通過土地整治,加深排水溝深度,可以有效延遲耕地根區(qū)土壤鹽分積累到障礙水平的時間。當(dāng)排水溝深度為2.2 m時,在未來10 a 內(nèi)耕地土壤鹽分基本保持不變,并小于1.7 dS/m,在此條件下的耕地可以保證區(qū)內(nèi)玉米正常生長;當(dāng)排水溝深為2.5 m 時,在未來10 a 內(nèi)耕地土壤鹽分甚至出現(xiàn)緩慢下降的趨勢,耕地土壤電導(dǎo)率值均小于0.5 dS/m。但是,綜合考慮建造排水溝工程量的大小及作物正常生長的需要,在本研究中推薦的排水溝深度為2.2 m。

    圖10 排水溝深度對季節(jié)1 耕地土壤電導(dǎo)率的影響 Fig.10 Effect of drainage ditch depth on soil electrical conductivity of cultivated land in season 1

    3 討 論

    經(jīng)過率定和驗證,SahysMod 模型可以綜合考慮土壤理化參數(shù)的空間變異性,能夠反映土壤實際狀況,可以用來對灌區(qū)的土壤鹽分、灌水量、地下水埋深等模擬預(yù)測。這和以往的研究一致,Sarangi 等[45]通過比較BP 神經(jīng)網(wǎng)絡(luò)、RBF 神經(jīng)網(wǎng)絡(luò)、GRNN 神經(jīng)網(wǎng)絡(luò)及SaltMod 模型,認(rèn)為SaltMod 模型在根區(qū)土壤鹽分的預(yù)測精度要高于人工神經(jīng)網(wǎng)絡(luò);Yao 等[18]已經(jīng)驗證SahysMod 模型能夠預(yù)測中國濱海農(nóng)田的土壤鹽分時空變化;Akram 等[46]基于SahysMod 模型設(shè)置不同情景(如不同的導(dǎo)水率、初始鹽分含量、灌溉水鹽分、初始地下水埋深等)模擬生物排鹽在解決鹽漬化問題的可行性,認(rèn)為SahysMod 模型可以用來模擬及預(yù)測干旱半干旱地區(qū)的土壤鹽漬化問題。但是,SahysMod 模型在率定Flr、Flx、Flq 及Kaq 參數(shù)時,實測值與模擬值仍存在一定的誤差,其可能原因包括收集到的灌區(qū)地下水動態(tài)觀測資料、含水層性質(zhì)和氣象資料有限;土壤鹽分存在較大的空間變異,SahysMod 模型沒有考慮植物體吸鹽量以及施用化肥引入的鹽量;最后,農(nóng)戶對于種植作物的選擇是基于市場經(jīng)濟效益,當(dāng)一種作物的經(jīng)濟價值降低時,農(nóng)戶可能會改變作物種植制度,進(jìn)而影響區(qū)域鹽分的空間變化。

    基于SahysMod 模型對研究區(qū)劃分了98 個內(nèi)部網(wǎng)格,但限于篇幅限制,只隨機選取了代表性的荒地多邊形網(wǎng)格1、26、60 及耕地網(wǎng)格7、61、76、87、89。所選的荒地及耕地網(wǎng)格基本反映了研究區(qū)典型的耕地與荒地的土壤水鹽動態(tài)變化。荒地主要分布在研究區(qū)的西南部,如圖2 所示,選取的網(wǎng)格1 位于整個荒地區(qū)域的下部,網(wǎng)格26 位于整個荒地區(qū)域的中部,網(wǎng)格60 位于整個荒地區(qū)域的上部,基本上可以代表整個荒地區(qū)域。同時,荒地土壤電導(dǎo)率范圍為1.07~5.50 dS/m,其中,網(wǎng)格1、26、60 的電導(dǎo)率分別為4.44、1.75、3.77 dS/m,基本上可以反映荒地土壤電導(dǎo)率范圍;耕地主要分布在研究區(qū)的北部,如圖2 所示,選取的網(wǎng)格7 位于整個耕地區(qū)域的下部,網(wǎng)格61、76 位于整個耕地區(qū)域的中部,網(wǎng)格87、89位于整個耕地區(qū)域的上部,基本上可以代表整個耕地區(qū)域。耕地土壤電導(dǎo)率范圍為0.14~1.07 dS/m,其中,網(wǎng)格7、61、76、87、89 的電導(dǎo)率分別為0.67、0.88、0.34、0.44、0.31 dS/m,基本上可以反映耕地土壤電導(dǎo)率范圍。另外,從荒地與耕地的位置關(guān)系考慮,耕地網(wǎng)格7 被荒地網(wǎng)格3、4、6 包圍,耕地網(wǎng)格61 緊鄰荒地網(wǎng)格70,耕地網(wǎng)格76 鄰近荒地網(wǎng)格68,耕地網(wǎng)格87、89 遠(yuǎn)離荒地(圖2)。因此,可以反映荒地與耕地的不同位置關(guān)系對灌排管理的影響。

    在當(dāng)前的灌排管理中,灌區(qū)灌水量不僅應(yīng)保證作物生長及排走灌溉本身帶來的鹽分,而且應(yīng)排走歷史上土壤的殘余鹽分。SahysMod 模型預(yù)測未來10 a 土壤鹽分變化時,發(fā)現(xiàn)隨灌水量的增加,“驅(qū)鹽”效果越好。這和以往的研究基本一致,王旭等[47]指出隨灌水量增加脫鹽效果越顯著,史海濱等[48]也認(rèn)為灌溉可以對土壤鹽分進(jìn)行有針對性的調(diào)控,進(jìn)而達(dá)到“驅(qū)鹽”效果。因此,加大灌水量是解決土壤鹽漬化的一個重要途徑。但是灌水量并不是越多越好,灌水量過大雖然可以增加鹽分的淋洗量,但同時也會帶來較強的土壤蒸發(fā),降低地下水埋深,帶動下層鹽分上移,使得作物根系活動層出現(xiàn)積鹽[49]。因此,在土地整治時灌區(qū)應(yīng)盡量選取適宜的灌水量。

    銀北灌區(qū)主要通過黃河水進(jìn)行灌溉,但是近幾年來為了兼顧全流域的農(nóng)業(yè)、生態(tài)、生產(chǎn)問題以及保證下游的社會發(fā)展,在銀北灌區(qū)的引黃量逐漸減少,采用微咸水灌溉是水資源約束區(qū)域的重要灌水方式之一。如劉娟等[50]在寧夏銀北地區(qū)采用微咸水灌溉白漿土鹽堿地,認(rèn)為采用礦化度為1 g/L 的微咸水滴灌,可以獲得較好的植物生長和較高產(chǎn)量;王詩景等[51]在寧夏銀北惠農(nóng)引黃灌區(qū)開展微咸水灌溉試驗,指出采用井渠1∶1 的混灌模式是春小麥的適宜微咸水灌溉利用模式。在水資源約束下,銀北灌區(qū)不少農(nóng)戶通過排水溝排出水直接進(jìn)行灌溉,但是研究區(qū)內(nèi)排水溝排出水含鹽量高(2016—2017 年期間對10 個排水溝中排出水的電導(dǎo)率進(jìn)行測定,其電導(dǎo)率均值為2.10 dS/m)。如圖10 所示,采用排水溝排出的水灌溉耕地,耕地土壤電導(dǎo)率逐年增加,到預(yù)測后期耕地已不再適用于種植作物。直接采用排水溝水進(jìn)行灌溉不僅增加土壤鹽分、破壞土壤結(jié)構(gòu),還將影響作物產(chǎn)量[52]。

    陳艷梅等[53]認(rèn)為排水溝深度達(dá)3.0 m 時,可有效減少以河套灌區(qū)沙壕渠灌域耕地土壤鹽分;Bah?eci 等[52]指出排水溝深度在1.2 m 時,可以保證土耳其Konya–C?umra平原的作物生長,并認(rèn)為在原有1.5 m 的排水溝深度上應(yīng)減小排水溝深。本研究認(rèn)為通過土地整治的排水溝深度達(dá)到2.2 m 以上時,在未來10 a 內(nèi)耕地土壤鹽分增加趨勢相對較小,并小于1.7 dS/m,可以延遲耕地根區(qū)土壤鹽分積累到障礙水平的時間,在此條件下的耕地可以保證區(qū)內(nèi)玉米正常生長。但是,這比Bah?eci 等[52]設(shè)定的排水溝深度大的多,原因可能是研究區(qū)灌溉水的電導(dǎo)率(1.05 dS/m)高于 Bah?eci 研究中設(shè)置的電導(dǎo)率(0.4 dS/m),同時耕地地下水礦化度較高(耕地地下水電導(dǎo)率為4.2 dS/m),灌溉及地下水垂直補給給耕地土壤帶來的鹽分較多。這又比陳艷梅等[53](2012)的排水溝深度小,可能是由于研究區(qū)耕地土壤初始鹽分相較于陳艷梅的研究區(qū)中的初始鹽分小,進(jìn)而影響了排水溝深度的設(shè)置。值得注意的是,通過土地整治加大排水溝深度,雖然可以增加排水排鹽的能力,可有效控制耕地土壤鹽分的增加,但是排水溝深度的增加會加大土地整治的工程量,同時經(jīng)濟投入量也相應(yīng)增加。另外,也有研究指出過分加深排水溝將導(dǎo)致水面蒸發(fā)量減少,間接引起耕地向荒地的干排鹽量減少,即削弱了“干排鹽”的作用,將導(dǎo)致耕地積鹽量上升[54]。因此,在灌區(qū)土地整治時應(yīng)合理地設(shè)計排水溝深度,以便同時實現(xiàn)土壤鹽分降低和經(jīng)濟效益最大化。

    4 結(jié) 論

    基于SahysMod 模型研究在土地整治過程中不同灌排管理下未來10 a 內(nèi)土壤水鹽動態(tài)變化。通過率定及驗證,認(rèn)為SahysMod 模型綜合考慮土壤理化參數(shù)的空間變異性,可以對土地整治過程中不同灌排管理情景下土壤水鹽動態(tài)變化進(jìn)行模擬及預(yù)測。不同灌水量、灌溉水水質(zhì)及排水溝深度中,耕地及荒地土壤鹽漬化程度不同。其主要結(jié)論如下:

    1)根區(qū)淋洗效率對土壤鹽分比較敏感,對地下水埋深無影響。過渡層淋洗效率及含水層淋洗效率對地下水埋深具有低的敏感性,含水層導(dǎo)水率對地下水埋深具有較高的敏感性。

    2)現(xiàn)有灌排管理模式下(即灌水量為670 mm,灌溉水電導(dǎo)率為1.05 dS/m,排水溝深1.5 m),荒地土壤鹽分在2017—2022 年逐年升高,在2023—2027 年變化平緩;耕地土壤鹽分在2017—2022 年變化緩慢,在2023—2027 年逐年緩慢增加。

    3)加大灌水量是解決土壤鹽漬化的一個重要途徑,可以延遲耕地根區(qū)土壤鹽分積累到障礙水平的時間,但是需要同時兼顧銀北灌區(qū)的水資源約束問題。

    4)在現(xiàn)有灌排管理模式下,2024 年以后作物生長就會受到鹽害脅迫。當(dāng)灌溉水電導(dǎo)率減小為0.60 dS/m 時,耕地土壤鹽分在未來10 a 內(nèi)基本保持不變,土壤電導(dǎo)率最大值約0.5 dS/m,可以滿足作物的正常生長需求;當(dāng)灌溉水電導(dǎo)率增加時,作物生長受到鹽害脅迫的時間相應(yīng)提前。由于黃河來水量的減少,部分農(nóng)民利用排水溝水進(jìn)行灌溉,這會加快鹽漬化脅迫的出現(xiàn)時間

    5)排水溝深度影響耕地土壤鹽分水平,通過土地整治加深排水溝深度,可以有效延遲土壤電導(dǎo)率達(dá)到障礙水平的時間。當(dāng)排水溝深為1.5 m 時,耕地根區(qū)土壤鹽分增加最快,并且在2024 年耕地土壤鹽分積累到障礙水平;當(dāng)排水溝深為2.2 m 時,在未來10 a 內(nèi)耕地土壤鹽分基本保持不變,并且土壤電導(dǎo)率值均小于1.7 dS/m,可以保證研究區(qū)內(nèi)玉米正常生長。

    猜你喜歡
    根區(qū)荒地排水溝
    獨登南山
    熱風(fēng)管道加溫下日光溫室根區(qū)溫度場的CFD模擬
    桉樹人工幼齡林根區(qū)和非根區(qū)土壤屬性特征分析
    Thalidomide for refractory gastrointestinal bleeding from vascular malformations in patients with significant comorbidities
    荒 地
    中國詩歌(2018年6期)2018-11-14 13:24:12
    LED補光和根區(qū)加溫對日光溫室起壟內(nèi)嵌式基質(zhì)栽培甜椒生長及產(chǎn)量的影響*
    對外發(fā)包的荒地為何被判無效
    動詞“Get”的用法
    做好種植屋面排水,防止植物爛根
    樹盤施肥區(qū)域大小對 15N吸收利用及桃幼樹生長的影響
    欧美日韩亚洲高清精品| 亚洲av日韩精品久久久久久密| av免费在线观看网站| 久久午夜综合久久蜜桃| 69av精品久久久久久 | 一区二区三区乱码不卡18| 欧美少妇被猛烈插入视频| 亚洲一区二区三区欧美精品| 免费在线观看完整版高清| 伊人久久大香线蕉亚洲五| av国产精品久久久久影院| 亚洲av男天堂| av网站免费在线观看视频| 精品人妻在线不人妻| 中文字幕制服av| 高清欧美精品videossex| 亚洲中文日韩欧美视频| 亚洲精品国产色婷婷电影| 欧美日韩视频精品一区| 日韩熟女老妇一区二区性免费视频| 亚洲少妇的诱惑av| a 毛片基地| 国精品久久久久久国模美| 各种免费的搞黄视频| 亚洲国产欧美一区二区综合| 亚洲视频免费观看视频| 国产成人欧美| 手机成人av网站| 精品高清国产在线一区| 欧美中文综合在线视频| 日韩一卡2卡3卡4卡2021年| 色精品久久人妻99蜜桃| 岛国毛片在线播放| 欧美 日韩 精品 国产| 精品卡一卡二卡四卡免费| 国产高清国产精品国产三级| 精品国产一区二区久久| 欧美精品人与动牲交sv欧美| 午夜福利在线免费观看网站| 国产免费av片在线观看野外av| 狠狠精品人妻久久久久久综合| 亚洲自偷自拍图片 自拍| 久久影院123| 国产精品秋霞免费鲁丝片| 日韩有码中文字幕| 久久久久国内视频| 这个男人来自地球电影免费观看| 久久国产精品大桥未久av| 午夜福利一区二区在线看| 涩涩av久久男人的天堂| 国产视频一区二区在线看| 老司机在亚洲福利影院| 亚洲专区字幕在线| 久久久久久久精品精品| 免费久久久久久久精品成人欧美视频| 汤姆久久久久久久影院中文字幕| 精品久久久久久久毛片微露脸 | 国产精品二区激情视频| 亚洲少妇的诱惑av| 热re99久久精品国产66热6| 国产亚洲精品一区二区www | av欧美777| 在线亚洲精品国产二区图片欧美| 91九色精品人成在线观看| 50天的宝宝边吃奶边哭怎么回事| 国产精品秋霞免费鲁丝片| 这个男人来自地球电影免费观看| 午夜福利在线免费观看网站| a级片在线免费高清观看视频| 午夜视频精品福利| 曰老女人黄片| 亚洲中文日韩欧美视频| 久久热在线av| 欧美日韩黄片免| 亚洲精品国产一区二区精华液| 久久精品亚洲av国产电影网| 亚洲精品中文字幕一二三四区 | 美女中出高潮动态图| 另类精品久久| 超碰97精品在线观看| 中文字幕精品免费在线观看视频| 精品第一国产精品| 99热国产这里只有精品6| 精品福利永久在线观看| 亚洲av电影在线进入| 在线永久观看黄色视频| 丝袜喷水一区| 老汉色∧v一级毛片| 美女国产高潮福利片在线看| 日本av手机在线免费观看| 国产亚洲午夜精品一区二区久久| 热99re8久久精品国产| 啪啪无遮挡十八禁网站| 午夜福利免费观看在线| 国产不卡av网站在线观看| 亚洲国产毛片av蜜桃av| 国产欧美日韩一区二区三区在线| 亚洲中文字幕日韩| 一本久久精品| 亚洲色图综合在线观看| 高清av免费在线| 国产一区二区在线观看av| 日韩中文字幕视频在线看片| 亚洲国产av影院在线观看| 免费人妻精品一区二区三区视频| 丝袜美腿诱惑在线| 亚洲精品美女久久av网站| 久久狼人影院| 久久精品熟女亚洲av麻豆精品| av免费在线观看网站| av不卡在线播放| 老司机影院成人| 免费一级毛片在线播放高清视频 | 一区二区日韩欧美中文字幕| 好男人电影高清在线观看| 麻豆乱淫一区二区| 99国产精品一区二区三区| 窝窝影院91人妻| 天堂8中文在线网| 免费久久久久久久精品成人欧美视频| 免费不卡黄色视频| 国产欧美日韩综合在线一区二区| 99热国产这里只有精品6| tocl精华| 久9热在线精品视频| 欧美精品一区二区大全| 一个人免费在线观看的高清视频 | 国产一区二区激情短视频 | 69精品国产乱码久久久| 国产黄色免费在线视频| 永久免费av网站大全| 久久久久精品人妻al黑| 久久精品国产亚洲av香蕉五月 | 少妇人妻久久综合中文| 99国产综合亚洲精品| 亚洲av日韩精品久久久久久密| 久久久久精品人妻al黑| 岛国毛片在线播放| 在线天堂中文资源库| 亚洲男人天堂网一区| 美女午夜性视频免费| 无遮挡黄片免费观看| 天天躁狠狠躁夜夜躁狠狠躁| 日韩制服丝袜自拍偷拍| 999精品在线视频| 亚洲三区欧美一区| 黄网站色视频无遮挡免费观看| av免费在线观看网站| 我要看黄色一级片免费的| 热re99久久国产66热| 丁香六月欧美| 久久 成人 亚洲| 日日夜夜操网爽| 亚洲 欧美一区二区三区| 淫妇啪啪啪对白视频 | 欧美国产精品一级二级三级| 老司机午夜福利在线观看视频 | 国产精品.久久久| 久久香蕉激情| 午夜福利,免费看| 他把我摸到了高潮在线观看 | 美女中出高潮动态图| 亚洲国产看品久久| 国产淫语在线视频| 国产主播在线观看一区二区| 国产淫语在线视频| 日本精品一区二区三区蜜桃| 热99久久久久精品小说推荐| 成年动漫av网址| 亚洲av欧美aⅴ国产| 久久久久久亚洲精品国产蜜桃av| 国产成人av激情在线播放| 一本久久精品| 久久久久国产精品人妻一区二区| 国产一级毛片在线| 日本91视频免费播放| 国产在线视频一区二区| 五月开心婷婷网| 大型av网站在线播放| 欧美日韩成人在线一区二区| 精品国产一区二区久久| 亚洲国产欧美在线一区| 成人影院久久| netflix在线观看网站| 亚洲va日本ⅴa欧美va伊人久久 | 少妇精品久久久久久久| 男人添女人高潮全过程视频| 91成人精品电影| 欧美黄色片欧美黄色片| 在线天堂中文资源库| videosex国产| 丝袜美腿诱惑在线| 国产精品 国内视频| 国产一区二区三区av在线| 国产一区二区三区av在线| 操出白浆在线播放| 国产男人的电影天堂91| 久久ye,这里只有精品| 欧美 亚洲 国产 日韩一| 亚洲精品成人av观看孕妇| 又紧又爽又黄一区二区| 日韩中文字幕欧美一区二区| 亚洲三区欧美一区| 亚洲精品国产av蜜桃| 国产亚洲av片在线观看秒播厂| 午夜福利,免费看| 51午夜福利影视在线观看| 男女边摸边吃奶| 麻豆国产av国片精品| 国产男人的电影天堂91| 亚洲欧洲精品一区二区精品久久久| 高清在线国产一区| 日日夜夜操网爽| 亚洲欧美清纯卡通| av免费在线观看网站| 亚洲伊人色综图| 宅男免费午夜| 亚洲av电影在线进入| 51午夜福利影视在线观看| 亚洲精品日韩在线中文字幕| 日韩有码中文字幕| 各种免费的搞黄视频| 久久久精品94久久精品| 人人妻人人添人人爽欧美一区卜| 精品少妇久久久久久888优播| 亚洲欧美清纯卡通| 秋霞在线观看毛片| 99久久人妻综合| 国产又色又爽无遮挡免| 欧美激情极品国产一区二区三区| 精品福利永久在线观看| 国产精品成人在线| 婷婷丁香在线五月| 久久国产精品人妻蜜桃| 欧美日韩国产mv在线观看视频| 一级a爱视频在线免费观看| 老汉色∧v一级毛片| 十八禁高潮呻吟视频| 精品视频人人做人人爽| 午夜福利一区二区在线看| 国产精品国产av在线观看| 国产真人三级小视频在线观看| 国产一卡二卡三卡精品| 国产有黄有色有爽视频| 最近最新中文字幕大全免费视频| 成人手机av| 亚洲精品国产一区二区精华液| 精品卡一卡二卡四卡免费| 狂野欧美激情性bbbbbb| 黄网站色视频无遮挡免费观看| 91字幕亚洲| 亚洲伊人色综图| 极品少妇高潮喷水抽搐| 在线观看一区二区三区激情| 午夜免费观看性视频| 日韩,欧美,国产一区二区三区| 欧美xxⅹ黑人| 免费高清在线观看视频在线观看| 美女主播在线视频| 亚洲免费av在线视频| 国产人伦9x9x在线观看| 久久精品亚洲熟妇少妇任你| 另类精品久久| 亚洲国产欧美日韩在线播放| 欧美日韩亚洲高清精品| 精品乱码久久久久久99久播| 国产男女超爽视频在线观看| 亚洲综合色网址| 国产亚洲一区二区精品| 精品人妻在线不人妻| 最近最新免费中文字幕在线| av一本久久久久| 午夜激情av网站| 捣出白浆h1v1| 国产成人欧美在线观看 | 免费高清在线观看视频在线观看| 国产成人av激情在线播放| 一级毛片女人18水好多| 久久香蕉激情| 99热国产这里只有精品6| 最黄视频免费看| 夜夜夜夜夜久久久久| 亚洲欧美清纯卡通| 亚洲,欧美精品.| 亚洲国产av新网站| www日本在线高清视频| 亚洲av美国av| 亚洲一区二区三区欧美精品| 色综合欧美亚洲国产小说| 精品国产国语对白av| 自线自在国产av| 成年人免费黄色播放视频| 91国产中文字幕| 久久久久久久久免费视频了| 啦啦啦啦在线视频资源| 亚洲成av片中文字幕在线观看| 男女边摸边吃奶| 中文精品一卡2卡3卡4更新| 国产高清国产精品国产三级| 91精品国产国语对白视频| www日本在线高清视频| 国产高清视频在线播放一区 | 精品一区二区三区四区五区乱码| 最近中文字幕2019免费版| 精品福利永久在线观看| 亚洲精品久久久久久婷婷小说| 日韩有码中文字幕| www.精华液| 欧美 日韩 精品 国产| 亚洲人成77777在线视频| 大片电影免费在线观看免费| 老汉色∧v一级毛片| 老司机影院毛片| 男女国产视频网站| 丝袜在线中文字幕| 欧美日韩亚洲综合一区二区三区_| 亚洲av电影在线观看一区二区三区| 日韩一区二区三区影片| 亚洲成av片中文字幕在线观看| 一区二区av电影网| 国产亚洲精品一区二区www | 亚洲色图 男人天堂 中文字幕| 久久久久视频综合| 99久久国产精品久久久| 精品亚洲成国产av| tube8黄色片| 午夜福利乱码中文字幕| 久久精品人人爽人人爽视色| 纵有疾风起免费观看全集完整版| e午夜精品久久久久久久| 亚洲精品国产色婷婷电影| av不卡在线播放| 午夜福利乱码中文字幕| 又大又爽又粗| 国产高清国产精品国产三级| 亚洲熟女毛片儿| 在线永久观看黄色视频| 欧美精品人与动牲交sv欧美| 一本—道久久a久久精品蜜桃钙片| 日日爽夜夜爽网站| 999精品在线视频| 最近中文字幕2019免费版| 老汉色∧v一级毛片| 日韩欧美一区二区三区在线观看 | 18禁黄网站禁片午夜丰满| 免费在线观看视频国产中文字幕亚洲 | 制服人妻中文乱码| 国产精品国产三级国产专区5o| 亚洲av日韩精品久久久久久密| av片东京热男人的天堂| 一级a爱视频在线免费观看| 国产在线免费精品| av网站在线播放免费| 亚洲成人免费电影在线观看| 中文字幕人妻丝袜一区二区| 中文字幕人妻熟女乱码| 少妇猛男粗大的猛烈进出视频| 在线观看一区二区三区激情| 亚洲熟女精品中文字幕| 免费在线观看日本一区| 性色av一级| 老司机影院毛片| 如日韩欧美国产精品一区二区三区| 亚洲av电影在线观看一区二区三区| 免费日韩欧美在线观看| av天堂在线播放| 大香蕉久久网| 免费观看a级毛片全部| 精品少妇久久久久久888优播| 韩国高清视频一区二区三区| 久久久欧美国产精品| 欧美精品一区二区大全| 成人18禁高潮啪啪吃奶动态图| 狠狠精品人妻久久久久久综合| 叶爱在线成人免费视频播放| 欧美一级毛片孕妇| 久久人人爽人人片av| 国产欧美日韩一区二区三 | 成年女人毛片免费观看观看9 | 99国产综合亚洲精品| 成人影院久久| 天天躁夜夜躁狠狠躁躁| 天天操日日干夜夜撸| 中国国产av一级| 亚洲全国av大片| 少妇裸体淫交视频免费看高清 | 狂野欧美激情性bbbbbb| 国产欧美日韩一区二区三 | 国产男女内射视频| a 毛片基地| 国产一卡二卡三卡精品| 国产熟女午夜一区二区三区| 免费久久久久久久精品成人欧美视频| 一区二区三区乱码不卡18| 波多野结衣av一区二区av| 人人妻人人澡人人爽人人夜夜| 亚洲第一欧美日韩一区二区三区 | 久久精品成人免费网站| 中文精品一卡2卡3卡4更新| 亚洲av片天天在线观看| 丝袜喷水一区| 久久免费观看电影| 美女扒开内裤让男人捅视频| 国产精品一二三区在线看| 大香蕉久久网| 两人在一起打扑克的视频| a 毛片基地| 久久中文看片网| 午夜精品久久久久久毛片777| 日日摸夜夜添夜夜添小说| 乱人伦中国视频| 青草久久国产| 欧美一级毛片孕妇| 精品久久久精品久久久| 黑人巨大精品欧美一区二区蜜桃| 人妻一区二区av| 后天国语完整版免费观看| 正在播放国产对白刺激| www.自偷自拍.com| 欧美日本中文国产一区发布| 亚洲天堂av无毛| 青春草亚洲视频在线观看| 国产精品久久久久久精品电影小说| 日本91视频免费播放| 精品久久久久久久毛片微露脸 | www.精华液| 亚洲欧美激情在线| 考比视频在线观看| 久久久精品国产亚洲av高清涩受| 天天躁日日躁夜夜躁夜夜| 欧美变态另类bdsm刘玥| 久久九九热精品免费| 亚洲av成人一区二区三| 午夜福利视频在线观看免费| 欧美在线一区亚洲| 成人18禁高潮啪啪吃奶动态图| 国产欧美日韩一区二区三区在线| 日韩熟女老妇一区二区性免费视频| 精品国产一区二区久久| 老司机影院成人| 欧美老熟妇乱子伦牲交| 久久久久国内视频| 国产成人精品久久二区二区免费| videos熟女内射| 一级,二级,三级黄色视频| 日日摸夜夜添夜夜添小说| 脱女人内裤的视频| 岛国在线观看网站| 久久久欧美国产精品| 99精品久久久久人妻精品| 五月天丁香电影| 高清视频免费观看一区二区| 国产精品亚洲av一区麻豆| 黄色视频在线播放观看不卡| 欧美激情 高清一区二区三区| 日韩制服骚丝袜av| 在线 av 中文字幕| 久久中文看片网| 欧美黄色淫秽网站| 啦啦啦啦在线视频资源| 日韩视频一区二区在线观看| 爱豆传媒免费全集在线观看| 久久久精品区二区三区| 性高湖久久久久久久久免费观看| 高清av免费在线| 飞空精品影院首页| 黄频高清免费视频| 国产日韩一区二区三区精品不卡| 丝瓜视频免费看黄片| 欧美精品人与动牲交sv欧美| 搡老熟女国产l中国老女人| 韩国精品一区二区三区| 高清黄色对白视频在线免费看| 亚洲欧洲精品一区二区精品久久久| 18禁观看日本| 一级片免费观看大全| 精品少妇内射三级| 一区二区日韩欧美中文字幕| 视频区图区小说| 五月天丁香电影| 女性生殖器流出的白浆| 久久久国产一区二区| 欧美日韩亚洲高清精品| av超薄肉色丝袜交足视频| 王馨瑶露胸无遮挡在线观看| 国产不卡av网站在线观看| 国产成人欧美在线观看 | 国产一区二区三区在线臀色熟女 | 国产精品欧美亚洲77777| 韩国精品一区二区三区| 成年人黄色毛片网站| 一区二区三区四区激情视频| 国产高清视频在线播放一区 | 亚洲情色 制服丝袜| 伊人久久大香线蕉亚洲五| 欧美激情高清一区二区三区| 国产成人啪精品午夜网站| 久久人人爽人人片av| 国产一区二区三区综合在线观看| 一区二区av电影网| 曰老女人黄片| 在线天堂中文资源库| 色婷婷av一区二区三区视频| 国产精品99久久99久久久不卡| 久久亚洲国产成人精品v| 欧美日韩福利视频一区二区| 超碰成人久久| 亚洲精品一二三| 欧美另类一区| 男女午夜视频在线观看| 欧美精品一区二区免费开放| www日本在线高清视频| 亚洲精品国产精品久久久不卡| 男女国产视频网站| 97在线人人人人妻| 久久精品亚洲av国产电影网| 欧美人与性动交α欧美软件| 国产一区二区三区综合在线观看| 91麻豆av在线| 999久久久精品免费观看国产| 高清黄色对白视频在线免费看| 91精品伊人久久大香线蕉| 高清视频免费观看一区二区| 国产男女内射视频| 国产日韩欧美亚洲二区| 亚洲色图 男人天堂 中文字幕| 色婷婷久久久亚洲欧美| 性高湖久久久久久久久免费观看| 久久热在线av| 美女午夜性视频免费| 精品国产乱码久久久久久小说| 秋霞在线观看毛片| 国产av精品麻豆| 亚洲精品久久午夜乱码| 他把我摸到了高潮在线观看 | 亚洲精品日韩在线中文字幕| 国产成+人综合+亚洲专区| 深夜精品福利| 一级毛片电影观看| 亚洲国产毛片av蜜桃av| 国产成人一区二区三区免费视频网站| 亚洲成人免费电影在线观看| 亚洲国产精品一区二区三区在线| 99久久综合免费| 91精品三级在线观看| 正在播放国产对白刺激| 久久国产精品男人的天堂亚洲| 99久久国产精品久久久| 亚洲精品中文字幕在线视频| 最近最新免费中文字幕在线| 国产在视频线精品| 亚洲中文av在线| 人人澡人人妻人| 考比视频在线观看| 老司机福利观看| 啦啦啦免费观看视频1| 两人在一起打扑克的视频| 欧美亚洲日本最大视频资源| 一区二区三区精品91| a级毛片黄视频| 91av网站免费观看| 亚洲av成人一区二区三| 女性生殖器流出的白浆| 精品国内亚洲2022精品成人 | 电影成人av| 高清欧美精品videossex| 老司机影院毛片| 国产成人a∨麻豆精品| 日韩,欧美,国产一区二区三区| 国产精品 欧美亚洲| 午夜福利视频在线观看免费| 菩萨蛮人人尽说江南好唐韦庄| 午夜视频精品福利| tocl精华| 美女视频免费永久观看网站| 中文欧美无线码| 亚洲精品第二区| 亚洲精品日韩在线中文字幕| 一边摸一边抽搐一进一出视频| 中文字幕人妻丝袜制服| 啦啦啦视频在线资源免费观看| 一本一本久久a久久精品综合妖精| 夫妻午夜视频| 少妇被粗大的猛进出69影院| 午夜久久久在线观看| 久久狼人影院| www.精华液| 久久亚洲精品不卡| 亚洲精品国产一区二区精华液| 亚洲第一av免费看| 爱豆传媒免费全集在线观看| 亚洲欧美色中文字幕在线| 十分钟在线观看高清视频www| 精品人妻1区二区| 日韩制服骚丝袜av| 成人国产一区最新在线观看| 别揉我奶头~嗯~啊~动态视频 | 国产成人av教育| 老汉色∧v一级毛片| 国产黄频视频在线观看| 国产精品一二三区在线看| 最新的欧美精品一区二区| 啦啦啦啦在线视频资源| 国产成人影院久久av| 欧美日韩精品网址| 一级片免费观看大全| 色老头精品视频在线观看| 精品国产乱码久久久久久男人| 99久久综合免费| 人人妻人人澡人人爽人人夜夜| 91精品国产国语对白视频| 精品人妻熟女毛片av久久网站| 中文字幕最新亚洲高清| 亚洲精品久久成人aⅴ小说| 日本猛色少妇xxxxx猛交久久|