首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献，搜索用时 15 毫秒
1.
A physically based inverse method is developed using hybrid formulation and coordinate transform to simultaneously estimate hydraulic conductivity tensors, steady‐state flow field, and boundary conditions for a confined aquifer under ambient flow or pumping condition. Unlike existing indirect inversion techniques, the physically based method does not require forward simulations to assess model‐data misfits. It imposes continuity of hydraulic head and Darcy fluxes in the model domain while incorporating observations (hydraulic heads, Darcy fluxes, or well rates) at measurement locations. Given sufficient measurements, it yields a well‐posed inverse system of equations that can be solved efficiently with coarse grids and nonlinear optimization. When pumping and injection are active, well rates are used as measurements and flux sampling is not needed. The method is successfully tested on synthetic aquifer problems with regular and irregular geometries, different hydrofacies and flow patterns, and increasing conductivity anisotropy ratios. All problems yield stable inverse solutions under increasing head measurement errors. For a given set of observations, inversion accuracy is strongly affected by the conductivity anisotropy ratio. Conductivity estimation is also affected by flow pattern: within a hydrofacies, when Darcy flux component is very small, the corresponding directional conductivity perpendicular to streamlines becomes less identifiable. Finally, inversion is successful even if the location of aquifer boundaries is unknown. In this case, the inversion domain is defined by the location of the measurements.  相似文献   

2.
Simulating groundwater flow in a water‐table (unconfined) aquifer can be difficult because the saturated thickness available for flow depends on model‐calculated hydraulic heads. It is often possible to realize substantial time savings and still obtain accurate head and flow solutions by specifying an approximate saturated thickness a priori, thus linearizing this aspect of the model. This specified‐thickness approximation often relies on the use of the “confined” option in numerical models, which has led to confusion and criticism of the method. This article reviews the theoretical basis for the specified‐thickness approximation, derives an error analysis for relatively ideal problems, and illustrates the utility of the approximation with a complex test problem. In the transient version of our complex test problem, the specified‐thickness approximation produced maximum errors in computed drawdown of about 4% of initial aquifer saturated thickness even when maximum drawdowns were nearly 20% of initial saturated thickness. In the final steady‐state version, the approximation produced maximum errors in computed drawdown of about 20% of initial aquifer saturated thickness (mean errors of about 5%) when maximum drawdowns were about 35% of initial saturated thickness. In early phases of model development, such as during initial model calibration efforts, the specified‐thickness approximation can be a very effective tool to facilitate convergence. The reduced execution time and increased stability obtained through the approximation can be especially useful when many model runs are required, such as during inverse model calibration, sensitivity and uncertainty analyses, multimodel analysis, and development of optimal resource management scenarios.  相似文献   

3.
We present a geostatistically based inverse model for characterizing heterogeneity in parameters of unsaturated hydraulic conductivity for three-dimensional flow. Pressure and moisture content are related to perturbations in hydraulic parameters through cross-covariances, which are calculated to first-order. Sensitivities needed for covariance calculations are derived using the adjoint state sensitivity method. Approximations of the conditional mean parameter fields are then obtained from the cokriging estimator. Correlation between parameters and pressure – moisture content perturbations is seen to be strongly dependent on mean pressure or moisture content. High correlation between parameters and pressure data was obtained under saturated or near saturated flow conditions, providing accurate estimation of saturated hydraulic conductivity, while moisture content measurements provided accurate estimation of the pore size distribution parameter under unsaturated flow conditions.  相似文献   

4.
We present a geostatistically based inverse model for characterizing heterogeneity in parameters of unsaturated hydraulic conductivity for three-dimensional flow. Pressure and moisture content are related to perturbations in hydraulic parameters through cross-covariances, which are calculated to first-order. Sensitivities needed for covariance calculations are derived using the adjoint state sensitivity method. Approximations of the conditional mean parameter fields are then obtained from the cokriging estimator. Correlation between parameters and pressure – moisture content perturbations is seen to be strongly dependent on mean pressure or moisture content. High correlation between parameters and pressure data was obtained under saturated or near saturated flow conditions, providing accurate estimation of saturated hydraulic conductivity, while moisture content measurements provided accurate estimation of the pore size distribution parameter under unsaturated flow conditions.  相似文献   

