首页 | 本学科首页   官方微博 | 高级检索  
相似文献
 共查询到20条相似文献,搜索用时 31 毫秒
1.
We investigated the flux footprints of receptors at different heights in the convective boundary layer (CBL). The footprints were derived using a forward Lagrangian stochastic (LS) method coupled with the turbulent fields from a large-eddy simulation model. Crosswind-integrated flux footprints shown as a function of upstream distances and sensor heights in the CBL were derived and compared using two LS particle simulation methods: an instantaneous area release and a crosswind linear continuous release. We found that for almost all sensor heights in the CBL, a major positive flux footprint zone was located close to the sensor upstream, while a weak negative footprint zone was located further upstream, with the transition band in non-dimensional upwind distances −X between approximately 1.5 and 2.0. Two-dimensional (2D) flux footprints for a point sensor were also simulated. For a sensor height of 0.158 z i, where z i is the CBL depth, we found that a major positive flux footprint zone followed a weak negative zone in the upstream direction. Two even weaker positive zones were also present on either side of the footprint axis, where the latter was rotated slightly from the geostrophic wind direction. Using CBL scaling, the 2D footprint result was normalized to show the source areas and was applied to real parameters obtained using aircraft-based measurements. With a mean wind speed in the CBL of U = 5.1 m s−1, convective velocity of w * = 1.37 m s−1, CBL depth of z i = 1,000 m, and flight track height of 159 m above the surface, the total flux footprint contribution zone was estimated to range from about 0.1 to 4.5 km upstream, in the case where the wind was perpendicular to the flight track. When the wind was parallel to the flight track, the total footprint contribution zone covered approximately 0.5 km on one side and 0.8 km on the other side of the flight track.  相似文献   

2.
Three surface-layer flux footprint models have been evaluated with the results of an SF6 tracer release experiment specifically designed to test such models. They are a Lagrangian stochastic model, an analytical model, and a simplified derivative of the analytical model. Vertical SF6 fluxes were measured by eddy correlation at four distances downwind of a near-surface crosswind line source in an area of homogeneous sagebrush. The mean fluxes were calculated for 136 half-hour test periods and compared to the fluxes predicted by the footprint models. All three models gave similar predictions and good characterizations of the footprint over the stability range -0.01 < z 0/L < 0.005. The predictions of the three models were within the limits of the uncertainty of the experimental measurements in all but a few cases within this stability range. All three models are unconditionally recommended for determining the area defined by the footprint over short vegetative canopies in this range. They are also generally appropriate for estimating flux magnitudes within the limits of experimental uncertainties. Most of the mean differences observed between the measured and predicted fluxes at each of the four towers reflect a tendency for the measured fluxes to be greater than those predicted by the three models. Rigorous verification of the models in strongly stable conditions was complicated by the need to obtain very accurate measurements of small fluxes in only marginally stationary conditions. Verification in strongly unstable conditions was hampered by the limited number of appropriate data.  相似文献   

3.
Attenuation of Scalar Fluxes Measured with Spatially-displaced Sensors   总被引:1,自引:0,他引:1  
Observations from the Horizontal Array Turbulence Study (HATS) field program are used to examine the attenuation of measured scalar fluxes caused by spatial separation between the vertical velocity and scalar sensors. The HATS data show that flux attenuation for streamwise, crosswind, and vertical sensor displacements are each a function of a dimensionless, stability-dependent parameter n m multiplied by the ratio of sensor displacement to measurement height. The scalar flux decays more rapidly with crosswind displacements than for streamwise displacements and decays more rapidly for stable stratification than for unstable stratification. The cospectral flux attenuation model of Kristensen et al. agrees well with the HATS data for streamwise sensor displacements, although it is necessary to include a neglected quadrature spectrum term to explain the observation that flux attenuation is often less with the scalar sensor downwind of the anemometer than for the opposite configuration. A simpler exponential decay model provides good estimates for crosswind sensor displacements, as well as for streamwise sensor displacements with stable stratification. A model similar to that of Lee and Black correctly predicts flux attenuation for a combination of streamwise and crosswind displacements, i.e. as a function of wind direction relative to the sensor displacement. The HATS data for vertical sensor displacements extend the near-neutral results of Kristensen et al. to diabatic stratification and confirm their finding that flux attenuation is less with the scalar sensor located below the anemometer than if the scalar sensor is displaced an equal distance either horizontally or above the anemometer.  相似文献   

