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

    生物炭對東北草甸黑土水力特性影響的數值化研究

    2018-05-31 03:34:57王睿垠魏永霞張翼鵬張雨鳳馬德才
    農業(yè)機械學報 2018年5期
    關鍵詞:導水率吸力毛細管

    王睿垠 魏永霞 張翼鵬 張雨鳳 馬德才

    (1.東北農業(yè)大學理學院, 哈爾濱 150030; 2.東北農業(yè)大學水利與土木工程學院, 哈爾濱 150030;3.農業(yè)部農業(yè)水資源高效利用重點實驗室, 哈爾濱 150030; 4.中山大學中法核工程與技術學院, 珠海 519082)

    0 引言

    東北黑土區(qū)是世界三大黑土帶之一,也是國家商品糧生產基地,糧食產量對于保障國家安全具有重要意義。水分是農業(yè)生產的重要因素,生物炭作為土壤改良劑對土壤水分運動的影響尚不明確。目前針對東北黑土施加生物炭的研究多為試驗研究,少有生物炭對土壤水分影響的數值研究。土壤是一個十分復雜的系統(tǒng),是植物生長中水分和養(yǎng)分的主要來源,是由固、液、氣三相物質組成的體系。生物炭是生物有機質在無氧或缺氧情況下,經高溫裂解的產物[1]。研究表明,在某些情況下,土壤中添加生物炭,可以修復土壤[2]并改變土壤的持水、導水特性[3-4]。土壤的導水率是土壤粒徑、黏土含量、有機質含量以及粒子大小分布等的函數,是全面反映土壤導水特性的一個參量[5]。土壤具有較高的導水率時可以加快滲透和排水[6],土壤快速排水有利于減少徑流的發(fā)生[7],但滲透速度過快,也會減少水中營養(yǎng)物和農用化學品充分溶解過濾的機會[8]。黑土是中國重要的土壤資源,由于黑土表土疏松、底土黏重以及人為不合理的耕作,使得東北黑土持水能力降低,農田生產力下降。生物炭作為土壤改良劑,在改善土壤的持水和導水能力、提高作物的生產力方面與其他有機添加物相比,更加安全和持久[9]。研究證實了土壤添加生物炭可以影響作物的產量[10-12],其物理機制多歸因于土壤孔隙度的增加、土壤容重的降低[13]、土壤保水導水能力的提高以及土壤可利用養(yǎng)分的增加[14]。研究表明,土壤結構在時間上還會受到濕潤與干燥循環(huán)、生物活性以及農業(yè)作業(yè)的影響[15]。

    雖然添加生物炭對土壤水力學特性產生不同程度的影響,但有關預測生物炭定量改善土壤水分特性、生物炭決定水分分布的確切解的研究[16]鮮見報道。有些研究給出了土壤水動力的空間變化,主要是土壤粒子分布的不均勻性[17]。土壤中的水分運動,取決于固體顆粒的尺寸、形狀和分布[18]。土壤的多樣性以及水分運動的復雜性,使得試驗的數量和耗時巨大。本文推導添加生物炭土壤的毛細管模型[19],并在此模型下利用有限元法或差分法對土壤進行數值模擬。對比此模型與添加生物炭土壤土柱滲透的數值分析結果和試驗結果,以證實此模型預測不同條件下的土壤水分分布或者校驗試驗結果的正確性。

    1 材料與方法

    1.1 研究區(qū)概況

    取土地點位于黑龍江省農墾北安管理局紅星農場實驗基地(48°10′E,127°1′N)。該地區(qū)為典型的東北黑土區(qū),屬寒溫帶大陸性季風氣候,冬季寒冷干燥,夏季高溫多雨,降雨時間集中,年降水量多介于450~650 mm,其中7—9月累積降水量約402 mm,約占年平均降水量的75%,因此旱澇災害頻繁,水土流失嚴重。該區(qū)域的土壤多為壤土和粘壤土,以粗粉沙和粘粒最多,各占30%~40%,容重1.0~1.3 g/cm3,透水、持水、通氣性均較好,總孔隙度多在40%~60%,毛管孔隙可占20%~30%,通氣孔隙度約占20%。

    1.2 試驗材料與設備

    供試材料為玉米秸稈生物炭,購自遼寧金和福農業(yè)開發(fā)有限公司,制備方式為450℃無氧條件下高溫裂解?;纠砘再|:粒徑范圍1.5~2.5 mm,其中2 mm粒徑占60%以上。各指標如下:全碳質量分數72.21%,全氮質量分數1.56%,全磷質量分數0.72%,全鉀質量分數1.64%,pH值為9.35,經測量生物炭密度約596 kg/m3。

    脫水面半徑為9.8 cm的CR- 21G3型高速離心機,用于室內測定土壤水分特征曲線;自制直徑8 cm、長為80 cm的有機玻璃筒,距供水端10、20、30、40、50、60、70 cm處設有取土孔,用于測量不同位置的土壤含水率。

    1.3 試驗方法與內容

    土樣為田間試驗土樣。取土小區(qū)面積為20 m×5 m,坡度3°。設置5種生物炭用量添加比例,即不添加生物炭處理(CK),以及生物炭添加量為25、50、75、100 t/hm2的處理。試驗小區(qū)采用相同的水肥管理,每個處理重復3次。2017年5月中旬將生物炭人工均勻混入耕層黑土(0~20 cm),生物炭占土壤的體積比約為2%(C1)、4%(C2)、6%(C3)、8%(C4),取樣時間為2017年9月25日。用環(huán)刀采樣器取對照土樣(CK),取土深度為10~15 cm,環(huán)刀高5 cm、體積100 cm3,取土后將環(huán)刀用蒸餾水吸水48 h,至飽和后稱量。生物炭放置環(huán)刀中,用蒸餾水吸水并壓實48 h后稱量。用離心機法分別測定土壤和生物炭的水分特征曲線。離心機轉速設定為500~6 000 r/min,離心時間為100 min[20],離心結束后稱量,土壤與生物炭各重復3次。

    室內試驗用定水頭法分別測定生物炭和5種處理的飽和導水率。利用馬氏瓶供水,維持水頭恒定為5 cm,每隔10 min記錄出流水量,出流水量保持穩(wěn)定時結束試驗,計算飽和導水率,各處理重復3次。

    測定土柱的含水率。在5個處理小區(qū)取原狀土,取土深度10~15 cm,土樣經風干、破碎、篩分后,實地測量土壤干容重為1.1 g/cm3。將土樣按此容重每層3 cm分層回填入土柱,裝填好的土柱直徑8 cm、高80 cm。進水端采用馬氏瓶供水,視進水端的土壤為飽和,當濕潤鋒到達土柱約3/4處,在各取土孔處迅速取土,測定各處的含水率。各處理重復3次。

    1.4 數據處理方法

    采用Excel 2010進行數據整理,所用數據均為各重復測定的平均值,土壤水分特征曲線采用Retc軟件進行擬合,采用python語言建模編程,采用通用的有限元分析軟件Abaqus 6.13進行數值模擬計算。

    2 理論分析

    2.1 土壤水分特征曲線

    目前學者大多采用VAN- GENUCHTEN(VG)模型給定土壤水分特征曲線(SWCC)。模型由VAN GENUCHTEN[21]在1980年提出,適用于描述低含水率時的土壤水分特征曲線。其表達式為

    (1)

    其中

    m=1-1/n(n>1)

    式中θr——土壤殘余含水率,cm3/cm3

    θs——土壤飽和含水率,cm3/cm3

    h——土壤吸力,cm

    a、n、m——曲線擬合參數

    θ——土壤含水率

    本研究采用VAN- GENUCHTEN模型,利用試驗測得土壤含水率,擬合土壤飽和含水率和殘余含水率。

    2.2 添加生物炭土壤的飽和導水率

    2.2.1假設生物炭在土壤中呈層狀分布的導水率

    穩(wěn)定水流分別穿過土壤和生物炭時,設K1為土壤飽和導水率,K2為生物炭飽和導水率。當生物炭在土壤中呈層狀分布時,設生物炭與土壤的混合土層的總厚度為L,令混合土層中土壤體積比為α、生物炭體積比為β,則α+β=1,因此土壤層的厚度為αL,生物炭層厚度為βL,結構如圖1a所示。設土壤層與生物炭層水力梯度分別為i1、i2,整個混合土層有效導水率為Ke,整個混合土層水頭損失應等于兩層水頭損失之和,即

    iL=αi1L+βi2L

    (2)

    式中i——總的水力梯度

    JW=Kei=K1i1=K2i2

    (3)

    式中JW——水流通量

    結合式(2)可以得到

    (4)

    圖1 生物炭分布示意圖Fig.1 Sketches of media in different directions

    一般情況下,土壤體積比α遠大于生物炭體積比β,另外,生物炭的飽和導水率K2也遠大于草甸黑土的飽和導水率K1,這樣式(4)可以進一步簡化為

    (5)

    當生物炭在土壤中層狀分布時,混合土壤的飽和導水率可以近似地認為是土粒的飽和導水率除以土粒在土壤中的比例。

    2.2.2假設生物炭在土壤中呈條狀分布的導水率

    當生物炭在土壤中條狀分布時,土壤區(qū)的面積視為αA,生物炭區(qū)的面積為βA,結構如圖1b所示。整個混合土壤截面的有效導水率為Ke、水力梯度為i,整個混合土層的水流通量應等于兩區(qū)域水流通量之和,即

    KeAi=K1αAi+K2βAi

    (6)

    Ke=αK1+βK2

    (7)

    當生物炭在土壤中呈條狀分布時,混合土壤的飽和導水率可以近似地用式(7)估計。

    由于生物炭在土壤中是隨機均勻分布的,雖然式(5)和式(7)均不能給出飽和導水率的準確值,但通過這2個公式大致能夠估算出飽和導水率。

    2.3 毛細管模型推導添加生物炭土壤的導水率

    前人以試驗為基礎,提出很多測定土壤水分運動參數的方法。由于一些模型參數較難獲得[22],因此限制了水分運動的數值模擬計算。本文利用相對容易獲得的土壤水分特征曲線,推導了添加生物炭土壤的毛細管理論,給出了土壤水分特征曲線和非飽和導水率的關系。

    當水與毛細管管壁完全浸潤時,水在毛細管內上升的最大高度為

    (8)

    式中ρW——水的密度

    H——水在毛細管內上升的最大高度

    g——重力加速度

    R——毛細管的曲率半徑

    σ——水的表面張力系數

    假定添加生物炭的土壤可以用不同形態(tài)的毛細管束來描述,每個毛細管束的長度為Le,土壤的實際長度為L,土壤的橫截面積為A,土壤兩端的水頭差為ΔH。根據泊肅葉公式,對于曲率半徑為Rj的毛細管,其流量為

    (9)

    式中μ——液體粘滯系數

    當土壤中的所有毛細管充滿水時,水流通量為

    (10)

    其中

    式中Nj——半徑為Rj的毛細管的數量

    P——不同半徑毛細管的數量

    nj——單位面積內半徑為Rj的毛細管數量

    由于土壤的基質吸力與土壤的毛細管孔徑有關,因此土壤水分特征曲線也描述了毛細管模型的尺寸分布規(guī)律。通過未添加生物炭土壤(CK)的水分特征曲線以及生物炭的水分特征曲線,能夠給出添加生物炭后土壤的非飽和導水率的范圍。

    根據土壤的含水率將土壤水分特征曲線畫成M個相同寬度的區(qū)間,如圖2所示,每個區(qū)間寬度為Δθ,每一個區(qū)域對應一個基質吸力hj(即與該吸力對應的所有毛細管都具有相同孔徑)和土壤含水率θs-jΔθ,即

    hj=h(θs-jΔθ)

    (11)

    圖2 計算非飽和導水率的毛細管模型Fig.2 Capillary model for calculating unsaturated hydraulic conductivity

    另外,由于生物炭和土壤顆粒共同在一起,所以同一位置處吸力相同,即在hj吸力下,根據生物炭的水分特征曲線θW(h)得到生物炭中的含水率差值為

    ΔθWj=θW(hj-1)-θW(hj)

    (12)

    θ(h0)=θWS,θWS為生物炭的飽和含水率。根據式(8),基質吸力hj對應的毛細管孔徑為

    (13)

    對于土壤和生物炭,對應的j區(qū)的單位面積內毛細管數量分別為

    (14)

    (15)

    式中nk——j區(qū)生物炭內單位面積的毛細管數量

    ΔθWk——生物炭的含水率區(qū)間寬度

    Rk——生物炭內毛細管半徑

    當h=hj時,所有孔徑小于Rj的毛細管充滿了水,提供了土壤的水流通量,并且這些毛細管是由土壤顆粒和生物炭顆粒共同提供的。另外,一般認為[23]生物炭的導水率遠大于土壤的導水率,意味著生物炭中毛細管的密度大于土壤。如果假定生物炭中多于土粒中毛細管的倍數可以由飽和導水率的比值近似給出,再同時考慮式(10),得到混合土壤的水流通量表達式為

    (16)

    其中

    q=Kcs/Kss

    式中q——生物炭與土壤飽和導水率的比值

    Kcs——生物炭飽和導水率

    Kss——土壤飽和導水率

    再由式(13)~(15),得到

    (17)

    式中hk——區(qū)間k內毛細管的吸力

    N——生物炭水分特征曲線劃分的最大區(qū)間數

    (18)

    其中

    l=L/Le

    式中l(wèi)——毛細管彎曲度

    式(18)與達西定律相比,得到混合土壤的非飽和導水率為

    (19)

    當毛細管全部充滿水,就得到混合土壤的飽和導水率

    (20)

    結合式(19)和式(20)得到混合土壤的非飽和導水率

    (21)

    基于土壤和生物炭的水分特征曲線,再根據添加生物炭的比例以及土壤和生物炭的飽和導水率,可以用式(21)近似給出添加不同量生物炭后土壤的非飽和導水率。

    3 結果與分析

    3.1 土壤水分特征曲線

    用離心機法測定的土壤和生物炭的水分特征曲線如圖3所示。

    圖3 土壤與生物炭的水分特征曲線Fig.3 Water characteristic curves of soil and biochar

    采用VAN- GENUCHTEN模型,利用Retc軟件對離心機法的實測數據進行擬合,擬合的飽和含水率為0.418 1 cm3/cm3,殘余含水率0.009 8 cm3/cm3。由圖3可以看出,當土壤吸力小于30 cm時,生物炭與土壤幾乎平行,當土壤吸力介于30~2 200 cm,土壤的SWCC位于生物炭上方,當吸力大于2 200 cm時,生物炭的SWCC位于土壤上方,即生物炭在中高吸力段的含水率大于土壤,因此添加生物炭后土壤持水能力的增加,主要依靠的是生物炭在中高吸力段持水能力的提升。也有文獻表明,生物炭的添加減小了土壤容重、增加了土壤孔隙度[24],使得土壤持水能力增強。

    3.2 毛細管模型計算土壤的非飽和導水率

    由于毛細管模型綜合考慮了生物炭和土壤的性能,因此根據毛細管模型公式(21)以及測定的土壤和生物炭的水分特征曲線,可以計算不同添加量生物炭土壤的非飽和導水率。采用圖3數據計算毛細管模型的數據如表1所示。計算中采用的生物炭和土壤的飽和導水率實際測量值見表2。由表1可以看出,[K(θs-jΔθ)/Ks]1/3項可以近似看作有效飽和度。在低吸力段,有效飽和度隨著生物炭含量的增加從0.710 446增加到0.739 290,有效飽和度的增高可看作土壤導水性能得到提高[20];在35.82~1 590 cm吸力段,有效飽和度隨生物炭含量的增加而減小,土壤導水性能降低。導水性能的增加是由于生物炭在低吸力下的高導水率,而高吸力下導水性能的降低,是由于生物炭在高吸力下的低導水率。

    3.3 土壤飽和導水率測量結果

    表2為實際測量的CK、C1、C2、C3、C4的飽和導水率,以及采用式(5)和式(7)計算的估計值。

    表1 采用毛細管模型計算不同添加量生物炭土壤的非飽和導水率Tab.1 Calculated unsaturated hydraulic conductivity by using capillary model

    表2 土壤飽和導水率的測量值和估計值Tab.2 Tested and estimated values of soil saturatedhydraulic conductivity mm/h

    由表2可以看出,C1、C2、C3、C4飽和導水率的實際測量值,均在層狀分布和條狀分布的估計值之間,但非常接近于層狀分布,與層狀分布相比的最大偏差為8.9%。分析其原因,可能是當土壤中的生物炭施加量較少時,不足以讓生物炭粒子連接起來,所以應該看作層狀。另外,試驗數據表明[16],添加生物炭土壤中的水分流動,多為水在生物炭和土粒之間的空隙流動。由于層狀分布的估計值沒有考慮生物炭和土粒之間的空隙,當生物炭添加后,空隙間的水分流動會對飽和導水率產生影響,因而導致實際測量的飽和導水率比層狀估算值大。

    4 添加生物炭土壤滲水過程的數值分析

    4.1 添加生物炭后的土壤模型

    為研究相同滲透時間后,不同生物炭添加量土樣中的水分分布,首先要對添加生物炭的土壤進行數值模擬。由于實際土壤的復雜性,計算的模型擬采用二維形式。添加生物炭的土壤被分成土壤區(qū)和生物炭區(qū)。計算時,這兩個區(qū)具有各自的非飽和導水率和初始水分特征。

    由于生物炭是隨機分布在土壤中,采用Python語言編程的方法給定土壤截面。將生物炭隨機分布在土壤截面中的步驟如下:首先建立計算土壤的面積區(qū)域[25]。面積選擇過大,生物炭顆粒的數量隨之變大,導致計算量過大難以實現。因此,計算的面積選取10 mm×157 mm。其次確定土樣中生物炭的個數,生物炭的個數通過添加生物炭的比例來校對。本研究中,添加生物炭的體積比約為2%、4%、6%、8%。根據添加比例和試驗所用生物炭的平均粒徑給出生物炭的個數(生物炭粒徑范圍1.5~2.5 mm,其中2 mm粒徑約占60%以上,故取2 mm)。隨機產生生物炭所在位置的圓心坐標,畫圓代表生物炭。校驗是否有圓相交,如果有,刪除后再補充相應數量的圓。最后,對生成的所有生物炭顆粒求面積,核算添加的比例和實際摻炭比例是否相等。完成上述工作后賦予土壤區(qū)和生物炭區(qū)不同的水分運動參數。圖4為最終生成的隨機添加生物炭的土壤模型。

    圖4 隨機分布生物炭顆粒的土壤截面 Fig.4 Profiles of biochar particles distributed randomly in soil

    4.2 有限元采用的模型

    采用通用的有限元分析軟件Abaqus進行數值計算。 對于飽和土,滲透系數可以視為常量,但對于液體與氣體并存的非飽和土,飽和度的大小直接影響滲透的結果,因此滲透系數K是飽和度S的函數。Abaqus采用BROOKS等[26]提出的滲透系數和飽和度的關系。利用土壤水分特征曲線給出的3個參數,即土壤的進氣值、有效飽和度以及空隙尺寸分布指標λ,滲透系數函數與有效飽和度的關系為(下文中出現的飽和度均指有效飽和度)

    (22)

    其中

    式中Sr——有效飽和度

    δ——經驗常數,取3

    (ua-uw)——基質吸力

    (ua-uw)b——進氣值

    4.3 實際土樣與毛細管模型土樣的數值計算

    為研究不同土樣經過相同滲透時間后的水分分布,數值模擬部分首先計算添加生物炭土壤的飽和度,其次采用前面推導的毛細管模型等效的非飽和導水率計算摻炭土壤的飽和度,初始給定的土壤吸力是相同的。

    當區(qū)域為對照土樣(CK)時,此時α=1,采用的吸力和對應的飽和度如表1所示,實際上這是采用單一土壤的毛細管模型給出土壤的非飽和導水率,飽和導水率Ks采用試驗測得的值1.152 mm/h。邊界條件為:左邊界初始吸力為零表示進水,其他邊界為不透水邊界。全部區(qū)域的初始吸力和飽和度分別為200 cm、0.152 7,滲透時間5 000 s,計算出CK的飽和度隨位置的分布規(guī)律如圖5a所示。

    當區(qū)域為C4時,區(qū)域如圖4c所示,采用的吸力和飽和度如表1所示,土壤區(qū)為α=1,生物炭區(qū)為α=0。土壤和生物炭的飽和導水率分別為1.152、22.068 mm/h,邊界條件仍為左邊界初始吸力為零,土壤的初始吸力和飽和度為200 cm、0.152 7,生物炭的初始吸力和飽和度為為200 cm、0.060 7,滲透時間5 000 s,計算出的C4的飽和度隨位置的分布規(guī)律如圖5b所示。

    當區(qū)域采用C4的毛細管模型C4*時,采用的吸力和飽和度如表1所示。飽和導水率為試驗測得的1.371 6 mm/h,邊界條件仍為左邊界初始吸力為零,C4*初始吸力和飽和度為200 cm、0.086 5,滲透時間5 000 s,計算后C4*的飽和度隨位置的分布規(guī)律如圖5c所示。圖5表明,毛細管模型很好地描述了施加生物炭后,土壤飽和度隨位置的分布規(guī)律。

    計算出的CK、C4中的土壤、C4中的生物炭以及C4*的飽和度隨位置的分布規(guī)律如圖6a所示。

    圖6 滲透5 000 s后的土壤飽和度曲線Fig.6 Saturation of soil after infiltration of 5 000 s

    由圖6a可以看到,C4*的飽和度分布更類似于C4中生物炭的飽和度分布形式,這是由于采用毛細管模型(式(21))的推導結果。由于生物炭中毛細管的密度遠多于土壤,與同體積的土壤相比,能夠吸收更多的水分,因此對土壤的飽和度影響更大。式(21)中的q值決定了C4*飽和度曲線的形狀。圖6a中,C4*的飽和度曲線在水平距離上略微超前于C4中生物炭的飽和度曲線,說明計算C4*飽和度時,采用測量的飽和導水率偏大。這可能是由于C4的計算中,沒有考慮生物炭和土壤的空隙,而這是影響飽和導水率的重要因素。

    設定滲透時間5 000 s,按上述方法計算C1、C2、C1*、C2*的飽和度分布規(guī)律。C1*、C2*飽和導水率為測得的1.195 2、1.249 2 mm/h。毛細管模型采用的土壤初始吸力均為200 cm,C1*和C2*的飽和度分別為0.114 4和0.099 9。圖6b為滲透時間5 000 s后土樣的飽和度隨位置的分布曲線。C1*、C2*為毛細管模型模擬的土樣(由于圖線過于密集,C3沒有畫出)。

    由圖6可以看到,水分滲透一定時間后,C1*、C2*、C4*曲線均落在CK上方,即相同位置上添加生物炭土壤的飽和度都比CK大,這說明添加生物炭可以提高土壤的導水率。并且導水率隨著生物炭添加量的增大而增大。C4*曲線落在CK下方,是由于高吸力下生物炭導水性能的降低造成的。毛細管模型綜合考慮了生物炭和土壤的特性,能較好地反映添加生物炭后土壤的導水特性。

    4.4 土柱的數值計算和試驗數據的比較

    本試驗采用的土柱尺度較大,因此在有限元計算時,如果考慮實際的生物炭粒子,劃分的網格必須很密,這樣會導致計算量過大,所以一般無法對土柱直接計算,這時可以采用毛細管模型近似計算。計算時,由于C1*、C2*、C4*的飽和含水率和殘余含水率各不相同,因此C1*、C2*、C4*的初始飽和度是不同的。為簡化計算,給定C1*、C2*、C4*的初始飽和度,與試驗測得CK的初始飽和度相同。

    為驗證毛細管模型的正確性,設定水分滲透時間為24 h,實際測量了C1、C2、C4在80 cm土柱中的含水率,采用CK擬合的飽和含水率和殘余含水率值計算飽和度。計算的初始數據如下:實際測定的CK初始飽和度0.294,與之對應的土壤的吸力約為68.6 cm。毛細管模型在此飽和度下的吸力分別為51.9 cm(C1*)、43.4 cm(C2*)、40.4 cm(C4*)。圖7為水分滲透24 h后,飽和度隨土柱位置變化關系圖。

    圖7 滲透24 h后土柱飽和度的計算值和測量值Fig.7 Calculated and measured values of soil saturation after infiltration of 24 h

    由圖7的計算結果看出,土柱中水分滲透24 h后,對于同一生物炭添加比例,毛細管模型的飽和度與實際測量值相比略低。這是由模型中沒有考慮實際土壤中生物炭粒子和土壤之間的空隙造成的。生物炭體積比為2%土壤的誤差范圍為4.27%~7.24%,4%土壤的誤差范圍為3.06%~8.6%,8%土壤的誤差范圍為2.90%~10.11%。可見,隨著生物炭體積比的增大,數值計算結果與毛細管模型的計算結果偏差有所增加,這種偏差是因為在毛細管模型假設中,生物炭中的不同孔徑的毛細管數量是均勻分布的,總是用一個恒定的飽和導水率的比值近似給出,這與實際是有偏差的。由圖7結果看出,在距供水端較遠處,即土壤吸力較大處,土壤的實際導水率大于毛細管模型的值,說明生物炭中的細小毛細管的數量多于估計值。因此,采用毛細管模型的數值計算,對于中等及以下生物炭施用量,在土壤水力特征的定量分析上更加準確。毛細管模型可以簡化土壤的水力特性,可以用來對實際的大塊土地進行數值分析。

    5 討論

    本研究針對東北黑土區(qū)的草甸黑土,推導了施加生物炭后土壤的毛細管模型,利用此模型,研究了一定滲透時間后土壤中的水分分布,選擇飽和度等參數將添加不同量生物炭的土壤進行數值模擬,最后利用土柱的實際測定結果對數值模擬結果進行驗證,以此開展生物炭對土壤水分影響的數值化定量化研究。

    從土壤水分特征曲線看出,添加生物炭后土壤的持水能力增加,由于生物炭在中高吸力段的含水率大于參照組土壤,因此添加生物炭后土壤持水能力的增加,主要依靠的是生物炭在中高吸力段持水能力的提升。毛細管模型表明,施加生物炭后土壤性能的改善,是由于低吸力段生物炭的良好導水特性。利用毛細管模型計算的土壤非飽和導水率,依賴于土壤的水分特征曲線和飽和導水率,為得到精確的數值計算結果,必須精確測量土壤水分特征曲線。

    土柱一端加水滲透24 h后,數值模擬的土壤水分分布值比實際測量值略低,主要是由于毛細管模型中生物炭多于土壤的毛細管數量,總是用一個恒定的值近似給出,這與實際有偏差。另外是由于模型中生物炭與土壤的空隙分布不能與實際情況完全吻合造成的。

    雖然直接測定的結果更接近真實值,但測量過程繁瑣、消耗時間長,因此采用毛細管模型的有限元數值計算優(yōu)勢明顯。另外,對于只能通過試驗研究解決的問題,如一定范圍的生物炭施加量、不同的施加面積、地勢情況(如坡耕地等)、不同的滲透時間等情況,數值模擬能夠更好更快地給出近似答案。采用毛細管模型還可以對試驗數據進行校驗,找出試驗中有偏差的數據,有助于提高試驗結果的準確性,從而為研究土壤的水力特性提供更多的方法與保障。

    6 結論

    (1)試驗測得的水分特征曲線表明,當土壤吸力大于2 200 cm時,生物炭的殘余含水率高于對照組,因此添加生物炭后土壤持水能力的增加,主要依靠的是生物炭在中高吸力段持水能力的提升。

    (2)通過飽和導水率試驗數據和理論估計值的比較可以看出,實際測量的飽和導水率與生物炭層狀分布的飽和導水率非常接近,最大偏差為8.9%,說明添加生物炭土壤更應考慮為層狀分布。試驗數據和理論數據的差別可能是生物炭和土壤之間的空隙增加造成的。

    (3)毛細管模型結果顯示,添加生物炭土壤導水率的提高,主要是由生物炭的導水特性決定的,添加的生物炭主要改變了土壤在低吸力下的導水率。

    (4)數值分析顯示,施用生物炭能夠提高東北草甸黑土的導水率。添加生物炭的比例越高,導水率增加越大。

    (5)生物炭中的不同孔徑的毛細管數量不是均勻分布的,孔徑越小,毛細管的數量越多。

    (6)毛細管模型簡單有效,生物炭體積比小于8%時,數值分析的結果與試驗數據吻合良好,說明該模型可以為土壤水分運動參數的數值計算提供比較準確的方法。

    1 何緒生, 耿增超, 佘雕, 等. 生物炭生產與農用的意義及國內外動態(tài)[J]. 農業(yè)工程學報, 2011,27(2):1-7.

    HE Xusheng, GENG Zengchao, SHE Diao, et al. Implications of production and agricultural utilization of biochar and its international dynamic[J]. Transactions of the CSAE, 2011, 27(2):1-7. (in Chinese)

    2 LEHMANN J, JOSEPH S. Biochar for environmental management: an introduction [J]. Science and Technology, 2009, 25(1):15801-15811.

    3 GITHINJI L. Effect of biochar application rate on soil physical and hydraulic properties of a sandy loam[J]. Archives of Agronomy and Soil Science, 2013:1-14.

    4 STREUBEL J D, COLLINS H P, GARCIA P M, et al. Influence of contrasting biochar types on five soils at increasing rates of application[J]. Soil Biology & Biochemistry,2011,75:1402-1413.

    5 WEST L T, ABREW M A, BISHOP J P. Saturated hydraulic conductivity of soils in the Southern Piedmont of Georgia, USA: field evaluation and relation to horizon and landscape properties[J]. CATENA, 2008, 73(2): 174-179.

    6 ABEL S, PETERS A, TRINKS S, et al. Impact of biochar and hydrochar addition on water retention and water repellency of sandy soil[J]. Geoderma, 2013, 202-203:183-191.

    7 魏永霞, 劉志凱, 馮鼎銳, 等. 生物炭對草甸黑土物理性質及雨后水分動態(tài)變化的影響[J/OL]. 農業(yè)機械學報, 2016, 47(8): 201-207. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160825&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2016.08.025.

    WEI Yongxia, LIU Zhikai, FENG Dingrui, et al. Influences of the biochar on physical properties of the meadow black soil and dynamic changes of soil water after the individual rainfall[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(8):201-207.(in Chinese)

    8 LI Y, GAO R, YANG R, et al. Using a simple soil column method to evaluate soil phosphorus leaching risk[J]. CLEAN-Soil, Air, Water, 2013, 41(11):1100-1107.

    9 ZIMMERMAN A R. Abiotic and microbial oxidation of laboratory-produced black carbon (biochar)[J]. Environmental Science & Technology, 2010, 44(4): 1295-1301.

    10 CHAN K Y, VAN ZWIETEN L, MESZAROS L, et al. Agronomic values of greenwaste biochar as a soil amendment[J]. Australian Journal of Soil Research,2007, 45:629-634.

    11 FENG L, LI G T, LIN Q M, et al. Crop yield and soil properties in the first 3 years after biochar application to a Calcareous soil[J]. Journal of Integrative Agriculture,2014, 13(3): 525-532.

    12 勾芒芒, 屈忠義, 楊曉, 等. 生物炭對砂壤土節(jié)水保肥及番茄產量的影響研究[J/OL]. 農業(yè)機械學報, 2014, 45(1): 137-142.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140122&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.01.022.

    GOU Mangmang, QU Zhongyi, YANG Xiao, et al. Study on the effects of biochar on saving water, preserving fertility and tomato yield[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(1):137-142. (in Chinese)

    13 AKHTAR S S, LI G, ANDERSEN M N, et al. Biochar enhances yield and quality of tomato under reduced irrigation[J]. Agricultural Water Management, 2014, 138:37-44.

    14 ASAI H, SAMSON B K, STEPHAN H M, et al. Biochar amendment techniques for upland rice production in Northern Laos 1. Soil physical properties, leaf SPAD and grain yield[J]. Field Crops Research, 2009, 111(1): 81-84.

    15 BODNER G, SCHOLL P, KAUL H P. Field quantification of wetting-drying cycles to predict temporal changes of soil pore size distribution[J]. Soil & Tillage Research, 2013, 133(5):1-9.

    16 LIU Z, DUGAN B, MASIELLO C A, et al. Biochar particle size, shape, and porosity act together to influence soil water properties.[J]. Plos One, 2017, 12(6):e0179079.

    17 ALLETTO L, COQUET Y, VACHIER P, et al. Hydraulic conductivity, immobile water content, and exchange coefficient in three soil profiles[J]. Soil Science Society of America Journal, 2006, 70(4):1272-1280.

    18 SHE D, GAO X, WANG P, et al. Comparison of soil hydraulic properties with different levels of soil salinity and sodicity[J]. Arabian Journal of Geosciences, 2015, 8(8):5351-5360.

    19 ROGER E J, RICHARD G, DEXTER A R, et al. Integration of soil structure variations with time and space into models for crop management: a review[J]. Agronomy for Sustainable Development, 2009, 29(1):135-142.

    20 尚熳廷,馮杰,劉佩貴,等. SWCC 測定時吸力計算公式與最佳離心時間的探討[J]. 河海大學學報, 2009, 37(1): 12-15.

    SHANG Manting, FENG Jie, LIU Peigui, et al. On formulas for soil suction and optimum centrifugal time by use of SWCC[J]. Journal of Hohai University, 2009, 37(1): 12-15. (in Chinese)

    21 VAN GENUCHTEN M T H. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5):892-898.

    22 BARNES R T, GALLAGHER M E, MASIELLO C A, et al. Biochar-induced changes in soil hydraulic conductivity and dissolved nutrient fluxes constrained by laboratory experiments[J]. Plos One, 2014, 9(9):e108340.

    23 YARGICOGLU E N, SADASIVAM B Y, REDDY K R, et al. Physical and chemical characterization of waste wood derived biochars[J]. Waste Management, 2015, 36:256-268.

    24 MELO C M T, HOLANDA N M, MADARI B E, et al. Biochar increases plant-available water in a sandy loam soil under an aerobic rice crop system[J].Solid Earth,2014,5:939-952.

    25 WONGSTO A, LI S. Micromechanical FE analysis of UD fibre-reinforced composites with fibres distributed at random over the transverse cross-section[J]. Composites Part A: Applied Science & Manufacturing, 2005, 36(9):1246-1266.

    26 BROOKS R H, COREY A T. Hydraulic properties of porous media and their relation to drain design[J].Transactions of the ASAE, 1964,7(1):26-28.

    猜你喜歡
    導水率吸力毛細管
    陜北山地蘋果園土壤飽和導水率和植物導水率特征
    深水大型吸力錨測試技術
    基于線性源法與圖像處理的土壤飽和導水率快速測量方法
    ROV在海上吸力樁安裝場景的應用及安裝精度和風險控制
    化工管理(2022年11期)2022-06-03 07:08:24
    毛細管氣相色譜法測定3-氟-4-溴苯酚
    云南化工(2020年11期)2021-01-14 00:50:54
    深水吸力樁施工技術研究
    海城市土壤飽和導水率空間分布特征
    超聲萃取-毛細管電泳測定土壤中磺酰脲類除草劑
    雜草學報(2015年2期)2016-01-04 14:58:05
    毛細管氣相色譜法測定自釀葡萄酒中甲醇的含量
    中藥與臨床(2015年5期)2015-12-17 02:39:28
    超強吸力
    少年科學(2015年7期)2015-08-13 04:14:32
    一个人看视频在线观看www免费 | 国产一区二区三区视频了| 亚洲 国产 在线| 亚洲精品在线观看二区| 他把我摸到了高潮在线观看| 日韩人妻高清精品专区| 黄色女人牲交| 国产一区二区三区在线臀色熟女| 亚洲九九香蕉| 欧美色视频一区免费| 成人鲁丝片一二三区免费| 成熟少妇高潮喷水视频| 99精品在免费线老司机午夜| 欧美乱色亚洲激情| 欧美三级亚洲精品| 亚洲美女黄片视频| 亚洲熟妇中文字幕五十中出| av视频在线观看入口| 禁无遮挡网站| 久久久国产成人精品二区| 亚洲avbb在线观看| 国产成人精品久久二区二区免费| 成人特级黄色片久久久久久久| 欧美一区二区精品小视频在线| 神马国产精品三级电影在线观看| 在线a可以看的网站| 国产伦在线观看视频一区| 成人一区二区视频在线观看| 国产人伦9x9x在线观看| 亚洲欧美激情综合另类| 亚洲国产精品成人综合色| 波多野结衣高清作品| 午夜福利在线在线| 一区二区三区国产精品乱码| 99久久综合精品五月天人人| 欧美成人一区二区免费高清观看 | 久9热在线精品视频| 精品久久久久久久毛片微露脸| 日韩欧美免费精品| 亚洲中文av在线| 少妇裸体淫交视频免费看高清| 亚洲国产精品成人综合色| 青草久久国产| 亚洲五月婷婷丁香| 免费看十八禁软件| 亚洲精华国产精华精| 免费看十八禁软件| 黄片大片在线免费观看| 色播亚洲综合网| 欧美xxxx黑人xx丫x性爽| 在线免费观看不下载黄p国产 | 性欧美人与动物交配| 精品久久久久久久久久久久久| 国产精品野战在线观看| 国产精品国产高清国产av| 国产一区在线观看成人免费| 天天添夜夜摸| 嫁个100分男人电影在线观看| 两人在一起打扑克的视频| 高清在线国产一区| 老汉色∧v一级毛片| 国产精品久久久av美女十八| 午夜福利免费观看在线| 成人av在线播放网站| 熟女少妇亚洲综合色aaa.| 成人欧美大片| 亚洲 欧美 日韩 在线 免费| 亚洲乱码一区二区免费版| 国产主播在线观看一区二区| 十八禁人妻一区二区| 亚洲成a人片在线一区二区| 国产综合懂色| 一二三四在线观看免费中文在| 好看av亚洲va欧美ⅴa在| 亚洲自偷自拍图片 自拍| 午夜精品一区二区三区免费看| 99国产精品一区二区三区| 色哟哟哟哟哟哟| 丁香欧美五月| 亚洲性夜色夜夜综合| 男人和女人高潮做爰伦理| 性欧美人与动物交配| 欧美精品啪啪一区二区三区| 国产野战对白在线观看| 成人特级av手机在线观看| 在线观看美女被高潮喷水网站 | 超碰成人久久| 又粗又爽又猛毛片免费看| 日本精品一区二区三区蜜桃| 欧美色视频一区免费| 欧美国产日韩亚洲一区| 亚洲国产精品合色在线| 免费看a级黄色片| 国产高清激情床上av| АⅤ资源中文在线天堂| 久久久久久大精品| 色视频www国产| 麻豆一二三区av精品| 国产视频一区二区在线看| 亚洲成人精品中文字幕电影| 精品久久久久久久毛片微露脸| 我要搜黄色片| 美女高潮的动态| 国产探花在线观看一区二区| 不卡一级毛片| 亚洲国产欧美人成| 一级作爱视频免费观看| 亚洲精品乱码久久久v下载方式 | 99热这里只有是精品50| 一二三四在线观看免费中文在| 亚洲国产精品999在线| 亚洲av美国av| 久久久久性生活片| 久久久久久久久中文| 动漫黄色视频在线观看| 丰满的人妻完整版| 日本 欧美在线| 中文字幕精品亚洲无线码一区| 99久久综合精品五月天人人| 国产视频一区二区在线看| 久久精品aⅴ一区二区三区四区| 听说在线观看完整版免费高清| 欧美成狂野欧美在线观看| 亚洲av电影不卡..在线观看| 国产亚洲精品一区二区www| 无限看片的www在线观看| 又粗又爽又猛毛片免费看| 深夜精品福利| 成年人黄色毛片网站| 亚洲精品乱码久久久v下载方式 | 天堂av国产一区二区熟女人妻| 国产视频一区二区在线看| av欧美777| 国产精品女同一区二区软件 | 天天添夜夜摸| 国产亚洲欧美98| 看黄色毛片网站| 国产高潮美女av| 色综合亚洲欧美另类图片| 免费大片18禁| 女人被狂操c到高潮| 啦啦啦韩国在线观看视频| 小蜜桃在线观看免费完整版高清| 国产高清视频在线播放一区| 美女高潮喷水抽搐中文字幕| 午夜福利视频1000在线观看| 麻豆一二三区av精品| 特级一级黄色大片| 精品电影一区二区在线| 十八禁网站免费在线| 一边摸一边抽搐一进一小说| 色尼玛亚洲综合影院| 婷婷亚洲欧美| 久久亚洲精品不卡| 国产av不卡久久| 亚洲欧美日韩高清专用| 亚洲欧美日韩高清在线视频| 变态另类成人亚洲欧美熟女| 国产成年人精品一区二区| 最近最新免费中文字幕在线| 精品久久久久久久久久久久久| 搞女人的毛片| 国产精品久久久久久亚洲av鲁大| 嫩草影院精品99| 操出白浆在线播放| 麻豆一二三区av精品| 国产高清视频在线播放一区| 黄色视频,在线免费观看| 国产美女午夜福利| 美女扒开内裤让男人捅视频| 国产亚洲av高清不卡| 在线视频色国产色| 成人午夜高清在线视频| 免费看日本二区| 91九色精品人成在线观看| 国产精品综合久久久久久久免费| 九色成人免费人妻av| 男插女下体视频免费在线播放| 国产高清videossex| 国产精品一区二区免费欧美| 亚洲熟女毛片儿| 久久亚洲精品不卡| 啪啪无遮挡十八禁网站| 亚洲男人的天堂狠狠| 成人永久免费在线观看视频| 99re在线观看精品视频| 制服丝袜大香蕉在线| 岛国在线免费视频观看| 一卡2卡三卡四卡精品乱码亚洲| 欧美精品啪啪一区二区三区| 97人妻精品一区二区三区麻豆| 精品久久久久久久毛片微露脸| 精品欧美国产一区二区三| 两个人看的免费小视频| 久久香蕉精品热| 一进一出好大好爽视频| a在线观看视频网站| 欧美丝袜亚洲另类 | svipshipincom国产片| 欧美精品啪啪一区二区三区| 男女下面进入的视频免费午夜| 日本a在线网址| 色视频www国产| 国产伦精品一区二区三区四那| 久久精品91蜜桃| 久久99热这里只有精品18| 国产高清有码在线观看视频| 欧美一级毛片孕妇| 日本熟妇午夜| 99精品欧美一区二区三区四区| 午夜免费激情av| 久久人妻av系列| 老汉色∧v一级毛片| 国产伦人伦偷精品视频| 亚洲第一电影网av| 日韩欧美 国产精品| 天堂√8在线中文| 国产成人影院久久av| 国产在线精品亚洲第一网站| 国产精品久久久久久久电影 | АⅤ资源中文在线天堂| 别揉我奶头~嗯~啊~动态视频| 中文字幕久久专区| 午夜影院日韩av| 亚洲av日韩精品久久久久久密| 免费看十八禁软件| 99热这里只有是精品50| 久久久成人免费电影| 一边摸一边抽搐一进一小说| 看免费av毛片| 久久精品人妻少妇| 亚洲精品国产精品久久久不卡| 午夜激情福利司机影院| 窝窝影院91人妻| 男人舔女人的私密视频| 日韩三级视频一区二区三区| 桃色一区二区三区在线观看| 欧美不卡视频在线免费观看| 午夜a级毛片| 国产伦人伦偷精品视频| 久久午夜综合久久蜜桃| 日日夜夜操网爽| 日韩欧美国产一区二区入口| 岛国在线免费视频观看| 久久午夜综合久久蜜桃| 又紧又爽又黄一区二区| 亚洲电影在线观看av| 午夜福利在线在线| 精品久久久久久久毛片微露脸| 在线观看午夜福利视频| 日日摸夜夜添夜夜添小说| 97碰自拍视频| 岛国在线免费视频观看| 亚洲欧洲精品一区二区精品久久久| 精品不卡国产一区二区三区| 免费高清视频大片| xxxwww97欧美| 免费观看精品视频网站| 一个人看的www免费观看视频| 在线永久观看黄色视频| 久久午夜亚洲精品久久| 免费在线观看日本一区| 欧美激情久久久久久爽电影| 他把我摸到了高潮在线观看| 午夜福利视频1000在线观看| 久久久久久久久中文| 亚洲色图av天堂| 国产免费av片在线观看野外av| 中亚洲国语对白在线视频| 亚洲人成网站高清观看| 国内少妇人妻偷人精品xxx网站 | 久久中文字幕一级| 欧美极品一区二区三区四区| 白带黄色成豆腐渣| 亚洲精品粉嫩美女一区| 五月玫瑰六月丁香| 一个人观看的视频www高清免费观看 | 91老司机精品| 亚洲av五月六月丁香网| 精品人妻1区二区| 欧美黑人欧美精品刺激| 最近最新免费中文字幕在线| 欧美日韩一级在线毛片| 欧美一级毛片孕妇| 看免费av毛片| 国产精品亚洲美女久久久| 国产高清激情床上av| 桃红色精品国产亚洲av| 国产成人精品久久二区二区免费| 亚洲精品一卡2卡三卡4卡5卡| 国产成人系列免费观看| 中亚洲国语对白在线视频| 国产 一区 欧美 日韩| 成人18禁在线播放| av天堂中文字幕网| 黄色 视频免费看| 男女床上黄色一级片免费看| 在线观看午夜福利视频| 在线观看舔阴道视频| 波多野结衣高清作品| 黄色丝袜av网址大全| 精华霜和精华液先用哪个| 国产麻豆成人av免费视频| 午夜两性在线视频| 麻豆av在线久日| 天天一区二区日本电影三级| 久久中文字幕一级| 亚洲av成人一区二区三| 久久久久久久午夜电影| 可以在线观看毛片的网站| 1024香蕉在线观看| 亚洲国产精品久久男人天堂| 国产高清激情床上av| 亚洲国产色片| 国产在线精品亚洲第一网站| 国产av在哪里看| 美女高潮喷水抽搐中文字幕| x7x7x7水蜜桃| 亚洲av五月六月丁香网| 国产av不卡久久| 美女高潮的动态| 欧美最黄视频在线播放免费| www日本在线高清视频| 美女高潮喷水抽搐中文字幕| 免费大片18禁| 一本精品99久久精品77| 亚洲成人精品中文字幕电影| 真人一进一出gif抽搐免费| 国产精品亚洲美女久久久| 国产av不卡久久| 老司机午夜福利在线观看视频| 在线a可以看的网站| 亚洲av日韩精品久久久久久密| 婷婷丁香在线五月| 1000部很黄的大片| 欧美色视频一区免费| 在线十欧美十亚洲十日本专区| 免费电影在线观看免费观看| 色视频www国产| 国产亚洲av高清不卡| 久久久精品欧美日韩精品| 91麻豆av在线| 国产毛片a区久久久久| 国产精品久久久久久亚洲av鲁大| 亚洲欧美一区二区三区黑人| 91九色精品人成在线观看| 久久天堂一区二区三区四区| 在线观看日韩欧美| 成人国产综合亚洲| 久久久水蜜桃国产精品网| 法律面前人人平等表现在哪些方面| 我的老师免费观看完整版| 制服丝袜大香蕉在线| 99热这里只有是精品50| 国产精品香港三级国产av潘金莲| 国产成人影院久久av| 小蜜桃在线观看免费完整版高清| 18禁国产床啪视频网站| 99热这里只有精品一区 | 88av欧美| 麻豆国产av国片精品| 亚洲在线观看片| 日韩欧美精品v在线| 久久久久九九精品影院| 亚洲色图av天堂| 五月伊人婷婷丁香| 两性夫妻黄色片| 欧美日本视频| 精品国产美女av久久久久小说| 亚洲成人久久爱视频| 免费无遮挡裸体视频| 国产三级黄色录像| 蜜桃久久精品国产亚洲av| www.自偷自拍.com| 久久午夜亚洲精品久久| 亚洲av五月六月丁香网| 国产高清激情床上av| 日本三级黄在线观看| 婷婷精品国产亚洲av在线| 亚洲精品久久国产高清桃花| 搡老熟女国产l中国老女人| 国产私拍福利视频在线观看| 欧美黄色淫秽网站| 九九在线视频观看精品| 99热只有精品国产| 精品一区二区三区av网在线观看| 美女免费视频网站| 一进一出抽搐动态| 一二三四在线观看免费中文在| 嫩草影院精品99| 日本一本二区三区精品| 热99在线观看视频| 国产成人欧美在线观看| 国产99白浆流出| 国产精品女同一区二区软件 | 欧美中文综合在线视频| 欧美+亚洲+日韩+国产| 美女cb高潮喷水在线观看 | 一进一出好大好爽视频| 精品午夜福利视频在线观看一区| 女人高潮潮喷娇喘18禁视频| 嫁个100分男人电影在线观看| 中文字幕高清在线视频| 国产精品免费一区二区三区在线| 熟女少妇亚洲综合色aaa.| 国内少妇人妻偷人精品xxx网站 | 亚洲精品乱码久久久v下载方式 | 国产精品久久久人人做人人爽| 天天添夜夜摸| av片东京热男人的天堂| 午夜日韩欧美国产| 久久久久久人人人人人| 淫秽高清视频在线观看| 国产精华一区二区三区| svipshipincom国产片| 久久国产精品影院| 日韩精品青青久久久久久| 哪里可以看免费的av片| 亚洲国产精品成人综合色| aaaaa片日本免费| 亚洲av成人一区二区三| 免费看a级黄色片| 免费搜索国产男女视频| 国产亚洲精品av在线| 国产精品久久久av美女十八| 欧美丝袜亚洲另类 | 国产精品亚洲一级av第二区| 免费观看人在逋| 一区二区三区国产精品乱码| 国产精品免费一区二区三区在线| 国产成人啪精品午夜网站| 男人舔奶头视频| 日韩高清综合在线| 国产精品一区二区三区四区免费观看 | av欧美777| 国产视频内射| 男人舔女人的私密视频| 十八禁人妻一区二区| 手机成人av网站| 19禁男女啪啪无遮挡网站| 日韩欧美在线二视频| 午夜免费成人在线视频| 国产成人一区二区三区免费视频网站| 中文亚洲av片在线观看爽| 国产精品av久久久久免费| 99精品在免费线老司机午夜| 97超视频在线观看视频| 色尼玛亚洲综合影院| 国产真人三级小视频在线观看| 狂野欧美白嫩少妇大欣赏| 色老头精品视频在线观看| 九九热线精品视视频播放| 一本精品99久久精品77| 九九久久精品国产亚洲av麻豆 | 午夜视频精品福利| 欧美性猛交╳xxx乱大交人| 国产毛片a区久久久久| 俺也久久电影网| 18禁观看日本| 精品久久蜜臀av无| 日韩欧美三级三区| 亚洲av片天天在线观看| 99久久无色码亚洲精品果冻| 亚洲精品在线美女| 99热6这里只有精品| 美女大奶头视频| 日本在线视频免费播放| 19禁男女啪啪无遮挡网站| 欧美日韩国产亚洲二区| 亚洲国产欧美一区二区综合| 国产高潮美女av| 首页视频小说图片口味搜索| 91麻豆精品激情在线观看国产| 国内精品久久久久精免费| 一个人免费在线观看的高清视频| 欧美精品啪啪一区二区三区| 欧美乱妇无乱码| 精品人妻1区二区| 法律面前人人平等表现在哪些方面| 五月玫瑰六月丁香| 国产精品亚洲美女久久久| 国产精品久久视频播放| 大型黄色视频在线免费观看| 人妻久久中文字幕网| 久久久久久久久久黄片| 91九色精品人成在线观看| 夜夜爽天天搞| 久久久久久九九精品二区国产| 国产单亲对白刺激| 亚洲av美国av| 好看av亚洲va欧美ⅴa在| 午夜福利在线观看吧| 亚洲成人久久性| 免费看日本二区| 美女黄网站色视频| 亚洲av片天天在线观看| 精品久久久久久久末码| 又粗又爽又猛毛片免费看| 欧美另类亚洲清纯唯美| 日韩三级视频一区二区三区| 757午夜福利合集在线观看| 欧美成狂野欧美在线观看| 国产三级在线视频| 国产又黄又爽又无遮挡在线| 综合色av麻豆| h日本视频在线播放| 99精品久久久久人妻精品| 无遮挡黄片免费观看| 亚洲国产中文字幕在线视频| 长腿黑丝高跟| 久久久久国产精品人妻aⅴ院| 久久午夜亚洲精品久久| 国产成人系列免费观看| 夜夜躁狠狠躁天天躁| 极品教师在线免费播放| 一级黄色大片毛片| 老司机午夜十八禁免费视频| 免费av毛片视频| 真人做人爱边吃奶动态| 日本 欧美在线| ponron亚洲| 熟女少妇亚洲综合色aaa.| 亚洲精品美女久久av网站| 少妇熟女aⅴ在线视频| 真人一进一出gif抽搐免费| 国产精品久久久av美女十八| 亚洲欧美日韩高清在线视频| 91在线精品国自产拍蜜月 | 国产精品一区二区三区四区久久| 91在线观看av| 成人无遮挡网站| 嫩草影视91久久| www.精华液| 19禁男女啪啪无遮挡网站| 婷婷亚洲欧美| 一进一出抽搐动态| 国产蜜桃级精品一区二区三区| 中文字幕高清在线视频| 欧美黑人巨大hd| 男人舔奶头视频| 国产av一区在线观看免费| 亚洲av电影在线进入| 免费av不卡在线播放| bbb黄色大片| 国产蜜桃级精品一区二区三区| 久久久久久久久中文| 成人无遮挡网站| 最好的美女福利视频网| 亚洲人成网站高清观看| 亚洲七黄色美女视频| 国产一区二区三区视频了| 精品午夜福利视频在线观看一区| 国产成人福利小说| 国产精品1区2区在线观看.| 国产美女午夜福利| 色综合亚洲欧美另类图片| 9191精品国产免费久久| 丝袜人妻中文字幕| 无遮挡黄片免费观看| 亚洲欧美日韩高清在线视频| 欧美zozozo另类| 91字幕亚洲| 国内久久婷婷六月综合欲色啪| 久久久久久国产a免费观看| 亚洲欧美精品综合一区二区三区| 美女 人体艺术 gogo| 国产精品,欧美在线| 搞女人的毛片| 欧美日韩一级在线毛片| 99久久精品热视频| 欧美日韩国产亚洲二区| 午夜精品在线福利| 国产精品香港三级国产av潘金莲| 精品国内亚洲2022精品成人| 99久久无色码亚洲精品果冻| 人人妻,人人澡人人爽秒播| 中文在线观看免费www的网站| 制服丝袜大香蕉在线| 成人av一区二区三区在线看| 成人18禁在线播放| 精品国产乱码久久久久久男人| 精品国产超薄肉色丝袜足j| 国产成人一区二区三区免费视频网站| 搡老岳熟女国产| 麻豆一二三区av精品| 久久精品91无色码中文字幕| 亚洲真实伦在线观看| 黄色丝袜av网址大全| 日本免费a在线| ponron亚洲| 丰满人妻熟妇乱又伦精品不卡| 哪里可以看免费的av片| 婷婷精品国产亚洲av| 色视频www国产| 亚洲熟妇熟女久久| 日本熟妇午夜| 女警被强在线播放| 欧美中文综合在线视频| 亚洲精品粉嫩美女一区| 熟妇人妻久久中文字幕3abv| 18美女黄网站色大片免费观看| 精品不卡国产一区二区三区| 男人和女人高潮做爰伦理| 国产高潮美女av| 午夜精品久久久久久毛片777| 久久久久亚洲av毛片大全| 国产单亲对白刺激| 久久久成人免费电影| 国产aⅴ精品一区二区三区波| 十八禁网站免费在线| 91九色精品人成在线观看| 日日干狠狠操夜夜爽| 又紧又爽又黄一区二区| 在线观看日韩欧美| 欧美av亚洲av综合av国产av|