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

    Self-adaptive Green-Ampt inf i ltration parameters obtained from measured moisture processes

    2016-04-18 10:34:58WenwenLingYongshuZhuLiChenbZhongboYu
    Water Science and Engineering 2016年3期

    *,Wen-wen LingYong-shu ZhuLi Chenb,Zhong-bo Yu

    aState Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,Hohai University,Nanjing 210098,China

    bDesert Research Institute,Las Vegas,NV 89119,USA

    Self-adaptive Green-Ampt inf i ltration parameters obtained from measured moisture processes

    Long Xianga,*,Wen-wen Linga,Yong-shu Zhua,Li Chena,b,Zhong-bo Yua

    aState Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,Hohai University,Nanjing 210098,China

    bDesert Research Institute,Las Vegas,NV 89119,USA

    Available online 30 September 2016

    The Green-Ampt(G-A)inf i ltration model(i.e.,the G-A model)is often used to characterize the inf i ltration process in hydrology.The parameters of the G-A model are critical in applications for the prediction of inf i ltration and associated rainfall-runoff processes.Previous approaches to determining the G-A parameters have depended on pedotransfer functions(PTFs)or estimates from experimental results,usually without providing optimum values.In this study,rainfall simulators with soil moisture measurements were used to generate rainfall in various experimental plots.Observed runoff data and soil moisture dynamic data were jointly used to yield the inf i ltration processes,and an improved self-adaptive method was used to optimize the G-A parameters for various types of soil under different rainfall conditions.The two G-A parameters,i.e.,the effective hydraulic conductivity and the effective capillary drive at the wetting front,were determined simultaneously to describe the relationships between rainfall,runoff,and inf i ltration processes.Through a designed experiment,the method for determining the GA parameters was proved to be reliable in ref l ecting the effects of pedologic background in G-A type inf i ltration cases and deriving the optimum G-A parameters.Unlike PTF methods,this approach estimates the G-A parameters directly from inf i ltration curves obtained from rainfall simulation experiments so that it can be used to determine site-specif i c parameters.This study provides a self-adaptive method of optimizing the G-A parameters through designed field experiments.The parameters derived from field-measured rainfall-inf i ltration processes are more reliable and applicable to hydrological models.

    Green-Ampt model;Levenberg-Marquardt algorithm;Parameter optimization;Ungauged basin;Pedotransfer function

    1.Introduction

    Inf i ltration plays an important role in terrestrial hydrologic processes.It affects the runoff generation process and is dramatically inf l uenced by the soil hydraulic properties and soil porous texture(Sivakumar,2015).Most physically based hydrologic models have described the relationships between rainfall,runoff,and inf i ltration processes implicitly or explicitly(Lee et al.,2015).The Green-Ampt(G-A)inf i ltration model(Green and Ampt,1911)is one such model based on the soil porous media characteristics(Prevedello et al.,2009).It is widely used in the hydrologic field due to its reasonable physical mechanism and easy-to-solve solution(Govindaraju et al.,1996;Ma et al.,2010;O'Brienet al.,2009;Silburn and Connolly,1995;Wang et al.,2010). For example,both the FLO-2D model in Maricopa County and the Water Erosion Prediction Project(WEPP)model of the U.S.Department of Agriculture(USDA)(Dun et al., 2009;Risse et al.,1995)have adopted the G-A model for rainfall-runoff prediction.Meanwhile,many researchers have improved the G-A model in order to adapt it to more complex soil systems(Gowdish and Munoz-Carpena,2009).They have focused on solving the equations and estimating its parameters theoretically(Regalado et al.,2005;Verbist et al., 2010).The most common method is to make use of pedotransfer functions(PTFs)(Brooks and Corey,1966)to calculate the parameters of the G-A model(Rawls et al., 1983;Regalado et al.,2005).They also suggest that it is reasonable to use half of the saturation conductivity as the effective conductivity(Bouwer,1966).This hypothesis has been commonly accepted in practice.However,it has limitations when the G-A model is used in the field(van den Putte et al.,2013).In order to avoid the limitations,ring inf i ltrometers were used to measure the hydraulic conductivity by imposing ponded conditions in some field experimental methods(Angermann et al.,2002;Esteves et al.,2000; Galbiati and Savi,1995;Mohamoud,1991;Reynolds,2010; Suleiman and Swartzendruber,2003).However,this method is not f i t for non-ponded initial conditions when the G-A model is used.On the other hand,some researchers have used rainfall simulators for the parameterization of the G-A model (Esteves et al.,2000;Rawls et al.,1992;Suleiman and Swartzendruber,2003;Valiantzas,2010;Taskinen et al., 2008).In those studies,the measured runoff data(the runoff process is an indirect part of the inf i ltration process) were used to calibrate the G-A model and to generate best fi tting parameters for the model.However,the uncertainty of observed data series decreases the reliability of the parameters and at the same time prevents derivation of multiparameters to describe the relationships between the coupled equations put forward by Athira and Sudheer(2015). To solve these problems,we designed an inf i ltration equipment for small-scale rainfall experiments through a moisture survey and developed a parameter optimization algorithm to derive the G-A parameters using the measured data.

    In this study,rainfall simulators were set up in various experimental plots based on a typical pedologic background. Accumulated inf i ltration curves and moisture dynamic data were obtained in the experiments.High-intensity rainfall experiments were conducted in each plot and inf i ltration curves were yielded under the simulated rainfall conditions.Then,a new self-adaptive optimization approach based on the Levenberg-Marquardt(LM)algorithm(Marquardt,1963)was validated theoretically,in order to estimate the G-A parameters directly from these data,and a quantified relationship between soil textures and the G-A parameters was built through comparison with the parameters derived directly from the PTF method.Based on these data,the method for extraction of the G-A parameters based on the pedologic background was developed for high-resolution distributed hydrologic models.

    2.Materials and methods

    2.1.Optimization setup

    Since the G-A model was presented by Green and Ampt (1911),it has been modified by several researchers.Mein and Larson(1973)extended the model from ponded conditions to constant intensity conditions.Chu(1978)applied this model to unsteadyrainfallintensities.Inthese studies,the G-Amodel was treated as two parts:under the steady state,the inf i ltration rate equals the rainfall intensity before ponding;as the wetting front moves downwards with time,ponding occurs and the integrated version of the G-A model after ponding can be expressed as

    where I is the vertical accumulative inf i ltration depth;S is the soil capillary drive at the wetting front;θiand θsare initial and saturated water contents,respectively;Keis the effective hydraulic conductivity;t is time;the tpis the ponding time,and tp=Ip/P,in which P is the rainfall intensity(P>Ke)and Ipis the inf i ltration depth at tp,which is calculated as follows:

    tsis a virtual time def i ned as follows:

    To implicitly calculate I in Eq.(1),four parameters are needed:Ke,S,θs,and θi.In Eq.(1),S,θs,and θialways have an integrity form of S(θs-θi).This can be simplified to one parameter M,where M=S(θs-θi).The parameters to be optimized are then reduced to two:Keand M.

    The self-adaptive optimization algorithms can be categorized as local and global search methods.Depending on the hill-climbing strategy,search algorithms can be divided into direct and gradient-based methods.Gradient-based methods use the information about the gradient of the objective function and direct search methods use only the information about the objective function value.In this study,we chose a gradientbased method,the LM algorithm,as our optimization method. The general object of the LM algorithm is to minimize the sum of the square residuals(Eq.(4))by gradually changing the optimized parameters.Its objective function is assumed to be the nonlinear least squares problems as follows:

    where yiis the observed data series at time step i,f(xi,β)is the optimized series using the optimized parameter β,xiis a variable,and m is the maximum series number.f(xi,β)in this study was solved implicitly through Eq.(1)using the Newton method because the G-A model is an implicit function for accumulated inf i ltration series.Each time step is recorded as i(i=1,2,···,t)in correspondence to accumulated inf i ltration at time t.As shown in Eq.(1),the total errors between simulated and observed values for accumulated inf i ltration can be provided by Eq.(5)at different duration times related to the same parameters through the Newton method.In Eq. (5),we def i ne the objective function,which is the same as that required in the LM algorithm,but the change of parameters changes the whole inf i ltration process.The objective function in Eq.(5)focuses on the parameters'domains and their optimized climbing paths.Under the two parameter conditions,the objective function can be written as

    where I is the vertical accumulative inf i ltration depth,and its observed value in this equation is related to time t;and the ponding time is related to M and Kein various combinations. We rebuilt tsand Ipas the function,so the Jacobian of the objective function can be derived as follows:

    Using tp=Ip/P and Eq.(3),we obtain

    Similarly,using Eqs.(2),(3),and(5),we obtain the following equations:

    Using Eqs.(9)through(11),we obtain

    Thus,Eqs.(8)and(12)compose a Jacobian,which is used in the optimization computation through consideration of rainfall intensity and ponding time.The various steps can be calculated by Eq.(13)in each time interval:

    where Xtis a group of parameters at time t,Xt+1is the parameter's value at the time following time t,J is the Jacobian matrix mentioned above,F(X)is the target optimized function,such as in Eq.(5),λ is a constant coeff i cient,and D is a constant number.When λ is very small,the optimized search step approaches that of the Newton method;when λ is large enough,the optimized search step is close to that of the gradient descent method.Using the input Jacobian for the objective function,the LM algorithm can search the optimum parameters.

    2.2.Theoretical calibration of optimization method

    The structure of the solution is the basis for searching optimization parameters.The existence of a unique global solution should be examined.Then,the impact factors(e.g.,recharge intensity,duration,and ponded conditions)should be validated theoretically for the solution scheme.Fig.1 shows a straightforward method(the approaching method)of calculating the minimumrootmeansquareerror(RMSE)betweenthesolutions of the G-A model and the Richards equation for all possible Keand M parameter combinations.Different recharge intensities and durations were tested for three types of theoretical soils (sandy,sandy clay loam,and clay loam),as shown in Figs.1(a through c).The RMSE,as the objective function,can well describe the structure of solution for the G-A model.Numerical analysis also shows that the optimized parameters are superior to the parameters extracted from theoretical and semitheoretical parameter estimation methods(e.g.,yield from Brooks and Corey(1966)or van Genuchten(1980)).

    In this study,the general object of the LM algorithm was to minimize the residual sum of squares,which is similar to the value of the RMSE.The automatic search program was revised in our study.A problem in the LM algorithm is the result of non-unique or local solutions,which are derived from numerical errors or uncertainty of input data(More et al., 1980).In order to solve this problem,the MINPACK-1 package from the Argonne National Laboratory(More et al., 1980)can be used to extract the G-A parameters from the theoretical inf i ltration curves and numerical approximation curves in various soils by repeated iterations.All locally optimized solutions can be compared and a global optimum solution can be obtained by identifying the smallest value of the objective function.This approach can avoid random errors in calculations(Hristopulos,2015).For example,as shown in Fig.2,the initial value affects the f i nal optimum parameters due to uneven distribution of objective functions.However, multiple trials under the various initial conditions should generate a series of objective function values whose minimum error is the best location for solution.Thus,the modifiedmethod is capable of searching the parameters of the G-A model in complex inf i ltration processes.

    Fig.1.RMSE distributions in approaching method for selected soils.

    Fig.2.Multiple trials for avoiding initial value problems.

    2.3.Experimental framework

    A series of rainfall simulation tests was conducted in each representative mapping unit in the Rainbow Wash and White Tank watersheds located in Arizona,in the United States.The rainfall simulator(RFS)was used in the tests(Bhardwaj and Singh,1992;Munn and Huntington,1976;Mutchler and Moldenhauer,1963).The cover area was 61× 61 cm2. Water drops were produced on the needles by providing a constant gravity head directly beneath the rainfall simulator. The recharge rates were designed for recording precipitation at 54 mm/h.In this condition,the experimental frame was close to an actual rainfall event.Unlike the traditional ponding inf i ltration experiment,most of the rainfall characteristics were retained in this study.For pedologic background, approximately 200 g of soil were collected from a depth ranging from 0 to 10 cm at each of the sampling locations within each mapping unit,and the soil particle size distributions and bulk densities were analyzed.Tested soils were classified by particle size distribution,as shown in Fig.3. Statistical data of the particle size distribution of soil in the experimental plots are shown in Table 1.

    To obtain the accumulated water inf i ltration,we measured the water content in soil porous media beneath the RFS and the surface runoff at the down-sloping boundary.A probe with a length of 13 cm was installed at an angle of 30。with respect to the horizontal plane to measure the water content, which means that the effective depth for analysis was 6.5 cm. The collected data were recorded manually and continuously using a water content recorder(WCR)(model CS-616, Campbell Scientif i c Inc.(CSI)Logan,UT,U.S.A.).The accumulative inf i ltration curves were then obtained from these data.Only measurements taken in the first 10—30 min in each test were adopted in order to ensure that the wetting front did not occur below the probe.The measurements were based on the estimated wetting front position using the rainfall intensity and collected runoff volume.In the experiment,the time of runoff occurrence was recorded to optimize the ponding time.

    Because soil structure strongly inf l uences inf i ltration and runoff characteristics,a semi-quantitative estimation of the soil structure was conducted.Although the soil structure was diff i cult to quantify,its inf l uences could be examined in order to explain unexpected measured results in rainfall simulator tests.Three parallel experiments were conducted for each plot unit.

    Fig.3.Soil catalog by particle size distribution for test plots'soils.

    Table 1Statistic analysis of soil's PSD(particle size distribution)at each plot.

    3.Results and discussion

    3.1.Evaluation of inf i ltration curves

    The observed and f i tting inf i ltration curves for each plot, representing three stochastic locations at each experimental plot,are labeled A,B,and C,in order to avoid uncertainty (Fig.4).For most experimental cases,the inf i ltration curves were drawn by continuous data from the WCR sensor in field experiments.However,some plot tests(e.g.,Q2A and Q0A, et al.)were recorded manually because there were coarse rocks on the plots'surfaces.The sensor could not identify soil layers or heterogeneous surface flows,so the monitoring data interrupted the homogenous hypothesis in the G-A model.In these cases,the accumulated inf i ltrations were calculated indirectly by subtracting observed runoff from rainfall at specif i c times.All measured and optimized predicted inf i ltration curves using the LM algorithm are shown in Fig.4.In most cases,the inf i ltration depths obtained from the WCR sensor ref l ect the same time series behaviors predicted by the G-A model.For these cases,the optimized accumulated infi ltrations agree with the observed ones.However,for some tests,the measured curves,which exhibit a distinctly nontheoretical S shape(e.g.,Q1A,Q2C,and Q3C in Fig.4),are different from the theoretical accumulated inf i ltration curves predicted by the G-A model.These non-theoretical cases were observed on older soils.Here,the inf i ltration rate was not a monotonically decreasing function,which existed for a theoretical uniform soil.Although the experiments were conducted in the same mapping unit,the inf i ltration processes were different.

    In these cases,the inf i ltration rates were low in the early period,then increased rapidly toward the middle period of the test,and f i nally decreased to lower and relatively stable values. In the field,the causes leading to S-shaped inf i ltration curves can be complex.Layered soil(Ma et al.,2010),the surface seal(Damodhara et al.,2006),preferential flow paths(Lepore et al.,2009),and the surface microtopography(Vˊazquez et al., 2005)affect the results.In experiments,heterogeneous soils in test plots and soil structures are major factors.Initially,the upper layer and microtopography decide the recharge rates.If the conditions of the upper layer and microtopography are the same,the curves are close to theoretical ones.When the heterogeneous layers appear,the accumulative curves change at the wetting front between the soil layer interfaces.The abrupt changes of recorded curves shown in Fig.4 depend on the difference between the soil layers.A stable inf i ltration rate is decided by the minimum permeation of the soil layers at specified plots.

    Using a standard G-A model,these non-theoretical inf i ltration curves cannot be matched.For these non-theoretical cases,the accumulated inf i ltration is overestimated by the G-A model before the measured inf i ltration rates increase signif i cantly.The reasons for the phenomena are not known entirely.However,we have found that the layering,presence of a surface crust,and microtopography can inf l uence the initial inf i ltration rates and preferential flow,causing the wetting front to penetrate the designed depth(6.5 cm)ahead of the expected time.The records show that a portion of the initial precipitation will puddle and prevent all water from inf i ltrating when clay particles exist.Most non-theoretical cases occurred on soils with the presence of silt(Q1 and Q2).This silt layer,with the higher water holding capacity and lower Ks(saturated conductivity),could have prevented water from moving deeper into the soil,within the range of inf l uence ofthe sensor.The non-theoreticalcasesappearmore frequently in the same plots(Q1 and Q2),implying that soil structure may contribute to the formation of non-theoretical curves.The proposed model cannot handle such complex curves,due to the hypotheses of the G-A model(such as a constant inf i ltration rate and uniform soil).

    3.2.Optimized f i tting parameters

    Fig.4.Extraction of optimized G-A parameters by LM algorithm from RFS experiments.

    Table 2Optimized G-A parameters for field RFS experiments and Saxton and Rawls's PTF(2006).

    The automatic search results from the method mentioned above are shown in Table 2 and the predicted data are plotted in Fig.4.For most of the cases,the measured inf i ltration curvehas an overall shape similar to the theoretical curve.For example, before ponding,the gradient of the curve(e.g.,an inf i ltration rate approximately equal torainfallintensity)remainsconstant; after ponding,the inf i ltration rate gradually decreases until it approaches a constant value for the remaining time.In these conditions,the results obtained from our proposed search approach agree with the observed series.In other cases,the measured curve has a different shape(e.g.,some measured curves in this study were S-shaped),and the proposed approach willnotbefullysatisfactory;fewerobservedpointsleadtomore fi ttingerrorsinQ2Aaswell.However,theoptimizedresultscan still describe the complex soil structure and other factors.For sandy soil cases(Q4A and Q2C),the in fi ltration curves were short and there were few sampling points due to the heterogeneous surface and a high in fi ltration rate,indicating that the sampling points cannot represent the experimental hypothesis. Thus,these cases were disregarded during parameter optimization.From the Nash-Sutcliffe test in Table 2,we note that the proposed approach could take into account all uncertainties and provide a series of representative parameters from the experimental data set.

    3.3.Analysis and comparison with previous PTF methods

    Modeling soil hydrologic processes for different landscape elements is important for many studies on land-use planning problems.The PTF has been developed as a simplified method of assessing soil hydraulic parameters obtained from soil physical properties,and it is much less laborious and less expensive.In model application,PTF can more routinely measurethesoilhydraulicpropertiesinungaugedbasins.Asfor the PTF in the G-A model,a commonly used one is Saxton and Rawls's PTF.According to Rawls et al.(1983),the KevaluerelatedtothisPTFishalfofKs.BasedonBrakensiek(1977)and Rawls et al.(1983),the soil capillary drive S from Saxton and Rawls's PTF is calculated by the following equations:

    where γ is the Brooks-Corey pore size distribution parameter,and ψbis the bubbling pressure(mm),assumed to be one half of the air entry value(Bouwer,1966;Brakensiek, 1977).In their studies,however,Saxton and Rawls(2006) considered the bubbling pressure to be equal to the air entry value.The calculated parameters are listed in Table 2. Fig.5 shows the comparison of the optimized parameter results and the estimated parameters by Saxton and Rawls's PTF(Saxton and Rawls,2006).Saxton and Rawls's PTF predicts higher Kevalues,but for most of the cases these values are ranked in the same order of magnitude as the optimized results.The S values predicted by the PTF generally range from 200 to 450 mm,whereas the optimized S values have a much larger range.For most cases,the PTF predicted values are higher than the optimized S values. Considering both Keand S effects,the parameters predicted by Saxton and Rawls's PTF tend to overestimate inf i ltration, and thus underestimate runoff when applied to hydrologic modeling.Nevertheless,the rainfall simulator tests were the closest to natural conditions,and the parameters derived from them represented the hydrologic units better than PTF-generated parameters.

    Fig.5.G-A parameters from RFS test and Saxton and Rawls's PTF in various plots.

    3.4.G-A parameters derived from rainfall-inf i ltration measurement

    When the set of PTF input parameters is def i ned,the PTF output can be obtained.In this study,multiple regression analysis was conducted to build the relationships between input and output parameters.An advantage of regression techniques is that most essential input parameters can be found automatically using stepwise regression.Multiple linear regressions is a parametric test,in which,for a given set of independent variables,the possible values for a dependent variable are assumed to be normally distributed and have a constant variance.Typically,the probability of Kein soils is logarithmically distributed,which requires log-transformed values of Keor lg Keto be the dependent variable in PTF. Correlations between the G-A parameters(Table 2)and soil texture are obtained from the optimization results of the RFS tests where soil information are also available,especially including the percentages of gravel,sand,silt,and clay,as well as bulk density(Brakensiek,1977;Brakensiek and Onstad, 1977).Effective hydraulic conductivity and soil capillary drive prediction using RFS tests are:

    where Gr,Sa,and Clare the percentages of gravel,sand,and clay,respectively,and Bdis the bulk density(g/cm3).Correlation coeff i cients for Keand S are 0.69 and 0.74,respectively. The relatively low values of the correlation coeff i cients in these regression equations ref l ect the inherent heterogeneity and uncertainty of natural soils.In addition,the limited number of successful experiments reduces the applicability of these equations.Although the specif i c form of these equations should be verified with a larger number of additional tests due to the limited RFS experiments,that form can ref l ect the integrated information for the specif i c locations.

    4.Conclusions

    The RFS test has the advantage of simulating the natural impact of rainfall and consequently eliminates the disadvantage associated with the traditional methods of obtaining hydrologic parameters in ungauged areas.It offers most of the integrated information for the inf i ltration process in specif i cexperimental locations.In this study,an optimized method of estimating the G-A parameters has been developed based on the RFS test.Unlike previous PTF methods,this method directly estimates the G-A parameters from in fi ltration curves. Through a parameter optimization procedure,the proposed method provides optimized results as to the effective hydraulic conductivity and the soil capillary drive.Further analyses show that the optimized parameters reproduce the measured in fi ltration curves.Field-observed non-theoretical in fi ltration curves(e.g.,S-shaped in fi ltration curves),usually due to complex soil structure or microtopography,are not adequately fi tted by the optimized parameters found by this method. However,the optimized results can provide the best approximation.Compared with the theoretical PTF approach,the optimized method for the G-A model is proper for estimating the soil hydraulic parameters for a speci fi c site.This study shows that the optimization approach with the RFS test can rapidly address potential error related to complex plot surface and soil structure,and,combined with the RFS test,it contributes to hydrologic modeling studies with improved parameters in ungauged watersheds.

    Angermann,T.,Wallender,W.W.,Wilson,B.W.,Werner,I.,Hinton,D.E., Oliver,M.N.,Zalom,F.G.,Henderson,J.D.,Oliveira,G.H.,Deanovic,L.A., et al.,2002.Runoff from orchard f l oors:Micro-plot field experiments and modeling.J.Hydrol.265(1—4),178—194.http://dx.doi.org/10.1016/S0022-1694(02)00109-9.

    Athira,P.,Sudheer,K.P.,2015.A method to reduce the computational requirement while assessing uncertainty of complex hydrological models. Stoch.Environ.Res.Risk Assess.29(3),847—859.http://dx.doi.org/ 10.1007/s00477-014-0958-4.

    Bhardwaj,A.,Singh,R.,1992.Development of a portable rainfall simulator inf i ltrometerforinf i ltration,runoffanderosionstudies.Agric.WaterManag. 22(3),235—248.http://dx.doi.org/10.1016/0378-3774(92)90028-U.

    Bouwer,H.,1966.Rapid field measurement of air entry value and hydraulic conductivity of soil as signif i cant parameters in flow system analysis.Water Resour.Res.2(4),729—738.http://dx.doi.org/10.1029/WR002i004p00729.

    Brakensiek,D.L.,1977.Estimating the effective capillary pressure in the Green and Ampt inf i ltration equation.Water Resour.Res.13(3),680—682. http://dx.doi.org/10.1029/WR013i003p00680.

    Brakensiek,D.L.,Onstad,C.A.,1977.Parameter estimation of the Green and Ampt inf i ltration equation.Water Resour.Res.13(6),1009—1977.http:// dx.doi.org/10.1029/WR013i006p01009.

    Brooks,R.H.,Corey,A.T.,1966.Properties of porous media affecting fl uid flow.J.Irrig.Drain.Div.72(IR2),61—88.

    Chu,S.T.,1978.In fi ltration during unsteady rain.Water Resour.Res.14(3), 461—466.http://dx.doi.org/10.1029/WR014i003p00461.

    Damodhara,R.M.,Raghuwanshi,N.S.,Singh,R.,2006.Development of a physically based 1D-in fi ltration model for irrigated soils.Agric.Water Manag.85(1—2),165—174.http://dx.doi.org/10.1016/j.agwat.2006.04.009.

    Dun,S.,Wu,J.Q.,Elliot,W.J.,Robichaud,P.R.,Flanagan,D.C., Frankenberger,J.R.,Brown,R.E.,Xu,A.C.,2009.Adapting the water erosion prediction project(WEPP)model for forest applications.J.Hydrol. 366(1—4),46—54.http://dx.doi.org/10.1016/j.jhydrol.2008.12.019.

    Esteves,M.,Faucher,X.,Galle,S.,Vauclin,M.,2000.Overland flow and in fi ltration modelling for small plots during unsteady rain:Numerical results versus observed values.J.Hydrol.228(3—4),265—282.http:// dx.doi.org/10.1016/S0022-1694(00)00155-4.

    Galbiati,G.,Savi,F.,1995.Evaluation of the comparative in fl uence of soil hydraulic properties and roughness on overland flow at the local scale.J. Agric.Eng.Res.61(3),183—190.http://dx.doi.org/10.1006/jaer.1995.1045.

    Govindaraju,R.S.,Kavvas,M.L.,Jones,S.E.,Rolston,D.E.,1996.Use of Green-Ampt model for analyzing one-dimensional convective transport in unsaturated soils.J.Hydrol.178(1—4),337—350.http://dx.doi.org/ 10.1016/0022-1694(95)02796-3.

    Gowdish,L.,Munoz-Carpena,R.,2009.An improved Green-Ampt inf i ltration and redistribution method for uneven multistorm series.Vadose Zone J. 8(2),470—479.http://dx.doi.org/10.2136/vzj2008.0049.

    Green,W.H.,Ampt,G.A.,1911.Studies on soil physics,part 1:The flow of air and water through soils.J.Agric.Sci.4(1),1—24.

    Hristopulos,D.,2015.Covariance functions motivated by spatial random field models with local interactions.Stoch.Environ.Res.Risk Assess.29(3), 739—754.http://dx.doi.org/10.1007/s00477-014-0933-0.

    Lee,T.,Shin,J.,Park,T.,Lee,D.,2015.Basin rotation method for analyzing the directional inf l uence of moving storms on basin response.Stoch.Environ.Res.Risk.Assess.29(1),251—263.http://dx.doi.org/10.1007/ s00477-014-0870-y.

    Lepore,B.J.,Morgan,C.L.S.,Norman,J.M.,Molling,C.C.,2009.A mesopore and matrix inf i ltration model based on soil structure.Geoderma 152(3—4), 301—313.http://dx.doi.org/10.1016/j.geoderma.2009.06.016.

    Ma,Y.,Feng,S.Y.,Su,D.Y.,Gao,G.Y.,Huo,Z.L.,2010.Modeling water inf i ltration in a large layered soil column with a modified Green-Ampt model and HYDRUS-1D.Comput.Electron.Agric.71(s1),S40—S47. http://dx.doi.org/10.1016/j.compag.2009.07.006.

    Marquardt,D.W.,1963.An algorithm for least squares estimation of nonlinear parameters.J.Soc.Ind.Appl.Math.11(2),431—441.

    Mein,R.G.,Larson,C.L.,1973.Modelinginf i ltrationduringasteadyrain.Water Resour.Res.9(2),384—394.http://dx.doi.org/10.1029/WR009i002p00384.

    Mohamoud,Y.M.,1991.Evaluating the Green and Ampt inf i ltration parameter values for tilled and crusted soils.J.Hydrol.123(1—2),25—38.http:// dx.doi.org/10.1016/0022-1694(91)90066-Q.

    More,J.J.,Garbow,B.S.,Hillstrom,K.E.,1980.User Guide for Minpack-1. Argonne National Laboratory,Argonne.

    Munn,J.R.,Huntington,G.L.,1976.Aportablerainfallsimulatorforerodibility and inf i ltration measurements in rugged terrain.Soil Sci.Soc.Am.J.60, 622—624.http://dx.doi.org/10.2136/sssaj1976.03615995004000040046x.

    Mutchler,C.K.,Moldenhauer,W.C.,1963.Applicator for laboratory rainfall simulator.Trans.Am.Soc.Agric.Eng.6,220—222.http://dx.doi.org/ 10.13031/2013.40871.

    O'Brien,J.S.,Jorgensen,G.R.,Garcia,R.,2009.FLO-2D Software Version 2009.FLO-2D Software,Inc,Nutrioso,AZ.

    Prevedello,C.L.,Loyola,J.M.T.,Reichardt,K.,Nielsen,D.R.,2009.New analytic solution related to the Richards,Philip,and Green-Ampt equations for inf i ltration.Vadose Zone J.8(1),127—135.http://dx.doi.org/10.2136/ vzj2008.0091.

    Rawls,W.J.,Brakensiek,D.L.,Miller,N.,1983.Green-Ampt inf i ltration parameters from soils data.J.Hydraul.Eng.109(1),62—69.http://dx.doi.org/ 10.1061/(ASCE)0733-9429(1983)109:1(62).

    Rawls,W.J.,Ahuja,L.R.,Brakensiek,D.L.,1992.Estimating soil hydraulic properties from soils data.In:Proceedings of the International Workshop on Indirect Methods for Estimating Hydraulic Properties of Unsaturated Soils.University of California,Riverside,pp.329—340.

    Regalado,C.M.,Ritter,A.,Alvarez-Benedi,J.,Munoz-Carpena,R.,2005. Simplified method to estimate the Green-Ampt wetting front suction and soil sorptivity with the Philip-Dunne falling-head permeameter.Vadose Zone J.4(2),291—299.

    Reynolds,W.D.,2010.Measuring soil hydraulic properties using a cased borehole permeameter:Steady flow analyses.Vadose Zone J.9(3), 637—652.http://dx.doi.org/10.2136/vzj2009.0136.

    Risse,L.M.,Nearing,M.A.,Zhang,X.C.,1995.Variability in Green-Ampt effective hydraulic conductivity under fallow conditions.J.Hydrol. 169(1—4),1—24.http://dx.doi.org/10.1016/0022-1694(94)02676-3.

    Saxton,K.E.,Rawls,W.J.,2006.Soil water characteristic estimates by texture and organic matter for hydrologic solutions.Soil Sci.Soc.Am.J.70(5), 1569—1578.http://dx.doi.org/10.2136/sssaj2005.0117.

    Silburn,D.M.,Connolly,R.D.,1995.Distributed parameter hydrology model (ANSWERS)applied to a range of catchment scales using rainfall simulator data I:inf i ltration modelling and parameter measurement.J.Hydrol. 172(1—4),87—104.http://dx.doi.org/10.1016/0022-1694(95)02740-G.

    Sivakumar,B.,2015.Networks:A generic theory for hydrology?Stoch.Environ.Res.Risk Assess.29(3),761—771.http://dx.doi.org/10.1007/ s00477-014-0902-7.

    Suleiman,K.A.,Swartzendruber,D.,2003.Measurement of sated hydraulic conductivity of surface soil in the field with a small-plot sprinkling inf i ltrometer.J.Hydrol.272(1—4),203—212.http://dx.doi.org/10.1016/S0022-1694(02)00265-2.

    Taskinen,A.,Sirvi¨o,H.,Bruen,M.,2008.Modelling effects of spatial variability of saturated hydraulic conductivity on autocorrelated overland flow data:Linear mixed model approach.Stoch.Environ.Res.Risk Assess. 22(1),67—82.http://dx.doi.org/10.1007/s00477-006-0099-5.

    Valiantzas,J.D.,2010.New linearized two-parameter inf i ltration equation for direct determination of conductivity and sorptivity.J.Hydrol.384(1—2), 1—13.http://dx.doi.org/10.1016/j.jhydrol.2009.12.049.

    van den Putte,A.,Govers,G.,Leys,A.,Langhans,C.,Clymans,W.,Diels,J., 2013.Estimating the parameters of the Green-Ampt inf i ltration equation from rainfall simulation data:Why simpler is better.J.Hydrol.476, 332—344.http://dx.doi.org/10.1016/j.jhydrol.2012.10.051.

    vanGenuchten,M.T.,1980.Aclosed-formequationforpredictingthehydraulic conductivity of unsaturated soils.Soil Sci.Soc.Am.J.44(5),892—898.

    Vˊazquez,E.V.,Miranda,J.G.V.,Gonzˊalez,A.P.,2005.Characterizing anisotropy and heterogeneity of soil surface microtopography using fractal models.Ecol.Model.182(3—4),337—353.http://dx.doi.org/10.1016/ j.ecolmodel.2004.04.012.

    Verbist,K.,Torfs,S.,Cornelis,W.M.,Oyarzun,R.,Soto,G.,Gabriels,D., 2010.Comparison of single-and double-ring inf i ltrometer methods on stony soils.Vadose Zone J.9(2),462—475.http://dx.doi.org/10.2136/ vzj2009.0058.

    Wang,L.L.,Li,Z.J.,Bao,H.J.,2010.Development and comparison of Gridbased distributed hydrological models for excess-inf i ltration runoffs.J. HohaiUniv.Nat.Sci.38(2),123—128.http://dx.doi.org/10.3876/ j.issn.1000—1980.2010.02.001(in Chinese).

    Received 12 October 2015;accepted 9 May 2016

    This work was supported by the National Natural Science Foundation of China(Grants No.51309078 and 51349015),the National Technology Support Program in the 12th Five-Year Plan of China(Grant No.2012BAK10B04),the Fundamental Research Funds for the Central Universities,the Program of Dual Innovative Talents Plan and Innovative Research Team in Jiangsu Province, and the Research on Spatio-Temporal Variable Source Runoff Model Based on Geomorphic Hydrological Response Units and Demonstration Application (Grant No.SHZH-IWHR-73).

    *Corresponding author.

    E-mail address:xianglonghhu@gmail.com(Long Xiang).

    Peer review under responsibility of Hohai University.

    http://dx.doi.org/10.1016/j.wse.2016.05.001

    1674-2370/?2016 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).

    ?2016 Hohai University.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http:// creativecommons.org/licenses/by-nc-nd/4.0/).

    一进一出抽搐gif免费好疼| 国产精华一区二区三区| 欧美三级亚洲精品| 国产成人a区在线观看| 简卡轻食公司| 日本av手机在线免费观看| 国产精品日韩av在线免费观看| 精品久久久久久久久亚洲| 日本免费一区二区三区高清不卡| 国产伦一二天堂av在线观看| 一本精品99久久精品77| 亚洲婷婷狠狠爱综合网| 免费观看人在逋| 我要看日韩黄色一级片| 九九热线精品视视频播放| 天美传媒精品一区二区| av女优亚洲男人天堂| 国产精品伦人一区二区| 亚洲人成网站在线观看播放| 成人亚洲欧美一区二区av| 国产精品国产三级国产av玫瑰| eeuss影院久久| 亚洲内射少妇av| 狠狠狠狠99中文字幕| 国产老妇伦熟女老妇高清| 成年免费大片在线观看| 精品日产1卡2卡| 特大巨黑吊av在线直播| 日韩成人伦理影院| 国产老妇女一区| 国产乱人视频| 丰满的人妻完整版| 最近视频中文字幕2019在线8| 亚洲国产精品成人综合色| 最近手机中文字幕大全| 久久草成人影院| 成人亚洲精品av一区二区| 中文字幕人妻熟人妻熟丝袜美| 国产综合懂色| 搞女人的毛片| 午夜福利高清视频| www日本黄色视频网| 深夜精品福利| 一夜夜www| 午夜免费激情av| 最近2019中文字幕mv第一页| 中文在线观看免费www的网站| 免费一级毛片在线播放高清视频| 看免费成人av毛片| 99国产极品粉嫩在线观看| 高清日韩中文字幕在线| 久久6这里有精品| 国产淫片久久久久久久久| 男人狂女人下面高潮的视频| 国产午夜精品一二区理论片| 国产精品一区二区性色av| eeuss影院久久| 九九热线精品视视频播放| 能在线免费看毛片的网站| 国产视频内射| 亚洲人成网站在线观看播放| 国产爱豆传媒在线观看| 欧美精品国产亚洲| 插阴视频在线观看视频| 欧美高清性xxxxhd video| 18+在线观看网站| 你懂的网址亚洲精品在线观看 | 91精品一卡2卡3卡4卡| 黄色一级大片看看| 97人妻精品一区二区三区麻豆| 丰满人妻一区二区三区视频av| 亚洲精品国产成人久久av| 免费一级毛片在线播放高清视频| 大香蕉久久网| 久久精品久久久久久噜噜老黄 | 久久精品夜夜夜夜夜久久蜜豆| 国产国拍精品亚洲av在线观看| 99九九线精品视频在线观看视频| 国产精华一区二区三区| 国产精品久久电影中文字幕| 亚洲成av人片在线播放无| 成年女人看的毛片在线观看| 国产精品av视频在线免费观看| 亚洲精品粉嫩美女一区| 成人性生交大片免费视频hd| 精品久久久久久久久久久久久| 欧美高清成人免费视频www| 精品熟女少妇av免费看| 丰满人妻一区二区三区视频av| 综合色丁香网| 美女高潮的动态| 久久久久久久久中文| 亚洲经典国产精华液单| 久久人人精品亚洲av| 黑人高潮一二区| 欧美在线一区亚洲| 国产精品一二三区在线看| 中文在线观看免费www的网站| 欧美+日韩+精品| 白带黄色成豆腐渣| 少妇丰满av| 成人漫画全彩无遮挡| 亚洲国产精品成人久久小说 | 九九在线视频观看精品| 菩萨蛮人人尽说江南好唐韦庄 | 久久99蜜桃精品久久| av国产免费在线观看| 丰满乱子伦码专区| 亚洲av免费高清在线观看| 亚洲不卡免费看| 尾随美女入室| 最好的美女福利视频网| 丰满乱子伦码专区| 韩国av在线不卡| 中文字幕久久专区| 国产老妇女一区| avwww免费| 国产亚洲91精品色在线| 最后的刺客免费高清国语| 99热这里只有精品一区| 欧美+日韩+精品| 国产亚洲5aaaaa淫片| 麻豆国产97在线/欧美| 亚洲四区av| 国内久久婷婷六月综合欲色啪| 国产伦理片在线播放av一区 | 国产国拍精品亚洲av在线观看| 99热网站在线观看| 成年女人永久免费观看视频| 日日摸夜夜添夜夜爱| 在线观看午夜福利视频| 久久精品国产清高在天天线| 精品免费久久久久久久清纯| 国产精品爽爽va在线观看网站| 午夜老司机福利剧场| 成人二区视频| 国产精品不卡视频一区二区| 男的添女的下面高潮视频| 丰满人妻一区二区三区视频av| 日本av手机在线免费观看| 日韩一区二区三区影片| 日韩成人伦理影院| 国产av在哪里看| 亚洲精品色激情综合| 一边摸一边抽搐一进一小说| 国产精品一区www在线观看| 成人综合一区亚洲| 22中文网久久字幕| 亚洲精品国产成人久久av| 国产麻豆成人av免费视频| 中文在线观看免费www的网站| 欧美不卡视频在线免费观看| www日本黄色视频网| 有码 亚洲区| 99热6这里只有精品| 国产精品久久久久久久久免| 国产高清三级在线| 精品久久国产蜜桃| 亚洲五月天丁香| 国产亚洲欧美98| 午夜久久久久精精品| 国产午夜精品论理片| 亚洲真实伦在线观看| 国产视频首页在线观看| 久久99热6这里只有精品| 麻豆一二三区av精品| 国产午夜福利久久久久久| 久久精品国产亚洲av香蕉五月| 听说在线观看完整版免费高清| 禁无遮挡网站| 国产精品久久电影中文字幕| 国内精品一区二区在线观看| 亚洲综合色惰| 1024手机看黄色片| 亚洲成av人片在线播放无| 日产精品乱码卡一卡2卡三| 久久久久久大精品| 日韩 亚洲 欧美在线| av专区在线播放| 欧美精品一区二区大全| 深夜精品福利| 久久九九热精品免费| 久久精品国产亚洲av天美| 久久久久久久久久久免费av| 国产高清激情床上av| 一个人免费在线观看电影| 美女高潮的动态| 欧美激情国产日韩精品一区| 国产 一区精品| 亚洲熟妇中文字幕五十中出| 国产精品久久视频播放| 亚洲图色成人| av在线蜜桃| 国产亚洲91精品色在线| 亚洲成a人片在线一区二区| a级一级毛片免费在线观看| 两个人的视频大全免费| a级毛片免费高清观看在线播放| 99久国产av精品国产电影| 日本三级黄在线观看| 特级一级黄色大片| 国产蜜桃级精品一区二区三区| av在线播放精品| 老司机福利观看| 一个人看的www免费观看视频| 成人三级黄色视频| 狠狠狠狠99中文字幕| 热99在线观看视频| 变态另类丝袜制服| 国国产精品蜜臀av免费| 91精品国产九色| 少妇高潮的动态图| 日韩制服骚丝袜av| 日韩亚洲欧美综合| 日韩成人伦理影院| 精品国产三级普通话版| 欧美xxxx黑人xx丫x性爽| 又爽又黄a免费视频| 亚洲电影在线观看av| 亚洲四区av| 大又大粗又爽又黄少妇毛片口| 高清午夜精品一区二区三区 | 精品欧美国产一区二区三| 日本撒尿小便嘘嘘汇集6| 97超碰精品成人国产| 蜜臀久久99精品久久宅男| 免费观看人在逋| 精品99又大又爽又粗少妇毛片| 欧美潮喷喷水| 国产高清不卡午夜福利| 国产精品电影一区二区三区| 国产成人a区在线观看| 99久久人妻综合| 伦精品一区二区三区| 91午夜精品亚洲一区二区三区| 蜜臀久久99精品久久宅男| 亚洲最大成人av| 男插女下体视频免费在线播放| 久久国内精品自在自线图片| 桃色一区二区三区在线观看| 久久99热这里只有精品18| 久久99蜜桃精品久久| 国产亚洲91精品色在线| 少妇猛男粗大的猛烈进出视频 | 亚洲人成网站高清观看| 欧美人与善性xxx| 不卡一级毛片| 2021天堂中文幕一二区在线观| 丝袜美腿在线中文| 日本撒尿小便嘘嘘汇集6| 欧美日韩国产亚洲二区| 村上凉子中文字幕在线| 天堂√8在线中文| 毛片一级片免费看久久久久| 欧美丝袜亚洲另类| 精品久久久久久久人妻蜜臀av| 免费观看在线日韩| 色综合亚洲欧美另类图片| 日本av手机在线免费观看| 蜜桃亚洲精品一区二区三区| 黄色欧美视频在线观看| 熟妇人妻久久中文字幕3abv| 一级黄色大片毛片| 中文字幕av成人在线电影| 日韩中字成人| 亚洲欧洲日产国产| 亚洲精品国产av成人精品| 午夜精品国产一区二区电影 | 欧美最黄视频在线播放免费| av在线亚洲专区| 少妇高潮的动态图| 美女脱内裤让男人舔精品视频 | 亚洲第一区二区三区不卡| 国产成人精品久久久久久| 久久国产乱子免费精品| 国产成年人精品一区二区| 国产一区二区激情短视频| 一个人免费在线观看电影| 中文字幕制服av| 国产精品国产高清国产av| 你懂的网址亚洲精品在线观看 | 赤兔流量卡办理| 你懂的网址亚洲精品在线观看 | 日韩欧美国产在线观看| a级一级毛片免费在线观看| 中文亚洲av片在线观看爽| 亚洲欧美精品专区久久| 亚洲一级一片aⅴ在线观看| 偷拍熟女少妇极品色| 久久国内精品自在自线图片| 午夜福利成人在线免费观看| 国产大屁股一区二区在线视频| 黄色欧美视频在线观看| 日韩成人伦理影院| 成人毛片60女人毛片免费| www.色视频.com| 成人漫画全彩无遮挡| 国产日本99.免费观看| 久久精品久久久久久久性| 99热这里只有是精品在线观看| 校园春色视频在线观看| av在线老鸭窝| 人妻制服诱惑在线中文字幕| 大香蕉久久网| 日韩一本色道免费dvd| 一区福利在线观看| 久久久久久国产a免费观看| 深爱激情五月婷婷| 91久久精品国产一区二区成人| 欧美在线一区亚洲| 午夜激情欧美在线| 亚洲久久久久久中文字幕| 久久久精品大字幕| 国产av麻豆久久久久久久| 天堂影院成人在线观看| 国产久久久一区二区三区| 亚洲国产精品国产精品| 亚洲不卡免费看| 久久久久久久久久久丰满| 亚洲成av人片在线播放无| 久久精品91蜜桃| 欧美精品一区二区大全| 神马国产精品三级电影在线观看| 国产亚洲精品久久久久久毛片| 特级一级黄色大片| 美女脱内裤让男人舔精品视频 | 亚洲av一区综合| 一夜夜www| 亚洲国产精品成人久久小说 | 欧美成人免费av一区二区三区| 精品久久久久久久末码| 春色校园在线视频观看| 国产精品一二三区在线看| 成人欧美大片| 国产精品久久久久久精品电影小说 | 免费av毛片视频| 五月伊人婷婷丁香| 中文字幕熟女人妻在线| 狠狠狠狠99中文字幕| 免费人成视频x8x8入口观看| 久久韩国三级中文字幕| 狂野欧美白嫩少妇大欣赏| 欧美成人精品欧美一级黄| 大又大粗又爽又黄少妇毛片口| 久久草成人影院| 国产亚洲av片在线观看秒播厂 | 嫩草影院精品99| 日本爱情动作片www.在线观看| 亚洲av免费高清在线观看| 最近2019中文字幕mv第一页| 成人美女网站在线观看视频| 欧美不卡视频在线免费观看| 国产精品国产高清国产av| 国产午夜福利久久久久久| 亚洲国产高清在线一区二区三| 亚洲成人精品中文字幕电影| 免费观看人在逋| 精品人妻一区二区三区麻豆| 国产成人精品久久久久久| 亚洲国产欧洲综合997久久,| 男插女下体视频免费在线播放| 国产大屁股一区二区在线视频| 久久久久久久久久久免费av| 国产三级在线视频| 欧美性猛交╳xxx乱大交人| 人妻系列 视频| 亚洲国产精品合色在线| 国产极品天堂在线| 国产黄色小视频在线观看| 免费看日本二区| 国产精品.久久久| 日日撸夜夜添| 51国产日韩欧美| 国产探花极品一区二区| 2021天堂中文幕一二区在线观| 九色成人免费人妻av| 亚洲人与动物交配视频| 中国美白少妇内射xxxbb| 91午夜精品亚洲一区二区三区| 久久久久久国产a免费观看| 久久精品国产亚洲网站| 人人妻人人看人人澡| 丰满人妻一区二区三区视频av| 99热全是精品| 小说图片视频综合网站| 亚洲av第一区精品v没综合| 乱码一卡2卡4卡精品| 91av网一区二区| 午夜福利在线观看吧| 老女人水多毛片| 在线观看av片永久免费下载| 一级毛片电影观看 | www日本黄色视频网| 在线观看午夜福利视频| 九九爱精品视频在线观看| 亚洲经典国产精华液单| 美女黄网站色视频| 99久久无色码亚洲精品果冻| 校园春色视频在线观看| 国产伦理片在线播放av一区 | 免费观看人在逋| 亚洲av电影不卡..在线观看| 成人亚洲欧美一区二区av| 久久久久久久久久黄片| 国产黄色视频一区二区在线观看 | 午夜久久久久精精品| 老熟妇乱子伦视频在线观看| 草草在线视频免费看| 春色校园在线视频观看| 99久久精品一区二区三区| 亚洲精品日韩av片在线观看| 男人的好看免费观看在线视频| 国产成人精品一,二区 | 中国美女看黄片| 亚洲在线自拍视频| 天堂中文最新版在线下载 | 国产乱人视频| 噜噜噜噜噜久久久久久91| 桃色一区二区三区在线观看| 久久久久久久久久成人| 国产不卡一卡二| 日本av手机在线免费观看| 美女 人体艺术 gogo| 中文字幕久久专区| 黄色日韩在线| 国产女主播在线喷水免费视频网站 | 能在线免费看毛片的网站| 大又大粗又爽又黄少妇毛片口| 身体一侧抽搐| 哪里可以看免费的av片| 久久人人爽人人片av| 欧美+日韩+精品| 深夜精品福利| 午夜激情福利司机影院| 不卡视频在线观看欧美| 精品久久久久久久久av| 欧美极品一区二区三区四区| av.在线天堂| 日本免费a在线| 国产av一区在线观看免费| 国产免费一级a男人的天堂| 精品免费久久久久久久清纯| 精品国内亚洲2022精品成人| 欧美成人a在线观看| 美女高潮的动态| av福利片在线观看| 亚洲国产欧美人成| 中文字幕免费在线视频6| 成年女人看的毛片在线观看| 在线观看美女被高潮喷水网站| 91精品国产九色| 一级黄色大片毛片| 成人特级av手机在线观看| 国产乱人偷精品视频| 国产乱人视频| 毛片一级片免费看久久久久| 真实男女啪啪啪动态图| eeuss影院久久| 一区二区三区免费毛片| a级一级毛片免费在线观看| 22中文网久久字幕| 亚洲精品影视一区二区三区av| 在线播放国产精品三级| 日本-黄色视频高清免费观看| 国产黄a三级三级三级人| 国产v大片淫在线免费观看| 日本黄大片高清| 最近中文字幕高清免费大全6| 老女人水多毛片| 免费观看a级毛片全部| 一级毛片我不卡| a级一级毛片免费在线观看| 菩萨蛮人人尽说江南好唐韦庄 | 亚洲欧美成人综合另类久久久 | 综合色丁香网| 午夜久久久久精精品| av天堂在线播放| 两个人视频免费观看高清| 老女人水多毛片| 日韩精品有码人妻一区| 国产片特级美女逼逼视频| 免费无遮挡裸体视频| 高清日韩中文字幕在线| 桃色一区二区三区在线观看| 黄色日韩在线| 久久综合国产亚洲精品| 免费人成在线观看视频色| 一区福利在线观看| 日韩一区二区三区影片| 免费av不卡在线播放| 日韩强制内射视频| 亚洲在线自拍视频| 日本色播在线视频| 久久人人爽人人片av| 亚洲人成网站在线播放欧美日韩| 97在线视频观看| 国产久久久一区二区三区| 欧美日本亚洲视频在线播放| 亚洲va在线va天堂va国产| 伦精品一区二区三区| 久久久久久久午夜电影| 99久国产av精品| av在线观看视频网站免费| 亚洲四区av| 在线a可以看的网站| 嫩草影院精品99| av黄色大香蕉| 久久精品国产自在天天线| 久久人人爽人人片av| 国产精品一区二区性色av| 午夜福利在线观看吧| 日韩,欧美,国产一区二区三区 | 亚洲欧美精品专区久久| 日韩高清综合在线| 欧美人与善性xxx| a级毛片a级免费在线| 啦啦啦观看免费观看视频高清| 国产高清激情床上av| 久久精品久久久久久久性| 熟女电影av网| 好男人视频免费观看在线| 国产精品99久久久久久久久| 99热只有精品国产| 亚洲av电影不卡..在线观看| 波多野结衣高清无吗| 亚洲欧美清纯卡通| 国产探花极品一区二区| 成人漫画全彩无遮挡| 一级毛片aaaaaa免费看小| 日韩欧美一区二区三区在线观看| 麻豆乱淫一区二区| 又粗又爽又猛毛片免费看| 能在线免费观看的黄片| 小说图片视频综合网站| 少妇裸体淫交视频免费看高清| 久久这里有精品视频免费| 色综合亚洲欧美另类图片| 麻豆乱淫一区二区| 天天躁夜夜躁狠狠久久av| 男人的好看免费观看在线视频| 国产午夜精品论理片| 麻豆国产97在线/欧美| 国产大屁股一区二区在线视频| 男女下面进入的视频免费午夜| 美女 人体艺术 gogo| 夫妻性生交免费视频一级片| 观看美女的网站| 久久这里有精品视频免费| 欧美精品国产亚洲| 国产 一区 欧美 日韩| 久久婷婷人人爽人人干人人爱| 99热全是精品| 午夜福利高清视频| 欧美性猛交╳xxx乱大交人| 色5月婷婷丁香| 欧美性猛交╳xxx乱大交人| 亚洲av中文字字幕乱码综合| 国产国拍精品亚洲av在线观看| 直男gayav资源| 97人妻精品一区二区三区麻豆| 国产黄色视频一区二区在线观看 | 91精品一卡2卡3卡4卡| 97超视频在线观看视频| 高清在线视频一区二区三区 | 波多野结衣巨乳人妻| 黑人高潮一二区| 成人美女网站在线观看视频| 久久久久久久久久黄片| 亚洲国产精品国产精品| 日韩制服骚丝袜av| 亚洲国产精品国产精品| 国产精品一区二区三区四区久久| 亚洲在线观看片| 久久久久久国产a免费观看| 国产激情偷乱视频一区二区| 色哟哟·www| 国产av麻豆久久久久久久| 精品久久久久久久久av| 啦啦啦韩国在线观看视频| 男女做爰动态图高潮gif福利片| 久久久国产成人精品二区| 99国产极品粉嫩在线观看| 一进一出抽搐动态| 久久久久久九九精品二区国产| 99久久精品国产国产毛片| 身体一侧抽搐| 一区二区三区免费毛片| 老熟妇乱子伦视频在线观看| 综合色丁香网| 97超碰精品成人国产| 免费观看人在逋| 国产精品99久久久久久久久| 国产高清有码在线观看视频| 内射极品少妇av片p| 91在线精品国自产拍蜜月| h日本视频在线播放| 最近手机中文字幕大全| 悠悠久久av| 久久精品国产亚洲网站| 日韩欧美精品v在线| 国产精品久久电影中文字幕| 男人舔女人下体高潮全视频| 国产v大片淫在线免费观看| 亚洲美女视频黄频| 欧美zozozo另类| .国产精品久久| 久久国内精品自在自线图片| 两个人的视频大全免费| 夜夜夜夜夜久久久久| av免费观看日本| 亚洲不卡免费看| 日韩三级伦理在线观看| 欧美潮喷喷水| 免费观看人在逋| 国语自产精品视频在线第100页|