4.
A Lagrangian particle dispersion model (LPDM) driven by velocity fields from large-eddy simulations (LESs) is used to determine the mean and variability of plume dispersion in a highly convective planetary boundary layer (PBL). The total velocity of a “particle” is divided into resolved and unresolved or random (subfilter scale, SFS) velocities with the resolved component obtained from the LES and the SFS velocity from a Lagrangian stochastic model. This LPDM-LES model is used to obtain an ensemble of dispersion realizations for calculating the mean, root-mean-square (r.m.s.) deviation, and fluctuating fields of dispersion quantities. An ensemble of 30 realizations is generated for each of three source heights: surface, near-surface, and elevated. We compare the LPDM calculations with convection tank experiments and field observations to assess the realism of the results. The overall conclusion is that the LPDM-LES model produces a realistic range of dispersion realizations and statistical variability (i.e., r.m.s. deviations) that match observations in this highly convective PBL, while also matching the ensemble-mean properties. This is true for the plume height or trajectory, vertical dispersion, and the surface values of the crosswind-integrated concentration (CWIC), and their dependence on downstream distance. One exception is the crosswind dispersion for an elevated source, which is underestimated by the model. Other analyses that highlight important LPDM results include: (1) the plume meander and CWIC fluctuation intensity at the surface, (2) the applicability of a similarity theory for plume height from a surface source to only the very strong updraft plumes—not the mean height, and (3) the appropriate variation with distance of the mean surface CWIC and the lower bound of the CWIC realizations for a surface source.  相似文献   

5.
Vertical dispersion in the neutral surface layer is investigated using a Markov Chain simulation procedure. The conceptual basis of the procedure is discussed and computation procedures outlined. Wind and turbulence parameterizations appropriate to the neutral surface layer are considered with emphasis on the Lagrangian time scale. Computations for a surface release are compared with field data. Good agreement is found for the variation of surface concentration and cloud height to distances 500 m downwind of the source. The functional form of the vertical concentration profile is examined and an exponential with exponent ∼1.6 is found to give the best fit with simulations. For elevated releases, it is demonstrated that an initial dip of the mass mean height from the simulation can be normalized for various release heights using a non-dimensionalized downwind coordinate incorporating advective wind speed and wind shear. The vertical distribution standard deviation (σz), as employed in Gaussian models, shows a fair degree of independence with source height but close examination reveals an optimum source height for maximum σz at a given downwind distance,x. This source height increases with downwind distance. Also the simulations indicate that vertical wind shear is more important than vertical variation of Lagrangian time scale close to the source, with a reverse effect farther downwind.  相似文献   

6.
Based on the theoretical background of existing models for the crosswind-integrated footprint, a new model is presented, which, in contrast to the existing models, describes the normalized footprint by a closed analytical formula. This was made possible by using well-known power profiles for wind speed and eddy viscosity instead of Monin–Obukhov based profiles at a certain stage of model development. However, the major difference between the new model and the existing models is that the so-called shape parameter of vertical plume dispersion, a function of upwind distance in the existing models, is set constant in the new model in order to circumvent a formal inconsistency found in the derivation of the existing models. Due to this inconsistency, the existing models do not generally satisfy the fundamental condition that the cumulative normalized footprint must approach unity for the upwind distance tending towards infinity.  相似文献   

7.
It is frequently observed in field experiments that the eddy covariance heat fluxes are systematically underestimated as compared to the available energy. The flux imbalance problem is investigated using the NCAR’s large-eddy simulation (LES) model imbedded with an online scheme to calculate Reynolds-averaged fluxes. A top–down and a bottom–up tracer are implemented into the LES model to quantify the influence of entrainment and bottom–up diffusion processes on flux imbalance. The results show that the flux imbalance follows a set of universal functions that capture the exponential decreasing dependence on u */w *, where u * and w * are friction velocity and the convective velocity scale, respectively, and an elliptic relationship to z/z i , where z i is the mixing-layer height. The source location in the boundary layer is an important factor controlling the imbalance magnitude and its horizontal and vertical distributions. The flux imbalance of heat and the bottom–up tracer is tightly related to turbulent coherent structures, whereas for the top–down diffusion, such relations are weak to nonexistent. Our results are broadly consistent with previous studies on the flux imbalance problem, suggesting that the published results are robust and are not artefacts of numerical schemes.  相似文献   

