Diffusion Model-based Probabilistic Downscaling for 180-year East Asian Climate Reconstruction

As our planet is entering into the"global boiling"era, understanding regional climate change becomes imperative. Effective downscaling methods that provide localized insights are crucial for this target. Traditional approaches, including computationally-demanding regional dynamical models or statistical downscaling frameworks, are often susceptible to the influence of downscaling uncertainty. Here, we address these limitations by introducing a diffusion probabilistic downscaling model (DPDM) into the meteorological field. This model can efficiently transform data from 1{\deg} to 0.1{\deg} resolution. Compared with deterministic downscaling schemes, it not only has more accurate local details, but also can generate a large number of ensemble members based on probability distribution sampling to evaluate the uncertainty of downscaling. Additionally, we apply the model to generate a 180-year dataset of monthly surface variables in East Asia, offering a more detailed perspective for understanding local scale climate change over the past centuries.


INTRODUCTION
As the average global temperature continues to rise and extreme weather events become more frequent, the impact of climate change on our lives is increasingly apparent 1,2 .However, the rate of warming varies across regions, as does the frequency of such extreme events [3][4][5][6] .Therefore, reliable and accurate regional climate information is crucial to deal with local climate change and its impacts.While many international efforts have been put to develop high-resolution global climate models (GCMs) in recent decades 7,8 , it is important to note that when assessing century-scale climate changes, most GCMs still operate with grid spacings of 100 km or more due to limitations in computational resources [9][10][11] .Such coarse spatial resolutions have proved to be inadequate for evaluating climate changes at the scale relevant to local communities 12 .As a result, this challenge has given rise to the adoption of downscaling techniques, which are now widely employed to bridge the gap between global climate projections and the specific climate information needs of local communities 9,[13][14][15][16] .
These downscaling techniques, including both dynamical and statistical methods, have been developed to generate high-resolution climate data [13][14][15][16][17][18][19][20] .Among these techniques, Regional Climate Models (RCMs) stand out as dynamical models that utilize topography and circulation conditions from GCMs to generate the regional climate information [13][14][15] .However, it is important to note that RCM downscaling results are reliant on the large-scale circulation conditions provided by the GCMs, which leads to the propagation of systematic biases from the GCMs and an increase in uncertainties in the downscaled result 21 .Furthermore, dynamical models incur significant computational expenses and necessitate substantial data storage and processing costs [13][14][15] .In contrast, statistical downscaling methods can provide high-resolution outputs akin to dynamical downscaling but with significantly reduced resource and computational demands.Statistical downscaling methods leverage statistical relationships between low-resolution and high-resolution climate data.These statistical downscaling methods rely on mathematical techniques, including deep learning and traditional statistical approaches, to establish statistical relationships between lowresolution and high-resolution climate data, enabling the derivation of detailed localscale information [17][18][19][20]22 . Curently, both dynamical and statistical downscaling techniques find widespread use in studies related to climate change, climate variability, hydro-climate extremes, and impact assessments at regional scales, particularly within sectors such as agriculture, energy, and water resources [23][24][25][26] .
Unfortunately, both dynamical and statistical downscaling methods currently prioritize deterministic modelling, which often leads to the oversight of the inherent uncertainties within the data and the ill-posed problem of the downscaling process 27,28 .
Those uncertainties have become a growing concern, as precise estimations of climate change and robust assessment methods are crucial for a comprehensive understanding of climate change.Many are now exploring innovative technologies and approaches to tackle the challenges posed by climate change [29][30][31][32] .Here, we introduce the Diffusion Probabilistic Downscaling Model (DPDM), a data-driven approach designed to simulate the probability distribution function of high-resolution climate based on the corresponding low-resolution data.Furthermore, the DPDM incorporates a conditional probability method that accounts for the influence of external factors, such as topography information during the downscaling process.Compared with deterministic downscaling, DPDM derives probability distribution functions and generates a large number of ensemble members 33 , which not only obtains the accurate local details but also allows for more robust estimates of uncertainties.Additionally, we apply DPDM to the NOAA-20th century reanalysis 34 to offer a detailed information on surface climate over East Asia from 1836 to 2015, which significantly enhance our understanding of local-scale climate changes since the late 19th century.