5.
This work presents a rigorous numerical validation of analytical stochastic models of steady state unsaturated flow in heterogeneous porous media. It also provides a crucial link between stochastic theory based on simplifying assumptions and empirical field and simulation evidence of variably saturated flow in actual or realistic hypothetical heterogeneous porous media. Statistical properties of unsaturated hydraulic conductivity, soil water tension, and soil water flux in heterogeneous soils are investigated through high resolution Monte Carlo simulations of a wide range of steady state flow problems in a quasi-unbounded domain. In agreement with assumptions in analytical stochastic models of unsaturated flow, hydraulic conductivity and soil water tension are found to be lognormally and normally distributed, respectively. In contrast, simulations indicate that in moderate to strong variable conductivity fields, longitudinal flux is highly skewed. Transverse flux distributions are leptokurtic. the moments of the probability distributions obtained from Monte Carlo simulations are compared to modified first-order analytical models. Under moderate to strong heterogeneous soil flux conditions (σ2y≥1), analytical solutions overestimate variability in soil water tension by up to 40% as soil heterogeneity increases, and underestimate variability of both flux components by up to a factor 5. Theoretically predicted model (cross-)covariance agree well with the numerical sample (cross-)covarianaces. Statistical moments are shown to be consistent with observed physical characteristics of unsaturated flow in heterogeneous soils.©1998 Elsevier Science Limited. All rights reserved  相似文献   

6.
Modeling unsaturated flow in porous media requires constitutive relations that describe the soil water retention and soil hydraulic conductivity as a function of either potential or water content. Often, the hydraulic parameters that describe these relations are directly measured on small soil cores, and many cores are needed to upscale to the entire heterogeneous flow field. An alternative to the forward upscaling method using small samples are inverse upscaling methods that incorporate soft data from geophysical measurements observed directly on the larger flow field. In this paper, we demonstrate that the hydraulic parameters can be obtained from cross borehole ground penetrating radar by measuring the first arrival travel time of electromagnetic waves (represented by raypaths) from stationary antennae during a constant flux infiltration experiment. The formulation and coupling of the hydrological and geophysical models rely on a constant velocity wetting front that causes critical refraction at the edge of the front as it passes by the antennae. During this critical refraction period, the slope of the first arrival data can be used to calculate (1) the wetting velocity and (2) the hydraulic conductivity of the wet (or saturated) soil. If the soil is undersaturated during infiltration, then an estimate of the saturated water content is needed before calculating the saturated hydraulic conductivity. The hydraulic conductivity value is then used in a nonlinear global optimization scheme to estimate the remaining two parameters of a Broadbridge and White soil.  相似文献   

7.
Infiltration capacity of bank filtration systems depends on water extraction and hydraulic resistance of the bed sediments. Lakebed hydraulics may be especially affected by clogging, which is dependent on settlement of fine particles, redox potential, and other factors. In the field, most of these processes are difficult to quantify, and thus, when calculating response to pumping the water flux across the sediment surface is assumed to be linearly dependent on the hydraulic gradient. However, this assumption was not adequate to describe conditions at a bank filtration site located at Lake Tegel, Berlin, Germany. Hence, we first assumed the leakage coefficient (or leakance) is spatially distributed and also temporally variant. Furthermore, observations show that the leakance is considerably higher in shallow than in deeper areas; hence, leakance was assumed to be dependent on the existence and thickness of an unsaturated zone below the lake. The proposed explanation of spatial and temporal variability in leakance involves a hypothesis for redox dependent and reversible biogeochemical clogging, supported by geochemical observations in surface water and ground water. Four leakance approaches are implemented in the ground water flow code MODFLOW2000 and calibrated by inverse modeling using the parameter estimation software PEST. These concepts are evaluated by examining the fit to the hydraulic heads, to infiltration measurements, transport modeling results, and considering the degrees of freedom due to the number of calibration parameters. The leakage concept based on the assumption of the influence of an unsaturated zone on clogging processes best explains the field data.  相似文献   