8.
Vertical dispersion in the neutral surface layer is investigated using a Markov Chain simulation procedure. The conceptual basis of the procedure is discussed and computation procedures outlined. Wind and turbulence parameterizations appropriate to the neutral surface layer are considered with emphasis on the Lagrangian time scale. Computations for a surface release are compared with field data. Good agreement is found for the variation of surface concentration and cloud height to distances 500 m downwind of the source. The functional form of the vertical concentration profile is examined and an exponential with exponent 1.6 is found to give the best fit with simulations.For elevated releases, it is demonstrated that an initial dip of the mass mean height from the simulation can be normalized for various release heights using a non-dimensionalized downwind coordinate incorporating advective wind speed and wind shear. The vertical distribution standard deviation ( z ), as employed in Gaussian models, shows a fair degree of independence with source height but close examination reveals an optimum source height for maximum z at a given downwind distance,x. This source height increases with downwind distance. Also the simulations indicate that vertical wind shear is more important than vertical variation of Lagrangian time scale close to the source, with a reverse effect farther downwind.  相似文献   

9.
By non-dimensionalizing a trajectory-simulation (TS) model of turbulent dispersion, it is shown that the dimensionless concentration z 0cu*/kQ (cu */kQ) due to a continuous line (area) source of strength Q in the atmospheric surface layer depends only on z/z 0, x/z 0, z 0/L and z s/z0, where z s is the source height. The TS model is used to tabulate concentration profiles due to ground-level line and area sources. Concentration profiles generated by the TS model for elevated sources are shown to be inconsistent with the Reciprocal Theorems of Smith (1957) and it is suggested that this is because the flux-mean gradient closure scheme inherent in the Reciprocal Theorem is invalid for an elevated source.  相似文献   

10.
A Lagrangian statistical-trajectory model based on a Markov chain relation is used to investigate vertical dispersion from elevated sources into the neutral planetary boundary layer. The model is fully two-dimensional, in that both vertical and longitudinal velocity fluctuations, and their correlation, are simulated explicitly. The best observational information currently available is used to characterize the mean and turbulent structure of the neutral boundary layer. In particular, a realistic vertical profile of the Lagrangian integral time scale is proposed, based partly on a review of direct measurements and partly on a comparison of the model predictions with published diffusion data. The model predictions are shown to agree well with a variety of dispersion observations. The model is used to study vertical diffusion as a function of release height H, friction velocity u* and surface roughness z 0 for downwind distances up to 10 km from the source. The equivalent Gaussian dispersion parameter Σ z is shown to decrease slightly with an increase in H, and to increase with increases in z 0 or u*. It is demonstrated that relationships valid in a field of homogeneous turbulence can be applied to vertical dispersion in the atmosphere if the release occurs above the region of strongest gradients in the mean and turbulent parameters. Scaling in terms of the standard deviation in elevation angle of the wind at the release point leads to a universal curve which provides accurate estimates of Σ z over a wide range of values of H, z 0 and the meteorological parameters.  相似文献   

11.
We examine daily (morning–afternoon) transitions in the atmospheric boundary layer based on large-eddy simulations. Under consideration are the effects of the stratification at the top of the mixed layer and of the wind shear. The results describe the transitory behaviour of temperature and wind velocity, their second moments, the boundary-layer height Z m (defined by the maximum of the potential temperature gradient) and its standard deviation σ m , the mixed-layer height z i (defined by the minimum of the potential temperature flux), entrainment velocity W e, and the entrainment flux H i . The entrainment flux and the entrainment velocity are found to lag slightly in time with respect to the surface temperature flux. The simulations imply that the atmospheric values of velocity variances, measured at various instants during the daytime, and normalized in terms of the actual convective scale w*, are not expected to collapse to a single curve, but to produce a significant scatter of observational points. The measured values of the temperature variance, normalized in terms of the actual convective scale Θ*, are expected to form a single curve in the mixed layer, and to exhibit a considerable scatter in the interfacial layer.  相似文献   