Climate Downscaling via the DPDM
For monthly climate downscaling, one of the main limitations of deep learning methods is the lack of long-term high-resolution observational data.To mitigate this, we use the easily available ERA5-Land dataset 35 , which offers the surface meteorological variable data on a high resolution of 0.1°.Furthermore, for a more robust estimation of uncertainty in downscaling, we use the lower-resolution 1° from ERA5 dataset rather than relying on interpolation 36 .The ERA5 and ERA5-Land are produced by different versions of the Integrated Forecast System model, with ERA5-Land uses an improved land surface model that can better represent processes such as surface-atmosphere interactions, soil moisture dynamics, and vegetation processes compared to the one used in ERA5.In addition, in order to ensure sufficient training data, we assume that the mapping relationship between large-scale areas and smallscale details in downscaling is consistent at different time scales.Consequently, we utilize 6-hourly data spanning the period from 1961 to 2015 for model training (details in Table S1 and Methods).
Recognizing that many of the intricate details in high-resolution data are affected by topography information, we challenge to focus on the East Asia region (65°E-135°E, 15°N-55°N), which boasts some of the most complex and diverse terrain conditions.To comprehensively account for the impact of varying terrains, we introduce terrain and land-sea boundary data into the input dataset of DPDM.Furthermore, we employ a training strategy that divides an input image into a series of patches and randomly selects a patch as the training input.This strategy enhances the diversity of training data and reduce training complexity.Besides this, DPDM train with this strategy demonstrates robust generalization capabilities and can easily be adapted to regions beyond East Asia or to smaller areas through simple fine-tuning.In terms of our model's architecture, we are inspired by the well-established model (SR3) 37 .Additionally, to improve the efficiency of data generation process, we apply the sample process from Denoising Diffusion Implicit Model (DDIM 38 ), leveraging its flexibility to skip some steps in the denoising process and trade off quality for generation speed, as fewer steps inevitably lead to a loss of quality but 250 steps were deemed sufficient for our purposes (ref.39, Fig. S1 and Table S2).

Table 1 unequivocally highlights the superiority of DPDM in downscaling results
compared to other statistical methods.This includes deterministic models such as Enhanced Deep Residual Networks for Super-Resolution with Generative Adversarial framework (EDSR-GAN), known for their superiority over traditional statistical techniques [40][41][42][43] , as well as the widely-used linear interpolation methods (Lerp) in meteorology 44,45 .The consistent results across a variety of meteorological data reveal that the DPDM outperforms other models in terms of key metrics, including the Anomaly Correlation Coefficient (ACC), Peak Signal-to-Noise Ratio (PSNR), Structure Similarity Index Measure (SSIM), and Root Mean Square Error (RMSE).It is worth mentioning that this superiority is particularly evident in the case of precipitation and temperature downscaling, which are very important for social life in a warming climate.Furthermore, a distinct advantage of DPDM is the capacity to generate multiple ensemble members by sampling from the probability distribution.We conduct a comparative analysis of these metrics for varying membership sizes, including ensemble mean with 33 and 100 members.The results clearly demonstrate that increasing the membership size enhances downscaling capabilities proportionally.
It is worth noting that DPDM with multi-member ensembles will increase the inference time compared to other deterministic models (Fig. S2).Nevertheless, when compared to dynamical regional climate models (RCMs), the computational efficiency of DPDM remains acceptable.
To assess their performance in reproducing the spatial distribution, we compare the patterns of climatological mean precipitation and its variance in the test dataset based on the different downscaling methods (Fig. 1).It is evident that the Lerp yields overly smooth results, failing to accurately reproduce local details (Figs.1a4-d4).In contrast, both the deterministic and DPDM impressively capture the spatial details of the precipitation distributions and its temporal variability, as evidenced by improved spatial pattern correlation and reduced spatially averaged Normalized Root Mean Square Error (NRMSE, see Figs. 1 and S4).
However, we note that the deterministic method exhibits lower correlation in midhigh latitudes and plateau areas.This divergence can be attributed to the complex and variable land surface characteristics in such terrains, resulting in variations even under the same large-scale background conditions.In addition, methods based on the deterministic approach (EDSR-GAN) sometimes introduce false features 46,47 (Figs.1c2   and 1d2).Interestingly, DPDM, especially that with multiple members, is able to effectively address these issues.This conclusion is parallel with the findings in distribution histograms, where the Lerp struggles to accurately represent data variation ranges, and where EDSR-GAN tends to exhibit a potential for overestimation (see Fig. 1e and Fig. S2).For the downscaling of other meteorological variables, DPDM also exhibits similar advantages over the deterministic model (Fig. S5 and S6).
To evaluate whether the downscaled results effectively capture high-resolution details, we employ an objective assessment based on the R squared (R 2 ) and blurriness.
Blurriness serves as a crucial metric for objectively assessing whether a model exhibits excessive smoothing or sharpness, avoiding subjective visual assessments 20 .We utilize the absolute difference between the high-resolution data and the Lerp results to evaluate whether the high-resolution data contains sufficient detail information (shown in Fig. 2 c1-c5) because the Lerp results are smooth enough (Fig. S7).If the data predominantly falls on the left side of the diagonal, it suggests that the downscaled data's information content is less than that of true high-resolution data, indicating a bias toward smoothing.
Conversely, if the data leans toward the right side of the diagonal, it signifies excessive sharpness, introducing more information that deviates from reality.Notably, our DPDM with 100 members exhibits a higher degree of model fitting when compared to deterministic models.The DPDM effectively preserves true local details in downscaling results, although some smoothing is visible in precipitation.The evaluation highlights DPDM's capability in reconstructing high-resolution details, even with multiple ensemble members, achieving similar or superior performance compared to other methods (Fig. S8), while also mitigating issues such as excessive smoothing or sharpness observed in other models.

