• <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/).

    亚洲欧洲精品一区二区精品久久久| 国产激情偷乱视频一区二区| 成人亚洲精品av一区二区| 午夜精品在线福利| 色精品久久人妻99蜜桃| 精品久久久久久久毛片微露脸| 久久精品国产清高在天天线| 成人三级黄色视频| 久久久久久免费高清国产稀缺| 亚洲一区二区三区不卡视频| 制服人妻中文乱码| 熟女少妇亚洲综合色aaa.| 99精品久久久久人妻精品| a在线观看视频网站| 人成视频在线观看免费观看| 国产高清有码在线观看视频 | 久久久久久大精品| 国产高清视频在线播放一区| 好看av亚洲va欧美ⅴa在| 久久久久国产一级毛片高清牌| 丝袜人妻中文字幕| 欧美一区二区国产精品久久精品 | 欧美日韩中文字幕国产精品一区二区三区| 无遮挡黄片免费观看| ponron亚洲| 国产精品久久久久久久电影 | 一个人观看的视频www高清免费观看 | 亚洲国产精品成人综合色| 欧美日韩中文字幕国产精品一区二区三区| 国产精品,欧美在线| 国产男靠女视频免费网站| 一区福利在线观看| 亚洲国产精品999在线| 亚洲18禁久久av| 国产精品爽爽va在线观看网站| 久久久久久久精品吃奶| 男女之事视频高清在线观看| 国产熟女午夜一区二区三区| 校园春色视频在线观看| 99精品在免费线老司机午夜| 特级一级黄色大片| 色精品久久人妻99蜜桃| 麻豆久久精品国产亚洲av| 老司机靠b影院| 美女免费视频网站| 欧美性猛交黑人性爽| 成人午夜高清在线视频| 欧美乱妇无乱码| 久久久久性生活片| 午夜久久久久精精品| 国内精品久久久久精免费| 久久久久免费精品人妻一区二区| www日本在线高清视频| 国内少妇人妻偷人精品xxx网站 | 国产精品乱码一区二三区的特点| 欧美成人性av电影在线观看| 俄罗斯特黄特色一大片| 久久久久久久精品吃奶| 国产黄色小视频在线观看| 日本免费a在线| 岛国在线免费视频观看| 丝袜美腿诱惑在线| 制服丝袜大香蕉在线| 成在线人永久免费视频| 国产精品精品国产色婷婷| 亚洲av第一区精品v没综合| 欧美在线一区亚洲| 久久中文字幕人妻熟女| 国内毛片毛片毛片毛片毛片| 99热只有精品国产| 久久亚洲真实| 少妇的丰满在线观看| 亚洲av中文字字幕乱码综合| 亚洲人成网站在线播放欧美日韩| 18禁黄网站禁片免费观看直播| 亚洲av成人av| 国产精品 国内视频| 老司机午夜十八禁免费视频| 久久久久久久久免费视频了| 国产精品乱码一区二三区的特点| 嫩草影视91久久| 国产又色又爽无遮挡免费看| 亚洲熟妇熟女久久| 国产黄色小视频在线观看| 欧美成狂野欧美在线观看| 欧美黄色淫秽网站| 观看免费一级毛片| 亚洲专区字幕在线| 日本黄大片高清| 一a级毛片在线观看| 亚洲国产欧洲综合997久久,| 在线观看美女被高潮喷水网站 | av在线天堂中文字幕| 久久精品91无色码中文字幕| 非洲黑人性xxxx精品又粗又长| 国产片内射在线| 一本综合久久免费| 国产又色又爽无遮挡免费看| tocl精华| 国产免费男女视频| 日韩欧美三级三区| 国产亚洲欧美98| 在线十欧美十亚洲十日本专区| 久久这里只有精品中国| 免费在线观看黄色视频的| 亚洲国产欧美网| 久久久国产欧美日韩av| 日韩高清综合在线| 91九色精品人成在线观看| 国产成人系列免费观看| 成年人黄色毛片网站| 丁香六月欧美| av在线播放免费不卡| 99riav亚洲国产免费| 国产探花在线观看一区二区| 好男人电影高清在线观看| 丝袜美腿诱惑在线| 两个人看的免费小视频| 男女午夜视频在线观看| 日韩大码丰满熟妇| 精品高清国产在线一区| 国产精品av视频在线免费观看| 亚洲五月婷婷丁香| 一级毛片高清免费大全| 日本 av在线| 黄片大片在线免费观看| 波多野结衣巨乳人妻| 欧美av亚洲av综合av国产av| 成人手机av| 精品久久久久久久人妻蜜臀av| av中文乱码字幕在线| 欧美zozozo另类| 日韩欧美在线二视频| 亚洲 欧美 日韩 在线 免费| 制服丝袜大香蕉在线| 欧美中文综合在线视频| 又粗又爽又猛毛片免费看| 性色av乱码一区二区三区2| 国产精品 欧美亚洲| 欧美三级亚洲精品| 51午夜福利影视在线观看| svipshipincom国产片| 国产99白浆流出| 亚洲av电影在线进入| 亚洲欧美精品综合久久99| 日本成人三级电影网站| 极品教师在线免费播放| 日本在线视频免费播放| 精品人妻1区二区| 夜夜爽天天搞| 亚洲午夜理论影院| 黄片小视频在线播放| 日本一区二区免费在线视频| 老熟妇乱子伦视频在线观看| 午夜福利免费观看在线| 久久人妻av系列| www.自偷自拍.com| 欧美zozozo另类| 亚洲在线自拍视频| 日韩精品青青久久久久久| 欧美成人性av电影在线观看| 91老司机精品| 久久欧美精品欧美久久欧美| 一本久久中文字幕| 欧美精品啪啪一区二区三区| 国产三级黄色录像| 伦理电影免费视频| 国产av在哪里看| 亚洲五月天丁香| 一级毛片女人18水好多| 成人国产一区最新在线观看| 久久精品国产清高在天天线| 制服丝袜大香蕉在线| 国产欧美日韩一区二区精品| 精品日产1卡2卡| 夜夜看夜夜爽夜夜摸| 亚洲av电影在线进入| 国产激情久久老熟女| 老司机靠b影院| 男女之事视频高清在线观看| 99热这里只有是精品50| 国产成人系列免费观看| 好男人在线观看高清免费视频| 看黄色毛片网站| 一区福利在线观看| 欧美日韩一级在线毛片| 亚洲欧美精品综合一区二区三区| 男女床上黄色一级片免费看| 欧美3d第一页| 免费在线观看影片大全网站| av超薄肉色丝袜交足视频| 国产探花在线观看一区二区| 两个人免费观看高清视频| 国产真实乱freesex| 免费在线观看视频国产中文字幕亚洲| 国产欧美日韩精品亚洲av| 欧美中文日本在线观看视频| 成人av一区二区三区在线看| 久久精品成人免费网站| 欧美 亚洲 国产 日韩一| 精品乱码久久久久久99久播| 欧美黄色淫秽网站| 亚洲精品在线美女| av超薄肉色丝袜交足视频| 看片在线看免费视频| 老鸭窝网址在线观看| 欧美最黄视频在线播放免费| 国产精品av视频在线免费观看| 女警被强在线播放| 国产精品98久久久久久宅男小说| 免费观看人在逋| 国产精品一区二区免费欧美| 国内揄拍国产精品人妻在线| av片东京热男人的天堂| 老司机靠b影院| 亚洲精品中文字幕一二三四区| 亚洲人成网站高清观看| 男女做爰动态图高潮gif福利片| 99久久国产精品久久久| 色哟哟哟哟哟哟| 亚洲av熟女| 日日爽夜夜爽网站| 亚洲av美国av| 精品久久蜜臀av无| 妹子高潮喷水视频| 国产成人精品无人区| 免费一级毛片在线播放高清视频| 美女高潮喷水抽搐中文字幕| 亚洲男人的天堂狠狠| 欧美一区二区精品小视频在线| 亚洲国产精品成人综合色| 成人精品一区二区免费| 亚洲国产欧洲综合997久久,| 丰满人妻熟妇乱又伦精品不卡| 亚洲真实伦在线观看| 午夜福利视频1000在线观看| 亚洲一区中文字幕在线| 亚洲九九香蕉| 精品高清国产在线一区| 免费高清视频大片| 欧美大码av| 午夜福利18| 亚洲欧美日韩高清专用| 变态另类成人亚洲欧美熟女| 欧美成人性av电影在线观看| 国产亚洲精品一区二区www| 丰满人妻一区二区三区视频av | av在线天堂中文字幕| 夜夜躁狠狠躁天天躁| 日韩 欧美 亚洲 中文字幕| 午夜精品久久久久久毛片777| 亚洲国产精品合色在线| 欧美性猛交╳xxx乱大交人| 国产在线精品亚洲第一网站| 国产欧美日韩一区二区三| 夜夜夜夜夜久久久久| 日本黄大片高清| 欧美黑人巨大hd| 国产亚洲精品一区二区www| 一二三四社区在线视频社区8| www.精华液| 国产精品影院久久| 国产麻豆成人av免费视频| 亚洲中文字幕日韩| 母亲3免费完整高清在线观看| 中文资源天堂在线| 亚洲av电影在线进入| 久久久精品大字幕| 丰满人妻一区二区三区视频av | 五月玫瑰六月丁香| 久久午夜综合久久蜜桃| 色综合亚洲欧美另类图片| 成人国语在线视频| 欧美日本亚洲视频在线播放| 中文字幕精品亚洲无线码一区| 成人特级黄色片久久久久久久| 18禁美女被吸乳视频| 欧美日韩中文字幕国产精品一区二区三区| 午夜久久久久精精品| 美女免费视频网站| 免费看日本二区| 午夜福利高清视频| 久久九九热精品免费| 国产黄片美女视频| 国产亚洲精品久久久久5区| 看片在线看免费视频| 亚洲av美国av| 久久热在线av| 久久精品aⅴ一区二区三区四区| 成年免费大片在线观看| 国产精品影院久久| x7x7x7水蜜桃| 每晚都被弄得嗷嗷叫到高潮| 一边摸一边抽搐一进一小说| 国产精品免费视频内射| 久久亚洲真实| 欧美精品亚洲一区二区| 亚洲五月天丁香| 亚洲欧美日韩东京热| 久久香蕉国产精品| 国产av麻豆久久久久久久| 久久久久久久久久黄片| 亚洲天堂国产精品一区在线| 精品国产美女av久久久久小说| 中出人妻视频一区二区| 精华霜和精华液先用哪个| 老司机午夜十八禁免费视频| 人人妻,人人澡人人爽秒播| av福利片在线| 午夜视频精品福利| 在线观看免费午夜福利视频| 欧美三级亚洲精品| 两性夫妻黄色片| 国产亚洲av高清不卡| 免费人成视频x8x8入口观看| 欧美日韩精品网址| 99久久无色码亚洲精品果冻| 国产黄片美女视频| 精品久久久久久久末码| 三级男女做爰猛烈吃奶摸视频| 久久久久久久精品吃奶| 老司机在亚洲福利影院| 最近最新免费中文字幕在线| av福利片在线观看| 中国美女看黄片| 国产精品国产高清国产av| 制服丝袜大香蕉在线| 欧美日本视频| 在线观看美女被高潮喷水网站 | 19禁男女啪啪无遮挡网站| 99热这里只有是精品50| 亚洲精品美女久久久久99蜜臀| 国产日本99.免费观看| 手机成人av网站| 久久精品国产亚洲av香蕉五月| 在线观看免费日韩欧美大片| 一本大道久久a久久精品| 嫩草影院精品99| 国产精品野战在线观看| 久久精品综合一区二区三区| 精品高清国产在线一区| 舔av片在线| 亚洲中文日韩欧美视频| 亚洲性夜色夜夜综合| 久久久久久久久中文| av有码第一页| 国产成人精品久久二区二区免费| 淫秽高清视频在线观看| 国产午夜精品论理片| 欧美性猛交╳xxx乱大交人| 亚洲av成人一区二区三| 99国产综合亚洲精品| 欧美日本视频| 一个人免费在线观看电影 | 精品久久久久久,| 97超级碰碰碰精品色视频在线观看| 神马国产精品三级电影在线观看 | av在线播放免费不卡| 91麻豆av在线| 亚洲精品久久成人aⅴ小说| 精品久久久久久成人av| 国产精品av视频在线免费观看| 蜜桃久久精品国产亚洲av| 午夜福利欧美成人| www.熟女人妻精品国产| 国产精品,欧美在线| 国产aⅴ精品一区二区三区波| 1024香蕉在线观看| 夜夜躁狠狠躁天天躁| 亚洲国产欧洲综合997久久,| 日日摸夜夜添夜夜添小说| 国产午夜精品久久久久久| 啦啦啦韩国在线观看视频| 国产黄a三级三级三级人| 欧美乱妇无乱码| 在线观看www视频免费| 久久性视频一级片| 国内少妇人妻偷人精品xxx网站 | 国产伦一二天堂av在线观看| av视频在线观看入口| 日韩精品中文字幕看吧| 一边摸一边抽搐一进一小说| 欧美日韩黄片免| 国产欧美日韩一区二区精品| 又粗又爽又猛毛片免费看| x7x7x7水蜜桃| 国产成人系列免费观看| 午夜成年电影在线免费观看| 18禁国产床啪视频网站| 中文字幕久久专区| 淫秽高清视频在线观看| 免费在线观看视频国产中文字幕亚洲| 色播亚洲综合网| 久热爱精品视频在线9| 久久精品成人免费网站| 我的老师免费观看完整版| 好男人在线观看高清免费视频| 12—13女人毛片做爰片一| 18禁黄网站禁片午夜丰满| 男人舔奶头视频| 日本撒尿小便嘘嘘汇集6| 女人高潮潮喷娇喘18禁视频| 人妻久久中文字幕网| 国产主播在线观看一区二区| 黑人欧美特级aaaaaa片| 国产亚洲欧美在线一区二区| 国产成+人综合+亚洲专区| 51午夜福利影视在线观看| 亚洲av五月六月丁香网| 嫁个100分男人电影在线观看| 国产一级毛片七仙女欲春2| 国产熟女午夜一区二区三区| 亚洲精品在线观看二区| 成熟少妇高潮喷水视频| 国产91精品成人一区二区三区| 91老司机精品| 国产精品爽爽va在线观看网站| 国产熟女xx| 最近最新中文字幕大全免费视频| 成人国产综合亚洲| 午夜免费激情av| 人妻夜夜爽99麻豆av| 日本一区二区免费在线视频| 欧美av亚洲av综合av国产av| 国产精品永久免费网站| 一二三四社区在线视频社区8| 天堂av国产一区二区熟女人妻 | 日日夜夜操网爽| 琪琪午夜伦伦电影理论片6080| 国产熟女xx| 99国产精品一区二区蜜桃av| 久久精品91蜜桃| 精品熟女少妇八av免费久了| 久久国产乱子伦精品免费另类| 熟女少妇亚洲综合色aaa.| 久久精品国产清高在天天线| 国产精品野战在线观看| 一个人免费在线观看电影 | 国产91精品成人一区二区三区| 精品国产乱子伦一区二区三区| 色综合站精品国产| 99国产精品一区二区三区| 最近最新免费中文字幕在线| 悠悠久久av| av片东京热男人的天堂| 亚洲美女视频黄频| 亚洲欧美日韩高清在线视频| 日韩三级视频一区二区三区| 九色国产91popny在线| 国产午夜精品久久久久久| 欧美乱码精品一区二区三区| 亚洲国产精品成人综合色| 国产亚洲精品久久久久久毛片| 99在线视频只有这里精品首页| 波多野结衣高清无吗| 老司机午夜福利在线观看视频| 国产免费av片在线观看野外av| 午夜精品一区二区三区免费看| or卡值多少钱| 老司机在亚洲福利影院| 麻豆国产av国片精品| 亚洲精品在线观看二区| 欧洲精品卡2卡3卡4卡5卡区| 天堂√8在线中文| 欧美日本视频| 色精品久久人妻99蜜桃| 午夜福利高清视频| 精品午夜福利视频在线观看一区| 精品乱码久久久久久99久播| netflix在线观看网站| 岛国在线免费视频观看| 国产99白浆流出| 久久久久性生活片| 99久久精品热视频| 黄色毛片三级朝国网站| 亚洲av熟女| 亚洲 欧美 日韩 在线 免费| 在线看三级毛片| avwww免费| 久久久精品大字幕| 免费观看人在逋| 小说图片视频综合网站| 国产麻豆成人av免费视频| 一级a爱片免费观看的视频| 中文在线观看免费www的网站 | 国内精品一区二区在线观看| 欧美乱妇无乱码| 男女午夜视频在线观看| 正在播放国产对白刺激| 91成年电影在线观看| 手机成人av网站| 国产欧美日韩一区二区精品| 国产成人系列免费观看| 久久久久久久久免费视频了| 国内毛片毛片毛片毛片毛片| 国产亚洲精品久久久久久毛片| 日本一区二区免费在线视频| av在线播放免费不卡| 亚洲欧美日韩高清在线视频| 午夜亚洲福利在线播放| 动漫黄色视频在线观看| 亚洲免费av在线视频| 欧美午夜高清在线| 亚洲午夜理论影院| 在线观看午夜福利视频| 久久久久亚洲av毛片大全| 欧美中文综合在线视频| 女生性感内裤真人,穿戴方法视频| 天天一区二区日本电影三级| 香蕉国产在线看| 无遮挡黄片免费观看| 91大片在线观看| 国产精品爽爽va在线观看网站| 手机成人av网站| 别揉我奶头~嗯~啊~动态视频| 欧美av亚洲av综合av国产av| 久久人妻av系列| 欧美性猛交╳xxx乱大交人| 日本成人三级电影网站| 在线看三级毛片| 久久天堂一区二区三区四区| 91在线观看av| 日韩高清综合在线| 2021天堂中文幕一二区在线观| 级片在线观看| 俄罗斯特黄特色一大片| 一区二区三区激情视频| 亚洲一区高清亚洲精品| 免费在线观看完整版高清| 91在线观看av| 黄片大片在线免费观看| 亚洲成人精品中文字幕电影| 亚洲av成人精品一区久久| 99国产综合亚洲精品| 亚洲成人免费电影在线观看| 日本撒尿小便嘘嘘汇集6| 91麻豆精品激情在线观看国产| 亚洲午夜理论影院| 久久久久九九精品影院| 人人妻人人澡欧美一区二区| 可以免费在线观看a视频的电影网站| 97人妻精品一区二区三区麻豆| 精品欧美一区二区三区在线| 国产三级在线视频| xxxwww97欧美| av中文乱码字幕在线| 日韩欧美国产在线观看| 久久中文字幕一级| 亚洲精品中文字幕在线视频| 色综合欧美亚洲国产小说| 国产黄a三级三级三级人| 亚洲av五月六月丁香网| 午夜老司机福利片| 久久精品国产清高在天天线| 国产伦一二天堂av在线观看| 国产av又大| 色噜噜av男人的天堂激情| 亚洲精品av麻豆狂野| 国产精品久久久久久人妻精品电影| 日本免费一区二区三区高清不卡| 色尼玛亚洲综合影院| 亚洲狠狠婷婷综合久久图片| 久久香蕉精品热| 国产成人欧美在线观看| 免费高清视频大片| 老司机在亚洲福利影院| 国内精品久久久久久久电影| 精品乱码久久久久久99久播| 国产精品乱码一区二三区的特点| 欧美日韩一级在线毛片| 国产精品久久久人人做人人爽| 在线观看舔阴道视频| 在线观看免费日韩欧美大片| 五月伊人婷婷丁香| 国产1区2区3区精品| 色精品久久人妻99蜜桃| 国产亚洲精品第一综合不卡| 国内精品久久久久精免费| 色播亚洲综合网| 婷婷精品国产亚洲av在线| 国产又黄又爽又无遮挡在线| 国产亚洲精品久久久久5区| 成人av一区二区三区在线看| 久久伊人香网站| 黄色视频,在线免费观看| 精品国产亚洲在线| 一卡2卡三卡四卡精品乱码亚洲| 亚洲无线在线观看| 精品国产亚洲在线| 校园春色视频在线观看| 久久久精品欧美日韩精品| 不卡一级毛片| 国产精品 欧美亚洲| 免费看美女性在线毛片视频| 少妇粗大呻吟视频| 高清在线国产一区| 国产伦人伦偷精品视频| 18禁观看日本| 欧美大码av| av视频在线观看入口| 欧美又色又爽又黄视频| 淫妇啪啪啪对白视频| 啪啪无遮挡十八禁网站| 两个人看的免费小视频| av视频在线观看入口| 久久欧美精品欧美久久欧美| 精品一区二区三区av网在线观看| 亚洲黑人精品在线| 免费看a级黄色片| 啦啦啦观看免费观看视频高清| 亚洲精华国产精华精| 最好的美女福利视频网| 日日夜夜操网爽| 高清毛片免费观看视频网站| 色噜噜av男人的天堂激情|