12.
Turf-grass lawns are ubiquitous in the United States. However direct measurements of land–atmosphere fluxes using the eddy-covariance method above lawn ecosystems are challenging due to the typically small dimensions of lawns and the heterogeneity of land use in an urbanised landscape. Given their typically small patch sizes, there is the potential that CO2 fluxes measured above turf-grass lawns may be influenced by nearby CO2 sources such as passing traffic. In this study, we report on two years of eddy-covariance flux measurements above a 1.5 ha turf-grass lawn in which we assess the contribution of nearby traffic emissions to the measured CO2 flux. We use winter data when the vegetation was dormant to develop an empirical estimate of the traffic effect on the measured CO2 fluxes, based on a parametrised version of a three-dimensional Lagrangian footprint model and continuous traffic count data. The CO2 budget of the ecosystem was adjusted by 135gCm−2 in 2007 and by 134gCm−2 in 2008 to determine the natural flux, even though the road crossed the footprint only at its far edge. We show that bottom-up flux estimates based on CO2 emission factors of the passing vehicles, combined with the crosswind-integrated footprint at the distance of the road, agreed very well with the empirical estimate of the traffic contribution that we derived from the eddy-covariance measurements. The approach we developed may be useful for other sites where investigators plan to make eddy-covariance measurements on small patches within heterogeneous landscapes where there are significant contrasts in flux rates. However, we caution that the modelling approach is empirical and will need to be adapted individually to each site.  相似文献   

13.
Summary This paper investigates the influence of the planetary boundary-layer (PBL) parameterization and the vertical distribution of model layers on simulations of an Alpine foehn case that was observed during the Mesoscale Alpine Programme (MAP) in autumn 1999. The study is based on the PSU/NCAR MM5 modelling system and combines five different PBL schemes with three model layer settings, which mainly differ in the height above ground of the lowest model level (z 1). Specifically, z 1 takes values of about 7 m, 22 m and 36 m, and the experiments with z 1 = 7 m are set up such that the second model level is located at z = 36 m. To assess if the different model setups have a systematic impact on the model performance, the simulation results are compared against wind lidar, radiosonde and surface measurements gathered along the Austrian Wipp Valley. Moreover, the dependence of the simulated wind and temperature fields at a given height (36 m above ground) on z 1 is examined for several different regions. Our validation results show that at least over the Wipp Valley, the dependence of the model skill on z 1 tends to be larger and more systematic than the impact of the PBL scheme. The agreement of the simulated wind field with observations tends to benefit from moving the lowest model layer closer to the ground, which appears to be related to the dependence of lee-side flow separation on z 1. However, the simulated 2 m-temperatures are closest to observations for the intermediate z 1 of 22 m. This is mainly related to the fact that the simulated low-level temperatures decrease systematically with decreasing z 1 for all PBL schemes, turning a positive bias at z 1 = 36 m into a negative bias at z 1 = 7 m. The systematic z 1-dependence is also observed for the temperatures at a fixed height of 36 m, indicating a deficiency in the self-consistency of the model results that is not related to a specific PBL formulation. Possible reasons for this deficiency are discussed in the paper. On the other hand, a systematic z 1-dependence of the 36-m wind speed is encountered only for one out of the five PBL schemes. This turns out to be related to an unrealistic profile of the vertical mixing coefficient. Correspondence: Günther Z?ngl, Meteorologisches Institut der Universitat München, 80333 München, Germany  相似文献   

14.
A meandering plume model that explicitly incorporates the effects of small-scale structure in the instantaneous plume has been formulated. The model requires the specification of two physically based input parameters; namely, the meander ratio,M, which is dependent on the ratio of the meandering plume dispersion to the instantaneous relative plume dispersion and, a relative in-plume fluctuation measure,k, that is related inversely to the fluctuation intensity in relative coordinates. Simple analytical expressions for crosswind profiles of the higher moments (including the important shape parameters such as fluctuation intensity, skewness, and kurtosis) and for the concentration pdf have been derived from the model. The model has been tested against some field data sets, indicating that it can reproduce many key aspects of the observed behavior of concentration fluctuations, particularly with respect to modeling the change in shape of the concentration pdf in the crosswind direction.List of Symbols C Mean concentration in absolute coordinates - C r Mean concentration in relative coordinates - C0 Centerline mean concentration in absolute coordinates - C r,0 Centerline mean concentration in relative coordinates - f Probability density function of concentration in absolute coordinates - f c Probability density function of plume centroid position - f r Probability density function of concentration in relative coordinates - i Absolute concentration fluctuation intensity (standard deviation to mean ratio) - i r Relative concentration fluctuation intensity (standard deviation to mean ratio) - k Relative in-plume fluctuation measure:k=1/i r 2 - K Concentration fluctuation kurtosis - M Meander ratio of meandering plume variance to relative plume variance - S Concentration fluctuation skewness - x Downwind distance from source - y Crosswind distance from mean-plume centerline - z Vertical distance above ground - Instantaneous (random) concentration - Crosswind dispersion ofnth concentration moment about zero - ny Mean-plume crosswind (absolute) dispersion - y Plume centroid (meandering) dispersion in crosswind direction - y,c Instantaneous plume crosswind (relative) dispersion - Normalized mean concentration in absolute coordinates:C/C 0 - Particular value taken on by instantaneous concentration,   相似文献   