8.
This paper describes a stochastic analysis of steady state flow in a bounded, partially saturated heterogeneous porous medium subject to distributed infiltration. The presence of boundary conditions leads to non-uniformity in the mean unsaturated flow, which in turn causes non-stationarity in the statistics of velocity fields. Motivated by this, our aim is to investigate the impact of boundary conditions on the behavior of field-scale unsaturated flow. Within the framework of spectral theory based on Fourier–Stieltjes representations for the perturbed quantities, the general expressions for the pressure head variance, variance of log unsaturated hydraulic conductivity and variance of the specific discharge are presented in the wave number domain. Closed-form expressions are developed for the simplified case of statistical isotropy of the log hydraulic conductivity field with a constant soil pore-size distribution parameter. These expressions allow us to investigate the impact of the boundary conditions, namely the vertical infiltration from the soil surface and a prescribed pressure head at a certain depth below the soil surface. It is found that the boundary conditions are critical in predicting uncertainty in bounded unsaturated flow. Our analytical expression for the pressure head variance in a one-dimensional, heterogeneous flow domain, developed using a nonstationary spectral representation approach [Li S-G, McLaughlin D. A nonstationary spectral method for solving stochastic groundwater problems: unconditional analysis. Water Resour Res 1991;27(7):1589–605; Li S-G, McLaughlin D. Using the nonstationary spectral method to analyze flow through heterogeneous trending media. Water Resour Res 1995; 31(3):541–51], is precisely equivalent to the published result of Lu et al. [Lu Z, Zhang D. Analytical solutions to steady state unsaturated flow in layered, randomly heterogeneous soils via Kirchhoff transformation. Adv Water Resour 2004;27:775–84].  相似文献   

9.
Cross-borehole flowmeter tests have been proposed as an efficient method to investigate preferential flowpaths in heterogeneous aquifers, which is a major task in the characterization of fractured aquifers. Cross-borehole flowmeter tests are based on the idea that changing the pumping conditions in a given aquifer will modify the hydraulic head distribution in large-scale flowpaths, producing measurable changes in the vertical flow profiles in observation boreholes. However, inversion of flow measurements to derive flowpath geometry and connectivity and to characterize their hydraulic properties is still a subject of research. In this study, we propose a framework for cross-borehole flowmeter test interpretation that is based on a two-scale conceptual model: discrete fractures at the borehole scale and zones of interconnected fractures at the aquifer scale. We propose that the two problems may be solved independently. The first inverse problem consists of estimating the hydraulic head variations that drive the transient borehole flow observed in the cross-borehole flowmeter experiments. The second inverse problem is related to estimating the geometry and hydraulic properties of large-scale flowpaths in the region between pumping and observation wells that are compatible with the head variations deduced from the first problem. To solve the borehole-scale problem, we treat the transient flow data as a series of quasi-steady flow conditions and solve for the hydraulic head changes in individual fractures required to produce these data. The consistency of the method is verified using field experiments performed in a fractured-rock aquifer.  相似文献   

10.
The strong vertical gradient in soil and subsoil saturated hydraulic conductivity is characteristic feature of the hydrology of catchments. Despite the potential importance of these strong gradients, they have proven difficult to model using robust physically based schemes. This has hampered the testing of hypotheses about the implications of such vertical gradients for subsurface flow paths, residence times and transit time distribution. Here we present a general semi‐analytical solution for the simulation of 2D steady‐state saturated‐unsaturated flow in hillslopes with saturated hydraulic conductivity that declines exponentially with depth. The grid‐free solution satisfies mass balance exactly over the entire saturated and unsaturated zones. The new method provides continuous solutions for head, flow and velocity in both saturated and unsaturated zones without any interpolation process as is common in discrete numerical schemes. This solution efficiently generates flow pathlines and transit time distributions in hillslopes with the assumption of depth‐varying saturated hydraulic conductivity. The model outputs reveal the pronounced effect that changing the strength of the exponential decline in saturated hydraulic conductivity has on the flow pathlines, residence time and transit time distribution. This new steady‐state model may be useful to others for posing hypotheses about how different depth functions for hydraulic conductivity influence catchment hydrological response. Copyright © 2016 John Wiley & Sons, Ltd.  相似文献   