Uncertainties of DPDM
The DDPM can generate multi-member downscaling results through probability distribution sampling.Therefore, we evaluate in detail how the number of ensemble members affects the downscaling performance and why they can be used to evaluate uncertainty (Fig. 3).
As the number of members is increased, the ability of DPDM to capture true details improves in all aspects 48 .An interesting feature is that the improvement in downscaling capability from 30 to 100 members is relatively small.This is reminiscent of the conclusion reached by dynamical models using large members for ensemble predictions.
Part of the reason for this phenomenon is that approximately 30 members are sufficient to represent most of the high-resolution details and to cover the uncertainty space in the downscaling.The ensemble scheme of the dynamical model is similar to the ensemble scheme of DPDM in that both the methods involve sampling multiple results, although the sampling methods are different.Dynamical models utilize randomly perturbed parameterization or randomly perturbed physical processes to generate multiple results, whereas DPDM samples multi-members from the probability distribution.
We find significant varieties among different members, as shown in the spatial pattern of precipitation in July 2020 in Fig. S9.Fig. 3b highlights the considerable uncertainty in the model's downscaling output.But it is worth noting that multi-member average of the model output will better improve the downscaling capability.
Furthermore, we conduct a specific assessment on extreme rainfall events.In July 2020, severe floods occurred in the Yangtze River Basin in East Asia, and drought occurred in southern China (Fig. 3c and 3d).In the two extreme cases, we find that the results obtained through multi-member mean ensemble scheme are always closer to the true values than the Lerp and single-member results, no matter whether dealing with sparse rainfall during droughts or heavy rainfall during floods.And the distribution of multimember results closely approximates a normal distribution.The normal distribution may be better able to measure the uncertainty in atmospheric processes that determine the local detail conditions.Similarly, downscaling of other atmospheric variables also show similar results (Fig. S10).Additionally, we have evaluated the spread of DPDM with 100 members (Fig. S11 and Fig. S12).Unlike many other downscaling models that tend to exhibit underestimated spread, DPDM demonstrates the capability to generate sufficient ensemble spread, even overestimating it for each variable.Furthermore, the overestimation of uncertainty suggests that DPDM can accommodate additional sources of uncertainty beyond solely the downscaling process itself.