15.
A simple mixed layer model is used to derive the following expressions for the maximum (daily) convective velocity scale. w *m = AQ m 1/2; w *m = Bz im The variables A and B are shown to vary within narrow limits thus allowing them to be treated as constants. This is very useful for routine computation of w *m , an important variable for dispersion under unstable conditions, from estimates of either the kinematic surface heat flux Q m (m-1) or the maximum mixed layer height z im .Analysis of observations made during the Minnesota boundary layer experiment shows that there is ample justification for assigning typical values to A and B in estimating w *m.  相似文献   

16.
We present a numerical study aimed at quantifying the effects of concentration-dependent density on the spread of a seeping plume of CO2 into the atmosphere such as could arise from a leaking geologic carbon sequestration site. Results of numerical models can be used to supplement field monitoring estimates of CO2 seepage flux by modelling transport and dispersion between the source emission and concentration-measurement points. We focus on modelling CO2 seepage dispersion over relatively short distances where density effects are likely to be important. We model dense gas dispersion using the steady-state Reynolds-averaged Navier-Stokes equations with density dependence in the gravity term. Results for a two-dimensional system show that a density dependence emerges at higher fluxes than prior estimates. A universal scaling relation is derived that allows estimation of the flux from concentrations measured downwind and vice versa.  相似文献   

17.
A numerical model of airflow in the lowest 50–100 m of the atmosphere above changes in surface roughness and temperature or heat flux has been developed based on boundary layer approximations, the Businger-Dyer hypotheses for the non-dimensional wind shear and heat flux and a mixing length hypothesis.Results have been obtained for several situations, in particular, airflow with neutral upstream conditions encountering a step change in surface temperature or heat flux with no roughness change. In these cases large increases in shear stress at the outer edge of the internal boundary layer are predicted. The case of unstable upstream flow encountering a step change to zero heat flux is also considered.Two situations that may be encountered near the shores of the Great Lakes are considered.Notation B Businger-Dyer constant (= 16.0) in form for M, H - c p Specific heat at constant pressure - g Acceleration due to gravity - H Upward vertical heat flux - H 0 , H 1 Surface heat fluxes for x < 0, x 0 - k von Kármán's constant ( = 0.4) - l Mixing length - L Monin-Obukhov length - L 0 Upstream value of L - m Ratio of roughness lengths (= z 1/z 0) - RL * Non-dimensional parameter, see Equations (20, 22 and 24) - RL 1 * Same as RL * but with z 1 scaling (= mRL *) - T Scaled temperature - T 0 (z) Upstream temperature profile - u 0, u 1(x) Surface friction velocities for x < 0, x 0 - U, W Horizontal and vertical mean velocities - U 0 (z) Upstream velocity profile - x, z Horizontal and vertical coordinates - z i Local roughness length  相似文献   