11.
In traditional applications in soil physics it is convention to scale porous media properties, such as hydraulic conductivity, soil water diffusivity, and capillary head, with the gravitational acceleration. In addition, the Richards equation for water flux in partially saturated porous media also contains a gravity term. With the plans to develop plant habitats in space, such as in the International Space Station, it becomes necessary to evaluate these properties and this equation under conditions of microgravitational acceleration. This article develops models for microgravity steady state two-phase flow, as found in irrigation systems, that addresses critical design issues. Conventional dimensionless groups in two-phase mathematical models are scaled with gravity, which must be assigned a value of zero for microgravity modeling. The use of these conventional solutions in microgravity, therefore, is not possible. This article therefore introduces new dimensionless groups for two-phase models. The microgravity models introduced here determined that in addition to porous media properties, important design factors for microgravity systems include applied water potential and the ratio of inner to outer radii for cylindrical and spherical porous media systems.  相似文献   

12.
Lu C  Chen Y  Luo J 《Ground water》2012,50(3):386-393
Prevention of sea water intrusion in coastal aquifers subject to groundwater withdrawal requires optimization of well pumping rates to maximize the water supply while avoiding sea water intrusion. Boundary conditions and the aquifer domain size have significant influences on simulating flow and concentration fields and estimating maximum pumping rates. In this study, an analytical solution is derived based on the potential-flow theory for evaluating maximum groundwater pumping rates in a domain with a constant hydraulic head landward boundary. An empirical correction factor, which was introduced by Pool and Carrera (2011) to account for mixing in the case with a constant recharge rate boundary condition, is found also applicable for the case with a constant hydraulic head boundary condition, and therefore greatly improves the usefulness of the sharp-interface analytical solution. Comparing with the solution for a constant recharge rate boundary, we find that a constant hydraulic head boundary often yields larger estimations of the maximum pumping rate and when the domain size is five times greater than the distance between the well and the coastline, the effect of setting different landward boundary conditions becomes insignificant with a relative difference between two solutions less than 2.5%. These findings can serve as a preliminary guidance for conducting numerical simulations and designing tank-scale laboratory experiments for studying groundwater withdrawal problems in coastal aquifers with minimized boundary condition effects.  相似文献   

13.
Measurement uncertainty is a key hindrance to the quantification of water fluxes at all scales of investigation. Predictions of soil‐water flux rely on accurate or representative measurements of hydraulic gradients and field‐state hydraulic conductivity. We quantified the potential magnitude of errors associated with the parameters and variables used directly and indirectly within the Darcy – Buckingham soil‐water‐flux equation. These potential errors were applied to a field hydrometric data set collected from a forested hillslope in central Singapore, and their effect on flow pathway predictions was assessed. Potential errors in the hydraulic gradient calculations were small, approximately one order of magnitude less than the absolute magnitude of the hydraulic gradients. However, errors associated with field‐state hydraulic conductivity derivation were very large. Borehole (Guelph permeameter) and core‐based (Talsma ring permeameter) techniques were used to measure field‐saturated hydraulic conductivity. Measurements using these two approaches differed by up to 3\9 orders of magnitude, with the difference becoming increasingly marked within the B horizon. The sensitivity of the shape of the predicted unsaturated hydraulic conductivity curve to ±5% moisture content error on the moisture release curve was also assessed. Applied moisture release curve error resulted in hydraulic conductivity predictions of less than ±0\2 orders of magnitude deviation from the apparent conductivity. The flow pathways derived from the borehole saturated hydraulic conductivity approach suggested a dominant near‐surface flow pathway, whereas pathways calculated from the core‐based measurements indicated vertical percolation to depth. Direct tracer evidence supported the latter flow pathway, although tracer velocities were approximately two orders of magnitude smaller than the Darcy predictions. We conclude that saturated hydraulic conductivity is the critical hillslope hydrological parameter, and there is an urgent need to address the issues regarding its measurement further. Copyright © 2000 John Wiley & Sons, Ltd.  相似文献   

