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

    截滲墻設(shè)計(jì)對(duì)水庫(kù)海岸帶預(yù)防海水入侵的影響

    2023-05-30 10:48:04馬軒季智靈孔俊
    河北工業(yè)科技 2023年2期
    關(guān)鍵詞:河口

    馬軒 季智靈 孔俊

    摘 要:為探究海岸水庫(kù)處于低水位時(shí),不同深度與位置工況的截滲墻對(duì)濱海潛水層海水入侵的防治效果,在物模實(shí)驗(yàn)驗(yàn)證的基礎(chǔ)上采用SUTRA模型建立了具有代表性的二維海岸模型,對(duì)有、無(wú)截滲墻、截滲墻底部到潛水層底部距離分別為14,12,10 m和截滲墻到海洋邊界分別為70,65,60 m工況下的潛水層中水鹽分布、水流運(yùn)動(dòng)特征、物質(zhì)輸移規(guī)律和地下水入侵動(dòng)力情況進(jìn)行模擬。結(jié)果表明:1)截滲墻能夠顯著減少海水入侵距離,且對(duì)于有海岸水庫(kù)的濱海潛水層的海水入侵防治效果更加明顯,減少幅度為19.5%;2)截滲墻越深入潛水層或越靠近海洋邊界,對(duì)海水入侵的防治效果和潛水層淡水儲(chǔ)量的增加效果越明顯,內(nèi)陸物質(zhì)輸移越慢,海側(cè)物質(zhì)輸移越快;3)截滲墻能有效阻止潮汐的波動(dòng)向內(nèi)陸的傳播,且截滲墻越靠近海洋、越深入潛水層,效果越明顯。研究進(jìn)一步明確了海岸水庫(kù)地區(qū)截滲墻作用下海水入侵的特征與規(guī)律,其成果可為優(yōu)化截滲墻空間設(shè)計(jì)提供模型參考和理論指導(dǎo),并為海岸水庫(kù)建設(shè)和科學(xué)管理提供科學(xué)依據(jù)。

    關(guān)鍵詞:河口、海岸學(xué);海岸水庫(kù);截滲墻;海水入侵;鹽水楔;物質(zhì)輸運(yùn)

    中圖分類號(hào):P731.2

    文獻(xiàn)標(biāo)識(shí)碼:A

    DOI: 10.7535/hbgykj.2023yx02001

    Influence of cutoff wall design on prevention of seawater intrusion of reservoir coastal zone

    MA Xuan, JI Zhiling, KONG Jun

    (College of Harbour, Coastal and Offshore Engineering, Hohai University, Nanjing, Jiangsu 210098, China)

    Abstract: In order to explore the prevention and control effect of cutoff walls with different depths and location conditions on seawater intrusion in the coastal diving layer when coastal reservoir is at low water level, a representative two-dimensional coastal model was established by SUTRA model on the basis of model experiment verification. The water and salt distribution, water flow movement characteristics, material transport rules and groundwater intrusion dynamics in the diving layer were simulated under the following conditions: whether there was the cutoff wall, the distance from the bottom of the cutoff wall to the bottom of the diving layer was 14, 12 and 10 m, and the distance from the bottom of the cutoff wall to the ocean boundary was 70, 65 and 60 m.The results show that: 1) The cutoff wall can significantly reduce the seawater intrusion distance, and the prevention effect of coastal diving layer with coastal reservoirs is more obvious, with a reduction of 19.5%; 2) The deeper the cutoff wall is in the diving layer or the closer it is to the marine boundary, the more obvious the prevention effect of seawater intrusion and the increase effect of fresh water reserves in the diving layer are; The slower the inland material transport is, the faster the marine material transport is; 3) The cutoff wall can effectively prevent the tidal wave from spreading inland, and the closer the cutoff wall is to the ocean, the deeper to the diving layer, the more obvious the effect is. This study further clarifies the characteristics and laws of seawater intrusion under the action of the cutoff wall in the coastal reservoir area, and the results can provide model reference and theoretical guidance for optimizing the spatial design of the cutoff wall, and provide a scientific basis for the construction and scientific management of coastal reservoir.

    Keywords:estuarine and coastal science;coastal reservoir;cutoff wall;seawater intrusion;salt wedge;material transport

    近年來(lái),在濱海地區(qū),經(jīng)濟(jì)的快速發(fā)展加劇了水資源的短缺,而濱海地區(qū)的海水入侵問(wèn)題常年廣泛存在。海水入侵會(huì)導(dǎo)致土壤鹽漬化、地下水水質(zhì)惡化等問(wèn)題,加劇了濱海地區(qū)的環(huán)境壓力,成為制約濱海地區(qū)生產(chǎn)生活和經(jīng)濟(jì)發(fā)展的重要因素1-2,也對(duì)水資源的開(kāi)發(fā)利用提出了更高的要求,防治海水入侵日趨成為研究熱點(diǎn)。常見(jiàn)的防治工程有淡水帷幕和地下壩[3,此外,在濱海地區(qū)大量興建海岸水庫(kù),利用堤壩抵擋咸水入侵,并在恰當(dāng)?shù)臅r(shí)間段蓄積淡水,提供大量?jī)?yōu)質(zhì)淡水資源4,改變?yōu)I海潛水層的水文環(huán)境,對(duì)濱海地區(qū)潛水層的海水入侵有著重大影響5。因此,研究海岸水庫(kù)工程和截滲墻對(duì)濱海潛水層海水入侵的影響規(guī)律,對(duì)掌握濱海地區(qū)的海水入侵狀況及制定防治對(duì)策具有重大意義。

    對(duì)海水入侵的大量研究主要聚焦于咸淡水界面形狀、運(yùn)移機(jī)理和規(guī)律,研究所用數(shù)學(xué)模型主要為突變界面模型6-8和過(guò)渡帶模型等概化模型[9-11。GLOVER[12基于突變界面理論在咸淡水界面位置方程的基礎(chǔ)上提出了Glover解,該方程內(nèi)陸邊界設(shè)置為恒定流量,很好地描述了淡水通過(guò)排泄通道排入到海水的過(guò)程。文獻(xiàn)[13]提到,Herry以承壓潛水層為穩(wěn)態(tài)條件,求解了潛水層的水流運(yùn)動(dòng)特征以及對(duì)應(yīng)位置的濃度分布特征。PINDER等14在Herry研究的基礎(chǔ)上,提出了基于數(shù)值模型的過(guò)渡帶模型,進(jìn)行了數(shù)值求解。HUYAKORN等15建立了三維有限元數(shù)值模型,研究了水流運(yùn)動(dòng)特征和溶質(zhì)運(yùn)移。VOSS等[16基于飽和、非飽和潛水層的地下水流運(yùn)動(dòng)特征和溶質(zhì)運(yùn)移特征,利用過(guò)渡帶模型建立了SUTRA

    (saturated-unsaturated transport)

    模型來(lái)進(jìn)行對(duì)海水入侵的數(shù)值模擬研究。ROBINSON等17在過(guò)渡帶模型中詳細(xì)闡述了潮汐波動(dòng)效應(yīng)下潛水層的水鹽變化。在海岸水庫(kù)和截滲墻對(duì)潛水層海水入侵方面,MO等18分析了海岸水庫(kù)底部深水池對(duì)潛水層的鹽分分布和海水入侵距離的影響,得出海岸水庫(kù)深水池加劇了海水入侵程度。MEHDIZADEH等

    [19研究了壩后蓄積的淡水對(duì)鹽水楔的形狀和推移距離的影響。ANWAR[20考慮了截滲墻對(duì)潛水層鹽水楔的作用,提出了鹽水楔位置方程。LUYUN等21在其研究中對(duì)該解析解進(jìn)行了科學(xué)驗(yàn)證,指出截滲墻越靠近海岸和潛水層底部,其對(duì)海水入侵程度的減弱作用越小。

    潮汐是濱海潛水層水流運(yùn)動(dòng)波動(dòng)和溶質(zhì)運(yùn)移的主要?jiǎng)恿?,?huì)加快海水與淡水的循環(huán)混合作用[22。目前有關(guān)海岸水庫(kù)和截滲墻對(duì)潛水層海水入侵綜合作用的研究很少,并且很多學(xué)者在研究中往往忽略潮汐動(dòng)力邊界。突變界面模型是基于潛水層水動(dòng)力靜態(tài)平衡提出的,而咸淡水的混合程度隨時(shí)間不斷發(fā)生變化,處于動(dòng)態(tài)平衡過(guò)程,當(dāng)過(guò)渡帶較寬時(shí),過(guò)渡帶模型更能有效地解決動(dòng)態(tài)平衡問(wèn)題,保證模擬結(jié)果的精度。因此,本文采用過(guò)渡帶模型SUTRA進(jìn)行模擬,建立具有代表性的二維海岸剖面模型,結(jié)合潮汐動(dòng)力邊界,設(shè)計(jì)不同的截滲墻深度和位置,對(duì)建有海岸水庫(kù)的濱海地區(qū)潛水層發(fā)生海水入侵后的水鹽運(yùn)移進(jìn)行數(shù)值模擬研究,以期為海岸水庫(kù)和截滲墻作用下潛水層的海水入侵防治和地下截滲墻的合理設(shè)置方式提供理論依據(jù)。

    1 地下水?dāng)?shù)值模型

    1.1 控制方程

    由美國(guó)地質(zhì)調(diào)查局(USGS)開(kāi)發(fā)出來(lái)的飽和帶/非飽和帶運(yùn)移(SUTRA)模型采用有限元和有限差分混合算法,在加權(quán)剩余方法體系中模擬地下水飽和及非飽和流動(dòng)、能量或溶質(zhì)在含水層多孔介質(zhì)中的運(yùn)移2個(gè)過(guò)程,來(lái)解決地下環(huán)境中水流溶質(zhì)和鹽度的飽和及不飽和運(yùn)輸問(wèn)題。當(dāng)對(duì)研究區(qū)域進(jìn)行適當(dāng)?shù)臅r(shí)空離散化后,這種方法的結(jié)果非常準(zhǔn)確。標(biāo)準(zhǔn)有限元逼近只適用于描述流體質(zhì)量、溶質(zhì)質(zhì)量和能量的平衡方程條件,而其他所有非通量條件都使用積分有限差分方法的有限元網(wǎng)格形式近似表示。SUTRA模型既保留了有限元方法的數(shù)學(xué)精確性和幾何靈活性,還具有有限差分方法的有效性。它既可以解決自由表面水體的運(yùn)動(dòng)問(wèn)題,也可以對(duì)特定濃度海水的入侵進(jìn)行模擬,還可以對(duì)復(fù)雜邊界下的溶質(zhì)傳輸問(wèn)題進(jìn)行模擬,例如在海水邊界加入潮汐信號(hào)。對(duì)于海水入侵問(wèn)題,一般采用地下水流動(dòng)方程和地下水溶質(zhì)傳輸方程16。SUTRA模型的二維有限元地下水流動(dòng)控制方程如式(1)所示,溶質(zhì)運(yùn)移方程如式(2)所示。

    式中:ρ為流體密度;φ為孔隙度;Sw為飽和度;t為時(shí)間;q→為達(dá)西流量;ρs為源流體密度;K(ψ)為滲透系數(shù);ψ為土壤負(fù)水頭;P為孔隙水壓力;g為重力加速度;z為水位高度。

    溶質(zhì)運(yùn)移方程又表示為式(3):

    式中:C為流體鹽度;D為水動(dòng)力彌散常量;C*為源流體鹽度;Q為源流量。

    根據(jù)van GENUCHTEN[23提出的經(jīng)驗(yàn)公式,非飽和土壤的水力參數(shù)由式(4)—式(6)確定。

    式中:α和n為形狀系數(shù);Ks為飽和系數(shù);Swres為孔隙水殘余飽和度。

    1.2 物理模型構(gòu)建及驗(yàn)證

    1.2.1 實(shí)驗(yàn)室小尺度物理模型

    為驗(yàn)證SUTRA模型的準(zhǔn)確性,建立了小尺度基本物理模型,如圖1所示。參考KUAN等[24制作的實(shí)驗(yàn)裝置如圖2所示,其中砂槽由3個(gè)部分組成,中間填有2 cm厚、0.70 m高的石英砂,左邊是模擬坡度為1∶3的灘面,長(zhǎng)度為1.35 m的濱海潛水層,右邊模擬的是長(zhǎng)度為2.65 m的內(nèi)陸潛水層。砂槽的右側(cè)為水位控制箱,控制內(nèi)陸水位為0.59 m。砂槽左側(cè)是潮波生成器,下面連接鹽水蓄水箱,無(wú)潮汐情況下,鹽水從下面的蓄水箱中通過(guò)泵輸送到潮波生成器中,并通過(guò)溢流平臺(tái)排出;有潮汐情況下通過(guò)自動(dòng)調(diào)節(jié)水位的溢流平臺(tái)來(lái)生成預(yù)設(shè)的波動(dòng)水位。本次實(shí)驗(yàn)?zāi)M采用正弦信號(hào)波動(dòng)水位,如式(7)所示。鹽水濃度設(shè)定為3.5%(質(zhì)量分?jǐn)?shù),下同),密度為1 025 kg/m3。為了將海水入侵過(guò)程直觀化、可視化,將34 g的NaCl和1 g紅色食用色素溶于1 L淡水中來(lái)制備鹽水蓄水箱的鹽水,淡水密度為1 000 kg/m3,淡水鹽度為0。本次實(shí)驗(yàn)主要模擬有潮汐情況,右側(cè)為波動(dòng)水位,振幅為0.038 m,周期為62 s。

    式中:h(t)為t時(shí)刻的壓力水頭,即潮汐水位;A為振幅; MSL為平均水面;T為波動(dòng)周期;t為經(jīng)過(guò)時(shí)間。

    1.2.2 邊界條件及模型參數(shù)

    為了進(jìn)一步分析濱海潛水層的鹽度分布特征和水流運(yùn)動(dòng)規(guī)律,進(jìn)行了基于實(shí)驗(yàn)小尺度的數(shù)值模擬,運(yùn)用SUTRA數(shù)值模型來(lái)模擬有潮汐情況下潛水層的鹽水入侵。

    模型的邊界條件如圖3所示,本模型底部邊界為不透水邊界,沒(méi)有水流和溶質(zhì)交換;泥沙區(qū)水平頂部設(shè)為無(wú)水流和溶質(zhì)交換邊界;右側(cè)為地下水流邊界,地下水位為0.59 m,設(shè)其流動(dòng)為均勻流動(dòng),溶質(zhì)濃度為0;AB,AC和CD段為海岸邊界,其溶液濃度設(shè)為3.5%,其壓力水頭隨時(shí)間的變化規(guī)律由式(7)表示。

    在前處理軟件Argus One中選取二維常規(guī)網(wǎng)格格式進(jìn)行網(wǎng)格劃分,并在運(yùn)行之前進(jìn)行多次調(diào)試,設(shè)置合適的時(shí)長(zhǎng)與步數(shù),使其在正式的模擬中可以很快地實(shí)現(xiàn)收斂,并獲得足夠的數(shù)據(jù)用于后處理。模型網(wǎng)格劃分遵循式(8)16。

    式中:Pem為佩克萊數(shù);ΔLL為沿著流線方向的網(wǎng)格單元長(zhǎng)度;αL為縱向彌散度。

    潛水層被劃分為800×80的64 000個(gè)網(wǎng)格,對(duì)應(yīng)節(jié)點(diǎn)數(shù)64 881個(gè)。此模型為二維模型,并忽略了海岸帶橫向鹽分交換過(guò)程,假定多孔介質(zhì)呈現(xiàn)各向同性,依照van GENUCHTEN方程[23]對(duì)相關(guān)參數(shù)進(jìn)行設(shè)定。同時(shí),在前處理軟件Argus One中對(duì)各相關(guān)參數(shù)進(jìn)行校正,最后選取的參數(shù)如表1所示。另外,為了進(jìn)行網(wǎng)格無(wú)關(guān)性驗(yàn)證,將網(wǎng)格進(jìn)行加密處理,除了采用800×80的網(wǎng)格,還分別劃分出1 600×160和1 200×120的網(wǎng)格進(jìn)行計(jì)算驗(yàn)證。

    1.2.3 驗(yàn)證

    小尺度數(shù)學(xué)模型模擬結(jié)果如圖4所示,物理模型實(shí)驗(yàn)結(jié)果如圖5所示。可以看出,數(shù)值模擬的結(jié)果和實(shí)驗(yàn)結(jié)果基本吻合,鹽水入侵距離的數(shù)值模擬結(jié)果為2.07 m,實(shí)驗(yàn)結(jié)果為2.05 m,僅相差了2 cm,鹽水楔的形態(tài)和上部混合區(qū)的形態(tài)也基本一致。另外采取的2種加密網(wǎng)格計(jì)算出來(lái)的海水入侵距離分別為2.065 m和2.068 m,誤差在1%以內(nèi),驗(yàn)證了網(wǎng)格的敏感性,符合網(wǎng)格無(wú)關(guān)性要求。數(shù)值模擬和物理模型對(duì)比結(jié)果表明SUTRA數(shù)值模擬的結(jié)果是可靠的,可對(duì)海岸水庫(kù)作用下濱海潛水層的海水入侵進(jìn)行研究,包括對(duì)其鹽度分布和水流特征進(jìn)行定量和定性的分析。

    1.3 大尺度概念模型

    為研究截滲墻對(duì)海水入侵造成的影響,本文構(gòu)建了一個(gè)概化的海灘剖面進(jìn)行研究。如圖6所示,在低潮位處設(shè)立壩體,壩后蓄淡水,形成海岸水庫(kù)。建立模型時(shí)將平均海平面與海岸的交點(diǎn)設(shè)為O(150,20),在F(155,19)處設(shè)置壩體,壩高4 m,壩厚1 m,AB=30 m,BE=200 m,AC=100 m,DE=10 m。潮汐振幅A=1 m,周期T=12 h。GH為設(shè)置的截滲墻,截滲墻的寬度為1 m,其與海邊界的距離為XL,與潛水層底部邊界的距離為ZL,截滲墻整體設(shè)置為低滲透區(qū),滲透系數(shù)為0.5 m/d,是潛水層土壤滲透率的1/20,庫(kù)區(qū)水位取海水入侵嚴(yán)重時(shí)的低水位情況,為0.5 m。

    1.3.1 網(wǎng)格及參數(shù)設(shè)置

    在SUTRA模型中將計(jì)算模型剖分為了24 000個(gè)單元,24 461個(gè)節(jié)點(diǎn),生成四邊形網(wǎng)格,網(wǎng)格最小尺度為0.5 m,在模型中假定為各向同性介質(zhì),計(jì)算所采用的具體參數(shù)參照了XIN等25在模擬中使用的參數(shù),見(jiàn)表2。

    1.3.2 邊界條件設(shè)置

    如圖6所示,潛水層水平頂部AC段設(shè)為無(wú)水流和溶質(zhì)交換邊界,溶質(zhì)濃度為0。模型左側(cè)邊界為內(nèi)陸邊界,設(shè)為定流量邊界,無(wú)潮汐狀態(tài)下,給定內(nèi)陸水頭為21 m時(shí),達(dá)到穩(wěn)態(tài)對(duì)應(yīng)的單寬內(nèi)陸邊界流量為557 kg/d,將該值設(shè)為流量邊界流量值。右側(cè)邊界為海洋邊界,底部邊界BE設(shè)為不透水邊界,無(wú)水流和溶質(zhì)交換,模型灘面位置的表面邊界考慮3種工況:

    1)無(wú)水庫(kù)、有截滲墻和無(wú)截滲墻。CDE段設(shè)為潮汐動(dòng)力邊界,壓力水頭隨時(shí)間變化規(guī)律如式(9)所示,鹽水濃度為3.5%。

    式中:YMSL為平均海水位20 m;AT為潮汐振幅1 m;T為潮汐周期12 h;t為當(dāng)前時(shí)刻值。

    2)有水庫(kù)、有截滲墻和無(wú)截滲墻。壩體為不透水邊界,無(wú)水流和溶質(zhì)交換,F(xiàn)DE為潮汐動(dòng)力邊界,潮汐水位如式(9)設(shè)置,鹽水濃度為3.5%;水庫(kù)水位以下FG段為水庫(kù)水位邊界,設(shè)置為靜水壓強(qiáng),其變化規(guī)律如式(10)所示,水庫(kù)水位以上CG段邊界條件在模型中進(jìn)行動(dòng)態(tài)調(diào)整。當(dāng)上一時(shí)步對(duì)應(yīng)的節(jié)點(diǎn)為飽和節(jié)點(diǎn),則此節(jié)點(diǎn)設(shè)置為滲透面,壓強(qiáng)為0;當(dāng)上一時(shí)步對(duì)應(yīng)的節(jié)點(diǎn)為不飽和節(jié)點(diǎn),則此節(jié)點(diǎn)設(shè)置為無(wú)水流邊界。CF段考慮庫(kù)區(qū)水位為0.5 m,庫(kù)區(qū)溶質(zhì)濃度邊界設(shè)置為入流濃度為0,出流濃度梯度為0。

    Pr=ρg(hr+19-y) ,(10)

    式中:Pr為庫(kù)區(qū)邊界靜水壓強(qiáng);hr為庫(kù)區(qū)水深。

    3)有水庫(kù)、有截滲墻。設(shè)置3種截滲墻深度(ZL=14,12,10 m)和3種截滲墻位置(XL=70,65,60 m),其他同工況2)。

    通過(guò)SUTRA運(yùn)行計(jì)算模型,得出每個(gè)節(jié)點(diǎn)的鹽度值、流速,繪制鹽度流場(chǎng)圖,計(jì)算整個(gè)模型區(qū)域的鹽量,結(jié)合濱海潛水層的水流流場(chǎng)、海水入侵距離、單寬鹽量、單寬淡水儲(chǔ)量、鹽分質(zhì)心位置、粒子運(yùn)移時(shí)間與粒子運(yùn)移軌跡等參數(shù),分析計(jì)算結(jié)果,總結(jié)規(guī)律。以此分析截滲墻及其不同深度、位置的設(shè)計(jì)對(duì)有海岸水庫(kù)的濱海潛水層防治海水入侵的效果進(jìn)一步探究截滲墻對(duì)地下水入侵的動(dòng)力影響機(jī)制。

    2 結(jié)果與分析

    2.1 截滲墻對(duì)潛水層水鹽運(yùn)移的影響

    2.1.1 無(wú)海岸水庫(kù)工況下的影響

    首先對(duì)比有、無(wú)截滲墻時(shí)濱海地區(qū)潛水層的咸淡水分布情況,以無(wú)潮穩(wěn)定狀態(tài)作為初始狀態(tài),模擬出有潮汐作用下地下含水層的時(shí)均穩(wěn)態(tài)鹽度分布。模型中內(nèi)陸水位邊界為21 m,以周期平均的方法繪制一個(gè)周期內(nèi)含水層的鹽度分布,并采用拉格朗日質(zhì)點(diǎn)法研究含水層物質(zhì)輸移特征,在灘面和內(nèi)陸邊界分別釋放一定數(shù)量的粒子,繪制其運(yùn)動(dòng)軌跡,結(jié)果如圖7所示,白線上的數(shù)字表示粒子從釋放點(diǎn)到排出點(diǎn)的歷時(shí)。圖中紅色代表海水,黑線代表50%海水鹽度等值線(17.5)。將咸淡水界面(50%海水鹽度等值線,即17.5)與底部邊界的交點(diǎn)稱為海水入侵點(diǎn),其與海向邊界的水平距離稱為海水入侵距離。

    對(duì)比圖7中的a)和b)可以看出,修建截滲墻后鹽水楔產(chǎn)生了部分后退,其中海水入侵距離從52.5 m減少到50.5 m,減少了2 m,減少幅度為3.8%,單寬鹽量從5.025×104 kg/m減少到4.980×104 kg/m,減小了450 kg/m,減少幅度為9.0‰。這是因?yàn)榻貪B墻擋住了潛水層上部水流,淡水排泄通道在截滲墻的作用下變窄,增加了向海側(cè)的水力梯度,逼退鹽水楔,從而減少了海水入侵。觀察水流特征,可以看出內(nèi)陸釋放的粒子先沿水平方向運(yùn)動(dòng),經(jīng)過(guò)截滲墻時(shí)在截滲墻的擋流作用下,粒子向下運(yùn)動(dòng),再沿著咸淡水界面上升至灘面排出潛水層。因?yàn)樾纬闪溯^窄的淡水排泄通道,此處流速加快,使得內(nèi)陸粒子的運(yùn)移時(shí)間變少。如從模型內(nèi)坐標(biāo)為(75,10)釋放的粒子在修建截滲墻后,運(yùn)移時(shí)間從635.4 d減少到604.9 d,減少了30.5 d,減少幅度為4.8%。而海側(cè)粒子在鹽水楔中先沿著潛水層向下運(yùn)動(dòng),到達(dá)咸淡水界面后沿著過(guò)渡帶上升至灘面排出,由于截滲墻對(duì)鹽水楔的逼退作用,減少了海側(cè)粒子的運(yùn)移路徑,因此海側(cè)粒子的運(yùn)移時(shí)間也相應(yīng)減少了。如從(190,12)釋放的粒子在截滲墻的作用下運(yùn)移時(shí)間減少了314.6 d,減少幅度為8.7%。

    2.1.2 有海岸水庫(kù)工況下的影響

    為了研究截滲墻對(duì)海岸水庫(kù)作用下濱海潛水層水鹽運(yùn)移的影響,在有潮建庫(kù)穩(wěn)定工況的基礎(chǔ)上修建截滲墻。其中,截滲墻的寬度為1 m,其與海邊界的距離XL為60 m,與潛水層底部邊界的距離ZL為10 m,庫(kù)區(qū)水位為0.5 m,內(nèi)陸邊界為定流量557 kg/d。將有潮建庫(kù)穩(wěn)定狀態(tài)作為初始狀態(tài),模型模擬出的鹽度分布和粒子路徑如圖8所示,其中a)和b)分別為無(wú)截滲墻和有截滲墻的工況。

    對(duì)比圖8中的a)和b)可以看出,在截滲墻的作用下,鹽水楔的形態(tài)發(fā)生了很大改變。鹽水楔上邊緣不變,下邊緣(即楔腳處)后退,反映了海水入侵距離、單寬鹽量的減少。模型中增設(shè)截滲墻后,海水入侵距離減少了16 m,減少幅度為19.5%,單寬鹽量減少了5.181×103kg/m,減少幅度為8.0%;可以看出,截滲墻對(duì)海水入侵的距離影響較大,對(duì)潛水層整體鹽量的影響相對(duì)較小。其原因是截滲墻的截滲作用使得潛水層下部的水流向下運(yùn)動(dòng),逼退鹽水楔下部,而靠近過(guò)渡帶的上方內(nèi)陸淡水水平流由于截滲墻的存在而減弱,上部淡水對(duì)鹽水楔上部的逼退作用減弱,所以在海水密度流的推動(dòng)下,鹽水楔上部向內(nèi)陸移動(dòng)。下方退咸,上方進(jìn)咸,雙重作用下,模型的整體鹽量雖然在減少,但幅度較海水入侵距離的變化幅度更小。通過(guò)計(jì)算模型得出,在截滲墻的作用下,單寬淡水儲(chǔ)量從4 038 m3/m增加到4 246 m3/m,增加幅度為5.1%;可以看出,截滲墻可增加海岸水庫(kù)濱海潛水層的淡水儲(chǔ)量,可用來(lái)緩解甚至解決該區(qū)域淡水資源短缺的問(wèn)題,為濱海地區(qū)的工農(nóng)業(yè)用水及人類生活用水提供了解決方案。

    進(jìn)一步觀察圖8所示的水流特征可以看出,截滲墻阻礙了內(nèi)陸水流運(yùn)動(dòng)。較無(wú)截滲墻情況,內(nèi)陸粒子的運(yùn)移路徑延長(zhǎng),導(dǎo)致運(yùn)移時(shí)間增加。例如:添加截滲墻前后,從模型內(nèi)坐標(biāo)為(75,10)處釋放的粒子的運(yùn)移時(shí)間從495.1 d增加到560.8 d,增加了65.7 d,增加幅度為13.3%;海側(cè)粒子在鹽水楔中先沿著潛水層向下運(yùn)動(dòng),到達(dá)咸淡水界面后沿著過(guò)渡帶上升至灘面排出,由于截滲墻對(duì)鹽水楔的逼退作用,海側(cè)粒子的運(yùn)移路徑減少,粒子的運(yùn)移時(shí)間也相應(yīng)減少。例如:添加截滲墻前后,從模型內(nèi)坐標(biāo)為(190,12)處釋放的粒子的運(yùn)移時(shí)間從675.7 d減少到575.1 d,減少了100.6 d,減少幅度為14.9%。說(shuō)明截滲墻減慢了海岸水庫(kù)作用下的濱海潛水層的內(nèi)陸物質(zhì)輸移,加快了海側(cè)物質(zhì)輸移,且這一效果比沒(méi)有海岸水庫(kù)存在時(shí)更加明顯。

    2.2 截滲墻不同設(shè)計(jì)參數(shù)對(duì)海岸水庫(kù)作用下潛水層水鹽運(yùn)移的影響

    2.2.1 截滲墻深度的影響

    保持截滲墻距海邊界距離XL為60 m不變,改變截滲墻底部到潛水層底部的距離ZL,分別設(shè)定為14,12,10 m,海岸水庫(kù)的庫(kù)區(qū)水位為0.5 m。將有潮建庫(kù)穩(wěn)定狀態(tài)作為初始狀態(tài),模擬出3種不同深度截滲墻作用下潛水層內(nèi)的鹽度分布和粒子路徑,如圖9所示。計(jì)算出3種不同截滲墻深度下的海水入侵距離、單寬鹽量、單寬淡水儲(chǔ)量以及質(zhì)心位置,如表3所示。

    從表3中海水入侵距離和單寬鹽量這2個(gè)衡量海水入侵程度的指標(biāo)可以看出,隨著截滲墻不斷深入潛水層,其在海岸水庫(kù)作用下的濱海潛水層中對(duì)海水入侵的防治效果愈加明顯。分析鹽分在x,z方向上質(zhì)心的變化容易得出,隨著截滲墻不斷深入潛水層,鹽水楔質(zhì)心在水平方向上向海側(cè)移動(dòng),在垂直方向向上移動(dòng)。從截滲墻深度對(duì)海岸水庫(kù)作用下濱海潛水層單寬淡水儲(chǔ)量的影響方面可以看出,隨著截滲墻的不斷深入,潛水層的淡水儲(chǔ)量逐漸增加,且越深入,淡水儲(chǔ)量的增加幅度越明顯。進(jìn)一步觀察圖9中的水流特征可以看出,隨著截滲墻不斷深入潛水層,內(nèi)陸釋放的粒子向下的運(yùn)動(dòng)路徑增加,使得運(yùn)移時(shí)間增加。如從(75,10)釋放粒子由ZL=14 m到ZL=10 m的運(yùn)移時(shí)間增加了46.7 d,增加幅度為9.0%;在海側(cè),截滲墻的設(shè)置對(duì)鹽水楔存在逼退作用,且截滲墻越靠近潛水層底部,逼退效果越明顯;截滲墻的存在減少了海側(cè)粒子的運(yùn)移路徑,使得運(yùn)移時(shí)間減少。如從(190,12)釋放的粒子由ZL=14 m到ZL=10 m,運(yùn)移時(shí)間從1 574.7 d減少到1 360.1 d,減少幅度為15.3%。即截滲墻設(shè)置越深,阻礙水流運(yùn)動(dòng)作用越強(qiáng),內(nèi)陸物質(zhì)輸移越慢,海側(cè)物質(zhì)輸移越快。

    2.2.2 截滲墻位置的影響

    保持截滲墻距潛水層底部的距離ZL為10 m不變,改變截滲墻底部到海洋邊界的距離XL,分別設(shè)定為70,65,60 m,庫(kù)區(qū)水位取0.5 m。將有潮建庫(kù)穩(wěn)定狀態(tài)作為初始狀態(tài),計(jì)算出修建截滲墻后潛水層的鹽度分布和粒子路徑,如圖10所示。同時(shí)計(jì)算出3種工況下的海水入侵距離、單寬鹽量、單寬淡水儲(chǔ)量以及質(zhì)心位置,如表4所示。

    由表4可知,隨著截滲墻不斷向海側(cè)移動(dòng),其在海岸水庫(kù)作用下的濱海潛水層中防治海水入侵的效果愈加明顯。綜合圖10和表4可以看出,隨著截滲墻不斷向海側(cè)移動(dòng),鹽水楔形狀發(fā)生改變,鹽水楔質(zhì)心的位置在水平方向上向海側(cè)移動(dòng),在垂直方向上向上移動(dòng),截滲墻的位置對(duì)鹽水楔的質(zhì)心在x軸上移動(dòng)距離影響更大。分析截滲墻的位置對(duì)海岸水庫(kù)作用下濱海潛水層的單寬淡水儲(chǔ)量的影響容易得出,隨著截滲墻不斷向海側(cè)移動(dòng),相應(yīng)的潛水層淡水儲(chǔ)量增加,且截滲墻越靠近海側(cè),其對(duì)淡水儲(chǔ)量的增加效果越明顯。

    觀察圖10所示水流特征,截滲墻位置的不斷向海側(cè)移動(dòng),讓內(nèi)陸釋放的粒子向下運(yùn)動(dòng)路徑增加,粒子運(yùn)移時(shí)間增加,如從(75,10)釋放粒子的運(yùn)移時(shí)間由XL=70 m時(shí)的506.0 d,增加到XL=60 m時(shí)的563.4 d,增加幅度為11.3%。而在海側(cè),截滲墻的設(shè)置對(duì)鹽水楔起逼退作用,且截滲墻越靠近海洋邊界,逼退效果越明顯,減少了海側(cè)粒子的運(yùn)移路徑和運(yùn)移時(shí)間。即截滲墻越靠近海洋邊界,內(nèi)陸物質(zhì)輸移越慢,海側(cè)物質(zhì)輸移越快。

    2.3 截滲墻設(shè)計(jì)對(duì)潛水層地下水入侵的動(dòng)力影響

    2.3.1 截滲墻深度的動(dòng)力影響

    前面主要討論了有潮汐條件下設(shè)立截滲墻并達(dá)到穩(wěn)態(tài)時(shí),周期平均情況下濱海潛水層的水鹽特性。為了更好地研究截滲墻作用下潛水層的水流運(yùn)動(dòng)機(jī)理,在有潮汐作用、海水入侵達(dá)到穩(wěn)態(tài)的基礎(chǔ)上繼續(xù)分析周期過(guò)程中地下水位的變化過(guò)程,以明確不同截滲墻設(shè)計(jì)參數(shù)對(duì)潛水層中地下水波動(dòng)特性的影響。

    保持截滲墻距海邊界的距離XL為60 m不變,改變截滲墻底部到潛水層底部的距離ZL,分別設(shè)定為14,12,10 m,水庫(kù)水位取0.5 m。選取2處位置,分別是截滲墻靠海一側(cè)距離海岸邊界55 m,截滲墻靠?jī)?nèi)陸一側(cè)距離海岸邊界75 m。計(jì)算穩(wěn)定后一個(gè)周期(12 h)內(nèi)2個(gè)位置的水位變化,來(lái)研究不同截滲墻深度下2處位置對(duì)應(yīng)的地下水動(dòng)力特性。不同截滲墻深度下2處位置的地下水位隨時(shí)間變化情況如圖11所示。

    圖11 a)給出了截滲墻靠海一側(cè)距海岸邊界55 m處的水位變化特征,可以看出截滲墻靠海一側(cè)的地下水位隨著截滲墻深度的減小而有所抬高,水位波動(dòng)有所減小。當(dāng)截滲墻底部到潛水層底部的距離ZL從10 m增加到14 m,地下平均水位從19.87 m增加到19.88 m,增加了0.01 m,增加幅度為0.5‰,同時(shí)水位的波動(dòng)振幅從0.091 m降低到0.084 m,減少了0.007 m,減小幅度為7.7%。圖11 b)給出了在截滲墻靠?jī)?nèi)陸一側(cè)距海岸邊界75 m處的水位變化特征,可以看出隨著截滲墻深度的減小,地下水位在降低,水位波動(dòng)振幅在增加。截滲墻的深度降低了4 m,地下平均水位從20.32 m降低到20.25 m,降低了0.07 m,減小幅度為3.4‰,振幅從0.011 m變?yōu)?.013 m,增加了0.002 m,增加幅度為18.2%。這反映了截滲墻能夠阻礙潮汐向內(nèi)陸傳播,隨著截滲墻深度的減小,潮汐在地下水中引發(fā)的波動(dòng)效應(yīng)變得明顯。

    2.3.2 截滲墻位置的動(dòng)力影響

    保持截滲墻底部到潛水層底部的距離ZL為10 m不變,改變截滲墻距海邊界的距離XL,分別設(shè)定為60,65,70 m,水庫(kù)水位取0.5 m,選取2處位置,分別是截滲墻靠海一側(cè)距離海岸邊界55 m,截滲墻靠?jī)?nèi)陸一側(cè)距離海岸邊界75 m。計(jì)算穩(wěn)定后一個(gè)周期(12 h)內(nèi)2處位置的水位變化。不同截滲墻位置下2處位置的地下水位隨時(shí)間變化情況如圖12所示。

    截滲墻靠海一側(cè)距海岸邊界55 m處的地下水位,隨著截滲墻遠(yuǎn)離海洋逐漸抬高,水位波動(dòng)振幅有所降低。如截滲墻向內(nèi)陸移動(dòng)了10 m,地下平均水位從19.87 m增加到19.89 m,增加了0.02 m,增加幅度為1‰,同時(shí)水位波動(dòng)振幅從0.091 m降低到0.075 m,降低了0.016 m,減少幅度為17.6%。截滲墻靠?jī)?nèi)陸一側(cè)距海岸邊界75 m處的地下水位隨著截滲墻遠(yuǎn)離海洋逐漸降低,水位波動(dòng)振幅有所降低。截滲墻位置向內(nèi)陸移動(dòng)了10 m,地下平均水位從20.32 m降低到20.30 m,降低了0.02 m,減少幅度為1.0‰,振幅從0.011 m到0.009 m,降低了0.002 m,減少幅度為18.2%,這反映了截滲墻能夠阻礙潮汐向內(nèi)陸傳播,隨著截滲墻距海岸邊界距離的減小,潮汐在地下水中引發(fā)的波動(dòng)效應(yīng)變得明顯。另外從圖7可以看出,在截滲墻距離海岸邊界為65 m和70 m時(shí),2個(gè)計(jì)算點(diǎn)水位較接近,這是因?yàn)榻貪B墻設(shè)立在65 m和70 m處時(shí),其對(duì)潮汐波動(dòng)效應(yīng)的阻礙作用較低,此時(shí)潛水層的地下水位主要受潮汐控制。

    3 結(jié) 語(yǔ)

    本文采用SUTRA模型針對(duì)濱海水庫(kù)區(qū)域,構(gòu)建了帶有截滲墻濱海潛水層的數(shù)值模型,研究了截滲墻的深度及其位置對(duì)海岸水庫(kù)作用下的濱海潛水層的水鹽分布和水流運(yùn)動(dòng)特征的影響,結(jié)論如下。

    1)截滲墻能擋住潛水層上部水流,使得淡水排泄通道變窄,并增加向海側(cè)的水力梯度,逼退鹽水楔,從而減小海水入侵的距離;在具有海岸水庫(kù)的濱海潛水層中修建截滲墻,能夠顯著減少潛水層中的海水入侵程度,增加淡水儲(chǔ)量,減慢內(nèi)陸物質(zhì)向海側(cè)輸移速度,加快海側(cè)物質(zhì)向內(nèi)陸輸移速度。

    2)截滲墻的深度和位置設(shè)計(jì)參數(shù)顯著影響建有海岸水庫(kù)的濱海潛水層中的海水入侵程度,截滲墻越深入潛水層或越靠近海洋邊界,內(nèi)陸物質(zhì)輸移越慢,海側(cè)物質(zhì)輸移越快,對(duì)海水入侵的防治效果和增加潛水層中淡水儲(chǔ)量的效果越明顯。

    3)截滲墻阻止潮汐波動(dòng)向內(nèi)陸傳播的效果明顯,且截滲墻位置越靠近海洋、深度越深入潛水層,效果越明顯。為更好地防治海水入侵,工程應(yīng)用中可根據(jù)實(shí)際綜合水文地質(zhì)情況將截滲墻設(shè)置得盡量靠近海洋邊界和深入潛水層。截滲墻建設(shè)也將隨著新材料和新結(jié)構(gòu)(例如基于廢渣粉末的再生水泥基材料[26]和鋼-混凝土波形鋼腹板組合梁[27] 等)的不斷發(fā)展而降低成本成為防治海水入侵經(jīng)濟(jì)有效的方式。

    本文使用的研究模型未考慮降雨、蒸發(fā)、河流補(bǔ)給等源匯項(xiàng)的作用,且潮汐邊界為單一組分潮汐,在實(shí)際濱海地區(qū),潮汐多為復(fù)合潮汐與多組分潮汐;此外,本文采用的大尺度模型缺少實(shí)測(cè)數(shù)據(jù)的驗(yàn)證,在后續(xù)研究中需要進(jìn)一步對(duì)模型進(jìn)行調(diào)整與完善,結(jié)合實(shí)際工程進(jìn)行深化研究。

    參考文獻(xiàn)/References:

    [1] 李冬順,楊勁松,姚榮江.生態(tài)風(fēng)險(xiǎn)分析用于蘇北灘涂土壤鹽漬化風(fēng)險(xiǎn)評(píng)估研究[J].土壤學(xué)報(bào),2010,47(5):857-864.

    LI Dongshun,YANG Jinsong,YAO Rongjiang.Application of ecological risk analysis to soil salinization risk assessment of coastal tidal flat in North Jiangsu Province[J].Acta Pedologica Sinica,2010,47(5):857-864.

    [2] 姚榮江,楊勁松,陳小兵,等.蘇北海涂典型圍墾區(qū)土壤鹽漬化風(fēng)險(xiǎn)評(píng)估研究[J].中國(guó)生態(tài)農(nóng)業(yè)學(xué)報(bào),2010,18(5):1000-1006.

    YAO Rongjiang,YANG Jinsong,CHEN Xiaobing,et al.Evaluating soil salinization risk in typical coastal reclaimed regions in North Jiangsu Province[J].Chinese Journal of Eco-agriculture,2010,18(5):1000-1006.

    [3] 王輝,許學(xué)工.淡水帷幕防治海水入侵的國(guó)際經(jīng)驗(yàn)[J].中國(guó)農(nóng)村水利水電,2017(4):62-66.

    WANG Hui,XU Xuegong.International experience on hydraulic barrier serving as a seawater intrusion countermeasure[J].China Rural Water and Hydropower,2017(4):62-66.

    [4] GU Hanbin,GUO Qianli,LIN Pengzhi,et al.Feasibility study of coastal reservoirs in the Zhoushan Islands,China[J].Journal of Coastal Research,2019,35(4):835-841.

    [5] JIN Guangqiu,MO Yuminig,LI Mengdi,et al.Desalinization and salinization:A review of major challenges for coastal reservoirs[J].Journal of Coastal Research,2019,35(3):664-672.

    [6] JAKOVOVIC D,WERNER A D,DE LOUW P G B,et al.Saltwater upconing zone of influence[J].Advances in Water Resources,2016,94:75-86.

    [7] KARATZAS G P,DOKOU Z.Optimal management of saltwater intrusion in the coastal aquifer of Malia,Crete (Greece),using particle swarm optimization[J].Hydrogeology Journal,2015,23(6):1181-1194.

    [8] PARK N,SHI Lei.A comprehensive sharp-interface simulation-optimization model for fresh and saline groundwater management in coastal areas[J].Hydrogeology Journal,2015,23(6):1195-1204.

    [9] COBANER M,YURTAL R,DOGAN A,et al.Three dimensional simulation of seawater intrusion in coastal aquifers:A case study in the Goksu Deltaic Plain[J].Journal of Hydrology,2012,464/465:262-280.

    [10]XIAO Han,WANG Dingbao,HAGEN S C,et al.Assessing the impacts of sea-level rise and precipitation change on the surficial aquifer in the low-lying coastal alluvial plains and barrier islands,east-central Florida (USA)[J].Hydrogeology Journal,2016,24(7):1791-1806.

    [11]薛顯武,陳喜,魏玲娜,等.基于GIS的海水入侵可視化數(shù)值模擬系統(tǒng)[J].地下水,2007,29(1):35-39.

    XUE Xianwu,CHEN Xi,WEI Lingna,et al.Visualized numeric simulation system of seawater intrusion based on GIS[J].Ground Water,2007,29(1):35-39.

    [12]GLOVER R E.The pattern of fresh-water flow in a coastal aquifer[J].Journal of Geophysical Research,1959,64(4):457-459.

    [13]丁玲,李碧英,張樹(shù)深.海岸帶海水入侵的研究進(jìn)展[J].海洋通報(bào),2004,23(2):82-87.

    DING Ling,LI Biying,ZHANG Shushen.Development in researches on seawater intrusion of coastal areas[J].Marine Science Bulletin,2004,23(2):82-87.

    [14]PINDER G F,COOPER H H,Jr.A numerical technique for calculating the transient position of the saltwater front[J].Water Resources Research,1970,6(3):875-882.

    [15]HUYAKORN P S,ANDERSEN P F,MERCER J W,et al.Saltwater intrusion in aquifers:Development and testing of a three-dimensional finite element model[J].Water Resources Research,1987,23(2):293-312.

    [16]VOSS C I,PROVOST A M.SUTRA:A model for 2D or 3D saturated-unsaturated,variable-density ground-water flow with solute or energy transport[R].Reston, U.S: Geological Survey,2002.

    [17]ROBINSON C,LI L,BARRY D A.Effect of tidal forcing on a subterranean estuary[J].Advances in Water Resources,2007,30(4):851-865.

    [18]MO Yuming,JIN Guangqiu,TANG Hongwu,et al.Effects of deep pools on salinization in coastal reservoirs[C]//Proceedings of the 38th IAHR World Congress.Panama:IAHR,2019:3974-3985.

    [19]MEHDIZADEH S S,KETABCHI H,GHOROQI M,et al.Experimental and numerical assessment of saltwater recession in coastal aquifers by constructing check dams[J].Journal of Contaminant Hydrology,2020,231:103637.

    [20]ANWAR H O.The effect of a subsurface barrier on the conservation of freshwater in coastal aquifers[J].Water Research,1983,17(10):1257-1265.

    [21]LUYUN R J,MOMII K,NAKAGAWA K.Effects of recharge wells and flow barriers on seawater intrusion[J].Ground Water,2011,49(2):239-249.

    [22]NIELSEN P.Tidal dynamics of the water table in beaches[J].Water Resources Research,1990,26(9):2127-2134.

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

    [24]KUAN W K,JIN Guangqiu,XIN Pei,et al.Tidal influence on seawater intrusion in unconfined coastal aquifers[J].Water Resources Research,2012,48(2):W02502.

    [25]XIN Pei,ROBINSON C,LI Ling,et al.Effects of wave forcing on a subterranean estuary[J].Water Resources Research,2010,46(12):W12505.

    [26]李運(yùn)生,李孟彪,王亞萍,等.鋼-混凝土波形鋼腹板組合梁彎扭性能研究[J].河北科技大學(xué)學(xué)報(bào),2022,43(1):99-109.

    LI Yunsheng,LI Mengbiao,WANG Yaping,et al.Research on bending-torsional behaviors of steel-concrete composite beam with corrugated steel webs[J].Journal of Hebei University of Science and Technology,2022,43(1):99-109.

    [27]白衛(wèi)峰,亢勇,呂亞軍,等.基于廢渣粉末的再生水泥基材料性能試驗(yàn)研究[J].河北工業(yè)科技,2022,39(3):245-252.

    BAI Weifeng,KANG Yong,LYU Yajun,et al.Experimental study on properties of recycled cement-based materials based on waste slag powder[J].Hebei Journal of Industrial Science and Technology,2022,39(3):245-252.

    猜你喜歡
    河口
    A New Species of Cyrtodactylus (Squamata:Gekkonidae) from the Karst Forests of Daweishan National Nature Reserve,Yunnan,China
    河口沖淤對(duì)圍墾的響應(yīng)
    河口正陷入危機(jī)
    他們?yōu)槭裁催x擇河口
    河口,我們的家
    特殊的河口水
    河口的類型
    河口的樣子
    河口
    河口正陷入危機(jī)
    天天操日日干夜夜撸| 黑人操中国人逼视频| 久热爱精品视频在线9| 老熟女久久久| av有码第一页| 午夜精品久久久久久毛片777| 亚洲成国产人片在线观看| 黄色a级毛片大全视频| 在线十欧美十亚洲十日本专区| 欧美成人午夜精品| 国产精品久久久人人做人人爽| 欧美激情极品国产一区二区三区| 久久性视频一级片| 成人18禁高潮啪啪吃奶动态图| 夜夜骑夜夜射夜夜干| 免费在线观看影片大全网站| 99国产精品免费福利视频| 电影成人av| 中文字幕色久视频| 欧美日韩国产mv在线观看视频| 一级黄色大片毛片| 80岁老熟妇乱子伦牲交| 亚洲av日韩在线播放| 国产精品秋霞免费鲁丝片| 黄色a级毛片大全视频| 一边摸一边做爽爽视频免费| 亚洲成人国产一区在线观看| 青草久久国产| 50天的宝宝边吃奶边哭怎么回事| 亚洲三区欧美一区| 欧美老熟妇乱子伦牲交| 免费看十八禁软件| 国产精品欧美亚洲77777| 免费在线观看视频国产中文字幕亚洲 | 19禁男女啪啪无遮挡网站| 国产精品成人在线| 成人18禁高潮啪啪吃奶动态图| av不卡在线播放| 亚洲专区中文字幕在线| 人人妻人人澡人人爽人人夜夜| 欧美日韩成人在线一区二区| 两个人看的免费小视频| 国产精品久久久久久精品古装| 成人av一区二区三区在线看 | 久久天躁狠狠躁夜夜2o2o| 中文精品一卡2卡3卡4更新| 久久精品国产综合久久久| 欧美日韩亚洲高清精品| 19禁男女啪啪无遮挡网站| 国产熟女午夜一区二区三区| 国精品久久久久久国模美| 亚洲精品在线美女| 国产欧美亚洲国产| 国产一区二区在线观看av| 99国产精品99久久久久| 精品亚洲成a人片在线观看| 久久久国产精品麻豆| 波多野结衣一区麻豆| 亚洲av电影在线观看一区二区三区| 99国产精品一区二区三区| 亚洲少妇的诱惑av| 免费观看人在逋| 久久久精品国产亚洲av高清涩受| 色视频在线一区二区三区| 侵犯人妻中文字幕一二三四区| 搡老岳熟女国产| 亚洲人成77777在线视频| 国产熟女午夜一区二区三区| 中文字幕高清在线视频| 精品人妻一区二区三区麻豆| 叶爱在线成人免费视频播放| 精品国产一区二区久久| 大片免费播放器 马上看| 久久久久久久国产电影| 国产成人精品在线电影| 中文字幕色久视频| 九色亚洲精品在线播放| 久久亚洲国产成人精品v| a级毛片在线看网站| www.精华液| 亚洲一区中文字幕在线| 国产不卡av网站在线观看| 最近最新免费中文字幕在线| 午夜免费鲁丝| 亚洲av成人不卡在线观看播放网 | 亚洲精品国产av成人精品| 国产深夜福利视频在线观看| 男女免费视频国产| 久久ye,这里只有精品| 日韩 亚洲 欧美在线| 国产精品久久久人人做人人爽| 天堂中文最新版在线下载| 爱豆传媒免费全集在线观看| 男人爽女人下面视频在线观看| av天堂在线播放| 无遮挡黄片免费观看| 亚洲av电影在线进入| 国产精品香港三级国产av潘金莲| 老司机影院毛片| 麻豆国产av国片精品| 久久99热这里只频精品6学生| 黄网站色视频无遮挡免费观看| 人妻 亚洲 视频| 亚洲黑人精品在线| 国产精品一区二区免费欧美 | 亚洲天堂av无毛| 精品第一国产精品| 亚洲专区中文字幕在线| 男女之事视频高清在线观看| 久久久精品94久久精品| 悠悠久久av| 91精品三级在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 啦啦啦中文免费视频观看日本| 国产精品一二三区在线看| 精品人妻在线不人妻| 大香蕉久久网| 18禁裸乳无遮挡动漫免费视频| 美国免费a级毛片| 精品国产一区二区三区四区第35| 香蕉丝袜av| 国产亚洲欧美精品永久| 亚洲精品成人av观看孕妇| 亚洲,欧美精品.| 波多野结衣一区麻豆| 精品少妇内射三级| 久久精品熟女亚洲av麻豆精品| 国产福利在线免费观看视频| 国产男女内射视频| 麻豆乱淫一区二区| 国产精品一区二区在线观看99| 日本wwww免费看| 午夜福利影视在线免费观看| 欧美老熟妇乱子伦牲交| 别揉我奶头~嗯~啊~动态视频 | 国产精品.久久久| 国产成人一区二区三区免费视频网站| 午夜福利影视在线免费观看| 国产伦理片在线播放av一区| 日本vs欧美在线观看视频| 夜夜骑夜夜射夜夜干| 精品一品国产午夜福利视频| 久久香蕉激情| 免费日韩欧美在线观看| 中文字幕人妻丝袜一区二区| 日本vs欧美在线观看视频| 色婷婷久久久亚洲欧美| 美女高潮喷水抽搐中文字幕| 国产区一区二久久| 亚洲专区国产一区二区| 日韩熟女老妇一区二区性免费视频| 中亚洲国语对白在线视频| 黑人欧美特级aaaaaa片| 丁香六月天网| 久久精品成人免费网站| 丁香六月欧美| avwww免费| 成人手机av| 精品高清国产在线一区| 精品视频人人做人人爽| 国产欧美日韩综合在线一区二区| 亚洲国产欧美网| 日韩精品免费视频一区二区三区| 丝袜人妻中文字幕| 美国免费a级毛片| 少妇被粗大的猛进出69影院| 在线观看免费午夜福利视频| 国产精品1区2区在线观看. | 老司机靠b影院| 久久国产精品大桥未久av| 人妻久久中文字幕网| 国产精品一二三区在线看| 国产成人精品在线电影| 两个人看的免费小视频| av天堂在线播放| 男女边摸边吃奶| 久久99热这里只频精品6学生| 女人被躁到高潮嗷嗷叫费观| 亚洲国产精品一区二区三区在线| 国产高清国产精品国产三级| 国产亚洲av片在线观看秒播厂| 亚洲欧美精品自产自拍| 国产av国产精品国产| 精品亚洲成a人片在线观看| 大型av网站在线播放| 黑人欧美特级aaaaaa片| 国产精品久久久久久精品古装| 国产一区有黄有色的免费视频| 性色av乱码一区二区三区2| 一级毛片电影观看| av线在线观看网站| 女性生殖器流出的白浆| av又黄又爽大尺度在线免费看| 精品卡一卡二卡四卡免费| 国产精品久久久久久精品电影小说| 1024香蕉在线观看| 少妇的丰满在线观看| 久久人人爽人人片av| 大香蕉久久网| 国产一区二区 视频在线| 国产野战对白在线观看| 最黄视频免费看| 久久久久网色| 人人妻人人澡人人爽人人夜夜| 国产成人欧美在线观看 | 中文字幕高清在线视频| 欧美黄色片欧美黄色片| 熟女少妇亚洲综合色aaa.| 国产精品久久久久久精品古装| 国产亚洲欧美在线一区二区| 高清在线国产一区| 国产精品久久久av美女十八| 久久精品久久久久久噜噜老黄| 无限看片的www在线观看| 国产一卡二卡三卡精品| 自拍欧美九色日韩亚洲蝌蚪91| 亚洲精品第二区| 成人免费观看视频高清| av天堂久久9| 男女午夜视频在线观看| 麻豆国产av国片精品| 欧美激情久久久久久爽电影 | 丝袜人妻中文字幕| 五月天丁香电影| 十八禁网站网址无遮挡| 亚洲国产中文字幕在线视频| 精品亚洲成a人片在线观看| 汤姆久久久久久久影院中文字幕| 黑人操中国人逼视频| kizo精华| 免费不卡黄色视频| av免费在线观看网站| 美女高潮到喷水免费观看| 精品欧美一区二区三区在线| 十八禁高潮呻吟视频| 视频在线观看一区二区三区| 9191精品国产免费久久| 亚洲久久久国产精品| 可以免费在线观看a视频的电影网站| a 毛片基地| 90打野战视频偷拍视频| 国产日韩一区二区三区精品不卡| 欧美av亚洲av综合av国产av| 国产精品秋霞免费鲁丝片| 黑人巨大精品欧美一区二区mp4| av有码第一页| 色视频在线一区二区三区| 精品一区二区三卡| 色婷婷av一区二区三区视频| 午夜两性在线视频| 亚洲精华国产精华精| 黑人操中国人逼视频| 人人妻人人爽人人添夜夜欢视频| 国产成人免费观看mmmm| 欧美 日韩 精品 国产| 99国产精品一区二区蜜桃av | 美女扒开内裤让男人捅视频| 在线观看人妻少妇| av福利片在线| 国产成人欧美| 十八禁人妻一区二区| 亚洲精品国产av成人精品| 一二三四社区在线视频社区8| 一本—道久久a久久精品蜜桃钙片| 国产精品1区2区在线观看. | 日本91视频免费播放| 男人操女人黄网站| 黄网站色视频无遮挡免费观看| 精品国产一区二区三区久久久樱花| www.熟女人妻精品国产| 男女无遮挡免费网站观看| bbb黄色大片| 亚洲精品第二区| 黄频高清免费视频| 王馨瑶露胸无遮挡在线观看| 性高湖久久久久久久久免费观看| 国产激情久久老熟女| 精品亚洲乱码少妇综合久久| 两性夫妻黄色片| 欧美日韩精品网址| 久久精品国产亚洲av高清一级| 高潮久久久久久久久久久不卡| 波多野结衣av一区二区av| 亚洲一卡2卡3卡4卡5卡精品中文| 国产视频一区二区在线看| 亚洲七黄色美女视频| 啦啦啦在线免费观看视频4| 久久精品aⅴ一区二区三区四区| 亚洲人成电影观看| 狂野欧美激情性bbbbbb| 亚洲欧美精品自产自拍| 一本久久精品| 成在线人永久免费视频| 国产99久久九九免费精品| 国产精品国产av在线观看| 极品人妻少妇av视频| 国产成人欧美| 国产欧美日韩一区二区精品| 日韩一区二区三区影片| 欧美精品人与动牲交sv欧美| 热re99久久国产66热| 97在线人人人人妻| 国产精品av久久久久免费| 日本撒尿小便嘘嘘汇集6| 夫妻午夜视频| 高清欧美精品videossex| 999久久久国产精品视频| 久久天堂一区二区三区四区| 国产精品成人在线| 两个人看的免费小视频| 婷婷成人精品国产| 岛国毛片在线播放| 性少妇av在线| 久久精品亚洲熟妇少妇任你| 两个人免费观看高清视频| 亚洲美女黄色视频免费看| 91精品国产国语对白视频| 国产一区有黄有色的免费视频| 91av网站免费观看| 欧美精品av麻豆av| 美女视频免费永久观看网站| 亚洲av日韩精品久久久久久密| 女警被强在线播放| 免费在线观看视频国产中文字幕亚洲 | 亚洲精品粉嫩美女一区| 首页视频小说图片口味搜索| 侵犯人妻中文字幕一二三四区| 99国产精品免费福利视频| 黄色片一级片一级黄色片| 99久久99久久久精品蜜桃| 国产精品一区二区精品视频观看| 久久九九热精品免费| 9色porny在线观看| 免费在线观看影片大全网站| 国产淫语在线视频| 伊人亚洲综合成人网| 伦理电影免费视频| 三上悠亚av全集在线观看| 香蕉丝袜av| 丰满迷人的少妇在线观看| 日韩制服丝袜自拍偷拍| 国产成人精品无人区| 亚洲av电影在线观看一区二区三区| 国产无遮挡羞羞视频在线观看| 亚洲男人天堂网一区| 国产区一区二久久| 国产视频一区二区在线看| 91成年电影在线观看| 欧美日韩亚洲国产一区二区在线观看 | 岛国毛片在线播放| 国产99久久九九免费精品| 亚洲成人国产一区在线观看| 黄片播放在线免费| 国产91精品成人一区二区三区 | 91精品伊人久久大香线蕉| 亚洲国产日韩一区二区| 久久99热这里只频精品6学生| 久久天堂一区二区三区四区| 亚洲av电影在线观看一区二区三区| 免费少妇av软件| 国产97色在线日韩免费| 精品亚洲成国产av| 大码成人一级视频| 亚洲国产看品久久| 多毛熟女@视频| 国产又爽黄色视频| www.精华液| 午夜免费观看性视频| 大片电影免费在线观看免费| 久久亚洲精品不卡| 男女高潮啪啪啪动态图| a级片在线免费高清观看视频| 高清av免费在线| 午夜影院在线不卡| 国产男女内射视频| 91麻豆精品激情在线观看国产 | 亚洲成av片中文字幕在线观看| 精品一区在线观看国产| 18禁国产床啪视频网站| 久久久久久久久免费视频了| 亚洲精品美女久久av网站| 亚洲色图综合在线观看| 国产亚洲欧美在线一区二区| 亚洲欧美一区二区三区久久| 亚洲欧美一区二区三区黑人| 色精品久久人妻99蜜桃| 一边摸一边做爽爽视频免费| 日本av免费视频播放| 国产无遮挡羞羞视频在线观看| 中文字幕人妻丝袜制服| 老司机影院成人| 大型av网站在线播放| 9热在线视频观看99| 母亲3免费完整高清在线观看| av不卡在线播放| 男女免费视频国产| av免费在线观看网站| av超薄肉色丝袜交足视频| 狂野欧美激情性bbbbbb| 国产色视频综合| 热re99久久精品国产66热6| 免费看十八禁软件| 美女视频免费永久观看网站| 18禁黄网站禁片午夜丰满| 在线亚洲精品国产二区图片欧美| 国产色视频综合| 999久久久精品免费观看国产| 一区二区三区精品91| 999久久久精品免费观看国产| 久久精品亚洲av国产电影网| 丝袜在线中文字幕| 中文字幕人妻熟女乱码| 狠狠精品人妻久久久久久综合| 99精品久久久久人妻精品| 纯流量卡能插随身wifi吗| 99精品久久久久人妻精品| 日韩制服丝袜自拍偷拍| 黄色片一级片一级黄色片| 黄网站色视频无遮挡免费观看| 俄罗斯特黄特色一大片| 精品一区二区三卡| 十分钟在线观看高清视频www| 十八禁高潮呻吟视频| 爱豆传媒免费全集在线观看| 两个人免费观看高清视频| 久久久久精品国产欧美久久久 | 十分钟在线观看高清视频www| 女人高潮潮喷娇喘18禁视频| 大码成人一级视频| 久热这里只有精品99| 久久热在线av| 国产成人精品在线电影| 国产精品熟女久久久久浪| 中亚洲国语对白在线视频| 精品国产超薄肉色丝袜足j| 亚洲人成电影免费在线| 十八禁高潮呻吟视频| 自拍欧美九色日韩亚洲蝌蚪91| 在线观看人妻少妇| 成人三级做爰电影| 1024香蕉在线观看| 午夜福利视频精品| 久久国产精品大桥未久av| 欧美性长视频在线观看| 人成视频在线观看免费观看| 97精品久久久久久久久久精品| 亚洲av欧美aⅴ国产| 亚洲成人免费电影在线观看| 亚洲欧美成人综合另类久久久| 国产深夜福利视频在线观看| 大片电影免费在线观看免费| 亚洲av美国av| 国产免费一区二区三区四区乱码| 国产男女超爽视频在线观看| av电影中文网址| 美女大奶头黄色视频| 动漫黄色视频在线观看| 菩萨蛮人人尽说江南好唐韦庄| www.自偷自拍.com| 亚洲成国产人片在线观看| 国产av一区二区精品久久| 精品国产乱子伦一区二区三区 | 欧美人与性动交α欧美软件| 免费看十八禁软件| 天天躁夜夜躁狠狠躁躁| 成人国语在线视频| 国产av一区二区精品久久| 日本av手机在线免费观看| 亚洲久久久国产精品| 国产成人av激情在线播放| 亚洲精品美女久久av网站| 人妻人人澡人人爽人人| 午夜福利在线观看吧| 国产成人欧美| 久久久久视频综合| 亚洲国产精品999| 在线观看免费高清a一片| 久久人妻福利社区极品人妻图片| 欧美性长视频在线观看| 成人手机av| 97在线人人人人妻| 欧美另类一区| 19禁男女啪啪无遮挡网站| 国产精品久久久久成人av| 亚洲精华国产精华精| 欧美 日韩 精品 国产| 一本—道久久a久久精品蜜桃钙片| 波多野结衣av一区二区av| 亚洲欧洲精品一区二区精品久久久| 亚洲成av片中文字幕在线观看| 免费一级毛片在线播放高清视频 | 免费观看a级毛片全部| 成人av一区二区三区在线看 | 91成人精品电影| 久久久国产一区二区| 亚洲第一青青草原| 国产1区2区3区精品| 亚洲欧美一区二区三区黑人| 啦啦啦啦在线视频资源| 日本撒尿小便嘘嘘汇集6| 日韩精品免费视频一区二区三区| 欧美另类一区| 亚洲精品粉嫩美女一区| 在线十欧美十亚洲十日本专区| 国内毛片毛片毛片毛片毛片| 九色亚洲精品在线播放| 99香蕉大伊视频| 国产又爽黄色视频| 精品国产乱子伦一区二区三区 | 久久久精品94久久精品| 国产精品二区激情视频| 搡老熟女国产l中国老女人| av超薄肉色丝袜交足视频| 不卡一级毛片| 亚洲精品乱久久久久久| 久久毛片免费看一区二区三区| 自拍欧美九色日韩亚洲蝌蚪91| 老熟女久久久| 大片电影免费在线观看免费| 在线观看免费视频网站a站| 狂野欧美激情性xxxx| 无遮挡黄片免费观看| 国产日韩一区二区三区精品不卡| 国产一区二区三区综合在线观看| 国产av国产精品国产| 黄频高清免费视频| 欧美日韩黄片免| 女警被强在线播放| 日韩电影二区| 一区二区三区四区激情视频| 欧美xxⅹ黑人| 不卡一级毛片| 91麻豆av在线| 十八禁网站网址无遮挡| 岛国在线观看网站| 在线观看免费午夜福利视频| 欧美97在线视频| 久久久久久久精品精品| 俄罗斯特黄特色一大片| 韩国高清视频一区二区三区| 国产极品粉嫩免费观看在线| 久久精品国产亚洲av香蕉五月 | 久久久久久亚洲精品国产蜜桃av| 丁香六月天网| 亚洲精品国产色婷婷电影| 人人妻人人澡人人爽人人夜夜| 久久 成人 亚洲| 亚洲一码二码三码区别大吗| 午夜视频精品福利| √禁漫天堂资源中文www| 美女高潮到喷水免费观看| 一边摸一边抽搐一进一出视频| 久久精品成人免费网站| 新久久久久国产一级毛片| 青春草亚洲视频在线观看| 91av网站免费观看| 国产一级毛片在线| 亚洲avbb在线观看| 国产精品九九99| 9热在线视频观看99| 三级毛片av免费| 国产欧美日韩一区二区三 | 亚洲av美国av| 国产在视频线精品| 久久性视频一级片| 啦啦啦免费观看视频1| 亚洲 欧美一区二区三区| 亚洲黑人精品在线| 久久精品国产a三级三级三级| 中文精品一卡2卡3卡4更新| 啦啦啦 在线观看视频| 欧美日韩精品网址| 久久久精品国产亚洲av高清涩受| avwww免费| 国产成人精品久久二区二区免费| 伊人亚洲综合成人网| 国产在线一区二区三区精| 中国国产av一级| 男女之事视频高清在线观看| 飞空精品影院首页| 午夜成年电影在线免费观看| 精品久久久久久久毛片微露脸 | 精品国产乱码久久久久久男人| 欧美成人午夜精品| 亚洲免费av在线视频| 十八禁网站免费在线| 精品福利观看| 老司机午夜十八禁免费视频| 亚洲男人天堂网一区| 99国产精品99久久久久| 18在线观看网站| 国产精品久久久av美女十八| 免费一级毛片在线播放高清视频 | 免费av中文字幕在线| 精品久久久久久久毛片微露脸 | 狠狠婷婷综合久久久久久88av| 亚洲成人免费电影在线观看| 大码成人一级视频| 成人免费观看视频高清| 亚洲中文字幕日韩| 婷婷丁香在线五月| 91成年电影在线观看| 久久久久久免费高清国产稀缺| 国产成人免费观看mmmm| 精品人妻熟女毛片av久久网站| 久久久久网色| 精品久久久久久电影网| 一本一本久久a久久精品综合妖精| 超碰成人久久| 人人妻人人澡人人看| 国产精品国产av在线观看| 老司机深夜福利视频在线观看 | 久久久国产欧美日韩av| 精品国内亚洲2022精品成人 | 午夜福利乱码中文字幕|