Bias adjustment to preserve changes in variability: the unbiased mapping of GCM changes

ABSTRACT Standard quantile mapping (QM) performs well, as a bias adjustment method, in removing historical climate biases, but it can significantly alter a global climate model (GCM) signal. Methods that do incorporate GCM changes commonly consider mean changes only. Quantile delta mapping (QDM) is an exception, as it explicitly preserves relative changes in the quantiles, but it might present biases in preserving GCMs changes in standard deviation. In this work we propose the unbiased quantile mapping (UQM) method, which by construction preserves GCM changes of the mean and the standard deviation. Synthetic experiments and four Chilean locations are used to compare the performance of UQM against QDM, QM, detrended quantile mapping, and scale distribution mapping. All the methods outperform QM, but a tradeoff exists between preserving the GCM relative changes in the quantiles (QDM is recommended in this case), or changes in the GCM moments (UQM is recommended in this case).


Introduction
Outputs from global climate models (GCMs) have biases and cannot be used directly at local scales to assess the impacts of climate change (Wood et al. 2004). In particular, using biased climate time series in hydrological modelling produces unrealistic simulations of the physical processes that are involved in the natural system (Maraun and Widmann 2018). Thus, a bias adjustment (BA) method becomes essential in the downscaling process to adjust the statistical properties of the raw GCM outputs using surface climatic data to be consistent with local climate conditions (Maraun 2013, Maraun andWidmann 2018). The performance of the BA depends on the method (i.e. theoretical basis, assumptions and simplifications), the geographical location, the available observations and the selected climate simulations (Chen et al. 2013a(Chen et al. , 2013b. A wide variety of BA methods are reported in the literature (e.g. Maraun et al. 2010, Teutschbein and Seibert 2012, Maraun 2016. Because of its simplicity and easy implementation, quantile mapping (QM) (Panofsky and Brier 1968) is one of the most widely used BA methods (Andres et al. 2014, Kim et al. 2016; it was originally used to adjust outputs of GCM for hydrological forecasts (Wood et al. 2002). The method is used not only for bias adjusting projected precipitation and temperature time series (Ines and Hansen 2006, Piani et al. 2010a, 2010b, Themeßl et al. 2012, Thrasher et al. 2012 but also for assessing climate change impacts (e.g. Maurer et al. 2009, Maraun et al. 2010, Gutmann et al. 2014, Walton et al. 2017, Shi et al. 2018, evaluating the response of hydrological models (e.g. Teutschbein and Seibert 2012, Chen et al. 2013b, Teng et al. 2015, and estimating changes in soil erosion (Eekhout and de Vente 2019), among other uses. This method takes percentiles from the cumulative distribution function (CDF) of the GCM output and corrects them using a transfer function built with the values of the historical GCM and historical observed period. When consistency between climatic variables such as precipitation and temperature is needed, there is also multivariate BA based on QM (Piani andHaerter 2012, Tootoonchi et al. 2022).
There are concerns regarding the capability of the QM method in reproducing projected changes in the outputs from raw GCMs (Hageman, 2011, Ehret et al. 2012, Maurer et al. 2013, Maurer and Pierce 2014, Cannon et al. 2015, Pierce et al. 2015, Hnilica et al. 2017, Chadwick et al. 2018. Some researchers consider the inability of the QM method to reproduce projected changes to actually be a beneficial part of the BA itself (Boberg and Christensen 2012, Bellprat et al. 2013, Gobiet et al. 2015, Ivanov et al. 2018. They argue that GCMs may project implausible climate change signals. In this case, the modification induced by the standard QM is unrealistic and not physically justified (Maraun et al. 2017). Other researchers consider that BA methods should preserve the GCM changes (e.g. Lanzante et al. 2020, Tegegne andMelesse 2021), particularly when the changes of the GCMs are credibly simulated (Maraun et al. 2017). This latter approach will be adopted in our study. For a more detailed discussion on this issue the reader is referred to Maraun et al. (2017) and Maraun and Widmann (2018).
The discussion regarding whether BA should retain stationarity assumptions to the future has also extended to multivariate BA, and researchers hold different positions regarding whether historical relationships between climatic variables (e.g. the correlation between precipitation and temperature) or spatial properties should be preserved or not in the future (Tootoonchi et al. 2022). In an overall comparison of different univariate and multivariate BA methods, François et al. (2020) concludes that before choosing a method one should know the relevant statistical properties to be evaluated in the study, and also analyse whether they can be properly obtained from climate projections.
A variety of quantile mapping methods have been proposed to preserve raw GCM changes, which have advantages and drawbacks as compared to the standard QM approach (note that we will use the acronym QM only when referring to the traditional quantile mapping bias adjustment method, but not to the concept of quantile mapping). For example, detrended quantile mapping (DQM) (Bürger et al. 2013) applies the QM BA to detrended climatic data. After this correction, the trend is reintroduced additively for temperature data or multiplicatively in the case of precipitation. DQM performs well as it is able to preserve long-term raw GCM changes of the mean. The DQM approach is similar to that applied to temperature by Wood et al. (2002Wood et al. ( , 2004 who detrended temperature data before applying the QM, and is also similar to other QM or BA methods that preserve the raw GCM changes (Walton et al. 2020). Quantile delta mapping (QDM) (Cannon et al. 2015) is similar to the DQM approach, but it was built specifically to preserve the raw GCM changes of the quantiles. However, QDM does not consider changes in other statistical moments like the standard deviation. More recently, Switanek et al. (2017) proposed the scaled distribution mapping (SDM) method, which works at a daily scale and is able to preserve the raw GCM changes in the mean and in higher order moments (i.e. standard deviation and skewness).
Nevertheless, none of the previous QM methods formally preserve the relative changes of the GCM standard deviation through time. Moreover, Switanek et al. (2017) studied the impact of the depicted methods on the relative changes of the standard deviation, as well as other moments, and found differences in performance among the methods. In fact, their analysis shows that GCM biases are not constant through time, an aspect that needs to be further studied. In historical records, long droughts might partially be explained by significant changes in the standard deviation of precipitation; an example of this is the megadrought in Central Chile since 2010 (Barría et al. 2021). All these issues are relevant as standard deviation plays a key role in shaping the distributions used to represent the simulated climate variables. With other downscaling methods, researchers have tried to preserve changes in the standard deviation (Chadwick et al. 2018); nonetheless, QM-based methods are not necessarily capable of retaining such changes.
The objective of this paper is to present the newly developed unbiased quantile mapping (UQM) method, a QM method for BA that, through a fairly simple implementation, is able to preserve, by construction, both the mean and the standard deviation changes of raw GCMs, for precipitation. The developed method is compared against existing QM methods (i.e. QM, DQM, QDM, SDM) using a synthetic example. Furthermore, for a real application example, we chose four locations with different climate in Chile. The study is organized as follows: section 2 summarizes the existing QM methods (i.e. QM, DQM, QDM, SDM) and develops the UQM method, while section 3 describes the study area and the climatic data used in the assessment. Results are given and discussed in section 4, and the main conclusions are presented in section 5.

Methodology
Four of the most common QM approaches are the standard QM, DQM, QDM, and SDM, which are summarized in the following subsection. Subsequently, the new proposed UQM approach is presented in detail. The equations for all the QM methods are developed for precipitation, and thus consider multiplicative factors. Finally, in this section, we describe the implementation of the BA and the assessment methods that are used.

Standard QM
The QM method (Panofsky andBrier 1968, Wood et al. 2002) corrects the GCM bias with CDFs of the historical period of the GCM (F m;h ) and the historical observed period (F o;h ), using the following equation: where x m t ð Þ is the GCM output (subscript m) at time t, and x m t ð Þ is the bias-adjusted precipitation or temperature data. The subscript h is used to denote the historical periods, while the subscript o denotes observed data.

Detrended quantile mapping
The DQM method (Bürger et al. 2013) starts by extracting the long-term trend from the GCM climate data, after which the QM is applied to the detrended GCM data to finally reimpose the trend by scaling and rescaling, using a methodology similar to the delta change or the delta factor in the final step. The following equation summarizes the procedure: where x m t ð Þ is the mean of the GCM climate output within the time window ending in year t, while x m;h is the mean of the GCM historical period. Note that x m t ð Þ changes over time, so the factor used to detrend x m t ð Þ also changes over time.

Quantile delta mapping
QDM was designed to explicitly preserve the relative changes in precipitation quantiles. QDM was developed by Cannon et al. (2015), and was built after the quantile perturbation, quantile delta change and DQM methods (Bürger et al. 2013). Cannon et al. (2015) also showed that the equidistant CDF matching algorithm of Li et al. (2010) is equivalent to an additive QDM method. QDM applies relative changes to the quantiles, but not explicitly to the moments (e.g. mean, standard deviation). The first step of QDM is to compute the non-exceedance probability of x m from GCM climate data (Prob m t ð Þ), which is obtained by using a time-dependent CDF (F t ð Þ m;p ), with the subscript p denoting the projected period, as shown in Equation (3): x m is bias adjusted by a QM method and a time-dependent quantile delta factor (Δ m t ð Þ), to which it is subsequently applied. Δ m t ð Þ is estimated as follows: Then, Prob m t ð Þ is used in the inverse CDF of the observed historical data (F À 1 o;h ), which is then scaled by multiplying by the delta factor obtained in Equation (4): SDM (Switanek et al. 2017) was design to preserve the raw GCM changes while performing a BA at a daily scale; it performs slightly better than QDM and DQM in this regard. Like QDM, SDM preserves relative changes in precipitation quantiles. Although SDM was originally designed to bias adjust daily climatic data, Casanueva et al. (2020) used the SDM implemented in the pyCAT package, which works at a daily or monthly scale.

Scaled distribution mapping
In the method, precipitation values below a 0.1 mm threshold are removed, to then obtain the expected number of rainy months (R M m ) for the bias-adjusted GCM (step 1). The following equation is used for this purpose: where RM m is the number of rainy months of the projected raw GCM; RM m;h and RM m;o are the rainy months of the historical raw GCM and observed precipitation, respectively; and TM o;m and TM o;h are the total number of months of the historical raw GCM and observed climate. Then a gamma distribution is fitted to each precipitation series (i.e. historical observed, historical raw GCM and projected raw GCM) using the maximum likelihood method (step 2).
Next, a scaling factor vector SF (vector elements are presented in bold to simplify notation) similar to the timedependent quantile delta factors (Δ m t ð Þ, Equation (4) is estimated (step 3). However, GCM projected precipitation values must be sorted in ascending order to estimate SF. Subsequently, the recurrence interval vectors (RI) are estimated for the three sorted precipitation vectors (x) (step 4): where F½ � is the CDF for each respective precipitation series (x). When the number of rainy months differs between precipitation series, the RI of the historical observed and GCM are linearly interpolated to obtain the same number of months as the projected GCM (RI m ).
A scaled recurrence interval vector (SRI) is estimated as follows (step 5): where RI o;h and RI m;h are the RI of the interpolated historical observed and GCM series, respectively (note that the operator "�" is the Hadamard product, which multiplies, element by element, two equally sized vectors or matrixes, preserving their dimensions). The SRI is used to obtain a probability scaled vector (Prob scaled ): The Prob scaled values are sorted in ascending order, and are then used to obtain x m , the bias-adjusted precipitation (step 6): Eventually, x m must be linearly interpolated to ensure that it has R M m rainy months. Finally, x m values are placed back into the modelled time series in the correct temporal location (step 7), process which is done in descending order. More details regarding the SDM method are given by Switanek et al. (2017) (temperature equations of SDM are presented in Text S1 in Supplementary material).
SDM was originally developed to correct a projected time window, while preserving the original trend in the bias-adjusted series. Here, the SDM was adapted to a continuously changing climate instead. For this purpose, the final bias-adjusted precipitation series x m ðtÞ, is obtained by applying the SDM several times, and keeping the last year of each time window t. This procedure ensures a fair comparison against the other QM methods, as the same climatic information is used in each case.

Proposed BA methods: UQM
UQM is a QM BA that by construction preserves the mean and standard deviation of raw GCM changes. This is done by extracting the changes in the mean and the standard deviation of the GCM output, and then applying a QM BA method which incorporates these changes. A similar approach is used in the weather generator-based (i.e. statistical models that are based on randomly generating synthetic data and are required to replicate historical climate, as well as incorporate GCM changes when used to simulate climatic projections) unbiased mapping of GCM changes to local stations proposed and used by Chadwick et al. (2018Chadwick et al. ( , 2019Chadwick et al. ( , 2021. That unbiased mapping approach is a two-step process in which annual long-term trend percentiles from a GCM group are extracted first, and then climate time series at a local station are generated around them. In UQM, for each GCM, the precipitation change, or delta factor of the mean, is obtained by estimating the change between the moving averages of the GCM climate output (x m t ð Þ), associated with a moving time window, and the average from the GCM historical period (x m;h ). This change of the mean is called ðΔ μ ðtÞÞ: An analogous process is used to extract the change in the standard deviation of the GCM climatic output (Δ σ t ð Þ): where c x m t ð Þ is the GCM climatic output moving standard deviation (i.e. the standard deviation is estimated several times for a moving time window), and d x m;h is the standard deviation from the historical period. Then, Δ μ t ð Þ and Δ σ t ð Þ are used to parameterize time-dependent expressions of the mean, μ � t ð Þ, and standard deviation, The bias-HT as follows: where Prob m t ð Þ is the non-exceedance probability of x m from the GCM, just as in QDM (Equation 3), and θ � o;h t ð Þ is the parameter set. In UQM the inverse CDF (which can be that of any distribution) of the observed historical data becomes . Although a parameter set is always used to evaluate the probability distributions associated with the other BA methods, the parameter set θ � o;h t ð Þ is explicitly noted in Equation (15), given the relevant fact that its time dependence is essential to the UQM method. θ � o;h t ð Þ is estimated with the method of moments using μ � t ð Þ and σ � t ð Þ, hence they change with time according to Δ μ t ð Þ and Δ σ t ð Þ (Equations 13 and 14).
The same process is followed for temperature but the change in the mean, Δ μ t ð Þ, and standard deviation Δ σ t ð Þ, are estimated as the absolute changes. Hence, μ � t ð Þ and σ � t ð Þ are obtained by adding Δ μ t ð Þ and Δ σ t ð Þ to the historical observed mean (μ) and standard deviation (σ) of temperature, respectively. Text S2 (in the Supplementary material) shows the development of the corresponding expressions for the BA of temperatures. For a schematic representation of the UQM method we refer the reader to a short video in the Supplementary material in which the methodology is explained.

Implementation and assessments of the UQM method
In section 2.3 we will explain how the UQM is implemented, and how the UQM assessment and evaluation exercises are performed. The assessment exercises are divided into a synthetic exercise (section 2.3.1) and a real case exercise, which is performed in several Chilean basins (section 2.3.2). Note that the methodology assessment methods described in sections 2.3.1 and 2.3.2, primarily focus on the ability of the UQM in preserving raw GCM changes. The reason for this is that UQM, as well as other trend-preserving methods (e.g. DQM and QDM), assumes there is no trend, signal, or change (i.e. Δ μ and Δ σ ) in the historical period, hence the UQM becomes the standard QM in that period.
The UQM is implemented at a monthly and annual level in this study. These scales are sufficient to evaluate whether the UQM is an improvement compared to other QM-based methods in terms of preserving raw GCM changes (i.e. Δ μ and Δ σ ). Further studies should analyse possible adaptations of UQM to a daily scale.
In the synthetic exercise a two-parameter gamma distribution is used in Equations (1), (2), (5), (10), and (15), which is fitted with the method of moments. On the other hand, for the real case we used the coding package climQMBC (for more details, see the Code availability section at the end of the article), which has the newly developed UQM, as well as the other QM methods used in this study (QM, DQM, QDM, and SDM). The climQMBC uses strictly nonnegative distributions in precipitation (i.e. two-parameter gamma, log normal, and log Pearson Type III). The lowest Kolmogorov-Smirnov test is used to choose between these nonnegative distributions (further detail available in the documentation of climQMBC, at https://github.com/saedoquililongo/climQMBC).
Fitting probability distribution functions in locations with only a couple of months with positive precipitation might be challenging. For example, for a given month, 80 to 90% of them within a record period of several years might have zero precipitation (e.g. northern Chile). Then, if a monthly BA is applied using a 30 years record, maybe a few years (e.g., 3-5 years) with non-zero precipitation values for a given month may have to be used, to fit a censored shift distribution or mixed distribution. Because of this challenge, the climQMBC package in the UQM method uses an approach we call "random near-zero replacement" to handle no precipitation months, while fitting distributions. Each no precipitation month is replaced with a random number that is barely positive (i.e. by default it uses a random uniform number in the range 0; 1=100 ½ � mm), but for practical purposes it is physically zero. As a final step, after the BA, one eliminates the precipitations under a threshold (by default the threshold is 1 mm per month). Even in the case where all the months have zero precipitation values, this approach allows us to fit a distribution, which for practical purposes will give dry months that at the end will be replaced with zeros. The "random near-zero replacement" proposed here is similar to the "singularity stochastic removal" proposed by Vrac et al. (2016) and the jittering algorithm used in the copula approach (Pappadà et al. 2017) (the difference between "singularity stochastic removal" and "random nearzero replacement" is explained in Text S3).

Synthetic exercise
The synthetic example follows the one developed by Maurer and Pierce (2014) and Cannon et al. (2015): we used twoparameter gamma distributions to generate observed climate and GCM outputs. The gamma probability density function with shape parameter β and scale parameter α is given by: where Γ � ð Þ is the gamma function. Using the method of moments, the parameters are estimated from the mean μ and standard deviation σ using the relationships μ ¼ αβ and σ ¼ ffi ffi ffi β p α. The synthetic exercise uses the gamma distribution to generate three sets of climate data: (1) observed historical, (2) GCM historical, and (3) GCM future. Then each BA is calibrated and applied, to measure its performance in terms of its ability to preserve changes in the mean, standard deviation and quantiles of raw GCM changes.
In particular, four synthetic exercises are performed: (1) an illustrative exercise, which consists in having historical observed data and GCM data that has bias in the standard deviation, but not in the mean, and a GCM future where the mean increases while the standard deviation decreases.
(2) A more realistic exercise that starts with the same historical observed and GCM data as the "illustrative exercise," but in this case the GCM future has an increasing mean and standard deviation. (3) A transient time-dependent exercise, which has the same historical observed and GCM values, but linearly changes in time and after 100 periods of time arrives at the future realistic GCM values. (4) A sensitivity analysis having the same historical observed data and a GCM that increases its mean in the future, while in the standard deviation it presents the sensitivity analysis having a GCM (in historical and future period) with a standard deviation lower than, the same as or higher than the historical observed.
The illustrative exercise is developed to test the BAs under extreme raw GCM changes, while the more realistic exercise is a case with a more common change. The transient time-dependent exercise was developed to test the capability of each method in preserving changes in time. Finally, the sensitivity analysis was performed to understand how the performance of each method is affected under different scenarios. The results of these exercises are presented in section 4.1.

Real case exercise
For the real case exercise, we have used four Chilean basins: Limarí, Maipo, Maule and Valdivia ( Fig. 1; a detailed description of each basin is given in section 3). The basin average climate from the CAMELS-CL (Catchment Attributes and Meteorology for Large Sample Studies, Chile Dataset) database (Alvarez-Garreton et al. 2018) is used. Note that for this real case exercise average basin monthly precipitation values, with basin sizes ranging from 6600 to 21 000 km 2 (Table 1), have been used to avoid the variability inflation issues of performing a QM at a basin scale (Maraun 2013). The time period considered is 1979-2014, which is restricted by the initial year of CAMELS-CL (1979), and the final year (2014) of the historical GCMs outputs from the sixth Coupled Model Intercomparison Project (CMIP6) (Eyring et al. 2016).
The GCM climate projections for each basin are estimated first by making a spatial average of the gridded data that each GCM has over the shape file of each basin from Fig. 1. Once we had the GCM projection time series per basin, the previously described BA methods are applied. Finally, the changes in the mean and standard deviation, from raw GCM data, as well as for the bias-adjusted GCM data were estimated with Equation (11) and (12). Given the great number of climate change projections, three representative percentiles (i.e. 25 th , 50 th and 75 th percentile from all the GCMs projection) were estimated for each scenario (i.e. SSP2-4.5 and SSP5-8.5). Since each of the GCMs is assumed to have the same weight in the results, the inverse of the number of realizations of each GCM is used as a weight factor for each one of the GCM projections to estimate the percentiles presented in the results section (e.g. given that for ACCESS-ESM1-5 we are using three realizations (see Table S1), each projection has a weighting factor of 1/3), as suggested by Chadwick et al. (2018).
For the real case exercise, we have tested the capacity of each BA method to preserve changes in the mean and standard deviation. We performed two relevant analyses: (1) A qualitative visual comparison of the capability of each BA method in preserving changes in the mean and the standard deviation is performed.
(2) A quantitative comparison is performed using the mean absolute error (MAE) and the percentage of time the absolute error (AE) exceeds the 0.05 threshold of each BA method capability of preserving raw GCM changes in the mean and the standard deviation. The MAE is obtained by averaging AE: where y i is the raw GCM change in the mean or the standard deviation in the case i, while x i is a specific BA change in the mean or standard deviation in case i. Then the MAE is: where the AE is averaged over the total n cases (detailed results are presented in Tables S2-S5).
On the other hand, the percentage of time in which the AE of the changes in the mean or standard deviation exceeds the 0.05 threshold is: with T i being a count of each time the AE of a specific BA change in the mean or standard deviation exceeds the 0.05 difference against the raw GCM changes for case i. Then T i is used to estimate the percentage of time the threshold is exceeded (Prc Time AE > 0:05) (see detailed results in Tables S2 and S6-S8).

Description of Chilean basins
In addition to an assessment using the synthetic case, we have chosen four basins in Chile for a real case application. In this application the UQM approach was applied and compared against the other QM methods in central and southern Chile (between 30 and 40°S), as depicted in Fig. 1. According to the Köppen-Geiger climate classification for Chile (Sarricolea et al. 2017), the Limarí basin in the northern part of central Chile has a semi-arid climate, with a mean annual precipitation of 240.8 mm. The Maipo and Maule basins have a Mediterranean climate; the former has a mean annual precipitation of 623.3 mm, whereas the latter has a mean annual precipitation of 1294.3 mm. Finally, the Valdivia basin located in southern Chile has a marine west coast climate, with mean annual precipitation of 2570.7 mm. All these basins have the Pacific Ocean to the west and the Andes Mountains to the east, with ~150 km of separation between the coast and the Andes Mountains. Annual temperature in these basins tends to decrease with latitude, while the annual rainfall tends to increase (Quintana and Aceituno 2012) and its variability tends to decrease. Basic information about the study basins is given in Table 1. The Limarí, Maipo and Maule basins have strong seasonal precipitation regimes due to the influence of the South East Pacific Anticyclone and the inter-annual variability of El Niño Southern Oscillation (ENSO) (Saavedra and Foppiano, 2002, Montecinos and Aceituno 2003, Falvey and Garreaud 2007. These forces generate a high-variability anomaly for precipitation, as studied by Jurković and Pasarić (2013), which explains the high coefficient of variation (CV) values shown in Table 1. In the winter season (June, July and August) 68%, 60% and 55% of the annual rainfall takes place in the Limarí, Maipo and Maule basins, respectively, while the summer season (January, February and March) is drier (Fig. 2). On the other hand, the Valdivia basin is located in a transitional region where precipitation patterns, which are more stable through the year, are mainly controlled by the influence of both Antarctic Oscillation (AAO) and ENSO climate forcing (Montecinos and Aceituno 2003, Aravena and Luckman 2009, González-Reyes and Muñoz 2013. Winter precipitation in the Valdivia basin corresponds to ~45% of the annual total, with a significant portion of it occurring during the summer. Furthermore, there is a lower thermic oscillation in the basin (Quintana and Aceituno 2012), as shown in Fig. 2. By 2100, and under the RCP 8.5 scenario, the average temperature is expected to increase by around 4.2°C, 3.8°C and 3.6°C in the Limarí, Maipo, and Maule basins respectively, while annual precipitation is expected to significantly decrease, with average reductions of approximately 24%, 28% and 29%, respectively (Chadwick et al. 2018). Results were obtained using more than 45 GCM projections and a weather generation-type downscaling method (more detail in Chadwick et al. 2018). These precipitation reductions are so severe that by year 2090, the probability of observing three consecutive years with precipitation below the historical mean is expected to increase from 23% to 49% in Limarí, 24% to 54% in Maipo, and 20% to 59% in Maule basins. More details about the expected changes in temperature and precipitation for the three study basins are provided by Chadwick et al. (2018). Regarding the Valdivia River basin, Araya-Osses et al. (2020) report climate change projections from six different CMIP5 GCMs under the RCP 8.5 scenario. According to these projections, the winter season will be warmer and drier by the end of the century (2081-2100), while precipitation might decrease by up to ~55%. These projections agree with those reported in a previous study that used CMIP3 GCMs under the A2 scenario (Garreaud 2011).

Results and discussion
In this section we compare the ability of the UQM method to preserve the raw GCM changes (i.e. changes in the mean and in the standard deviation) against that of other QM approaches. For this purpose, the section is divided into two parts. First, we adopt a simple synthetic example to identify and assess the pros and cons of each method. Second, all the QM methodologies are applied to the study basins in Chile, for which their ability to preserve raw GCM changes under different climates and periods of time are shown.

Synthetic example
The synthetic illustrative and realistic examples generated observed historical data using values of μ ¼ 30 and σ ¼ 16 for the gamma distribution ( Fig. 3(a) and (c)); the historical data from the GCM follow a gamma distribution with μ ¼ 30 and σ ¼ 10, while for the future GCM data we developed an illustrative extreme example with μ ¼ 42 and σ ¼ 2:5 ( Fig. 3(a)) and a more realistic example with μ ¼ 40 and σ ¼ 12 (Fig. 3(c)). Note that in the illustrative example quite different σ values were used between the projected and historical GCM data (i.e. 2.5 vs. 10) to evaluate the ability and robustness of the different QM methods in preserving changes in σ. Figure 3(b) and (d) presents the relative changes from Fig. 3(a) and (c), respectively. The third exercise consists in the transient time-dependent synthetic example that was developed in Fig. 3(e). This example starts with 30 years of stationary climate with a historical observed Δ σ and σ ¼ 16, a historical GCM with μ ¼ 30 and σ ¼ 10. After these historical stationary years, the timedependent exercise has a future GCM where the μ and σ linearly change to the values of μ ¼ 40 and σ ¼ 12 by the year 130 (i.e. the transient period is developed between the years 31 and 130). Note that Fig. 3(e) shows the transient time-dependent exercise with the values from Fig. 3(c). The objective in Fig. 3(e) is to obtain the observed historical in the stationary part (years 1-30), and the historical with the change from the raw GCM (blue line) in the transient future period (years 31-130).
The raw GCM changes of the synthetic illustrative example (Fig. 3(a)) are the percentage of change between the projected and the historical values of the mean and standard deviation of the GCM data (i.e. Δ μ ¼ 40% and Δ σ ¼ À 75%). If the objective of the QM method is not only to preserve quantiles, but also to maintain Δ μ and Δ σ , then μ ¼ 42 and μ ¼ 42 after applying the QM. For the more realistic synthetic example (Fig. 3c) the changes in the GCM are Δ μ ¼ 33:3% and Δ σ ¼ 20%, hence the expected μ and σ after applying the BA should be 40 and 40 , respectively. Note that this is the equivalent of having a QM BA method and a delta change or delta factor method able to preserve the mean and standard deviation at the same time.
After applying the different QM methods, we find the UQM is the only method able to obtain the exact values of μ ¼ 42 and σ ¼ 4 for the illustrative example (UQM F in Fig. 3(a)), and μ ¼ 40 and σ ¼ 19:2 in the more realistic example (UQM F in Fig. 3(c)). This is expected, given that by construction UQM was built to preserve Δ μ and Δ σ . The second closest is the DQM, which has values of μ ¼ 40:6 and σ ¼ 3:9 in the illustrative example (Fig. 3(a)) and σ ¼ 39:7 and σ ¼ 19:1 in the realistic example (Fig. 3(c)). The QDM in the illustrative example ( Fig. 3(a)) has a mean of 39.8, which is quite close to 42, but a standard deviation of 10.7, which is far from 4. In contrast, in the realistic example (Fig. 3(c)) QDM has values of σ ¼ 19:1 and σ ¼ 19:9, which is close to μ and somewhat overshooting σ. The SDM presents results quite similar to those of the QDM, with these being the two QM methods with the worst performance in preserving Δ σ . Note that QDM and SDM are similar because the scaling factor (SF) of SDM is built after the quantile delta factor (Δ m t ð Þ) of QDM. QM, as reported before Pierce 2014, Cannon et al. 2015), has trouble preserving Δ μ , with a mean value of 49.2 in the illustrative example ( Fig. 3(a)) and 46.6 in the more realistic one (Fig. 3(c)). QM outperforms QDM and SDM in preserving Δ σ in the illustrative example ( Fig. 3(a)), as the value of σ ¼ 4:4 is quite close to 4, but gives a worse result with a value of σ ¼ 21:2 relative to the 19.2 objective in the more realistic example (Fig. 3(c)). Figure 3(b) and (d) compares the capability of different QM methodologies in preserving the relative change in the 25 th , 50 th , 75 th , 95 th and 99 th quantiles between the historical and future periods. The QM methods capable of preserving this type of changes are QDM and SDM. Both the UQM and the DQM perform similarly, but with some departures from the 1:1 line. Finally, the standard QM comes last, being farthest away from the objective, presenting a quite large error in the more realistic example (Fig. 3(d)).
In the transient synthetic example (Fig. 3(e)), the transient changes in time are evaluated, with the QM grossly overshooting the objective (blue line). All the trendpreserving QMs, including the proposed method (i.e. DQM, QDM, SDM, and UQM), tend to slightly overshoot. This occurs because the trend-preserving QM methods (DQM, QDM, SDM and UQM) assume that a moving window is reasonably stationary in the future. For example, QDM and UQM estimates Prob m t ð Þ in Equation (3), using a moving window that is assumed to be reasonably stationary, to use a CDF in the estimation of probabilities. However, because the climate is changing in the moving window used to estimate Prob m t ð Þ, these probabilities are slightly biased. We call this the "pseudo stationary bias", and it can be detected in a synthetic example but is quite challenging to detect in a real example.
A possible solution for the pseudo stationary bias, which can be applied in QDM and UQM, is to detrend the climate being used to estimate Prob m t ð Þ in Equation (3). We have called these the QDM* and UQM*, with the results from UQM* in the transient synthetic exercise presented in Fig. 3(e). UQM* is able to obtain the desired result in the transient synthetic exercise and avoid the pseudo stationary bias, but it assumes that the raw GCM trend is linear, which in several cases is not true. Because the assumption of linear GCM trend might be unrealistic, we still recommend UQM and QDM over UQM* and QDM*.
An additional comparison between the QM methods using the synthetic approach, the sensitivity analysis, is presented in Fig. 4. In this case the observed historical data follow a gamma distribution with μ ¼ 30 and σ ¼ 16. On the other hand, μ ¼ 30 and μ ¼ 42 for the historical and projected GCM data, respectively, while σ is set not to change in time (i.e. Δ σ ¼ 0%) and vary, due to a sensitivity analysis in intervals of 10%, between −60% and +60% around a reference value of 16 (i.e. from 6.4 to 25.6).
As expected, when comparing the ability of QM methods to preserve changes in Δ μ (Fig. 4(a)) and Δ σ (Fig. 4(b)), the best method is UQM, which by construction preserves both changes. The second-best method is DQM, followed by QDM and SDM (these two with identical results), with the standard QM coming in last. Note that QDM and SDM have trouble with the change in standard deviation (Fig. 4(b)), although the historical and projected GCM series have the same standard deviation (Δ σ ¼ 0%). On the other hand, in terms of the ability to preserve the changes in the different quantiles and the median (Fig. 4(c-f)), QDM and SDM are the best QM methods, followed by UQM and DQM. The standard QM is clearly the one with the worst performance.
Interestingly, none of the analysed QM methods can preserve jointly both Δ μ and Δ σ as well as the relative change for the different quantiles. Because of its construction, the UQM method is able to preserve the moments, but cannot preserve the relative changes of the quantiles; the opposite happens with the QDM and SDM methods. Hence, the researcher would have to choose between preserving one or the other, and could not preserve both. We suggest to focus the QM on Δ μ and Δ σ , for two reasons. First, the moments are commonly used to describe any distribution, hence focusing on getting the mean and the standard deviation right helps to obtain the future probability density function. Second, given they are obtained from the mean and the standard deviation of all the values of a time span, the values of Δ μ and Δ μ should be more stable, and thus better estimated than the relative change in extreme quantiles.

Comparison of the BA performance in a real exercise
Precipitation data for the four study basins (i.e. Limarí, Maipo, Maule, and Valdivia) were bias adjusted using the five QM methods (i.e. QM, DQM, QDM, SDM and UQM) at monthly and annual scales, under SSP2-4.5 and SSP5-8.5, for each GCM from Table S1. In this section we present, compare, and assess the ability of each QM method to preserve Δ μ and Δ μ , for different time periods, SSP-RCPs and locations. Only winter months (i.e. July, August and September) are considered in the analysis, as zero precipitation has been recorded in other seasons for several months in the Limarí and the Maipo river basins, and hence winter months are better suited to test the capability of UQM in preserving raw GCM changes. Figure 5 shows Δ μ (first row) and Δ σ (second row) of each QM method against the GCM-projected change of the 50 th percentile of the GCMs, under SSP5-8.5 for annual precipitation. The bias of each QM methodology is presented as the difference in Δ μ (Fig. 5, third row) and Δ σ (Fig. 5, fourth row) between the BA method and the raw GCM changes.
In Fig. 5 the QM method presents the worst performance in preserving Δ μ , with biases that exceed 10% for Limarí and Maipo (Fig. 5, third row), while in Maule and Valdivia, the five BA methods (QM, DQM, QDM, SDM and UQM) present similar performance. The SDM has trouble preserving Δ σ at Maule and Valdivia, with biases greater than 40% in Valdivia (Fig. 5, second and fourth row). The performance of all the QM methods is constant neither in time nor in location. All the methods present some bias in Δ μ between 2060 and 2070 in Limarí, bias that disappears by the year 2100 (Fig. 5, first row and column). Limarí and Maipo present greater bias in Δ μ than Maule and Valdivia (Fig. 5, first row). The first row of Fig. 5 shows a precipitation reduction (Δ μ ) in Valdivia by the end of the century by ~30% under SSP5-8.5, which is much less severe than the ~55% reduction in annual precipitation under RCP 8.5 previously reported by Araya-Osses et al. (2020), or the ~60-70% under the A2 scenario reported by Garreaud (2011). The former study considered only six GCMs and an analogue method as downscaling (Zorita and Von Storch 1999), while the latter considered the regional climate model (RCM) PRECIS (Providing REgional Climates for Impacts Studies), which only considers the Hadley Center GCM. The difference in Valdivia's projected precipitation illustrates the relevance of the GCM and downscaling selection, and motivates the careful analysis and assessment of the different uncertainty sources associated with future climate estimation at local scales.
The same results as in Fig. 5 are presented in Figs 6-8, but at a monthly scale for July, August and September, respectively. These figures were built using a monthly BA. Again, the QM method has more trouble preserving Δ μ and Δ σ , while the other methods present similar results. Overall, all the QM methods have their worst performance in the Limarí basin, an intermediate performance in the Maipo basin, and the best in Maule and Valdivia (further analysis regarding the performance of BA in different climates is presented below in the discussion of Table S2). Furthermore, no clear worsening or improvement in their performance is observed over time (Figs  6-8). Both Δ μ and Δ σ decrease with time, with a signal that is Figure 5. Changes in annual precipitation measured in Δ μ (first row) and Δ σ (second row) of the 50 th percentile of GCMs under SSP5-8.5, for the raw GCM, QM, DQM, QDM, SDM and UQM, and the four study basins. Differences between each type of bias adjustment (BA) and the GCM projected change of Δ μ (third row) and Δ σ (fourth row) of the 50 th percentile of GCMs under SSP5-8.5 in each basin. stronger for the former (Figs 6-8, first and second row). This has important implications as the climate in Chile is projected to be not only drier but also less variable, with more unusual extremely wet years or months.
Similar results to those shown in Fig. 5, but for the 25 th and 75 th percentile of the GCMs under SSP5-8.5, are presented in the Supplementary material (Figs S1 and S2). For the winter months of July, August and September, the same results as those reported in Figs 6-8 are presented for the 25 th and 75 th percentiles in Figs S3-S8 (Supplementary material). Finally, the Supplementary material also contains the same results as Figs 5 to 8, but for SSP2-4.5 (Figs S9-S12). Overall, SDM has some trouble preserving Δ σ (e.g. Figs S1 and S2 in Valdivia), QM presents more issues preserving Δ μ , and all the methods present more trouble preserving Δ μ in Limarí and their best performance in the Maule and Valdivia.
The overall performance of all the BA methods presented in the previous figures is measured with the average MAE and average percentage of time in which the AE exceeds the 0.05 threshold (Table 2). These quantities were calculated for Δ μ and Δ σ of precipitation under SSP5-8.5, and correspond to the averages considering all the study basins and winter months (July, August, and September). The results for the 25 th , 50 th , and 75 th percentiles of the GCMs are presented, as well as the overall average of these three percentiles. The overall average shows that DQM, QDM, and UQM have similar errors in both measurements, and clearly outperform QM (i.e. a 35-40% reduction in the MAE of Δ μ and at least a 10% reduction in the MAE of Δ μ ). The best preservation of Δ μ and Δ σ measured with the average MAE is obtained by the UQM method, with a slight improvement compared to DQM and QDM.
When measuring the capability of each method in preserving Δ μ and Δ σ , the best performance measured with AE > 0.05 Figure 6. Changes in July's precipitation measured in Δ μ (first row) and Δ σ (second row) of the 50 th percentile of GCMs under SSP5-8.5, for the raw GCM, QM, DQM, QDM, SDM and UQM, and the four study basins. Differences between each type of BA and the GCM projected change of Δ μ (third row) and Δ σ (fourth row) of the 50 th percentile of GCMs under SSP5-8.5 in each basin.
is obtained with UQM (Table 2). UQM only presents 2% of the time in which AE > 0.05 for Δ μ , while it has a 9% for Δ σ . The second-best method is DQM, with 4% and 9% for Δ μ and Δ σ , respectively, which is followed closely by QDM with 4% and 11% for Δ μ and Δ σ , respectively. This shows, there is some gain in using UQM, over other trend-preserving BA methods, if one wants to preserve Δ μ and Δ σ . The QM comes last in preserving Δ μ , with 12% of the time in which AE > 0.05, and SDM presents some issues with Δ σ , having 20% of the time in which AE > 0.05.
To better understand the effect of the climatic conditions on the BA, Table S1 presents, for each study basin, the average MAE and average percentage of time AE exceeds the 0.05 threshold for Δ μ and Δ σ of precipitation under . In this case, these quantities correspond to the averages over the winter months (July, August, and September) and the different BA methods (QM, QDM, DQM, SDM, and UQM). These errors are presented for the 25 th , 50 th , and 75 th percentile of the GCMs, and also for the overall average of these percentiles. The BA methods have more trouble preserving Δ μ in the Limarí, followed by Maipo and Maule, and they present a really good performance in Valdivia. For example, the average MAE values of Δ μ go from 0.025, to 0.005 from Limarí to Valdivia, while the AE > 0.05 values, go from 15% to 0%, for the same basins. While in preserving Δ σ the worse basin is Limarí, followed by Maipo, Valdivia, and Maule last. The average MAE values of Δ σ go from 0.046 to 0.015 from Limarí to Maule, while the AE > 0.05 values go from 26% to 10%, for the Limarí to Maule. Although the objective of this study is to present the UQM method, knowing which locations present a greater challenge when applying a BA method is Figure 7. Changes in August's precipitation measured in Δ μ (first row) and Δ σ (second row) of the 50 th percentile of GCMs under SSP5-8.5, for the raw GCM, QM, DQM, QDM, SDM and UQM and the four study basins. Differences between each type of BA and the GCM projected change of Δ μ (third row) and Δ σ (fourth row) of the 50 th percentile of GCMs under SSP5-8.5 in each basin. Figure 8. Changes in September's precipitation measured in Δ μ (first row) and Δ σ (second row) of the 50 th percentile of GCMs under SSP5-8.5, for the raw GCM, QM, DQM, QDM, SDM and UQM and the four study basins. Differences between each type of BA and the GCM projected change of Δ μ (third row) and Δ σ (fourth row) of the 50 th percentile of GCMs under SSP5-8.5 in each basin.  (Table 1). Less variable records are expected to better characterize the local precipitation regime, which in turn is more stable, making the BA more robust. Nevertheless, this hypothesis should be further tested in future studies using other synthetic experiments, as well as more locations with different precipitation regimes. More details about the statistical bias for each BA method and month is presented in the Supplementary material (Tables S2-S7). In the future, increasing our knowledge regarding the relationship between climates and the performance of BA methods will help researchers understand and develop new BAs methods, and will give significant insights regarding how they should be tested.

Conclusions
BA methods are used to adjust the statistical properties of raw GCM outputs to be consistent with local climate conditions. All the BA methodologies focus on eliminating the GCM precipitation output bias of the mean and the standard deviation of the historical period. Some trend-preserving biasadjusting methods also focus on the mean changes from the raw GCM outputs, and only a few available methods centre their efforts in the entire distribution. Examples of the latter are quantile delta mapping (QDM, Cannon et al. 2015) and the scale distribution method (SDM, Switanek et al. 2017), which explicitly preserve the relative changes in precipitation quantiles, but do not necessarily preserve changes in the standard deviation. This paper proposes the unbiased quantile mapping (UQM) method, a new QM method that preserves the changes in the mean (Δ μ ) and standard deviation (Δ σ ) of future versus historical GCM outputs (i.e. precipitation). The method is compared with other state-of-the-art QM methods using both a synthetic case and real study basins (i.e. Limarí, Maipo, Maule, and Valdivia) located in Chile. Our main conclusions are: • Based on the synthetic experiment, the UQM method appears to be the best in preserving Δ μ and Δ σ , closely followed by detrended quantile mapping (DQM), QDM and SDM, with the standard quantile mapping (QM) being the last. Furthermore, while comparing the performance in preserving Δ μ and Δ σ in the four study basins, UQM is slightly better than the other methods, with a similar performance to DQM and QDM, clearly outperforming the QM method and slightly outperforming SDM. • There is a tradeoff between preserving the relative changes in precipitation quantiles and the changes of the moments (Δ μ and Δ σ ) from the GCM outputs. If preserving the changes in the moments is the objective, UQM appears to be the best method. On the other hand, QDM and SDM are the best in preserving the relative changes in precipitation quantiles. Preserving Δ μ and Δ σ should be a desirable goal because these values are obtained from the mean and standard deviation of all the values of a time span, and thus they are more stable and better estimated than the relative changes in extreme quantiles. • QDM and SDM use the relative changes (or deltas) from the quantiles, which not only are time dependent, but also depend on the specific quantile at a given time. UQM explicitly uses Δ μ and Δ σ , which are only time dependent, making its implementation and use slightly easier. Because DQM does not use a delta but detrends the series before applying the QM to add the trend afterwards, its ease of use is similar to that of UQM. SDM is the most complex and complicated method, as it requires more steps in its implementation. • SDM sometimes presents big issues in preserving Δ σ in precipitation, reaching biases up to 40%, in a couple of cases. • The overall performance of the BA methods and the local coefficient of variation (CV) of the climate variable appear to be related. The smaller the CV value, the better the performance of the BA, regardless of the method being used. This conclusion must be further studied in other locations and also with synthetic experiments, to discard any particularity based on the four study basins. • Some locations that have registered Δ μ in historical records have shown important changes in Δ σ as well. These have happened in the historical record of the four Chilean study basins due to a megadrought since 2010 (Boisier et al. 2016, Garreaud et al. 2017, which has impacted both the historical mean and the standard deviation of precipitation (Barría et al. 2021). Further studies are needed to understand the implications of the projected Δ σ under climate change, and their possible impact over droughts and extreme precipitation. UQM is a simple BA alternative specifically developed to study Δ σ .
Overall, in real applications different BA methods should be tested, and the final choice should consider the specific application and objectives (Tong et al. 2020). UQM is now another available method that could be taken into account, which outperforms other QM methods in preserving raw GCM changes. UQM should be further tested in other locations and climates different to the ones used here. Future studies should also analyse possible adaptations of UQM to daily precipitation, and the performance of UQM when used with temperature. Furthermore, note that weather generation downscaling may be better suited to simulate certain characteristics of climate under changing conditions, especially when the GCM trend is not credibly simulated (Maraun and Widmann 2018).
Meteorológica de Chile (DMC, https://climatologia.meteochile.gob.cl/) for the availability of the data. Finally, we give special thanks to all the reviewers; their insight allowed a significant improvement of the paper.

Code availability
In this work we developed a coding package called climQMBC (climatic variables Quantile Mapping Bias Correction), written in R, Python, and MATLAB, which implements five quantile mapping methods: QM, DQM, QDM, UQM, and SDM, and produces a report that compares the methods. The climQMBC package can be downloaded at: https:// github.com/saedoquililongo/climQMBC.

Disclosure statement
No potential conflict of interest was reported by the authors.