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

    Bias Correction and Ensemble Projections of Temperature Changes over Ten Subregions in CORDEX East Asia

    2020-10-15 10:09:28ChenweiSHENQingyunDUANChiyuanMIAOChangXINGXueweiFANYiWUandJingyaHAN
    Advances in Atmospheric Sciences 2020年11期

    Chenwei SHEN, Qingyun DUAN, Chiyuan MIAO, Chang XING, Xuewei FAN, Yi WU, and Jingya HAN

    State Key Laboratory of Earth Surface Processes and Resource Ecology, Faculty of Geographical Science,Beijing Normal University, Beijing 100875, China

    ABSTRACT Regional climate models (RCMs) participating in the Coordinated Regional Downscaling Experiment (CORDEX)have been widely used for providing detailed climate change information for specific regions under different emissions scenarios. This study assesses the effects of three common bias correction methods and two multi-model averaging methods in calibrating historical (1980?2005) temperature simulations over East Asia. Future (2006?49) temperature trends under the Representative Concentration Pathway (RCP) 4.5 and 8.5 scenarios are projected based on the optimal bias correction and ensemble averaging method. Results show the following: (1) The driving global climate model and RCMs can capture the spatial pattern of annual average temperature but with cold biases over most regions, especially in the Tibetan Plateau region. (2) All bias correction methods can significantly reduce the simulation biases. The quantile mapping method outperforms other bias correction methods in all RCMs, with a maximum relative decrease in root-mean-square error for five RCMs reaching 59.8% (HadGEM3-RA), 63.2% (MM5), 51.3% (RegCM), 80.7% (YSU-RCM) and 62.0% (WRF). (3)The Bayesian model averaging (BMA) method outperforms the simple multi-model averaging (SMA) method in narrowing the uncertainty of bias-corrected results. For the spatial correlation coefficient, the improvement rate of the BMA method ranges from 2% to 31% over the 10 subregions, when compared with individual RCMs. (4) For temperature projections, the warming is significant, ranging from 1.2°C to 3.5°C across the whole domain under the RCP8.5 scenario. (5) The quantile mapping method reduces the uncertainty over all subregions by between 66% and 94%.

    Key words: CORDEX-EA, bias correction, BMA, temperature projection

    1. Introduction

    Climate change has attracted much attention as its effects on human society and ecological environments have grown over past decades (Grimm et al., 2013; Sun et al.,2015, 2016, 2020; Huang et al., 2016). According to the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report, several independent datasets show that the global average surface temperature increased by 0.65°C?1.06°C from 1880 to 2012 (IPCC, 2013). A set of global climate model (GCM) simulations shows that, relative to current climate (1986?2015), by the end of the 21st century(2081?2100), global average surface temperature further increases by from 0.3°C to 4.8°C under different scenarios[from Representative Concentration Pathway (RCP) 2.6 to RCP8.5] (IPCC, 2014). Accompanied by further increases in surface temperature, obvious impacts on natural and anthropogenic systems around the world have been reported, such as diminishing levels of snow and ice, increasing sea level and more occurrences of climate extremes (e.g., heat waves,droughts and floods) (Mann and Gleick, 2015; Schlaepfer et al., 2017; Kang and Eltahir, 2018). To evaluate climate change impacts and mitigate risks of climate changes, reli-able projections of future climate change are crucial for decision makers in government, nongovernmental organizations and the general public, especially for regions vulnerable to the adverse effects of climate change. For example,with precipitation and temperature being the main driving variables for hydrological processes, climate change will affect the hydrological cycle (Luo et al., 2018). Since water is an essential resource, hydrological changes will influence flood control policies, hydropower production management and agricultural irrigation adjustment. To investigate the changes in the hydrological cycle, reliable climate change information from regional climate models (RCMs)and GCMs is required to drive the hydrological models(Teutschbein and Seibert, 2010; Mearns et al., 2015).

    In past decades, GCMs have been widely used to investigate the mechanisms of climate changes (Gulizia and Camilloni, 2015; Eyring et al., 2016; Soden et al., 2018). Several studies show that GCMs can satisfactorily reproduce the spatial pattern and variability of historical temperature (Miao et al., 2014; McSweeney et al., 2015; Tang et al., 2016; Ashfaq et al., 2017; Ruan et al., 2019). However, GCMs are less skillful in regional climate factor simulations due to their coarse resolution (ranging from 100 to 200 km) (Bao et al.,2015; Li et al., 2018a). The dynamical downscaling of outputs from GCMs by RCMs can reproduce regional climate information that is more detailed and reliable. Compared with GCMs, RCMs can provide higher-resolution, moredetailed local information such as topography, land-use categories and soil moisture data, which are important in simulating regional climate information (Fulakeza et al., 2002;Salzmann et al., 2007; Akhtar et al., 2009; Jones and Brunsell, 2009; Halder et al., 2016). Although RCMs can alleviate the deficiencies in GCMs, the simulations of RCMs also greatly depend on (1) the quality of initial conditions, (2) lateral and boundary conditions provided by GCMs used to drive the RCMs and (3) simplified physical parameterization schemes for subgrid-scale physical processes (Au-Yeung and Chan, 2012; Rocheta et al., 2017). As a result,various techniques, such as bias correction and ensemble post-processing methods, have been developed to remove the systemic biases in RCM simulations and give better projections of future climate. The Coordinated Regional Downscaling Experiment (CORDEX) is a program sponsored by the World Climate Research Programme (Giorgi et al.,2009). It was established to provide an international coordinated downscaling framework for advancing RCM development, evaluation and applications (Gutowski et al., 2016).Within CORDEX, RCM ensembles have been created for multiple regions throughout the world (through dynamic downscaling driven by CMIP5 GCMs). Since the influence of climate change varies across regions and time scales,regional climate studies are necessary for correctly detecting climate change signals. East Asia is one of the most vulnerable regions to climate change, since it is a large domain that comprises diverse terrestrial features and complex climate systems (Li et al., 2018b; Miao et al., 2019; Zheng et al., 2019). The main climate system of East Asia is a monsoonal system, and the monsoons are always accompanied by extreme events such as heat waves, typhoons, droughts and floods (Ding and Chan, 2005; Chang et al., 2012; Lee et al., 2017). Observations show that since the middle of the 20th century, the average surface temperatures and frequencies of heat waves have increased for most regions of East Asia (Hijioka et al., 2014; Zhou et al., 2016).

    Focusing on the CORDEX East Asia domain (CORDEX-EA), five RCMs have been employed to provide ensemble simulations of regional climate information. Several studies have examined the performance of RCMs in simulating temperature and the projections of future temperature over CORDEX-EA. Kim et al. (2016) analyzed the spatial pattern of projected temperature data from RCMs participating in studies of the CORDEX-EA domain. The results showed that there will be a warming of 1°C?3°C over the whole domain by the year 2050 and the temperature will increase more at high latitude. Park et al. (2016) focused on the performance of RCMs participating in CORDEX-EA research in simulating summer temperature means and extremes.They found that, compared to the Asian Precipitation-Highly-Resolved Observational Data Integration Towards Evaluation dataset, the RCMs show systematic biases in seasonal means and the simulations of temperature means are more accurate than those of the extremes. Gu et al. (2018)showed that the RCMs have improved model performances as compared to the raw GCM outputs, and the projected trends of the RCM temperatures are increasing from the ensemble mean by around 1°C yr?1over the entire domain.Further, there have also been studies focused on the bias correction and ensemble calibration of RCM simulations over CORDEX-EA. Ngai et al. (2017) showed that RCM bias is comparable to the bias of the driving GCMs, and the bias correction methods can substantially reduce the bias in simulating historical temperature. Kim and Suh (2013) used a Bayesian model averaging (BMA) method to calibrate the probabilistic predictions of temperature. They found that BMA outperforms the equal-weighted method and other ensemble calibration methods in calibrating seasonal means and distributions of the simulations. In addition, the BMA forecasts outperform the single-RCM forecasts.

    Bias correction methods and ensemble calibration methods have been used in the aforementioned studies to improve the simulation of historical temperature and improve the reliability of temperature projections. However,little attention has been paid to the performance of different bias correction methods in improving the accuracy of RCM simulations, and comparison of their ability in narrowing future projections over the CORDEX-EA region is insufficient as well.

    In this study, we evaluated the effects of three common bias correction methods on improving historical temperature simulations, and then used the most effective bias correction method together with an ensemble method to narrow the uncertainties in temperature projections. The changes in future climate were analyzed based on the most reliable projection of temperatures. The paper is organized as follows:Section 2 describes the observational data, the climate model data, and the bias correction and ensemble averaging methods. Section 3 presents the design of the bias correction experiment, including the selection of calibration and validation cases and the evaluation of bias correction methods,as well as the results of ensemble calibration. Finally, the future temperature changes under the RCP4.5 and RCP8.5 scenarios are also discussed, based on the projection period(2030?49) and reference period (1980?99). A summary and conclusions are provided in section 4.

    2. Data and methods

    2.1. Observations and model setup

    The reference data of monthly temperature over CORDEX-EA were taken from the Climatic Research Unit Time-Series (CRU TS) 4.03 dataset developed by the University of East Anglia (available at http://data.ceda.ac.uk/badc/cru/data/cru_ts/cru_ts_4.03/data/tmp), with a resolution of 0.5° (Harris et al., 2014) and a time frame of 1901-2018. To conduct bias correction, only the data from 1980 to 2005 were used in the study, to coincide with the historic period of CORDEX-EA.

    Five RCMs were used in the CORDEX-EA experiment: the Hadley Centre Global Environmental Model, version 3, with Regional Atmosphere configurations(HadGEM3-RA); the Fifth-generation Pennsylvania State-National Center for Atmospheric Research Mesoscale Model (MM5); the Weather Research and Forecasting(WRF) model; the Regional Climate Model, version 4(RegCM4); and the Yonsei University Regional Climate Model (YSU-RCM). The selected RCMs include three nonhydrostatic models (HadGEM3-RA, MM5 and WRF) and two hydrostatic models (RegCM4 and YSU-RCM) (von Storch et al., 2000; Cha et al., 2008; Giorgi et al., 2012;Baek et al., 2013; Wang et al., 2013). Table 1 lists the detailed configurations of the five RCMs, including dynamics processes, physical parameterization schemes, and spectral nudging (Gu et al., 2018). The CORDEX-EA domain covers East Asia, India, South Asia, and the northern part of Australia (Fig. 1), and the spatial resolution is 50 km (except HadGEM3-RA, whose resolution is 0.44°). The historical experiment (1980-2005) and the future projections under the RCP4.5 and RCP8.5 scenarios (2006?49) are driven by the outputs of HadGEM2-AO (Hadley Centre Global Environmental Model, version 2, with Atmosphere and Ocean and sea ice configurations), whose horizontal resolution is 1.875° × 1.25°. Several studies have confirmed the good performance of HadGEM2-AO for simulating East Asia’s climatology (Martin et al., 2011; Baek et al., 2013; Sperber et al.,2013).

    To further correct the biases at smaller spatial scales,10 subregions were selected: northwestern China (NW;36°?43°N, 75°?103°E); the Tibetan Plateau (TP; 28°?35°N, 75°?103°E); northeastern China (NE; 30°?42°N,104°?121°E), northern China (NC; 42°?55°N, 113°?132°E);southern China (SC; 18°?30°N, 104°?122°E); the Korean Peninsula and Japan (KJ; 43°?51°N, 91°?112°E); Mongolia (MG; 30°?42°N, 125°?141°E); India (5°?27°N, 69°?91°E), Indochina (InC; 8°?28°N, 92°?110°E), and Southeast Asia (SEA; 10°S?6°N, 95°?151°E) (Zou and Zhou,2016; Zhou et al., 2016; Li et al., 2018a; Tang et al., 2018).The reference period (1980?99) and projected period(2030?49) are analyzed to explore future climate change.

    2.2. Bias correction methods

    As mentioned above, due to the deficiencies in RCMs,large biases can be found in simulations when compared to the observations. Thus, bias correction methods based on the observations can be implemented to improve the performance of RCMs. Considering the relatively short time sequence of data, here three common stationary temperature bias correction methods are used from the R packages“hyfo” and “DownscaleR”; namely, variance scaling, additive scaling, and quantile mapping based on empirical distribution (Wilcke et al., 2013). To facilitate bias correction and ensemble calibration, temperatures from RCMs were interpolated to a common grid of 0.5° × 0.5° latitude/longitude, following CRU, using bilinear interpolation.

    Table 1. Configurations of the five RCMs in the CORDEX-EA region (after Gu et al., 2018).

    Fig. 1. Simulation domain and topology of CORDEX-EA and the 10 selected subregions:northwestern China (NW), Tibetan Plateau (TP), northeastern China (NE), northern China(NC), southern China (SC), Korean Peninsula and Japan (KJ), Mongolia (MG), India,Indochina (InC), and Southeast Asia (SEA).

    2.2.1. Additive scaling

    The additive scaling method performs bias correction based on the bias between the average simulation and average observation during the calibration period. It is expressed as

    whereTvalandTbcdenote the validation temperature sequence before and after bias correction, respectively, andandare the average of the simulation and observation temperatures in the calibration, respectively.

    2.2.2. Variance scaling

    Building on the additive scaling method, the variance scaling method further corrects the variance of the temperature (Terink et al., 2010) using the following procedure.First, we apply the additive scaling method, which corrects the average of the temperature:

    Then, we convert the average temperature to 0, using

    Next, we scale the variance of temperature according to the ratio of the temperature variance between the calibration and validation periods:

    And finally, we add the temperature in step 3 to the temperature in step 1:

    2.2.3. Quantile mapping

    The quantile mapping (QM) method derives from the empirical transformation developed by Panofsky and Brier(1968). It has been widely used in the bias correction of both GCMs and RCMs (Wilcke et al., 2013; Miao et al.,2016). By contrast with the methods mentioned above, the QM method focuses not only on the mean of the distribution but also on correcting the quantiles of the distribution.The QM method estimates the cumulative distribution function (CDF) from the simulation data in the calibration period and then finds the corresponding percentile values of the model projections. The corrected projections can be derived through inverse CDFs of the observations. The transfer function is shown as follows:

    where the subscripts O and mc denote the observation and model calibration periods, respectively.Fmcis the CDF from simulation data in the calibration period, andis the inverse CDF of observations.

    2.3. Simple multi-model averaging method

    Simple multiple-model averaging (SMA), which gives each member in the ensemble equal weight, is the most common ensemble post-processing method. For an ensemble of n members, the weight of each member is 1 /n(Arsenault et al., 2015).

    2.4. BMA method

    The BMA method was introduced by Raftery et al.(2005) to combine different model forecasts into an ensemble and calibrate the under-dispersion during the ensemble forecasts (Duan and Phillips, 2010; Miao et al.,2013). BMA can be viewed as a post-processing method for producing the forecast probability density function (PDF) of output variables, which is a weighted average of the bias-corrected PDF of each individual ensemble member. The weights reflect the relative performance of each member in the ensemble during the training period. Following the notation in Raftery et al. (2005), the BMA-weighted forecast PDF of variable y is wkσ2

    The BMA weights and variance are estimated through the maximum likelihood method (Raftery et al.,2005). The expectation maximization algorithm is adopted in this study to calculate the BMA weights and variance.For a detailed description of the BMA method, see Raftery et al. (2005).

    2.5. Model evaluation in spatial simulation

    In this study, we conducted the bias correction for each subregion and for each month of the year. In addition, considering the time dependence of the model biases and the relatively short time sequence of historical simulations (only 26 years, from 1980 to 2006), cross validation (Miao et al.,2016) was used to calibrate and validate the performance of the RCMs. For calibration, 20 years were randomly selected out of the 26 years and then three bias correction methods were applied to get the bias correction factors. The remaining 6 years were used for validation, and the bias correction factors were applied. The sampling method used in cross validation was simple random sampling without replacement, which guarantees that each member has an equal opportunity to be selected. Validation was conducted by comparing the corrected temperatures with CRU data. We repeated the whole cross-validation process 30 times to overcome the limitation of insufficient sample sizes and enhance the robustness of the validation results. The performance of different bias correction methods was also evaluated based on the cross-validation results.

    While the most effective bias correction method was being obtained, the historical temperature simulations(1980-2005) were corrected for all subregions. Then, the ensemble averaging methods were applied to the bias-corrected temperatures. To validate the performance of ensemble averaging methods, monthly temperature distributions, interannual variability and Taylor diagrams were used. The Taylor diagram is especially useful in evaluating multiple aspects of complex models (IPCC, 2001). It incorporates three evaluation terms?spatial correlation, centered RMSE and standard deviations?and graphically measures how closely the simulation patterns match the observations(Taylor, 2001).

    3. Results

    3.1. Bias correction evaluation

    Fig. 2. Spatial distribution of annual average temperatures according to (a) CRU and (b) the HadGEM2-AO GCM. (c?h) Annual temperature biases of (c) the HadGEM2-AO GCM and (d?h) five RCMs during the years 1980?2005.

    Figure 2 illustrates the spatial pattern of annual average (1980-2005) temperatures of CRU and the driving GCM (HadGEM2-AO) as well as the annual average temperature biases of the driving GCM and the five RCMs. In addition, the annual average temperatures of the five RCMs are shown in Fig. S1 in the electronic supplementary material(ESM). We can see that both the GCM and the RCMs captured the spatial pattern of the observations with a decreasing south to north temperature gradient. However, both the driving GCM and the RCMs generally underestimated the annual average temperature in most regions, especially in the TP region, where the greatest bias exceeded ?8°C. The cold biases in the simulation came mainly from December?January?February (DJF), and the biases in June-July-August (JJA) were small (Fig. S2). Several studies have also reported similar annual average and seasonal temperature bias patterns (Ham et al., 2016; Guo et al., 2018; Li et al., 2018a; Hui et al., 2019). The cold biases in DJF were remarkable in most RCMs, with the largest bias exceeding?8°C in most RCMs. The only exception was RegCM,which had warm biases exceeding 4°C at high latitudes in CORDEX-EA, which is consistent with previous studies(Gao and Giorgi, 2017). The RCMs’ performance also varied in JJA temperature simulations. For example, the YSURCM model presented a large cold bias in NE and the TP(exceeding ?8°C), while the bias was small in other RCMs.

    In addition, compared to the driving GCM, some RCMs had improved their skill in simulating the temperature in several regions. For example, WRF and MM5 reduced the biases over the NC, SC and KJ regions, and RegCM and YSU-RCM improved their performance in the India region. However, for the other regions (e.g., TP and NW), the improvement of the RCMs was less and the cold biases were even larger than those of the driving GCM.Increased resolution does not always lead to improvement in simulations (Pr?mmel et al., 2010; Gu et al., 2018), and there are several reasons for this phenomenon. For instance,due to their simplified physical parameterization schemes,longwave radiation is underestimated in RCMs, which limits the heating of the low-level atmosphere (Hui et al.,2019). Also, the overestimation of albedo in RCMs (in the lower boundary conditions) may also lead to cold biases(Meng et al., 2018).

    Fig. 3. Comparisons of different bias correction methods at the 25th, 50th and 75th quantiles.

    Given the large biases in simulating annual and seasonal average temperatures, bias correction is necessary to improve the RCMs’ performance. We calculated the percentage of effective pixels for each method over 30 rounds of validation (Fig. 3). We found that all the methods can effectively improve RCM performance. The percentage of effective pixels at the 25% quantile was higher than at the 50%and 75% quantiles. The QM method outperformed the other two methods for the 25%, 50% and 75% quantiles. This may be due to the fact that, although all bias correction methods corrected the biases in the RCMs, their foci are different. For example, the additive scaling method focuses on the mean difference between the calibration model and observation data (Teutschbein and Seibert, 2012; Ezéchiel et al.,2016), while the variance scaling method focuses on adjustment of variance (Luo et al., 2018). The two methods can adjust the monthly mean values, but they neglect the cumulative distribution of the temperature, and thus they cannot adjust the quantiles of the simulation. However, in contrast with these two methods, the QM method constructs CDFs of RCMs and adjusts the distributions according to the corresponding distribution of observation data (Bennett et al.,2014; Singh et al., 2017; Ayugi et al., 2020). Therefore, the QM method can adjust not only the mean value of temperature but also the quantile values.

    We calculated the RMSE for both raw and bias-corrected temperature data, then the relative decrease in the RMSE for the 10 subregions and five RCMs was calculated based on the 30 rounds of cross validation (Fig. 4). Results showed that all the bias correction methods effectively reduced biases for all RCMs. The QM method was most effective among the methods, with the maximum relative decrease in the RMSE reaching 59.8% (HadGEM3-RA),63.2% (MM5), 51.3% (RegCM), 80.7% (YSU-RCM) and 62.0% (WRF). For subregions, although all the bias correction methods significantly reduced the biases over most subregions, results varied. For example, in the SEA region,MM5, bias-corrected by the additive scaling method, was worse than the raw simulation. We analyzed the results(Fig. S3) and found that this may be due to the fact that the additive scaling method can only adjust the mean difference between the model values and observations, while for extreme values, the additive scaling method failed (Fang et al., 2015). In the SEA region, for MM5, the additive scaling method narrowed the biases for high temperature but amplified the bias for low temperature when compared to the raw simulations. However, the QM-corrected results almost perfectly fitted the CRU distribution. Several previous studies have also shown that bias correction methods are regiondependent, and the scaling method may have adverse effects in some regions (Berg et al., 2012; Ayugi et al., 2020).Moreover, for almost all subregions (except the MG region as simulated by WRF and MM5), the QM method outperformed the other two methods. We also found that, for the TP and NW regions, where the bias was largest among the subregions, the bias reduction was significant, with the maximum relative decrease in the RMSE reaching 61.5% and 80.7%, respectively (both for YSU-RCM).

    Fig. 4. Relative decrease in RMSE (%) of three bias correction methods for the 10 subregions for the five RCMs: (a)HadGEM3-RA, (b) RegCM, (c) MM5, (d) WRF, and (e) YSU-RCM.

    Figure 5 gives the spatial distribution of the relative decrease in the MAE for the 10 subregions (the QM method performed best, so only the results of QM are shown here).The spatial distribution results are consistent with results in Fig. 4. The bias correction was effective in most regions, especially in the TP, NW, SC and NC subregions, where the MAE was reduced more than 50% compared to the MAE of the raw model output. For the YSU-RCM model, which had the largest biases in the simulation, the bias reductions were remarkable: the relative decrease in the MAE was more than 60% for most regions and more than 70% in most parts of China. The spatial distribution results for winter (DJF)and summer (JJA) (Fig. S4) also showed that the bias correction method effectively reduced the bias in seasonal simulations, especially in winter. The correction for summer focused mainly on the southern part of the CORDEX region, including India, InC and SC. In summary, all bias correction methods effectively reduced the biases in the simulations, and the QM method was most effective for almost all subregions and all RCMs. Therefore, the QM method was chosen as the most suitable method to correct the temperature simulations in the following text. Figure S5 provides the annual average temperature error distribution of the bias-corrected data using the QM method. The results show that, compared with the raw RCM simulations, the errors were removed remarkably well by the bias correction process, especially for the NE, NC, KJ, SC, INC, SEA and India regions,where the errors were small. However, for the NW, TP and MG regions, although bias correction did decrease the cold biases, some cold biases remained. Furthermore, due to the remarkable cold biases in the NW, TP and MG regions, the bias correction tended to overcorrect cold biases in these regions, especially for the TP and NW regions in the YSURCM model.

    3.2. Multi-model averaging based on bias correction

    Fig. 5. Spatial distribution of relative decrease in MAE (%) for annual temperatures using the QM method. Panels (a?e)show results from each of the five RCMs. See Fig. S4 in the ESM for spatial distribution of relative decrease in MAE (%) for seasons.

    In addition to the bias correction method, BMA and SMA were also used to further narrow the uncertainty in the corrected RCMs. Figure 6 gives the seasonal distribution of CRU, raw simulated, bias-corrected, SMA-weighted [bias-corrected (BC)] and BMA-weighted (BC) temperatures. Results show that, compared to the CRU data, the raw simulation had significant cold biases in winter (DJF) and biases were small in summer (JJA) over the 10 subregions. This is consistent with the results in Fig. S2 and several previous studies, where the cold biases in winter were large and in summer were small (Ham et al., 2016; Hui et al., 2019).Moreover, RCM performance was region-dependent; more specifically, the RCMs performed well for the regions at high latitudes (e.g., the NE, KJ and MG regions) and badly for regions at low latitudes (e.g., the SEA, TP, NW and InC regions). Several previous studies have analyzed the possible causes of the biases in these regions and suggested that careful configuration of the RCM parameterization schemes could help. For example, according to Hui et al. (2019), the radiation parameterization schemes in RCMs underestimate the monthly longwave upward and downward fluxes throughout the year, especially in cold months over subtropical regions, which leads to significant cold biases in subtropical regions. For the cold biases in the TP region, they may be attributable to the overestimation of upward shortwave radiation and the corresponding overestimation in albedo(Tangang et al., 2015; Chen et al., 2017; Hui et al., 2019;Yin et al., 2020). When the bias corrections of SMA and BMA were applied, the seasonal cycles of the raw models were adjusted and fitted the CRU data well. The bias corrections of SMA and BMA significantly improved the performance of the RCMs.

    The model’s ability to capture the real interannual variability is another important performance measure. Here, we used the variance among 26 historical years as the indicator of interannual variability. Figure 7 gives the results of interannual variability of the SMA (BC), BMA (BC), SMA (raw)and CRU temperature values. The results show that the RCMs’ ability to capture the real interannual variability varied among subregions. For the NW, TP, NE, MG and SEA regions, the interannual variability of bias-corrected data was closer to the real interannual variability. But for other regions, the bias correction narrowed the variability. The results may be due to the fact that the bias correction methods mainly focus on the mean and trend of the data, with less focus on the variance (Ayugi et al., 2020). Although, for most subregions, the performances of SMA and SMA (BC)were similar, but for the MG region, where the interannual variability was greater, the variabilities of SMA (BC) were closer to the variability of CRU. Thus, SMA (BC) was considered as a better method. In addition, the results based on SMA (BC) were better than those based on BMA (BC)when compared to the variability of CRU. This is due to the fact that the objective function of BMA only considers the minimum bias without adjusting for variance (Raftery et al.,2005; Fragoso et al., 2018). In future studies, we will pay more attention to the variance and consider multi-objective optimization.

    The spatial variability statistics of the models in reproducing annual average temperature are shown using Taylor diagrams in Fig. 8. The Taylor diagrams show that the bias correction improved the performance of the RCMs, contributing to a higher spatial correlation and lower normalized standard deviation. Furthermore, the BMA and SMA ensemble results both reduced the uncertainties in simulation with a closer distance to the observation. For some subregions, the performances of BMA and SMA were similar (e.g., in NE,NC, TP and SC). However, for other regions, such as NW,MG and SEA, the BMA method performed better. This is reasonable because the BMA weights are estimated according to the RCMs’ performance in the training period (Duan et al., 2007). A previous study also showed that the BMA method outperformed the SMA method when applied to the CORDEX-EA data (Kim and Suh, 2013). For the correlation coefficient, the improvement rate of the BMA method was between 2% and 31% when compared with individual RCMs. Although the SMA performed better with respect to interannual variability, we focus mainly on the mean and trend in the projection. Thus, we chose the BMA method for the projection to narrow the uncertainties. Although the BMA-related improvement was not substantial, the amount of improvement was reasonable considering the temperature had already been corrected by bias correction methods.

    Fig. 6. Observed (CRU), raw simulated, bias-corrected, SMA [bias-corrected (BC)] and BMA [bias-corrected (BC)] monthly temperatures of 10 subregions in the validation period.

    Fig. 7. Interannual variability for SMA [bias-corrected (BC)], BMA [biascorrected (BC)], SMA (raw model), and CRU temperature values among the 26 years of the historical period.

    3.3. Temperature projections

    Based on the most effective bias correction method(QM) and BMA weights derived from the training period(1980?2005), the bias-corrected and BMA-weighted temperature projections for the 10 subregions under the two scenarios (RCP4.5 and RCP8.5) were generated. Figure 9 shows the average temperature projections for the BMA ensemble under the RCP4.5 and RCP8.5 scenarios (results of five RCMs and the driving GCM are shown in Figs. S6 and S7 under the RCP4.5 and RCP8.5 scenarios, respectively). Similar warming trends were detected over the 10 subregions for the 2030?2049 period under both scenarios but with a more obvious warming trend under the RCP8.5 scenario.The warming trend was more remarkable in the northern part of CORDEX-EA than in the southern part, especially for the TP, NW, MG and NE regions, where the warming was over 3°C in the 2030?49 period. Moreover, analysis of seasonal warming results indicated that the warming was more remarkable in winter than in summer (not shown). Similar warming patterns have also been detected in previous studies (Ham et al., 2016; Gu et al., 2018), although these studies focused mainly on the China region. The BMA results indicate a clear increase in average temperature under the RCP4.5 and RCP8.5 scenarios for all subregions (Table 2).However, for a given subregion, the warming varied among RCMs. For example, the annual temperature increase over NE ranged from 0.5°C to 3.5°C under the RCP4.5 scenario.Figure 10 also illustrates the obvious warming trend over the 10 subregions from 2006 to 2049 under both scenarios.Note that the warming trends in the TP, NW, NE and MG regions [reaching 0.6°C (10 yr)?1?0.7°C (10 yr)?1] were more remarkable than in the other regions [0.3°C (10 yr)?1?0.5°C (10 yr)?1] under the RCP8.5 scenario, which is consistent with results in Fig. 9.

    The future temperature changes varied among subregions and months. Figure 11 illustrates the changes of monthly temperature in the 10 subregions under the RCP8.5 scenario (changes under RCP4.5 were similar but with smaller amplitude; not shown). The BMA projection also indicates that monthly temperature increased more notably in the northerly regions of CORDEX-EA, especially the NE and MG regions, where the most rapid increases were in November (more than 4.5°C). Furthermore, we also found that the increases in monthly temperature varied by latitude. For example, the MG and NE regions exhibited similar increasing patterns; the same was true for the KJ, NC, NW and TP regions, as well as for the SEA, India and InC regions. The SC region was dissimilar to the other regions. According to Hui et al. (2019), this may be due to the cloud cover, but determining detailed reasons for how the pattern of increases in temperature varies with latitude needs further study. Finally, the monthly warming pattern in the subregions of China under the RCP4.5 scenario ranged from 0.8°C to 4.2°C. These values are similar to but larger than findings in a previous study (0.3°C?2.2°C) (Gu et al.,2018), which was based on the raw temperature and simple multi-model ensemble averaging.

    Fig. 8. Taylor diagrams evaluating the model skill in simulating the annual temperature and bias correction effects over 10 subregions. CRU observation data used as reference. The x- and y-axes refer to the standard deviations (normalized) and the azimuthal axis refers to the spatial pattern correlation between two fields.

    The uncertainty in projections needs to be considered when using them in applications (such as driving hydrological models). Here, we took the standard deviation among five RCMs as an uncertainty indicator (Nordhaus, 2018).Figure 12 gives the uncertainty of the raw model outputs and the bias-corrected results for both the RCP4.5 and RCP8.5 scenarios. The results show that the uncertainty had no apparent relationship with time or solar radiation forcing(RCP4.5 and RCP8.5). The results were different from uncertainties projected by GCMs in Asia, where the uncertainty increases with time (Miao et al., 2016). This may be due to the fact that the RCMs in CORDEX-EA use the same driving GCM, so there is only internal variability (Chen et al.,2019). The uncertainty was greatest for the TP and NW regions, exceeding 2.5°C both for the RCP4.5 and RCP8.5 scenarios throughout the projection period. For the SEA, KJ and India subregions, the uncertainty was lowest?less than 0.7°C for both scenarios. This indicates that there was more uncertainty in the high-latitude subregions. Several previous studies have also shown that the uncertainty was greater at high latitude (Deser et al., 2012; Miao et al., 2016;Woldemeskel et al., 2016). Because the QM method narrowed the differences among RCMs, the uncertainty was reduced for most subregions and both scenarios (except for the NE region under the RCP4.5 scenario). The reductions were more remarkable for the RCP8.5 scenario, ranging from 66% to 94% across all subregions. The lower uncertainty in the RCP8.5 scenario indicates a consistent warming trend under the scenario for all subregions.

    Fig. 9. Spatial distribution of annual temperature changes [RCP 8.5 (4.5) minus baseline] projected by five RCMs and the BMA method. See supplementary material for results of five RCMs and the driving GCM under RCP4.5 (Fig. S6) and RCP8.5 (Fig. S7).

    Table 2. Annual temperature changes (future minus baseline) in °C projected by five RCMs and the BMA method for 10 subregions under the RCP4.5 and RCP8.5 scenarios.

    4. Summary and conclusions

    In this study, five RCMs that participated in CORDEXEA, with lateral and boundary forcing from the HadGEM2-AO model, were used to derive temperature projections over 10 subregions in CORDEX-EA under the RCP4.5 and RCP8.5 scenarios. To remove the biases in RCMs and narrow uncertainty in the projections, bias correction methods and ensemble calibration methods were used. At the same time, we analyzed the performance of bias correction and ensemble calibration methods over different subregions and RCMs. The major findings of the study can be summarized as follows:

    Fig. 10. Results from five RCMs and BMA for annual temperature change under two scenarios, RCP4.5 and 8.5. The straight lines, marked with uneven dashes, represent the trends of the temperature changes.

    (1) The RCMs all revealed reasonable representation of annual average temperature when compared with the reference data (CRU TS 4.03). All RCMs presented cold biases in most subregions, especially in regions with complex topography (e.g., the TP and NW regions). Moreover, the seasonal analysis shows that the cold bias mainly came from DJF, while the biases were relatively small in JJA. Therefore, there are obvious biases in the RCM dynamic downscaling results, and it is necessary to conduct bias correction.

    (2) The validation results of the bias correction show that all bias correction methods could reduce the simulation biases, and the QM method was most effective in all subregions and RCMs. Compared to the additive scaling and vari-ance scaling methods, the QM method corrected not only the mean but also the quantiles of the distribution. We found that the QM method could reduce the simulation bias by more than 50% when compared to the raw model output.For the TP and NW regions, the relative decreases in RMSE reached 61.5% and 80.7%, respectively. The spatial distribution of relative decreases in MAE further confirmed the results for RMSE. Finally, the QM method was used to correct the temperatures in projections of the future.

    Fig. 11. Projected monthly temperature changes (RCP 8.5 minus baseline) for 10 subregions. The monthly temperature changes for the RCP4.5 scenario are similar to RCP8.5 and not shown.

    Fig. 12. Uncertainty of raw model data and bias-corrected data under two scenarios. We used the standard deviation as a measure of uncertainty.

    (3) The SMA and BMA methods could further narrow the uncertainty in the bias-corrected ensemble. Although the performances of the SMA and BMA methods were similar in some subregions, the BMA method performed better overall. Compared to the individual RCMs, the BMA method improved the correlation coefficient from 2% to 32% over the 10 subregions. Note that the BMA method may produce poorer results for interannual variability of the temperature,since it focuses only on the bias and trend of the data.

    (4) For temperature projections (2030?49), both single RCMs and the BMA results showed consistent warming over all subregions under both scenarios (RCP4.5 and 8.5).BMA results indicated that the warming ranged from 1.2°C to 3.5°C over the 10 subregions under the RCP8.5 scenario(and from 1.0°C to 2.6°C for RCP4.5). The warming was more pronounced in the northern part of the CORDEX-EA domain. The monthly temperature changes seem to vary by latitude, and the detailed reasons need further study.

    (5) Furthermore, the QM method reduced uncertainty in temperature projections for all subregions. The reduction was more notable for the RCP8.5 scenario for all subregions. The QM method decreased the uncertainty among the 10 subregions by 66%?94% for the RCP8.5 scenario.

    Acknowledgements. This work was supported by the Stra-tegic Priority Research Program of the Chinese Academy of Sciences (Grant No. XDA20060401) and the National Natural Science Foundation of China (Grant Nos. 41622101 and 41877155).We acknowledge each of the CORDEX-EA modeling groups for making their simulations available for analysis, including the National Institute of Meteorological Research, three universities in the Republic of Korea (Seoul National University, Yonsei University and Kongju National University), and the World Climate Research Programme’s Working Group on the Coordinated Regional Climate Downscaling Experiment, for making the CORDEX data set available (http://cordex-ea.climate.go.kr/cordex/treeP-age.do).

    Electronic supplementary material:Supplementary material is available in the online version of this article at https://doi.org/10.1007/s00376-020-0026-6.

    成人18禁高潮啪啪吃奶动态图| 狠狠狠狠99中文字幕| 一边摸一边做爽爽视频免费| 天堂av国产一区二区熟女人妻 | 欧美午夜高清在线| 一卡2卡三卡四卡精品乱码亚洲| 午夜免费观看网址| av天堂在线播放| 亚洲最大成人中文| 啦啦啦观看免费观看视频高清| xxxwww97欧美| 亚洲九九香蕉| 久久 成人 亚洲| 日本撒尿小便嘘嘘汇集6| 我要搜黄色片| 久久久久久久久久黄片| 精品人妻1区二区| 曰老女人黄片| 级片在线观看| 久久午夜亚洲精品久久| 舔av片在线| 久久亚洲精品不卡| 悠悠久久av| 无限看片的www在线观看| 久久伊人香网站| 制服诱惑二区| 欧美日韩乱码在线| 国产91精品成人一区二区三区| www.自偷自拍.com| 欧美 亚洲 国产 日韩一| 99热这里只有精品一区 | 少妇的丰满在线观看| 欧美日韩乱码在线| 亚洲精品美女久久av网站| 亚洲自拍偷在线| 精品少妇一区二区三区视频日本电影| 久久久久久久久中文| 欧美日韩亚洲国产一区二区在线观看| 国模一区二区三区四区视频 | 午夜日韩欧美国产| 午夜福利成人在线免费观看| 国产av又大| 久久久久久久精品吃奶| 香蕉av资源在线| 亚洲全国av大片| 亚洲av五月六月丁香网| 亚洲国产欧美网| 久久国产精品人妻蜜桃| 欧美中文日本在线观看视频| 午夜福利在线观看吧| 久久久久久久午夜电影| 婷婷六月久久综合丁香| 午夜免费成人在线视频| 人人妻人人看人人澡| 人妻久久中文字幕网| a在线观看视频网站| 俄罗斯特黄特色一大片| 美女免费视频网站| 亚洲国产高清在线一区二区三| 国产在线精品亚洲第一网站| 欧美成人免费av一区二区三区| 国产精华一区二区三区| 免费看美女性在线毛片视频| 国产91精品成人一区二区三区| 男人舔奶头视频| 欧美性长视频在线观看| 国产探花在线观看一区二区| 欧美丝袜亚洲另类 | 一二三四社区在线视频社区8| 久久精品人妻少妇| 国产精品99久久99久久久不卡| 看片在线看免费视频| 色av中文字幕| 操出白浆在线播放| 国产精品久久视频播放| 国产精品自产拍在线观看55亚洲| 99re在线观看精品视频| 色哟哟哟哟哟哟| 白带黄色成豆腐渣| 真人做人爱边吃奶动态| 最新美女视频免费是黄的| 香蕉国产在线看| 香蕉丝袜av| 19禁男女啪啪无遮挡网站| 亚洲成人国产一区在线观看| 老司机在亚洲福利影院| 一级毛片女人18水好多| 一级片免费观看大全| 舔av片在线| 亚洲成人国产一区在线观看| 欧美av亚洲av综合av国产av| 色哟哟哟哟哟哟| 日本成人三级电影网站| 国产伦在线观看视频一区| 国产又黄又爽又无遮挡在线| 久久草成人影院| 国产黄色小视频在线观看| 男人舔奶头视频| 老熟妇乱子伦视频在线观看| 啦啦啦观看免费观看视频高清| 老司机午夜福利在线观看视频| 午夜免费激情av| 黄频高清免费视频| АⅤ资源中文在线天堂| 99re在线观看精品视频| 人人妻人人看人人澡| 午夜精品一区二区三区免费看| 欧美一级毛片孕妇| 最新在线观看一区二区三区| 久9热在线精品视频| 国产成人精品无人区| 最近视频中文字幕2019在线8| 欧美午夜高清在线| 最好的美女福利视频网| 成人国产综合亚洲| 久久人妻av系列| av片东京热男人的天堂| 好男人在线观看高清免费视频| a级毛片在线看网站| 精品免费久久久久久久清纯| 久久久水蜜桃国产精品网| 亚洲片人在线观看| 大型黄色视频在线免费观看| 亚洲最大成人中文| 麻豆一二三区av精品| 精品国内亚洲2022精品成人| 伦理电影免费视频| 成人三级做爰电影| 99久久精品国产亚洲精品| 两性夫妻黄色片| 青草久久国产| 夜夜爽天天搞| 亚洲一码二码三码区别大吗| 国产精品精品国产色婷婷| 深夜精品福利| 免费av毛片视频| 国产蜜桃级精品一区二区三区| 国产精品久久久久久亚洲av鲁大| 久久久久久九九精品二区国产 | 巨乳人妻的诱惑在线观看| 我的老师免费观看完整版| 色哟哟哟哟哟哟| 欧美不卡视频在线免费观看 | 国产精品影院久久| 白带黄色成豆腐渣| 成人国产综合亚洲| 国产精品一区二区三区四区免费观看 | 免费人成视频x8x8入口观看| 国产伦人伦偷精品视频| 19禁男女啪啪无遮挡网站| 久久久水蜜桃国产精品网| 黄色视频不卡| 久久久精品欧美日韩精品| 99精品欧美一区二区三区四区| 国产不卡一卡二| 久久精品人妻少妇| 麻豆久久精品国产亚洲av| 亚洲一区二区三区色噜噜| 国产精品久久久久久久电影 | 最好的美女福利视频网| 亚洲 欧美 日韩 在线 免费| 成年版毛片免费区| 国产免费男女视频| 亚洲av成人一区二区三| 最近最新中文字幕大全电影3| 黑人欧美特级aaaaaa片| 国产人伦9x9x在线观看| 中出人妻视频一区二区| 在线看三级毛片| 69av精品久久久久久| 91在线观看av| 香蕉久久夜色| 最近最新中文字幕大全电影3| 99热这里只有精品一区 | 一区二区三区激情视频| 激情在线观看视频在线高清| 国产伦在线观看视频一区| 不卡av一区二区三区| 免费在线观看视频国产中文字幕亚洲| 天天一区二区日本电影三级| 岛国在线观看网站| 国产精品1区2区在线观看.| 久久国产精品影院| 亚洲专区国产一区二区| 国产精品一区二区精品视频观看| 亚洲欧美日韩东京热| 国产激情偷乱视频一区二区| 窝窝影院91人妻| 久久热在线av| 久久午夜综合久久蜜桃| 黑人欧美特级aaaaaa片| 黑人欧美特级aaaaaa片| 国产99久久九九免费精品| 中文字幕久久专区| 成人一区二区视频在线观看| 黄色丝袜av网址大全| 免费人成视频x8x8入口观看| 国产精品亚洲美女久久久| 窝窝影院91人妻| 精品高清国产在线一区| 国产97色在线日韩免费| 国模一区二区三区四区视频 | 亚洲专区国产一区二区| 欧美日韩精品网址| 成人国产综合亚洲| 欧美日韩精品网址| 三级男女做爰猛烈吃奶摸视频| 脱女人内裤的视频| 精品高清国产在线一区| 男男h啪啪无遮挡| 在线播放国产精品三级| 欧洲精品卡2卡3卡4卡5卡区| 精品午夜福利视频在线观看一区| 国产精品亚洲av一区麻豆| 最新在线观看一区二区三区| 桃色一区二区三区在线观看| bbb黄色大片| 亚洲成a人片在线一区二区| 欧美黄色淫秽网站| 欧美成人免费av一区二区三区| 精品久久久久久成人av| 欧美另类亚洲清纯唯美| 欧美另类亚洲清纯唯美| 制服人妻中文乱码| 久久久久久久久免费视频了| 亚洲精品av麻豆狂野| 免费在线观看日本一区| 欧美久久黑人一区二区| 亚洲 欧美一区二区三区| 国产伦在线观看视频一区| 一本大道久久a久久精品| 老司机在亚洲福利影院| 欧美不卡视频在线免费观看 | 伦理电影免费视频| 日本在线视频免费播放| www.熟女人妻精品国产| 欧美久久黑人一区二区| 日日摸夜夜添夜夜添小说| 日日摸夜夜添夜夜添小说| 一区二区三区国产精品乱码| 欧美性长视频在线观看| 国产乱人伦免费视频| 精品一区二区三区视频在线观看免费| 男人舔女人的私密视频| 人人妻,人人澡人人爽秒播| 97人妻精品一区二区三区麻豆| 天堂动漫精品| 国产成人精品久久二区二区91| 亚洲一卡2卡3卡4卡5卡精品中文| 精品日产1卡2卡| 人人妻人人看人人澡| 免费无遮挡裸体视频| 国内精品一区二区在线观看| 国产精品1区2区在线观看.| 少妇粗大呻吟视频| av天堂在线播放| 成人三级黄色视频| 一进一出好大好爽视频| 老汉色av国产亚洲站长工具| 俄罗斯特黄特色一大片| 日韩高清综合在线| 免费在线观看黄色视频的| 亚洲欧美日韩高清在线视频| 久久香蕉精品热| 国产蜜桃级精品一区二区三区| 亚洲av成人一区二区三| 国产亚洲av高清不卡| av欧美777| 大型av网站在线播放| 国产单亲对白刺激| 九色国产91popny在线| 中文在线观看免费www的网站 | 精品久久久久久久久久免费视频| 欧美中文综合在线视频| 免费在线观看日本一区| 色综合婷婷激情| 婷婷亚洲欧美| 精品一区二区三区视频在线观看免费| 18禁观看日本| 成人国产一区最新在线观看| 久久久国产欧美日韩av| 久久亚洲精品不卡| 亚洲精品美女久久av网站| 亚洲精品色激情综合| av天堂在线播放| 日本精品一区二区三区蜜桃| 国产精品爽爽va在线观看网站| 亚洲aⅴ乱码一区二区在线播放 | 欧美大码av| 久久人妻av系列| 欧美高清成人免费视频www| 香蕉久久夜色| 女同久久另类99精品国产91| 男人舔女人的私密视频| 午夜福利18| 国产精品一区二区三区四区久久| 天天添夜夜摸| 亚洲七黄色美女视频| 好男人在线观看高清免费视频| 国产精品亚洲av一区麻豆| 天堂影院成人在线观看| 亚洲中文日韩欧美视频| 欧美精品亚洲一区二区| 日韩欧美 国产精品| 在线观看一区二区三区| 欧美黑人欧美精品刺激| 黄频高清免费视频| 亚洲成av人片免费观看| bbb黄色大片| 久久午夜亚洲精品久久| 国产日本99.免费观看| 人妻久久中文字幕网| 正在播放国产对白刺激| 国内精品久久久久精免费| 岛国在线观看网站| 午夜福利高清视频| 国产蜜桃级精品一区二区三区| 亚洲无线在线观看| 午夜福利在线观看吧| 久久久精品国产亚洲av高清涩受| 久久久久国产一级毛片高清牌| 午夜久久久久精精品| 亚洲中文日韩欧美视频| 麻豆成人av在线观看| 人人妻人人澡欧美一区二区| 久久精品影院6| 日本免费a在线| 麻豆国产av国片精品| 99热这里只有精品一区 | 国产精品一区二区精品视频观看| 亚洲一区二区三区色噜噜| 亚洲自拍偷在线| 国产av不卡久久| 中国美女看黄片| 久久午夜亚洲精品久久| 久久精品国产99精品国产亚洲性色| 亚洲在线自拍视频| 久久久久国产精品人妻aⅴ院| 亚洲欧美精品综合一区二区三区| 成人国产一区最新在线观看| 久久亚洲精品不卡| a级毛片a级免费在线| 午夜福利视频1000在线观看| 特大巨黑吊av在线直播| 看黄色毛片网站| 精品不卡国产一区二区三区| 日本免费一区二区三区高清不卡| av在线播放免费不卡| a级毛片a级免费在线| 欧美成狂野欧美在线观看| 国产区一区二久久| 久久久久久久久免费视频了| 亚洲一区中文字幕在线| 97碰自拍视频| 此物有八面人人有两片| 国产精品免费一区二区三区在线| 久久中文字幕一级| 欧美日韩一级在线毛片| 最新在线观看一区二区三区| 看免费av毛片| 一二三四在线观看免费中文在| 亚洲国产精品999在线| 巨乳人妻的诱惑在线观看| 黄频高清免费视频| 欧美日韩福利视频一区二区| 99久久综合精品五月天人人| 亚洲人成伊人成综合网2020| 久久热在线av| 桃色一区二区三区在线观看| 国产片内射在线| 又黄又粗又硬又大视频| 日韩欧美 国产精品| 少妇粗大呻吟视频| 国产精品久久久人人做人人爽| 99久久无色码亚洲精品果冻| 成人国语在线视频| 久久久精品大字幕| 国产亚洲av高清不卡| 美女免费视频网站| 午夜福利免费观看在线| 久久久国产成人免费| 九色国产91popny在线| 婷婷精品国产亚洲av| 五月伊人婷婷丁香| 国产三级黄色录像| 1024手机看黄色片| 日本撒尿小便嘘嘘汇集6| 久久人妻av系列| 又黄又爽又免费观看的视频| 少妇人妻一区二区三区视频| 视频区欧美日本亚洲| 18禁裸乳无遮挡免费网站照片| 亚洲av电影不卡..在线观看| 国产亚洲精品综合一区在线观看 | 亚洲国产高清在线一区二区三| 国产精品久久久人人做人人爽| 久久久久久久精品吃奶| x7x7x7水蜜桃| 精品国内亚洲2022精品成人| av福利片在线观看| 亚洲精品在线观看二区| 国产一区二区激情短视频| 久9热在线精品视频| 亚洲天堂国产精品一区在线| 777久久人妻少妇嫩草av网站| 在线观看一区二区三区| 我要搜黄色片| 久久久久久久精品吃奶| 狂野欧美激情性xxxx| 搡老熟女国产l中国老女人| 日韩欧美国产一区二区入口| 国语自产精品视频在线第100页| 日韩成人在线观看一区二区三区| 亚洲无线在线观看| 精品熟女少妇八av免费久了| 色av中文字幕| 成人高潮视频无遮挡免费网站| 欧美av亚洲av综合av国产av| 久久久水蜜桃国产精品网| 1024手机看黄色片| 日本一本二区三区精品| 国产精品一区二区三区四区免费观看 | 在线观看免费视频日本深夜| 99热6这里只有精品| 亚洲国产看品久久| 国产亚洲精品久久久久5区| 国产av又大| 少妇的丰满在线观看| 黄色a级毛片大全视频| 国产亚洲av高清不卡| 97人妻精品一区二区三区麻豆| 国产伦人伦偷精品视频| 日本在线视频免费播放| 在线观看免费视频日本深夜| 国产高清视频在线观看网站| www.www免费av| 欧美久久黑人一区二区| 亚洲 国产 在线| 国产高清有码在线观看视频 | 超碰成人久久| av欧美777| a级毛片在线看网站| 99精品在免费线老司机午夜| 成人三级黄色视频| 久久久久久久精品吃奶| 国产精品一区二区免费欧美| 女人高潮潮喷娇喘18禁视频| 国产成人精品无人区| 欧美久久黑人一区二区| 日本五十路高清| 母亲3免费完整高清在线观看| 日韩大尺度精品在线看网址| 草草在线视频免费看| 国产av又大| 欧美 亚洲 国产 日韩一| 久久久水蜜桃国产精品网| 国产激情久久老熟女| 欧美+亚洲+日韩+国产| 男女之事视频高清在线观看| 国产三级黄色录像| 人妻丰满熟妇av一区二区三区| 啦啦啦观看免费观看视频高清| 日韩欧美国产一区二区入口| 大型av网站在线播放| 男女那种视频在线观看| 欧美一区二区精品小视频在线| 久久中文字幕一级| 制服人妻中文乱码| 人人妻,人人澡人人爽秒播| 日韩av在线大香蕉| 精品国内亚洲2022精品成人| 国产精品久久久av美女十八| 给我免费播放毛片高清在线观看| 国产亚洲精品综合一区在线观看 | 99在线视频只有这里精品首页| 一本综合久久免费| 宅男免费午夜| 国产日本99.免费观看| 亚洲欧美精品综合一区二区三区| 亚洲精品国产精品久久久不卡| 成人亚洲精品av一区二区| 国产午夜福利久久久久久| 成在线人永久免费视频| 日本在线视频免费播放| 欧美又色又爽又黄视频| 亚洲精品久久国产高清桃花| 在线国产一区二区在线| 久久精品国产清高在天天线| 精品无人区乱码1区二区| av福利片在线观看| 视频区欧美日本亚洲| 国产片内射在线| 男女午夜视频在线观看| 国产av又大| av免费在线观看网站| 一个人免费在线观看的高清视频| 亚洲男人的天堂狠狠| 两个人的视频大全免费| 亚洲aⅴ乱码一区二区在线播放 | 91麻豆精品激情在线观看国产| 久久久久性生活片| 夜夜看夜夜爽夜夜摸| 18禁美女被吸乳视频| 啦啦啦韩国在线观看视频| 久久这里只有精品中国| 一进一出抽搐动态| 亚洲熟女毛片儿| 在线观看美女被高潮喷水网站 | 香蕉国产在线看| 色老头精品视频在线观看| av中文乱码字幕在线| 婷婷六月久久综合丁香| 欧美黑人精品巨大| 国产av不卡久久| а√天堂www在线а√下载| 一夜夜www| 丝袜美腿诱惑在线| 久久精品国产亚洲av香蕉五月| 好男人电影高清在线观看| 小说图片视频综合网站| 曰老女人黄片| 又黄又粗又硬又大视频| 免费在线观看完整版高清| 亚洲第一欧美日韩一区二区三区| 麻豆久久精品国产亚洲av| 亚洲欧美一区二区三区黑人| 波多野结衣巨乳人妻| 亚洲成人国产一区在线观看| 国产单亲对白刺激| 国产三级中文精品| 黄片大片在线免费观看| www日本在线高清视频| 琪琪午夜伦伦电影理论片6080| 中文字幕最新亚洲高清| 欧美在线黄色| 很黄的视频免费| 18禁观看日本| 桃红色精品国产亚洲av| 亚洲欧美日韩东京热| 麻豆国产97在线/欧美 | 亚洲专区字幕在线| 99在线人妻在线中文字幕| 国产精品 欧美亚洲| 欧美av亚洲av综合av国产av| 国产视频内射| 好男人电影高清在线观看| 成年人黄色毛片网站| 欧美精品亚洲一区二区| 午夜亚洲福利在线播放| 久久久久久久久中文| 欧美一区二区国产精品久久精品 | 手机成人av网站| 欧美日韩亚洲综合一区二区三区_| 丁香六月欧美| 亚洲激情在线av| 婷婷丁香在线五月| 国产男靠女视频免费网站| 国产成人av教育| 两个人的视频大全免费| 成年人黄色毛片网站| 夜夜看夜夜爽夜夜摸| 又大又爽又粗| 国产精品一区二区三区四区免费观看 | 精品一区二区三区视频在线观看免费| 亚洲国产欧美一区二区综合| 国产免费男女视频| 亚洲性夜色夜夜综合| 午夜日韩欧美国产| 国产午夜精品久久久久久| 高潮久久久久久久久久久不卡| 免费av毛片视频| 国产精品 欧美亚洲| 欧美黄色片欧美黄色片| 亚洲av成人精品一区久久| 久久99热这里只有精品18| 亚洲国产中文字幕在线视频| 国产不卡一卡二| 老司机午夜福利在线观看视频| av天堂在线播放| 黄色片一级片一级黄色片| 一区二区三区高清视频在线| 12—13女人毛片做爰片一| 中国美女看黄片| 日韩中文字幕欧美一区二区| 国产成人精品久久二区二区91| 国产午夜精品论理片| 久久中文字幕一级| 91九色精品人成在线观看| 欧美日韩中文字幕国产精品一区二区三区| 91老司机精品| 美女黄网站色视频| 国产精品一区二区免费欧美| 久久精品aⅴ一区二区三区四区| 亚洲一区高清亚洲精品| 无遮挡黄片免费观看| 久久久国产成人精品二区| 美女免费视频网站| 人妻久久中文字幕网| 啦啦啦免费观看视频1| 性欧美人与动物交配| 少妇人妻一区二区三区视频| 久久久精品国产亚洲av高清涩受| 亚洲成人久久性| av福利片在线| 99久久久亚洲精品蜜臀av| 最近在线观看免费完整版| 一进一出抽搐gif免费好疼| 操出白浆在线播放| 国产爱豆传媒在线观看 | 一边摸一边做爽爽视频免费| 在线观看午夜福利视频| av福利片在线观看| 免费看十八禁软件| 亚洲电影在线观看av| 制服诱惑二区|