Application for downscaling 180-year surface climate dataset
With the success of DPDM, we now apply the method to reconstruct highresolution historical climate data in East Asia for the past century.Additionally, we explore several potential application scenarios for this high-resolution dataset.These application scenarios not only demonstrate the model's capability for small-scale reconstruction but also address several climate change-related issues of broad academic interest.
We select the NOAA-20C dataset as the low-resolution data because it is the only global dataset covering nearly two centuries from 1836 to 2015, encompassing the entire industrialization period.It also provides different temporal resolutions and circulation data, which not only provides the basis for constructing high-resolution surface data at six-hourly intervals but also uses the reconstructed high-resolution data with circulation to explore more on mechanisms, especially the attribution of extreme events.
To ensure the reliability of reconstruction data, we perform evaluations against the widely utilized CRU-station dataset.Given inherent data credibility issues with reanalysis datasets and the fact that observation stations solely provide precipitation and temperature data, we conduct a relative error analysis between DPDM and the Lerp results for these two datasets.It may be a crude way to evaluate our results but we do not have any other choice in the absence of high-resolution observational data for such a long period from 1836 to 2015.Fig. 4a and 4b reveal a noticeable reduction in relative errors on most stations when employing DPDM.It shows that the reconstructed data are closer to the station-observed data.
To assess the changes in aridity over the past centuries by use of the highresolution reconstructed dataset, we have computed the Aridity Index (see Methods), a commonly applied metric with a threshold of 0.65.As shown in Fig. 4c High-resolution data can provide important details for assessing extreme events.Therefore, we also conduct a simple evaluation of extreme hot and dry compound events (Figs 4e and 4f, Methods).In North China, the high-resolution reconstructed data clearly provides greater detail, detecting more extreme events and offering a finer feature of areas that are prone to such events.Additionally, it augments the available samples for subsequent attribution and synthesis analyses.Note that, while wind speed data lacks observational records, the climate statistics indicate that DPDM results also yield more detailed characters (Fig. S14).In previous studies 49 , the trend in wind power change under the global warming is often estimated with the low-resolution data.
However, the high-resolution data notably reveals several regions undergoing faster changes (Figs.4g-i).
In summary, the high-resolution data generated by the DPDM not only exhibits a certain level of credibility but also enhances our understanding of climate details.This high-resolution dataset, covering the past centuries, may provide important details for improved understanding of the historical climate change.

DISCUSSION
In this study, we have introduced a novel probabilistic downscaling model, DPDM, for climate downscaling.We evaluate the downscaling capabilities of DPDM, which not only accurately simulates the probability distribution function of high-resolution data, but also generates a large number of members to quantify the uncertainty of the downscaled information.The latter is important since small-scale conditions under a large-scale background are never deterministic.In addition, the downscaling framework of DPDM has great potential in medium-term weather forecasts, climate predictions and future scenario projections.For instance, to generate an adequate number of ensemble members, it can be used to emulate traditional methods like the single model initial-condition large ensemble for identifying and robustly sampling extreme events 50,51 .
It is undeniable that DPDM still has more room for improvement.Introducing additional circulation conditions and external forcing could enhance the model ability with more physical constraints 52 .This new approach holds promise for applications in bias correction and downsizing of dynamical model predictions.As for high-resolution climate datasets over the past centuries, while it offers valuable insights and applications, the existing datasets are limited in terms of the number of available variables and their temporal resolution.With sufficient computing resources, there is potential for increasing the temporal resolution to 6-hour intervals and downscaling other surface or upper-atmosphere variables to enable more comprehensive analyses.
Additionally, we have noticed that NOAA-20CR provides corresponding ensemble spread data, which can effectively quantify the uncertainty in the initial conditions.In the future, we could sample from this ensemble spread to generate ensemble members that incorporate data uncertainty, and then apply the DPDM separately for each member.This approach would produce a larger ensemble distribution, incorporating not only the uncertainties from the downscaling process but also the uncertainties from the initial conditions.These may be reserved for future study.
In addition, the probabilistic essence and robust mathematical foundation of the diffusion model may open up a wealth of new possibilities for its practical applications in climate science as a promising tool.Its applicability extends far beyond downscaling; it holds potential for forecasting, assimilation, data reconstruction, model bias correction, sensitivity experiments, scientific inquiry, and even causal analysis.We believe that the time has come to explore the applications of the diffusion model, tackling intriguing scientific questions and contributing to the advancement of climate science.     .Here, we give a brief introduction to DPDM in forward and reverse process (more details in Ho et al 33,38,39,54 ).

Forward Diffusion process
The forward process of DPDM is similar to iteratively constructing a mapping relationship from a high-resolution distribution to pure gaussian distribution using the Markov chain with total length T.More specifically, given a data sampled from the real high resolution data distribution  0 ~(), the forward diffusion process, in which we gradually add Gaussian noise to  0 in total steps (T) with the noisy density for each step (t) being controlled by a variance schedule  1 …   , produces the sequence of highresolution data with noisy samples  1 …   .