14.
Surface soil hydraulic properties are key factors controlling the partition of rainfall and snowmelt into runoff and soil water storage, and their knowledge is needed for sound land management. The objective of this study was to evaluate the effects of three land uses (native grass, brome grass and cultivated) on surface soil hydraulic properties under near‐saturated conditions at the St Denis National Wildlife Area, Saskatchewan, Canada. For each land use, water infiltration rates were measured using double‐ring and tension infiltrometers at ?0·3, ?0·7, ?1·5 and ?2·2 kPa pressure heads. Macroporosity and unsaturated hydraulic properties of the surface soil were estimated. Mean field‐saturated hydraulic conductivity (Kfs), unsaturated hydraulic conductivity at ?0·3 kPa pressure head, inverse capillary length scale (α) and water‐conducting macroporosity were compared for different land uses. These parameters of the native grass and brome grass sites were significantly (p < 0·1) higher than that of the cultivated sites. At the ?0·3 kPa pressure head, hydraulic conductivity of grasslands was two to three times greater than that of cultivated lands. Values of α were about two times and values of Kfs about four times greater in grasslands than in cultivated fields. Water‐conducting macroporosity of grasslands and cultivated fields were 0·04% and 0·01% of the total soil volume, respectively. Over 90% of the total water flux at ?0·06 kPa pressure head was transmitted through pores > 1·36 × 10?4 m in diameter in the three land uses. Land use modified near‐saturated hydraulic properties of surface soil and consequently may alter the water balance of the area by changing the amount of surface runoff and soil water storage. Copyright © 2004 John Wiley & Sons, Ltd.  相似文献   

15.
Hydraulic properties of saturated and unsaturated stony soils were studied on a 3.35 m long column, 1.24 m in diameter, filled with alternating sand and boulder layers. The boulders averaged 6.2 × 15 × 20 cm in size and were laid down on their flat side. Tensiometers and a neutron probe access tube were placed in the column for measuring pressure heads and water contents, respectively. Saturated conditions were obtained by ponding the column. The resulting hydraulic conductivity K was 5.1 m/day. This value could also be calculated from the measured K for the sand alone on separate samples, using a simple equation that takes into account the void ratio of the sand alone and that of the boulder-sand mixture. Unsaturated K was determined by applying water at less than ponded infiltration rates. Resulting relations between the unsaturated K and water content or negative pressure head could also be estimated from the relation between unsaturated K and pressure head for the sand alone and the calculated saturated K of the boulder-sand mixture. The method of Millington and Quirk for calculating the relation between unsaturated K and water content also gave reasonable results. The dispersivity of the boulder-sand column was 18 times that of the sand alone. Pore velocity was accurately estimated as the Darcy velocity divided by the volumetric water content. Hydraulic properties of stony vadose zones are difficult to determine. This work shows that they can be estimated from K relations measured in the laboratory on samples of the soil between the rocks. Knowledge of hydraulic properties of vadose zones is important in predicting movement of water and pollutants to the underlying ground water.  相似文献   

16.
Using the first-order analysis, we investigate the spatial cross-correlation between hydraulic conductivity variation and specific discharge (flux) as well as its components measured in a borehole under steady-state flow conditions during cross-hole pumping tests in heterogeneous aquifers. These spatial correlation patterns are found to be quite different from that between the hydraulic conductivity variation and the hydraulic head measurement in the same borehole. This finding suggests that a specific discharge measurement carries non-redundant information about the spatial distribution of heterogeneity, even this measurement is collected from the same location where the head measurement is taken. As such, specific discharge observations should be included in the analysis of hydraulic tomography to increase the resolution of estimated aquifer heterogeneity. Using numerical experiments, we demonstrate the effectiveness of the joint interpretation of both hydraulic heads and fluxes for mapping fracture distributions in a hypothetic geologic medium.  相似文献   

17.
This study presents analytical solutions of the three‐dimensional groundwater flow to a well in leaky confined and leaky water table wedge‐shaped aquifers. Leaky wedge‐shaped aquifers with and without storage in the aquitard are considered, and both transient and steady‐state drawdown solutions are derived. Unlike the previous solutions of the wedge‐shaped aquifers, the leakages from aquitard are considered in these solutions and unlike similar previous work for leaky aquifers, leakage from aquitards and from the water table are treated as the lower and upper boundary conditions. A special form of finite Fourier transforms is used to transform the z‐coordinate in deriving the solutions. The leakage induced by a partially penetrating pumping well in a wedge‐shaped aquifer depends on aquitard hydraulic parameters, the wedge‐shaped aquifer parameters, as well as the pumping well parameters. We calculate lateral boundary dimensionless flux at a representative line and investigate its sensitivity to the aquitard hydraulic parameters. We also investigate the effects of wedge angle, partial penetration, screen location and piezometer location on the steady‐state dimensionless drawdown for different leakage parameters. Results of our study are presented in the form of dimensionless flux‐dimensionless time and dimensionless drawdown‐leakage parameter type curves. The results are useful for evaluating the relative role of lateral wedge boundaries and leakage source on flow in wedge‐shaped aquifers. This is very useful for water management problems and for assessing groundwater pollution. The presented analytical solutions can also be used in parameter identification and in calculating stream depletion rate and volume. Copyright © 2011 John Wiley & Sons, Ltd.  相似文献   

