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

    Forecasting Zonda Wind Occurrence with Vertical Sounding Data

    2022-01-15 07:22:18FedericoOTEROandDiegoARANEO
    Advances in Atmospheric Sciences 2022年1期

    Federico OTERO and Diego C. ARANEO

    Instituto Argentino de Nivología, Glaciología y Ciencias Ambientales (IANIGLA) CCT Mendoza - CONICET Av.Ruiz Leal s/n., Parque Gral. San Martín, Mendoza 5500, Argentina

    (Received 20 January 2021; revised 6 July 2021; accepted 12 August 2021)

    ABSTRACT Zonda wind is a typical downslope windstorm over the eastern slopes of the Central Andes in Argentina, which produces extremely warm and dry conditions and creates substantial socioeconomic impacts. The aim of this work is to obtain an index for predicting the probability of Zonda wind occurrence. The Principal Component Analysis (PCA) is applied to the vertical sounding data on both sides of the Andes. Through the use of a binary logistic regression, the PCA is applied to discriminate those soundings associated with Zonda wind events from those that are not, and a probabilistic forecasting tool for Zonda occurrence is obtained. This index is able to discriminate between Zonda and non-Zonda events with an effectiveness close to 91%. The best model consists of four variables from each side of the Andes. From an eventbased statistical perspective, the probability of detection of the mixed model is above 97% with a probability of false detection lower than 7% and a missing ratio below 1%. From an alarm-based perspective, models exhibit false alarm rate below 7%, a missing alarm ratio lower than 1.5% and higher than 93% for the correct alarm ratio. The zonal component of the wind on both sides of the Andes and the windward temperature are the key variables in class discrimination. The vertical structure of Zonda wind includes two wind maximums and an unstable lapse rate at midlevels on the lee side and a wind maximum at 700 hPa accompanied by a relatively stable layer near the mountain top.

    Key words: Zonda wind, foehn, downslope windstorm, forecasting

    1. Introduction

    Zonda wind (Argentinian foehn) is a strong, warm, and very dry wind associated with adiabatic compression upon descending over the eastern slopes of the Andes Cordillera and occurs most often in winter and spring, mainly in the provinces of Mendoza and San Juan (Norte, 2015; Otero and Norte, 2015). Despite the constant improvement of numerical weather prediction (NWP) models and the advances in the understanding of mountain meteorology dynamics over the past decades, downslope windstorms forecasting is still limited by several factors. These factors include the dependency and sensitivity to the model’s resolution (e.g., Reinecke and Durran, 2009b; Jackson et al.,2013), numerical schemes, vertical coordinates and diffusion parametrizations (e.g., Smith, 2007), physical formulations (the boundary layer especially, see Smith, 2007), and initial condition uncertainties (e.g., Reinecke and Durran,2009a). Thus, it is necessary to predict occurrence using statistical forecasting models specifically developed for these stations, which can improve or contribute to the NWP models.

    In the 1960s, Widmer developed a “foehn test” for the Altdorf, Switzerland foehn station that was refined by Courvoisier and Gutermann (1971). This test remains as the operational tool used today. Later, Dürr (2008) developed an automated method for identifying foehn (i.e., nowcasting foehn). His procedure is based on 10-min real-time data from the automated Swiss surface network. Quite recently,Drechsel and Mayr (2008) developed an objective, probabilistic forecasting method for foehn in the Wipp Valley (Innsbruck) based on the ECMWF model output. In Sprenger et al. (2017), a new objective method for foehn prediction based on a machine learning algorithm (called AdaBoost,short for adaptive boosting) is proposed to distinguish between foehn and non-foehn events. Further improvement will require the use of not only deterministic, but also statistical methods. In particular, ongoing work shows that model output statistics (MOS) are a promising tool for improving foehn forecasting.

    Considering the influence on the flow response when interacting with a topographic barrier, the development of the phenomenon is closely related to the vertical structure of the atmosphere, where mountain height, buoyancy frequency, and incident wind control the mountain waves activity and the dynamics of downslope wind (Durran, 1990;Smith and Skyllingstad, 2011; Damiens et al., 2018, among others). The application of statistical techniques, such as the study of Empirical Orthogonal Functions (EOF) or the Principal Component Analysis (PCA), allow the objective statistical characterization of vector and scalar variables or other physical variables, like temperature, humidity, wind, and stability. Otero et al. (2018) obtained an index for predicting the Zonda occurrence through the vertical sounding of the lee side of the Andes. A PCA is used to identify the patterns of the vertical structure of the atmosphere leading up to a Zonda wind event and used to construct the index model. A Zonda/non-Zonda index is calculable fromTandTdprofiles and is dependent on the climatological features of the region.

    In the present work, the methodology of Otero et al.(2018) is followed. Here, a substantial improvement is presented with respect to the previous work, where only two combined leeward variables were used (temperature and dewpoint temperature). In this new version, not only are windward soundings (Chilean side) added, but so are new and different combinations between variables and soundings, combining up to four variables on each side of the Andes. The newly added variables are wind (UandV) and the squared Brunt-V?is?l? frequency (N2), which are key factors in the atmospheric conditions for the development of this kind of downslope windstorm. Likewise, new metrics of the prediction model (not presented in Otero et al., 2018) are calculated such as the probability of detection (POD), probability of false detection (PODF), missing ratio (MR), false alarm rate (FAR), missing alarm ratio (MAR), and correct alarm ratio (CAR). In this case, the PCA is used to characterize the vertical structure on both sides of the Andes Mountains prior to the onset of a downslope windstorm (Zonda).A complete description of the vertical structure of wind, temperature, dewpoint, and stability for Zonda and non-Zonda events and the characteristic structure that discriminates between both classes are obtained.

    The structure of this paper is as follows: Section 2 describes the data and methodology, where the PCA and index model are described. Section 3 presents the results including mean and Zonda soundings, the index model efficiency and metrics, and the discriminant sounding. Finally,the discussion and conclusions and their applications are presented in section 4.

    2. Data and methodology

    2.1. Data

    For this study, available daily sounding data at 1200 UTC for the 1981-2019 period are taken for both sides of the Andes Mountains. The Santo Domingo surface station(33.65°S, 71.61°W, 75 m A. S. L.) on the Chilean side(Andes windward) and the Mendoza Airport surface station(32.83°S, 68.77°W, 704 m A. S. L.) on the Argentinean side(Andes lee side) are used (Fig. 1). The selected pressure levels for Mendoza's soundings are those standard levels between 850 hPa and 300 hPa (i.e., 850, 700, 500, 400, and 300 hPa) and for Santo Domingo those between 1000 hPa and 200 hPa (i.e., 1000, 850, 700, 500, 400, 300 and 200 hPa). The sectioned levels are due to the vertical resolution of the vertical soundings of each side, those on the windward side having the highest vertical (as well as temporal) resolution, as well as the altitude of the locations. The selectedvariables are temperature (T), dewpoint temperature (Td),zonal (U) and meridional (V) wind components, and the squared Brunt-V?is?l? frequency (N2). A data consistency is made by removing those inconsistent values. Missing data are detected and marked along with suspicious and out of range values, which is extremely challenging due to the event’s extreme nature.

    2.2. Methodology

    2.2.1.Zonda wind events classification and climatology

    The Zonda wind classification is made using hourly data from three surface stations on the lee side (Fig. 1, black dots) belonging to the National Weather Service of Argentina (SMN): Mendoza Airport, Mendoza Observatory(32.9°S, 68.866°W, 827 m A. S. L.), and San Juan Airport(31.56°S, 68.5°W, 598 m A. S. L.). The predictor variables are temperature (T), dewpoint temperature (Td), surface pressure (P), and 10-m wind speed (V). The surface stations report hourly and have data records of 35 years for Mendoza Airport, 24 years for San Juan Airport, and 13 years for Mendoza Observatory. This classification is made manually (subjective) according to an abrupt increase in temperature and decrease in dewpoint temperature in conjunction with an increase in wind speed.

    Figure 2a shows the Zonda climatology for the three surface stations and for the sounding data. Note that the Zonda wind presents higher frequencies of occurrence in winter and spring (mainly from June to October). The soundings used in this work must meet the requirements of section 2.1,in addition to having soundings on both sides of the Andes on the same day. This implies a total of 58 soundings (previous to the onset of the events) that presents the maximum frequency in October. The duration of Zonda events shows maximum frequencies between two to eight hours for Mendoza’s station and longer for San Juan Airport (Fig. 2b). The onset time presents highest frequencies in the afternoon, from 1500 UTC to 2100 UTC during winter, and a second maximum between 0000 UTC to 0300 UTC in spring(Fig. 2c).

    Fig. 2. (a) Zonda events distribution for the surface stations and for sounding data used. (b) Zonda duration annual frequency distribution. (c) Zonda onset time annual frequency distribution. The straight line corresponds to the sounding hour (1200 UTC).

    Fig. 3. (a) Mean (solid line) and Zonda (dotted line) vertical soundings at 1200 UTC for Mendoza Airport and (b)Santo Domingo (right). Wind barbs on the left correspond to the Zonda sounding, and wind barbs on the right correspond to the mean sounding.

    Once all Zonda events and their onset times have been found, the sounding data series is constructed. For this, the 1200 UTC sounding closest to the onset time of each event is considered (according to the onset time frequency, the 1200 UTC sounding is more likely to be closer to the event than the 0000 UTC sounding, and at the Mendoza Airport,no sounding for 0000 UTC is carried out, at least in a large part of the record). For example, if the event starts at 1500 UTC, the 1200 UTC sounding on the same day is considered (i.e., three hours after the event). If the event starts at 0900 UTC, the 1200 UTC sounding of the previous day is considered (i.e., 21 hours after the event). It is in this way that these soundings represent the characteristics of the atmosphere prior to the development of the Zonda wind.

    The selection of the days without Zonda is done in such a way that the closest Zonda day is at least five days away.In this way, it is ensured that the synoptic conditions are as different as possible from those of a Zonda day. Data from a total of 116 soundings associated with Zonda and non-Zonda events are chosen to perform the PCA. This dataset is the same used in Otero and Araneo (2021).

    2.2.2.Construction of the Zonda/non-Zonda index

    The methodology for the PCA follows that used in Araneo et al. (2011) and Otero et al. (2018). The Zonda/non-Zonda probability index is obtained from a logistic regression between the PCA loading components and a vector of 0 and 1, associated with Zonda/non-Zonda.To detect the patterns capable of discriminating between Zonda and non-Zonda cases, they must first be compared with random non-Zonda conditions. The selection of those dates is carefully made so that they are not close to any Zonda events. As a necessary condition, it is considered that these days are at least five days away from any Zonda events.

    For the calculation of the PCA, the sounding data are arranged forming a matrix [ X]5×116( [ X]5×116is a matrix of size 5 × 116) in which the rows contain the values corresponding to each pressure level for the variable to evaluate (T,Td,U,V, orN2), and each column represents the sounding for a given day. After that, the mean sounding is removed, obtaining the deviations matrixThis matrix is then standardized by columns. Of the 116 days, half correspond to Zonda soundings, and the other half correspond to non-Zonda soundings. In the case of Santo Domingo, the matrix has dimensions of [7 × 116] because more vertical levels are considered.

    The probability model (index) is built with the fit coefficients of a logistic regression for a binomial function. So,the logistic regression is transformed into probability values between 0 and 1. Then, for theith element, the probability will be given by:

    wherewiis theith element ofw=b0+b1f1+b2f2+...+bnfn, f1,f2,...,fnare the principal component loadings(columns ofF),b0,b1,......,bnare the regression coefficients (fitted with the maximum likelihood method), andnis the number of significant components retained. Then,︿ciis an estimator of theith coefficient ofc(known vector of zeros and ones, associated with Zonda/non-Zonda events) corresponding to theith day. This estimator represents the probability that day has of being classified as Zonda or non-Zonda based on a preset cutoff threshold of the index.

    The component loadings matrix can be written aswheremis the number of rows (days) ofis the matrix containing the standardized score components corresponding to the predictor variables used in the model (i.e., all those components related to the significant fit coefficientsbi) , b*is the vector matrix containing the coefficientsb0,b1,......,bn, and [1] is a column vector of elements equal to 1.

    In order to estimate prediction errors, the leave-one-out cross-validation method is implemented. In other words, for the construction of the Zonda index, all the dates except one are used in each step, which is used for its verification.Thus, in each step, the set formed byn-1 observations is considered the fitting set. The observation which is left out is then used to test the regression model obtained with the remaining ones. This procedure is performed with each date, obtaining a verification for each particular date.

    This procedure can be generalized using all the predictor variables in a single prediction index to improve the efficiency of the forecast index. Suppose the case with 2 variables, in principle with the same number of casesn,X1andX2, with dimensionstandr, respectively. For example, suppose that we havendays,X1is the matrix containing the Mendoza station’s soundings for each day’s temperatures attvertical levels, andX2is the matrix containing the Santo Domingo station’s soundings for each day’s squared Brunt-V?is?l? frequency atrvertical levels; then:

    Note that [A1]t×1and [A2]r×1have the same dimensions of the predictor variable (i.e., number of vertical levels of Mendoza and Santo Domingo’s soundings, respectively) and that the ︿cvalue depends on the dot product with the particular variable of the day to forecast. The value of︿cwill be closer to 1 as the larger the scalar products between the particular variables and the vectorsAwith a positive sign, and it will be closer to 0 as less the smaller the scalar products with a negative sign be. Therefore, theAvectors can be interpreted as a "discriminant" vertical sounding of Zonda or non-Zonda situations.

    2.2.3.Models metrics

    To compare among the models and the predictive skills of each one, different metrics, obtained from the confusion matrix, are used following the corrigendum of Barnes et al.(2007) (Tables 1 and 2).

    Table 1. Contingency table. TP (true positives), FP (false positives), FN (false negatives), and TN (true negatives).

    Table 2. Description of metrics used for the validation of the model.

    3. Results

    3.1. Mean soundings vs Zonda soundings

    The objective of this work is to detect those vertical profiles that may be able to detect the development of a Zonda event. For the statistical analysis (i.e., PCA), soundings anomalies with respect to the mean sounding of each station are used.

    Figure 3a shows Mendoza’s mean sounding (full line)and the mean sounding for Zonda events on a skew-T diagram. The mean sounding is a statically stable profile throughout the vertical. The vertical profile associated with Zonda is also stable, but in this case, it is less stable (than the mean sounding) due to the lower levels warming. Dry conditions throughout the profile are observed, with a dew point depression of approximately 10°C and higher at midlevels. Surface winds are weak up to 700 hPa from a northwest direction, rotating to the west with a maximum of 45 kt(knot; where 1 kt = 0.51 m s-1) at 300 hPa. The mean sounding associated with Zonda events shows positive temperature anomalies from 850 hPa up to 600 hPa in conjunction with negative dewpoint anomalies, with a maximum dewpoint depression at 700 hPa. Above 500 hPa (approximately the height of the orographic barrier), there are no significant differences with the mean sounding, other than small negative temperature anomalies in the Zonda soundings. Positive wind anomalies are observed from 700 hPa to upper levels. The observed speeds are generally greater than mean conditions, exceeding 60 kt at 300 hPa, while directions are almost the same as in the mean sounding.

    The Santo Domingo station’s mean sounding presents a vertical profile with a quasi-isothermal layer between 1000 hPa and 850 hPa, mainly associated with subsidence due to the presence of the Semipermanent Pacific Anticyclone(Fig. 3b). Near-surface winds are weak and rotating to the NW at 700 hPa with a maximum speed of 60 kt at 200 hPa.The vertical windward profile associated with Zonda events presents a relatively wetter and colder environment throughout the vertical. The surface layer is less stable and the winds are more intense, with values of 45 kt at 500 hPa and a maximum of 95 kt at 200 hPa. The windward mid and upper level winds are greater than those observed on the lee side, indicating the incidence of the mountain waves (Durran, 1990).

    3.2. Index model

    The PCA is primarily used as an exploratory tool in data analysis and for making predictive models. In this work, this analysis is used to characterize the vertical structure of the atmosphere on both sides of the Andes Mountains previous to a Zonda wind event. Also, by implementing a multiple logistic regression model for a binomial function between the response vectorcand the loadings componentsfj, a probability index of Zonda occurrence is defined(probabilistic predictive model) as shown in the methodology.

    For the PCA and the index construction, 116 soundings from the Mendoza and Santo Domingo stations at 1200 UTC are taken, constructing the anomaly matrix for the selected variables and levels. Half of these soundings correspond to Zonda wind events and the other half to random dates, where the presence of Zonda was not recorded (see their selection in methodologies). To obtain the index, an iteration is carried out in the PCA, with 115 soundings and leaving 1 out to evaluate the efficiency of the model for each cut-off value of the index. This procedure is carried out for each date, obtaining the Principal Components (F), the regression coefficients (b0,b1,......,bn, related to the predictors used in the regression), the discriminant sounding (Avector,see section 2.2.3), the standardized score components (Z*s),and the eigenvalues (D) and eigenvectors (Q).

    After performing the PCA, an index value is obtained for each date. Taking into account that the Zonda and non-Zonda events are known in advance, the obtained index values for each date are assigned to each class, and their distribution is analyzed. The model efficiency depends on the selected cut-off value of the index. Let’s firstly consider the model for the Mendoza soundings (lee side) for temperature.

    3.2.1.Model Efficiency

    The index distribution for Zonda (dark) and for non-Zonda (light) classes, using only the vertical profile of temperature at Mendoza, the total efficiency for each cut-off value and a boxplot that indicates the dispersion of the index for the two classes are shown in Figs. 4a, 4b and 4c respectively. An index value equal to one corresponds to the Zonda class and an index value equal to zero corresponds to the non-Zonda class, whiskers represent the 3%-97% interval).If, for example, an index cut-off value of 0.5 is chosen, the number of hits for Zonda events is 80.4%, while the non-Zonda events hits represent 58.4%, and the total efficiency is below 70%. The total efficiency (Fig. 4b) reveals that, for this particular model, the best cut-off point (i.e., the value associated to the maximum efficiency) is located at 0.41, giving a total efficiency of 72.23% (86.44% hits of Zonda events and 53.84% for non-Zonda). Thus, the maximum effectiveness of the index is 72.23%, with an error of 27.77%, divided into 6.77% probability of surprise (i.e.,Zonda cases that present an index value lower than the cutoff and therefore are predicted as non-Zonda) and 21% probability of false alarm (i.e., cases of non-Zonda that present an index value higher than the cut-off, so they are predicted as Zonda events). Those error values can be modified by changing the cut-off value of the index. If this value is closer to one, the non-Zonda surprises are reduced, but Zonda false alarms are augmented. Likewise, if the cut-off value is closer to zero, the Zonda false alarms are reduced and the non-Zonda surprises are augmented.

    Fig. 4. Leeside model for temperature. (a) Index values distribution for Zonda(dark) and for non-Zonda (light) events, (b) model efficiency according to each cut-off value, and (c) index boxplot.

    For example, a cut-off value of 0.78 (50.54% efficiency) yields 1.89% and 47.56% for false alarms and surprises, respectively. Furthermore, looking at the boxplot diagrams (Fig. 4c) and their dispersion, index values higher than 0.78 correspond to the right tail of the 3rd percentile for the non-Zonda cases, while index values lower than 0.2 correspond to the left tail of the 97th percentile for Zonda cases. Therefore, for a given new sounding for which the index is calculated, values greater than 0.78 would indicate a Zonda occurrence, with an error lower than 3%. Likewise,an index value lower than 0.2 would indicate a non-Zonda occurrence with the same error. The range between those index values represents an uncertainty interval in which the index fails to discriminate between classes (at least with the previous error rates).

    The same analysis is done for the windward side using the Santo Domingo vertical soundings and using both soundings at the same time (same variables). Figure 5 shows the index distribution for the zonal component of wind for the leeward (Fig. 5a), windward (Fig. 5b), and mixed (Fig. 5c) models. The windward model efficiency surpasses that for the lee side. However, the windward true positives (82.36%)are considerably lower than those obtained for the lee side(91.53%), but the true negatives and false alarms results are better. If both soundings are combined using the same variables, the model efficiency increases for some variables.The total efficiency (Fig. 5) reaches 83.57% for the zonal component of wind, slightly higher than the windward model. Again, the true positives are lower than for the leeside model, but higher than for the windward model, and false alarms are substantially reduced. Observe that the separation between classes (boxplot) is more evident in the mixed model. This indicates that moving the cut-off value of the index to the left leads to a large increase in true positives with relatively fewer surprises, but there is a higher false alarm rate and vice versa for moving the cut-off value to the right.

    If different variables for each side are combined, a new mixed model is obtained. Figure 6a shows the mixed model using temperature (T) for the lee side andU-wind for the windward side. Here, the model efficiency is maximum(86.23%), and separation between classes is greater. Note that an index value lower than 0.18 leads to 3% of Zonda surprises and a value greater than 0.9 leads to 3% of Zonda false alarms (whiskers are 3%-97%). The total efficiency for all possible combinations with one variable for each side is shown in Fig. 6b. The highest efficiency is achieved using temperature for the lee side and theU-wind for the windward side (a total efficiency of 86.31%). It is worth noting that the zonal component of the wind as a predictor is a key factor on both sides, followed by temperature on the lee side and by the meridional component on the windward side. Although this model increases the total efficiency (i.e.,decreases the total error) and decreases the false alarms, the true positives and surprises results are better when only considering the leeside soundings, and the true negatives results are better for the windward soundings.

    Fig. 5. Same as Fig. 4 but for the zonal component model. (a) Leeside sounding, (b) windward sounding, and (c) mixed model.

    Fig. 6. (a) 1-var mixed model for temperature for the lee side and zonal component of wind for the windward side and (b) the total efficiency for all possible combinations with one variable on each side.

    Table 3 shows the total efficiency, the Zonda hits (true positives), the non-Zonda hits (true negatives), and the total error separated into surprises and false alarms for all variables and models. These results show that theU-wind component is the best variable predictor for all models, with a maximum efficiency of 83.57% for the mixed model. Also,the true positives are maximized for the leeward model with 91.53% and a minimum of surprises of 4.23% (leeside model,U), and true negatives and a minimum of false alarms are best predicted for the mixed model (U).

    Table 3. Model metrics for the best cut-off index value. MDZ, Mendoza’s sounding, CH, Santo Domingo’s sounding and MDZ-CH,mixed model for same variable on each side.

    Table 4. Model metrics for the best-cut value. MDZ, Mendoza’s sounding, CH, Santo Domingo’s sounding, and mixed model for same variable on each side. POD (probability of detection), POFD (probability of false detection), MR (miss ratio), CAR (correct alarm ratio),FAR (false alarm ratio), and MAR (missed alarm ratio).

    3.2.2.Model metrics

    To compare models’ skills, their metrics are shown in Table 4. Following an event-based statistic, with a POD value of 95.58% the leewardUmodel is the best of all. This model also has the lower MR. The POFD result is lower for the mixed modelU. From an alarm-based perspective(which may be more important for a forecasting point of view), the models’ performances show that the mixed model result is better than the others, except for the MR.The False Alarm Ratio (FAR) is below 8%, and the CAR is above 92% for the mixed modelU.

    Table 5. Mixed model metrics for the best combination.

    Figure 7 shows all metrics used for the mixed model,where the dark colors represent an enhancement of the models. The true positives are mainly controlled by the lee side(Fig. 7a). The leeside model surpasses the mixed model forT,Td, andN2. For the mixed model, the true positives are increased usingU-wind, followed by temperature and dewpoint on the lee side. On the windward side, the squared Brunt-V?is?l? frequency, aside fromU, also increases the true positives. To the contrary, the true negatives (Fig. 7d)appear to be controlled by the windward side, in a similar case as the true positives. According to the definition of the metrics, the surprises (Fig. 7g), POD (Fig. 7e), and MR(Fig. 7c) result in similar patterns to those of the true positives and, false alarms (Fig. 7b), and FAR (Fig. 7f), and POFD (Fig. 7h) results in similar patterns to the true negatives. The MAR (Fig. 7i) is a combination of both patterns.

    So, it is clear that the leeside soundings are better predictors for the detection of Zonda events, while the windward side helps to improve false detections and alarms.

    3.2.3.Discriminant sounding

    The discriminant sounding (Avector, see section 2.2.3)between Zonda and non-Zonda classes obtained from the PCA is shown in Fig. 8. For the lee side (Fig. 8a), this vertical sounding is characterized by positive temperature anomalies between 850 hPa and 600 hPa and negative anomalies in upper levels. The low-level temperature increase is accompanied by positive anomalies of dewpoint near the surface and negative anomalies between 800 hPa and 500 hPa. The zonal component of wind shows negative anomalies at low levels with two relative maximums of positive anomalies at 700 hPa and 400 hPa, while the meridional component is the opposite of the previous one (not shown). The squared Brunt-V?is?l? frequency has positive anomalies at low and upper levels with negative anomalies at midlevels.

    Fig. 7. 1-var mixed model metrics: true positives (a) and negatives (d), surprises (g) and false alarms (b), POD (e), POFD (h), MR(c), FAR (f), and MAR (i).

    Fig. 8. (a) Discriminant sounding between Zonda and non-Zonda classes for the 1-var models for the lee side and (b) for the windward side. Since the discriminant soundings are anomalies with respect to the mean profile,temperature and dewpoint temperature variables are multiplied by a factor of five and added to the mean sounding to better appreciate the differences with the mean sounding. The zonal component (U) of the wind is multiplied by a factor of 10 as anomalous values (i.e., the mean sounding is not added). The meridional component is omitted due to its lesser relevance in the models.The squared Brunt-V?is?l? frequency presents the original values of the discriminant sounding.

    This leeside discriminant sounding represents a sounding that tends to become unstable at midlevels due to strong heating at low levels and cooling at midlevels, with greater drying at 700 hPa and moistening near the surface. Above this layer, the sounding tends to normalize in humidity and temperature, approaching the climatological one, but with positive anomalies of zonal wind and more static stability at 400 hPa. According to this profile, a stable layer near the surface is present and the windstorm probably is already present at 700 hPa (the soundings are previous to the event at land level) with a NW component of wind. Also, the upper-level jet streak could be present at this time.

    For the windward model (Fig. 8b), the vertical sounding is characterized by a temperature positive anomaly near the mountain top, between 500 hPa and 400 hPa, with a stable region above this layer. The layer below this region has negative anomalies of the squared Brunt-V?is?l? frequency, withU-wind anomaly maximums at 700 hPa and near the surface with opposite values. The meridional component acts only at low levels (not shown).

    Fig. 9. Total efficiency, true positives and true negatives, false alarms, POD (right axis), POFD, MR, CAR (right axis), FAR,and MAR for the lee side (a, d, g, and j), for the windward side (b, e, h, and k), and for mixed models (c, f, i, and l).

    Fig. 10. Best mixed models considering the total efficiency for two, three, and four variables. (a) 2-var mixed model for temperature and zonal component for each side. (b) 3-var mixed model for temperature, zonal component, and squared Brunt-V?is?l? for both sides. (c) 4-var mixed model for temperature, zonal and meridional components, and squared Brunt-V?is?l? for the lee side and temperature, dewpoint, zonal component, and squared Brunt-V?is?l? for the windward side.

    So, in accordance with the highest efficiency values to discriminate between classes, a vertical sounding that presents positiveU-wind (westerly) anomalies at 700 hPa on the windward side and two maximums of positive anomalies on the lee side are the key predictors for Zonda occurrence. The leesideU-wind anomalies correspond to the presence of Zonda at midlevels and to the jet streak in upper levels. This vertical structure could be accompanied by a temperature inversion (and a stable layer above) near the mountain top windward side and significant drying at midlevels on the lee side.

    3.3. Multi-var models

    3.3.1.Multi-var model efficiency and metrics

    If more variables are added to the model, an efficiency gain is observed, especially for the lee side. Figure 9 (a, d, g and j) shows the best result (i.e., best variables combinations) of the metrics for the leeside models by adding more variables. The total efficiency and the true negatives present a considerable increase from one variable to two, with little variation when adding more than two. The efficiency goes from 80% to 88%, and the true negatives go from 65% to 82%. The true positives seem not to change much by adding more variables and present a maximum for two variables.False alarms are also reduced, but surprises are not. The event-based statistic presents a considerable diminution in the POFD and a rise in the POD between one and two variables, but the MR remains almost constant. From an alarmbased perspective, the CAR and the FAR are improved, and the missing alarm ratio does not present any changes.

    For the windward side (Figs. 9b, e, h and k), the total efficiency and true negatives increase almost linearly by adding more variables, while true negatives sightly increase adding more variables. Surprises decrease 2 point up to four variables, and the false alarms highly improve going from one to two variables. The event-based statistic presents almost no change in the POFD (or even decrease), the POD linearly rises, and the Missing Ratio presents a diminution.From an alarm-based perspective, the False Alarm Ratio and the Missing Alarm Ratio decrease, improving the model. The Correct Alarm Ratio increases 1.5 times going from two to four variables.

    The mixed model (Figs. 9c, f, i, and l) presents its maximums of true negatives and positives for two and three variables, while the total efficiency reaches its maximum for four variables but almost with a small change for two, three,and four variables. In this model, the false alarms seem not to improve by adding variables, as observed for the previous models. From an event-based perspective, the POD increases 4 point (from 94.94% to 98.9%) from one to four variables and the MR also presents an improvement, but the POFD does not. Finally, the MAR is highly reduced, while the CAR and FAR reach their optimal values for three variables.

    It is clear that the mixed models turn out to be significantly better than the separate models. On the other hand,adding a greater number of variables to the mixed model does not always improve the metrics, as is the case of the true positives and negatives, POFD, CAR, and FAR, whose maximums are reached with two or three variables and not with four (Table 5).

    Figure 10 shows the best mixed models (considering the total efficiency) for two, three, and four variables. As can be noted, the separation between classes is more evident (than the 1-var mixed model), causing the uncertainty interval to decrease considerably. The Zonda 25th-75th percentile intervals are reduced to the index interval [0.8-1](the 1-var mixed model is [0.7-0.9]), while for non-Zonda cases, this interval is reduced to [0-0.2] (the 1-var mixed model is [0.0.5-0.45]). This value is greater than the 3rd percentile of the Zonda cases. For the 2-var mixed model(Fig. 10a), an index value lower than 0.19 leads to only 3%surprises for Zonda events. Unlike the 1-var mixed model,the 25th percentile, which was located at 0.7, is above 0.78,and the 75th percentile moves towards a value of 0.98 in the 4-var mixed model (Fig. 10c). In this way, surprises are considerably reduced. The total efficiencies for all two-variable mixed models (figure not shown) clearly show maximums by consideringT-U(89.66%) andT-N2(89.2%) for the lee side combined withT-Ufor the windward side.

    Table 6 shows the best combinations of variables for the leeside, windward, and mixed models with two, three,and four variables. It can be seen that the variablesU,N2,andTandTd(to a lesser extent) are the best classifiers between Zonda and non-Zonda soundings. In addition, note the little difference that exists by adding a greater number of variables. For example, for the mixed model, there is a difference of 1.02 and 1.11 points between the 2-var model and 3-var and 4-var models, respectively. These differences are greater when considering the leeward model, where the differences in efficiency between two and three variables is almost 2.5%, and between three and four variables the difference is 0.5%. The windward models do not seem to be very sensitive to the addition of variables.

    From Table 6, it is also concluded that, regardless of whether there is a single sounding or soundings for both(windward and leeward) for the forecasting day, the efficiency values of the models presented here can be high,with a difference of less than 5% (this is the case between the windward and mixed models). However, the presence of the leeside sounding is of greater importance in obtaining greater precision in the detection of Zonda events and in reducing false alarms.

    Table 6. Best 10 mixed models considering the total efficiency for two and three variables and best 5 mixed models for four variables.First column: leeside 2-, 3- and 4-var best 10 models, with second column for the windward side and third column for the mixed models.

    3.3.2.Discriminant sounding multi-var model

    The combination of soundings on both sides of the Andes produces a slight change in the obtained profiles due to the change in the fit coefficients of the new model. The discriminant sounding between Zonda and non-Zonda classes obtained from the PCA is shown in Fig. 11 for the 4-var mixed model.

    This model presents a leeside (Fig. 11a) vertical structure characterized by positive temperature anomalies between 850 hPa and 650 hPa and a minimum of negative anomalies at midlevels, which tends to represent an unstable profile. The vertical structure of dewpoint presents negative anomalies at 700 hPa and a maximum of positive anom-alies near the surface. The vertical structure ofUanomalies correlates well with the presence of vertically (backward)propagating waves, with an alternating pattern of negative and positive anomalies, with a maximum of negative anomalies at 850 hPa and positive anomalies at 400 hPa. The vertical sounding is accompanied by a stable layer at 400 hPa and a statically more unstable layer at midlevels.

    The windward (Fig. 11b) discriminant sounding presents positive anomalies at 850 hPa and negative anomalies in upper levels. Similarly to the lee side, this sounding becomes unstable. This vertical sounding does not present a temperature inversion like the 1-var model (Fig. 8b). Positive anomalies of dewpoint can be seen near the surface and at 500 hPa, while the negative anomalies are observed at 800 hPa and at 400 hPa. TheU-wind anomalies present strong maximums at 700 hPa and near the surface, with opposite values. Like for the lee side, positive and negative changes in the zonal component anomalies can be seen and the vertical sounding is accompanied by a stable layer above 450 hPa and a statically more unstable layer at 500 hPa.

    Fig. 11. Discriminant sounding between Zonda and non-Zonda classes for the 4-var model (a) for the lee side and (b) for windward side.

    4. Conclusions and discussion

    This paper represents an improvement over work previously carried out. In Otero et al. (2018), soundings were used forTandTdalong with surface values for the same variables only for the lee side. The development of the Zonda phenomenon is closely related to the vertical structure of the atmosphere, where the mountain height, the buoyancy frequency, and incident wind control the mountain wave activity and the dynamics of downslope wind, which the previous work did not consider. Here, the characteristics on both sides of the barrier of stability, wind, temperature, and humidity that are more favorable for this type of event could be known. Thus, in this new paper, progress in understanding the characteristics of the atmosphere, both windward and leeward, prior to the development of a Zonda event, is made.The efficiency values obtained in this new work are even higher, without the need of surface data. Statistical values(event-based and alarm-based model metrics) associated with the model were also recognized and can provide additional information for decision makers.

    The results show that this index model is able to discriminate between Zonda and non-Zonda events with an effectiveness close to 80% using only one variable on one side of the barrier, increasing to almost 84% by using one variable on each side of the barrier. These values of effectiveness in class separation increase to almost 91% using four variables on each side of the Andes. The number of Zonda hits is mainly captured by the leeside soundings, while the non-Zonda hits are mostly controlled by the windward soundings. The soundings combination in the mixed model produces better class separation and significantly reduces false alarms. The correlation values between the real index and the estimated index increase considerably when adding variables, reaching values of more than 0.84.

    It is found that the zonal component of the wind (Uwind) on both sides and the windward temperature are the key variables in class discrimination. The true positives and surprises are mainly controlled by the leeside sounding,while the POFD improves with the mixed model (i.e., by adding the windward sounding). Even though mixed model maximum efficiency, surprises, POD, MR, and MAR turn out to be with four variables on each side, the true positives and negatives, POFD, CAR, and FAR reach their maximums with two or three variables.

    The vertical profile that best succeeds to discriminate between Zonda and non-Zonda is characterized by a leeside sounding that tends to become unstable at midlevels and presents two wind maxima, possibly associated with the foehn effect at low levels and the presence of the jet streak at high levels. On the windward side, there is also a wind maximum at 700 hPa, accompanied by a relatively more stable layer near the top of the barrier.

    The methodology used in this work could be replicated in any region that has sounding data (preferable on both sides of the barrier). Although the results obtained may vary from region to region, it could represent a valuable tool for this type of phenomena that are so difficult to forecast.After determining the discriminant vectorAfrom the PCA and the multiple regression model, the probability index for Zonda/non-Zonda occurrence can be easily obtained to forecast the Zonda occurrence up to 24 hours in advance.

    This work ensures, in a preliminary way, that the PCA is a useful tool in the detection of the Zonda phenomenon.However, using more data (such as surface values or reanalysis) could significantly improve the efficiency of the model. In addition, the monthly anomalous soundings were tested, but the model did not significantly improve, so it was decided not to consider separate seasons. The seasonality could affect the frequency of events but not the atmospheric conditions (synoptic and local) under which they originate.

    久久这里有精品视频免费| 一区在线观看完整版| 丰满人妻一区二区三区视频av| 亚洲欧洲精品一区二区精品久久久 | 香蕉精品网在线| 免费看不卡的av| 欧美最新免费一区二区三区| 国产精品一区二区性色av| 夫妻性生交免费视频一级片| 日韩熟女老妇一区二区性免费视频| 黄色日韩在线| 国产有黄有色有爽视频| 如何舔出高潮| 国产高清三级在线| 精品久久久噜噜| 精品视频人人做人人爽| 国产黄色视频一区二区在线观看| 美女xxoo啪啪120秒动态图| 观看美女的网站| 日韩亚洲欧美综合| 最近2019中文字幕mv第一页| 国产亚洲av片在线观看秒播厂| 国产又色又爽无遮挡免| 黑丝袜美女国产一区| 亚洲av欧美aⅴ国产| 成人黄色视频免费在线看| 有码 亚洲区| 国产成人a∨麻豆精品| 精品卡一卡二卡四卡免费| 韩国av在线不卡| 一级av片app| 又黄又爽又刺激的免费视频.| 国产伦理片在线播放av一区| 国产成人一区二区在线| 亚洲av二区三区四区| 精品熟女少妇av免费看| 欧美日韩一区二区视频在线观看视频在线| 亚洲国产精品一区二区三区在线| 91久久精品电影网| 这个男人来自地球电影免费观看 | 女的被弄到高潮叫床怎么办| 久久精品国产a三级三级三级| 久久精品熟女亚洲av麻豆精品| 高清视频免费观看一区二区| 最近中文字幕2019免费版| 精品人妻偷拍中文字幕| 欧美日韩视频精品一区| 国产伦精品一区二区三区视频9| av又黄又爽大尺度在线免费看| 久久6这里有精品| 妹子高潮喷水视频| 大香蕉久久网| 亚洲av男天堂| 国语对白做爰xxxⅹ性视频网站| 免费黄网站久久成人精品| 人人妻人人爽人人添夜夜欢视频 | 人妻人人澡人人爽人人| 欧美精品高潮呻吟av久久| 丝袜在线中文字幕| 国产毛片在线视频| 国产熟女午夜一区二区三区 | 啦啦啦啦在线视频资源| 人妻一区二区av| 久久久亚洲精品成人影院| 国产男人的电影天堂91| 男男h啪啪无遮挡| 欧美日韩亚洲高清精品| 大香蕉久久网| 国产伦精品一区二区三区视频9| 国产精品久久久久久av不卡| 黑人巨大精品欧美一区二区蜜桃 | 亚洲人成网站在线观看播放| 国产精品无大码| 久久综合国产亚洲精品| 精品一区二区三区视频在线| 纯流量卡能插随身wifi吗| 成年av动漫网址| 内地一区二区视频在线| 日本欧美视频一区| 亚洲国产精品999| 老司机影院毛片| 亚洲一区二区三区欧美精品| 少妇人妻 视频| 男女边吃奶边做爰视频| 日日摸夜夜添夜夜爱| 下体分泌物呈黄色| 国产无遮挡羞羞视频在线观看| 日本91视频免费播放| 亚洲性久久影院| 国产成人精品久久久久久| 亚洲av.av天堂| 2021少妇久久久久久久久久久| 午夜影院在线不卡| 美女大奶头黄色视频| 国国产精品蜜臀av免费| 性色av一级| 中国美白少妇内射xxxbb| 日本与韩国留学比较| 日本av手机在线免费观看| 晚上一个人看的免费电影| 国产在线免费精品| 自拍欧美九色日韩亚洲蝌蚪91 | 国产精品偷伦视频观看了| 久久99热这里只频精品6学生| 纯流量卡能插随身wifi吗| av福利片在线观看| 日本猛色少妇xxxxx猛交久久| 国产美女午夜福利| 日韩成人伦理影院| 少妇被粗大猛烈的视频| av卡一久久| 精品亚洲乱码少妇综合久久| 91久久精品电影网| 观看美女的网站| 熟妇人妻不卡中文字幕| 亚洲国产精品一区二区三区在线| 一级二级三级毛片免费看| 免费看日本二区| 国产精品一区二区三区四区免费观看| 99久久中文字幕三级久久日本| 亚洲四区av| 国产欧美日韩精品一区二区| 日韩中字成人| 欧美 亚洲 国产 日韩一| 日韩电影二区| 日本午夜av视频| 亚洲婷婷狠狠爱综合网| 99热这里只有是精品50| 在线观看免费日韩欧美大片 | 欧美日韩亚洲高清精品| 高清黄色对白视频在线免费看 | 女人久久www免费人成看片| 色吧在线观看| 亚洲国产色片| 久久精品国产亚洲av涩爱| 又大又黄又爽视频免费| 久久久久精品性色| 亚洲久久久国产精品| 日本wwww免费看| 最近手机中文字幕大全| 99久久精品一区二区三区| 日本与韩国留学比较| av卡一久久| av在线老鸭窝| 国产真实伦视频高清在线观看| 大片免费播放器 马上看| 久久久国产欧美日韩av| 只有这里有精品99| 91精品伊人久久大香线蕉| 少妇精品久久久久久久| 91午夜精品亚洲一区二区三区| 18禁在线播放成人免费| 国产亚洲av片在线观看秒播厂| 国内少妇人妻偷人精品xxx网站| 久久影院123| 人妻系列 视频| 三级国产精品片| 欧美国产精品一级二级三级 | 天堂俺去俺来也www色官网| 久久久久久久国产电影| 麻豆成人av视频| 久久久亚洲精品成人影院| 午夜免费鲁丝| 乱码一卡2卡4卡精品| 国产女主播在线喷水免费视频网站| 美女xxoo啪啪120秒动态图| 亚洲av在线观看美女高潮| av.在线天堂| 久久午夜综合久久蜜桃| av在线老鸭窝| 五月伊人婷婷丁香| 男人和女人高潮做爰伦理| 91精品一卡2卡3卡4卡| 99热国产这里只有精品6| 婷婷色综合www| 免费看不卡的av| 精品久久久久久久久亚洲| 午夜福利影视在线免费观看| 麻豆成人av视频| 亚洲av男天堂| 亚洲精品中文字幕在线视频 | 亚洲精品久久午夜乱码| 丁香六月天网| 3wmmmm亚洲av在线观看| 少妇人妻精品综合一区二区| 夜夜看夜夜爽夜夜摸| av天堂久久9| 亚洲av在线观看美女高潮| 在线观看av片永久免费下载| 嫩草影院新地址| 亚洲欧美中文字幕日韩二区| 亚洲欧美精品专区久久| 亚洲伊人久久精品综合| 97超视频在线观看视频| 国产免费一区二区三区四区乱码| 国产爽快片一区二区三区| 另类亚洲欧美激情| 天天操日日干夜夜撸| 国产亚洲午夜精品一区二区久久| 久久精品国产亚洲av天美| 国产女主播在线喷水免费视频网站| 一级毛片久久久久久久久女| 男人添女人高潮全过程视频| 极品人妻少妇av视频| 亚洲欧美一区二区三区国产| 你懂的网址亚洲精品在线观看| 亚洲精品国产色婷婷电影| 免费黄网站久久成人精品| 国产有黄有色有爽视频| 亚州av有码| 99九九在线精品视频 | 亚洲中文av在线| 免费人妻精品一区二区三区视频| 午夜激情久久久久久久| 欧美 亚洲 国产 日韩一| 久久久久久久亚洲中文字幕| 插逼视频在线观看| 七月丁香在线播放| 最近手机中文字幕大全| 一个人免费看片子| 肉色欧美久久久久久久蜜桃| 亚州av有码| 女性被躁到高潮视频| 久热这里只有精品99| 日韩人妻高清精品专区| 搡老乐熟女国产| 亚洲精品中文字幕在线视频 | 毛片一级片免费看久久久久| 国产伦理片在线播放av一区| 久久人人爽人人爽人人片va| 亚洲图色成人| 亚洲国产成人一精品久久久| 亚洲av不卡在线观看| 蜜桃在线观看..| 男女边摸边吃奶| 极品人妻少妇av视频| 人人澡人人妻人| 少妇熟女欧美另类| 人妻 亚洲 视频| 制服丝袜香蕉在线| www.色视频.com| 亚洲高清免费不卡视频| 午夜精品国产一区二区电影| 久久久久久久久久成人| 婷婷色麻豆天堂久久| 国产亚洲午夜精品一区二区久久| 91精品国产国语对白视频| 蜜桃久久精品国产亚洲av| av有码第一页| 多毛熟女@视频| 国产精品久久久久久精品古装| 久久久国产精品麻豆| 国产精品免费大片| 国产乱人偷精品视频| 2022亚洲国产成人精品| 国产免费一区二区三区四区乱码| 欧美+日韩+精品| 十分钟在线观看高清视频www | 国产免费又黄又爽又色| 欧美bdsm另类| 麻豆乱淫一区二区| 亚洲第一区二区三区不卡| 女人精品久久久久毛片| 亚洲国产精品成人久久小说| 在线亚洲精品国产二区图片欧美 | 国产精品福利在线免费观看| 九色成人免费人妻av| av女优亚洲男人天堂| 久久狼人影院| 国产女主播在线喷水免费视频网站| 日本猛色少妇xxxxx猛交久久| 免费观看的影片在线观看| 久久久久久久久久久丰满| 人人妻人人爽人人添夜夜欢视频 | 亚洲精品久久久久久婷婷小说| 亚洲精品乱久久久久久| 2021少妇久久久久久久久久久| 精品久久久噜噜| 男人舔奶头视频| 一级毛片久久久久久久久女| 亚洲精品一区蜜桃| 天天操日日干夜夜撸| 99九九在线精品视频 | 中文字幕精品免费在线观看视频 | 亚洲国产精品999| 午夜激情久久久久久久| 一级爰片在线观看| 少妇被粗大猛烈的视频| 高清欧美精品videossex| 99久久精品一区二区三区| 午夜av观看不卡| 国精品久久久久久国模美| 久久精品国产亚洲av天美| 欧美xxⅹ黑人| 午夜激情久久久久久久| 三级国产精品片| 五月开心婷婷网| 国产欧美另类精品又又久久亚洲欧美| 狂野欧美白嫩少妇大欣赏| 成年人免费黄色播放视频 | 日韩亚洲欧美综合| 妹子高潮喷水视频| 狠狠精品人妻久久久久久综合| 在线看a的网站| 伦理电影免费视频| 亚洲精品日韩在线中文字幕| 日本av手机在线免费观看| 国产乱人偷精品视频| 婷婷色麻豆天堂久久| 亚洲婷婷狠狠爱综合网| 美女中出高潮动态图| 婷婷色麻豆天堂久久| 午夜福利视频精品| 免费看日本二区| 久久狼人影院| 人妻夜夜爽99麻豆av| 亚洲av成人精品一二三区| 性色avwww在线观看| 成年女人在线观看亚洲视频| 亚洲一区二区三区欧美精品| 国产成人a∨麻豆精品| 欧美97在线视频| 亚洲真实伦在线观看| 大片免费播放器 马上看| 久热这里只有精品99| 日韩不卡一区二区三区视频在线| 亚洲国产成人一精品久久久| 国产欧美日韩综合在线一区二区 | 亚洲,欧美,日韩| 亚洲欧洲精品一区二区精品久久久 | 少妇人妻精品综合一区二区| 亚洲精品乱码久久久久久按摩| 麻豆成人午夜福利视频| 欧美精品人与动牲交sv欧美| 国产成人aa在线观看| 久久99一区二区三区| 亚洲va在线va天堂va国产| 欧美最新免费一区二区三区| 在线观看免费高清a一片| 精品人妻熟女毛片av久久网站| av天堂久久9| 另类亚洲欧美激情| 中国国产av一级| 免费看不卡的av| 激情五月婷婷亚洲| kizo精华| 一级二级三级毛片免费看| 免费看不卡的av| 激情五月婷婷亚洲| 亚洲精华国产精华液的使用体验| 免费av中文字幕在线| 大又大粗又爽又黄少妇毛片口| 国产精品国产三级国产av玫瑰| 97在线人人人人妻| 七月丁香在线播放| 成人美女网站在线观看视频| 一本—道久久a久久精品蜜桃钙片| 十八禁高潮呻吟视频 | 精品一品国产午夜福利视频| 国产精品免费大片| 午夜免费男女啪啪视频观看| 国产视频首页在线观看| 中国美白少妇内射xxxbb| 久久影院123| 欧美老熟妇乱子伦牲交| 日本av手机在线免费观看| 一个人看视频在线观看www免费| 少妇裸体淫交视频免费看高清| 人人妻人人爽人人添夜夜欢视频 | 秋霞伦理黄片| 免费av中文字幕在线| 两个人免费观看高清视频 | 国产成人精品福利久久| 99久久综合免费| av福利片在线观看| 春色校园在线视频观看| 我的老师免费观看完整版| 国产精品嫩草影院av在线观看| 欧美日韩一区二区视频在线观看视频在线| 在现免费观看毛片| 国内少妇人妻偷人精品xxx网站| 熟女人妻精品中文字幕| 国产极品天堂在线| 嘟嘟电影网在线观看| 国产熟女午夜一区二区三区 | 性高湖久久久久久久久免费观看| 午夜91福利影院| 插逼视频在线观看| 男人爽女人下面视频在线观看| 成人毛片a级毛片在线播放| 亚洲欧美中文字幕日韩二区| 婷婷色麻豆天堂久久| 国产精品国产三级国产专区5o| 精品人妻偷拍中文字幕| 最近手机中文字幕大全| 婷婷色综合www| 国产日韩欧美在线精品| 国产在线一区二区三区精| 边亲边吃奶的免费视频| 老司机影院成人| 国产精品久久久久久久电影| 美女xxoo啪啪120秒动态图| 国产伦精品一区二区三区四那| 国产亚洲91精品色在线| 久久精品国产a三级三级三级| 国产熟女午夜一区二区三区 | av不卡在线播放| av天堂中文字幕网| 中国国产av一级| 成年女人在线观看亚洲视频| 亚洲国产精品专区欧美| 麻豆乱淫一区二区| 亚洲精品aⅴ在线观看| 啦啦啦视频在线资源免费观看| 国产伦精品一区二区三区四那| 久久久久精品久久久久真实原创| 国产精品一二三区在线看| 国产有黄有色有爽视频| 视频中文字幕在线观看| 久久国产精品男人的天堂亚洲 | 日本黄色日本黄色录像| 自拍欧美九色日韩亚洲蝌蚪91 | 亚洲精品视频女| 欧美日韩国产mv在线观看视频| 国产伦在线观看视频一区| 精品久久久久久久久av| 精品亚洲乱码少妇综合久久| av免费观看日本| 久久99热6这里只有精品| 国内少妇人妻偷人精品xxx网站| 99久久精品国产国产毛片| av天堂中文字幕网| 日日爽夜夜爽网站| 国产精品99久久久久久久久| 欧美国产精品一级二级三级 | 欧美国产精品一级二级三级 | 成人黄色视频免费在线看| 亚洲成人手机| 国产成人一区二区在线| 久久久精品94久久精品| 免费黄频网站在线观看国产| 国内揄拍国产精品人妻在线| 国产日韩欧美视频二区| 国内少妇人妻偷人精品xxx网站| 国内精品宾馆在线| 看十八女毛片水多多多| 欧美另类一区| 成人二区视频| 亚洲精品自拍成人| 精品国产国语对白av| 亚洲欧美成人精品一区二区| 久久久久视频综合| av在线app专区| 久久久久国产精品人妻一区二区| 久久99热这里只频精品6学生| a级一级毛片免费在线观看| 22中文网久久字幕| 亚洲国产精品专区欧美| 国产精品久久久久久久久免| 汤姆久久久久久久影院中文字幕| 国产老妇伦熟女老妇高清| 亚洲精品日韩av片在线观看| 有码 亚洲区| 成年av动漫网址| 观看av在线不卡| 水蜜桃什么品种好| 看十八女毛片水多多多| 亚洲欧洲精品一区二区精品久久久 | 天堂中文最新版在线下载| 高清不卡的av网站| 国产伦理片在线播放av一区| 国产精品免费大片| 精品久久久噜噜| 免费不卡的大黄色大毛片视频在线观看| 蜜桃久久精品国产亚洲av| 亚洲精品日韩av片在线观看| 韩国av在线不卡| 国产精品久久久久久久电影| 卡戴珊不雅视频在线播放| 久久毛片免费看一区二区三区| 一级a做视频免费观看| 日韩视频在线欧美| 精品国产一区二区久久| 美女视频免费永久观看网站| 精品人妻偷拍中文字幕| 美女cb高潮喷水在线观看| 人体艺术视频欧美日本| 久久99热这里只频精品6学生| 少妇的逼水好多| 国产精品国产三级国产av玫瑰| 如何舔出高潮| 亚洲av国产av综合av卡| 国产精品嫩草影院av在线观看| 色94色欧美一区二区| 一区二区三区四区激情视频| 国产伦在线观看视频一区| 色视频在线一区二区三区| 亚洲国产毛片av蜜桃av| 日日撸夜夜添| 熟妇人妻不卡中文字幕| 国产一区二区在线观看日韩| 久久99热6这里只有精品| 日韩av不卡免费在线播放| 在线 av 中文字幕| 日本-黄色视频高清免费观看| 嫩草影院入口| 女人久久www免费人成看片| 国精品久久久久久国模美| 人妻一区二区av| 丰满迷人的少妇在线观看| 午夜免费男女啪啪视频观看| a级毛片免费高清观看在线播放| 欧美97在线视频| 亚洲精品色激情综合| 久久精品国产鲁丝片午夜精品| 婷婷色综合大香蕉| 三级国产精品片| 男女免费视频国产| 久久97久久精品| 岛国毛片在线播放| 在线观看国产h片| 欧美97在线视频| 搡女人真爽免费视频火全软件| 婷婷色av中文字幕| 亚洲av中文av极速乱| 久久久久网色| 国产av码专区亚洲av| 免费观看无遮挡的男女| 久久这里有精品视频免费| www.色视频.com| 亚洲av在线观看美女高潮| 久久久久精品久久久久真实原创| 国产欧美日韩一区二区三区在线 | 男人添女人高潮全过程视频| 精品国产一区二区久久| 婷婷色av中文字幕| 91午夜精品亚洲一区二区三区| 黄色视频在线播放观看不卡| 国产精品一区二区性色av| 亚洲高清免费不卡视频| 亚洲在久久综合| 日韩强制内射视频| 久久婷婷青草| 大码成人一级视频| 男人添女人高潮全过程视频| 日韩大片免费观看网站| 美女中出高潮动态图| 九色成人免费人妻av| 国产一区有黄有色的免费视频| 精品国产一区二区久久| 乱人伦中国视频| 免费播放大片免费观看视频在线观看| 久久久久久人妻| www.色视频.com| 狠狠精品人妻久久久久久综合| 亚洲熟女精品中文字幕| 国产精品久久久久久久久免| 99热国产这里只有精品6| 久久精品国产亚洲网站| 99热这里只有是精品50| 日韩中文字幕视频在线看片| 亚洲三级黄色毛片| 在线观看三级黄色| 精品久久久噜噜| 国产成人精品婷婷| h视频一区二区三区| 桃花免费在线播放| 久久久久久久精品精品| 久久av网站| 我要看日韩黄色一级片| 精品人妻熟女毛片av久久网站| 亚洲色图综合在线观看| 女性生殖器流出的白浆| 日本av手机在线免费观看| 97超碰精品成人国产| √禁漫天堂资源中文www| 国产在线男女| 国产精品人妻久久久影院| 又黄又爽又刺激的免费视频.| 日韩中文字幕视频在线看片| 国产日韩欧美亚洲二区| 国产一区二区在线观看av| 国产精品麻豆人妻色哟哟久久| 啦啦啦在线观看免费高清www| 久久狼人影院| 亚洲av电影在线观看一区二区三区| av卡一久久| h日本视频在线播放| 熟女人妻精品中文字幕| 少妇人妻精品综合一区二区| 亚洲av在线观看美女高潮| 菩萨蛮人人尽说江南好唐韦庄| 午夜影院在线不卡| 两个人免费观看高清视频 | 春色校园在线视频观看| 永久网站在线| 亚洲av中文av极速乱| 免费看av在线观看网站| 精品久久久久久电影网| 亚洲三级黄色毛片| 成人漫画全彩无遮挡| 久久久国产精品麻豆| 免费观看的影片在线观看| 波野结衣二区三区在线| 国产综合精华液| 岛国毛片在线播放| 亚洲综合精品二区| 美女大奶头黄色视频| 在线观看www视频免费| 少妇丰满av| tube8黄色片| 国国产精品蜜臀av免费| 夜夜爽夜夜爽视频| 免费av中文字幕在线| 久久鲁丝午夜福利片|