𝑞(𝑥 𝑡 |𝑥 𝑡
If the magnitude   of the noise added at each step is small enough, and the total step T is large enough (in our experiments, T is set to 1000 steps), then   is equivalent to an isotropic Gaussian distribution   ∼ N (0, I).With the help of the properties of Markov chain, we successfully connect the high-resolution data distribution () with the Gaussian distribution N (0, I).Then, we can obtain the noise samples at step t in the forward process, using the following equation: In order to obtain   without iteration for fast training, we can further expand   with the help of reparameterization trick and additive nature of Gaussian distributions.
Let   = 1 −   and,  ̅  = ∏    =1 , we can get the   and (  | 0 ): Through Bayes' rule, we can get the following equation ( 6) and (7).Subsequently, based on the Markov assumption, equation ( 8) can be obtained from (7).Finally, by utilizing Gaussian distribution functions from the forward process, all components of equation ( 8) can be expressed, simplifying it to equation ( 9): Following the standard Gaussian distribution function, the mean and variance can be parameterized as follows: 1−  ~   (10)   According to equation ( 9) and (10), we can obtain the distribution function (11), and we can sample it to get the  −1 thereby realizing the reverse denoising process.With the ablation study (Fig. S15), we investigated the effects of adding terrain information to the patch input and explored the impact of different patch sizes (32,64,128,256) on model performance.Our experiments revealed that incorporating terrain information led to faster loss convergence, patch size of 128 can obtained optimal performance and patch size of 256 resulted in non-converging loss.
Through the patch strategy, the model can dynamically process low-resolution data to generate high-resolution data for any region and condition, instead of just fitting in a certain region, thereby improving the generalization ability of the model.For example, it can obtain high-resolution results outside the target region without the training, albeit with limited effectiveness, as seen in the Gulf of Mexico (Fig. S16) and northern Australia (Fig. S17).It should be noted that in the DPDM, we also need to use the above patch strategy for inference, rather using the entire low-resolution data.The input image is divided into different patches, each patch is processed individually, and then the results are combined to obtain the final output.To ensure continuity at the boundaries, we introduce overlapping between some patches.However, there still be some quality degradation at the boundaries of patches.Increasing the amount of overlap or employing advanced techniques for boundary restoration can yield more consistent results, but this approach requires sufficient computational resources.

Climate index and Compound events
An aridity index (AI 56 ) is a numerical indicator of the degree of dryness of the climate at a given location.
where PET is the potential evapotranspiration, which is calculate by the python package of Climate Indices (https://climate-indices.readthedocs.io/en/latest/)and the P is the annual average precipitation.
Compound hot-dry events are defined as the co-occurrence of high mean temperature anomaly (above the 90th percentile) and low mean precipitation anomaly (below the 10th percentile) values over the warm season 51 (i.e., the three consecutive months with highest mean temperature during 1836-2015) for each grid point, The centennial trends of the two variables are removed.
Wind energy 57 is a typical measure of wind energy potential, defined as follows: where ρ represents the air density, which is assumed to be a constant value of 1.213  *  −3 at standard atmospheric conditions, and  ℎ is approximately expressed by the wind speed at a height of 10 m.

Metric
Based on the downscaling results, we compute some metrics, i.e., Anomaly , the shaded area represents the temporal average of AI from 2005 to 2015.The blue line, comprised of both solid and dashed segments, depicts the changes in the 11-year running mean AI at a low resolution with a constant value of 0.65.In contrast, the green line illustrates the AI values of 0.65 from 2005 to 2015 at high resolution.The findings suggest that the low-resolution data underestimates the expansion of aridity in the mid-to high-latitudes of East Asia and Northern China in the context of global warming.Fig. 4d further quantifies the area proportion of drought regions (65°E-135°E, 33°N-55°N) on a decade-by-decade basis.Remarkably, even when accounting for the uncertainty among different ensemble members, the low-resolution data persistently underestimates drought areas by approximately 3%.Regarding the warming and humidification trend in the northwestern China during recent decades, we examine the precipitation changes in northwestern China (75°E-110°E, 30°N-50°N).The analysis reveals an actual increasing trend in precipitation since 1970, but the low-resolution data significantly overestimates precipitation intensity at a rate of 0.79 mm/day per 10 years.In contrast, the high-resolution data offers finer estimations and provides uncertainty estimates for assessing credibility, indicating a trend of 0.69 mm/day per 10 years within an uncertainty range of 0.65-0.78mm/day per 10 years.In addition, the precipitation distribution histogram in rainy season (MJJAS) found that our high-resolution data can correctly understand the differences in climate characteristics in different regions.It tends to obtained less precipitation in drought regions while more precipitation in humid regions than low resolution (Fig. S13).

Figure 1 .
Figure 1.The downscaling results for climatological mean total precipitation and

Figure 2 .
Figure 2. Evaluation of R squared (R 2 ) and blurriness for downscaling results in

Figure 3 .
Figure 3.The role of different members in the DPDM.(a) Impact of ensemble

Figure 4 .
Figure 4. Applicable scenarios for high-resolution datasets over the past 180 years

11 )
(  , ,  ̃)) , √  )    (  , ,  ̃)) + √   ~ (0, ) (Training details and model structure To obtain a sufficient volume of training data, we employ a dataset with a 6-hour time resolution instead of the monthly dataset.We established a connection between the ERA5 product (1° spatial resolution), and the high-resolution ERA-land dataset (0.1° spatial resolution).These reanalyses are provided by the Copernicus Climate Change Service at ECMWF, combining a large range of satellite-based and land-based observations with high-resolution model simulations through state-of-the-art data assimilation techniques, spanning the period from 1961 to 2021.In this study, we split the ERA-Land and ERA5 datasets into the training period of 1960-2015 and the test period of 2015-2021.In order to consider the effect of terrain, we use the topography, land-sea mask and low-resolution data information as the input of DPDM for training.We randomly crop 128*128 patch from the original lowresolution data and get the corresponding high-resolution data for training, instead of using the original low-and high-resolution data for training.The reason why we use this strategy is because, this can not only increase the diversity of data to learn the mapping of low-resolution data to high-resolution data, but also reduce the consumption of computing resources, such as graphics memory and computational cost.
For different surface variable, we adopt different normalization strategies to facilitate rapid loss convergence and training stability.Precipitation data undergoes a transformation by adding one and then taking the logarithm.Other variables are standardized and normalized to a range between 0 and 1.We also standardize the topography data and change the land-sea mask information to a matrix containing only 0 and 1.The DPDM model architecture is similar to SR3, which is a U-Net-like architecture (Fig.S1).At each time step t, the model's input comprises the concatenation of two datasets: the conditional data  ̃ and the noisy data   , both have the same dimensions as the high-resolution data  0 .Conditional data  ̃ includes interpolation result obtained by Lerp interpolating low-resolution data, topography, and land-sea mask information.Concatenation is a simple and effective method to add the conditional data to the model.Then, a convolutional layer with 64 kennels is used to extract the input data information.Downsampling modules are applied in DPDM consisting of several residual blocks and self-attention layers55 .Based on empirical experience, we only use 3 downsampling modules, and the data dimensions are reduced from 128 to 16.The upsampling modules is similar to the downsampling modules.All the convolutions in our model use the 3*3 convolution kernel size and the 1*1 stride.For encoding the timestep t, we use the sinusoidal positional encoder55 that contains two fully-connected layers and a sigmoid linear unit (SiLU) activation function between the two layers.Then we add the timestep feature encoded by above sinusoidal positional encoder to our intermediate feature maps after the group normalization operators in each residual block.

Table 1 .
Evaluation of the downscaling performance of five surface variables from 2016 to 2021using Root Mean Square Error (RMSE), Now we can conveniently sample   at an arbitrary timestep t in the forward process.At the timestep T, we can sample gaussion noise map   from the gaussion distribution  (0, ).Then we need to estimate ( −1 |  ) to remove the gaussion noise added in the forward process of DPDM.Unfortunately, we cannot easily estimate ( −1 |  ) because it requires to use the entire dataset.Therefore, we learn the DPDM model   conditioned on the state   at the timestep t and the conditional data  ̃ to approximate these conditional probabilities for getting the next state  −1 . ( −1 |  ,  ̃) = ( −1 |  ,  0 ) =  ( −1 ;   (  , ,  ̃), Σ  (  , )) Here,  and  denote the observed and the downscaling results, respectively.  ̅̅̅̅ and   ̅̅̅̅ denote the climatology of observed and the downscaling results in each calendar month m (from 1 to 12).The label  denotes the forecast target year.Finally,  and  denote the earliest (that is, 2016) and the latest year (that is, 2021) of the validation, respectively. denote the maximum value of the normalized data (that is, 1). , and  , represents the spatial means and standard deviations of the downscaling results, while  , and  , represents the spatial means and standard deviations of observation. 1 and  2 are constants to avoid computation instability when the denominator approaches zero.