18.
The winter-time arctic atmospheric boundary layer was investigated with micrometeorological and SF6 tracer measurements collected in Prudhoe Bay, Alaska. The flat, snow-covered tundra surface at this site generates a very small (0.03 cm) surface roughness. The relatively warm maritime air mass originating over the nearby, partially frozen Beaufort Sea is cooled at the tundra surface resulting in strong (4 to 30 °C · (100 m)-1) temperature inversions with light winds and a persistent weak (1 to 2 °C · (100 m)-1) surface inversion with wind speeds up to 17 m s-1. The absence of any diurnal atmospheric stability pattern during the study was due to the very limited solar insolation. Vertical profiles were measured with a multi-level mast from 1 to 17 m and with a Doppler acoustic sounder from 60 to 450 m. With high wind speeds, stable layers below 17 m and above 300 m were typically separated by a layer of neutral stability. Turbulence statistics and spectra calculated at a height of 33 m are similar to measurements reported for non-arctic, open terrain sites and indicate that the production of turbulence is primarily due to wind shear. The distribution of wind direction recorded at 1 Hz was frequently non-Gaussian for 1-hr periods but was always Gaussian for 5-min periods. We also observed non-Gaussian hourly averaged crosswind concentration profiles and assume that they can be modeled by calculating sequential short-term concentrations, using the 5-min standard deviation of horizontal wind direction fluctuations () to estimate a horizontal dispersion coefficient ( y ), and constructing hourly concentrations by averaging the short-term results. Non-Gaussian hourly crosswind distributions are not unique to the arctic and can be observed at most field sites. A weak correlation between horizontal ( v ) and vertical ( w ) turbulence observed for both 1-hr and 5-min periods indicates that a single stability classification method is not sufficient to determine both vertical and horizontal dispersion at this site. An estimate of the vertical dispersion coefficient, z , could be based on or a stability classification parameter which includes vertical thermal and wind shear effects (e.g., Monin-Obukhov length, L).  相似文献   

19.
One-dimensional Lagrangian dispersion models, frequently used to relate in-canopy source/sink distributions of energy, water and trace gases to vertical concentration profiles, require estimates of the standard deviation of the vertical wind speed, which can be measured, and the Lagrangian time scale, T L , which cannot. In this work we use non-linear parameter estimation to determine the vertical profile of the Lagrangian time scale that simultaneously optimises agreement between modelled and measured vertical profiles of temperature, water vapour and carbon dioxide concentrations within a 40-m tall temperate Eucalyptus forest in south-eastern Australia. Modelled temperature and concentration profiles are generated using Lagrangian dispersion theory combined with source/sink distributions of sensible heat, water vapour and CO2. These distributions are derived from a multilayer Soil-Vegetation-Atmospheric-Transfer model subject to multiple constraints: (1) daytime eddy flux measurements of sensible heat, latent heat, and CO2 above the canopy, (2) in-canopy lidar measurements of leaf area density distribution, and (3) chamber measurements of CO2 ground fluxes. The resulting estimate of Lagrangian time scale within the canopy under near-neutral conditions is about 1.7 times higher than previous estimates and decreases towards zero at the ground. It represents an advance over previous estimates of T L , which are largely unconstrained by measurements.  相似文献   

20.
A wind-tunnel study was conducted to investigate ventilation of scalars from urban-like geometries at neighbourhood scale by exploring two different geometries a uniform height roughness and a non-uniform height roughness, both with an equal plan and frontal density of λ p = λ f = 25%. In both configurations a sub-unit of the idealized urban surface was coated with a thin layer of naphthalene to represent area sources. The naphthalene sublimation method was used to measure directly total area-averaged transport of scalars out of the complex geometries. At the same time, naphthalene vapour concentrations controlled by the turbulent fluxes were detected using a fast Flame Ionisation Detection (FID) technique. This paper describes the novel use of a naphthalene coated surface as an area source in dispersion studies. Particular emphasis was also given to testing whether the concentration measurements were independent of Reynolds number. For low wind speeds, transfer from the naphthalene surface is determined by a combination of forced and natural convection. Compared with a propane point source release, a 25% higher free stream velocity was needed for the naphthalene area source to yield Reynolds-number-independent concentration fields. Ventilation transfer coefficients w T /U derived from the naphthalene sublimation method showed that, whilst there was enhanced vertical momentum exchange due to obstacle height variability, advection was reduced and dispersion from the source area was not enhanced. Thus, the height variability of a canopy is an important parameter when generalising urban dispersion. Fine resolution concentration measurements in the canopy showed the effect of height variability on dispersion at street scale. Rapid vertical transport in the wake of individual high-rise obstacles was found to generate elevated point-like sources. A Gaussian plume model was used to analyse differences in the downstream plumes. Intensified lateral and vertical plume spread and plume dilution with height was found for the non-uniform height roughness.  相似文献   

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

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