18.
Constant head borehole infiltration tests are widely used for the in situ evaluation of saturated hydraulic conductivities of unsaturated soils above the water table. The formulae employed in analysing the results of such tests disregard the fact that some of the infiltrating water may flow under unsaturated conditions. Instead, these formulae are based on various approximations of the classical free surface theory which treats the flow region as if it were fully saturated and enclosed within a distinct envelope, the so-called ‘free surface’. A finite element model capable of solving free surface problems is used to examine the mathematical accuracy of the borehole infiltration formulae. The results show that in the hypothetical case where unsaturated flow does not exist, the approximate formulae are reasonably accurate within·a practical range of borehole conditions. To see what happens under conditions closer to those actually encountered in the field, the effect of unsaturated flow on borehole infiltration is investigated by means of two different numerical models: a mixed explicit-implicit finite element model, and a mixed explicit-implicit integrated finite difference model. Both of these models give nearly identical results; however, the integrated finite difference model is considerably faster than the finite element model. The relatively low computational efficiency of the finite element scheme is attributed to the large number of operations required in order to re-evaluate the conductivity (stiffness) matrix at each iteration in this highly non-linear saturated-unsaturated flow problem. The saturated-unsaturated analysis demonstrates that the classical free surface approach provides a distorted picture of the flow pattern in the soil. Contrary to what one would expect on the basis of this theory, only a finite region of the soil in the immediate vicinity of the borehole is saturated, whereas a significant percentage of the flow takes place under unsaturated conditions. As a consequence of disregarding unsaturated flow, the available formulae may underestimate the saturated hydraulic conductivity of fine grained soils by a factor of two, three, or more. Our saturated-unsaturated analysis leads to an improved design of borehole infiltration tests and a more accurate method for interpreting the results of such tests. The analysis also shows how one can predict the steady state rate of infiltration from data collected during the early transient period of the test.  相似文献   

19.
Modern ground water characterization and remediation projects routinely require calibration and inverse analysis of large three-dimensional numerical models of complex hydrogeological systems. Hydrogeologic complexity can be prompted by various aquifer characteristics including complicated spatial hydrostratigraphy and aquifer recharge from infiltration through an unsaturated zone. To keep the numerical models computationally efficient, compromises are frequently made in the model development, particularly, about resolution of the computational grid and numerical representation of the governing flow equation. The compromise is required so that the model can be used in calibration, parameter estimation, performance assessment, and analysis of sensitivity and uncertainty in model predictions. However, grid properties and resolution as well as applied computational schemes can have large effects on forward-model predictions and on inverse parameter estimates. We investigate these effects for a series of one- and two-dimensional synthetic cases representing saturated and variably saturated flow problems. We show that "conformable" grids, despite neglecting terms in the numerical formulation, can lead to accurate solutions of problems with complex hydrostratigraphy. Our analysis also demonstrates that, despite slower computer run times and higher memory requirements for a given problem size, the control volume finite-element method showed an advantage over finite-difference techniques in accuracy of parameter estimation for a given grid resolution for most of the test problems.  相似文献   

20.
As pumping proceeds in an unconfined aquifer the water-table is lowered and slowly drains out of the initially saturated material giving a delayed yield from storage. Although considerable work has been carried out on delayed yield, a definitive solution is still lacking. The unsaturated zone was treated here as a significant factor in delayed yield analysis. The one-dimensional unsaturated flow equation was solved numerically to simulate the flow in the zone located between successive water levels in the unconfined aquifer for the case with zero flux at the upper boundary. The moisture content distribution curve was integrated numerically to determine the amount of water drained with time. The study shows that the numerical results are in close agreement with the monitored field data, which implies the important role of the unsaturated zone and the reliability of the numerical simulation in describing this physical phenomenon. The results show that the relationship between the flux and time is sensibly in exponential form, which agrees with the general applicability of Boulton's assumption.  相似文献   

设为首页 | 免责声明 | 关于勤云 | 加入收藏

Copyright©北京勤云科技发展有限公司  京ICP备09084417号