A Markov Chain-Based Bias Correction Method for Simulating the Temporal Sequence of Daily Precipitation

: Bias correction methods are routinely used to correct climate model outputs for hydrological and agricultural impact studies. Even though superior bias correction methods can correct the distribution of daily precipitation amounts, as well as the wet-day frequency, they usually fail to correct the temporal sequence or structure of precipitation occurrence. To solve this problem, we presented a hybrid bias correction method for simulating the temporal sequence of daily precipitation occurrence. We did this by combining a ﬁrst-order two-state Markov chain with a quantile-mapping (QM) based bias correction method. Speciﬁcally, a QM-based method was used to correct the distributional attributes of daily precipitation amounts and the wet-day frequency simulated by climate models. Then, the sequence of precipitation occurrence was simulated using the ﬁrst-order two-state Markov chain with its parameters adjusted based on linear relationships between QM-corrected mean monthly precipitation and the transition probabilities of precipitation occurrence. The proposed Markov chain-based bias correction (MCBC) method was compared with the QM-based method with respect to reproducing the temporal structure of precipitation occurrence over 10 meteorological stations across China. The results showed that the QM-based method was unable to correct the temporal sequence, with the cumulative frequency of wet- and dry-spell length being considerably underestimated for most stations. The MCBC method can could reproduce the temporal sequence of precipitation occurrence, with the generated cumulative frequency of wet-and dry-spell lengths ﬁtting that of the observation well. The proposed method also performed reasonably well with respect to reproducing the mean, standard deviation, and the longest length of observed wet- and dry-spells. Overall, the MCBC method can simulate the temporal sequence of precipitation occurrence, along with correcting the distributional attributes of precipitation amounts. This method can be used with crop and hydrological models in climate change impact studies at the ﬁeld and small watershed scales. writing—original draft preparation, H.L.; writing—review and editing, C.-Y.X.; funding acquisition, J.C. All authors


Introduction
As an essential part of the third working group report from the Intergovernmental Panel on Climate Change's (IPCC, 2018) sixth assessment report, food security under climate change conditions model in reproducing the precipitation series at five stations in Oklahoma, where they found that the sequence of the observed precipitation was well reproduced by this method for all stations. The Markov chain-based model is usually incorporated into weather generators to simulate the temporal sequence of precipitation occurrence. When using the weather generator as a downscaling tool, its parameters need to be adjusted based on climate change information between future and historical periods. However, the adjustment of weather generator parameters is not an easy task, especially for downscaling precipitation occurrence.
This study proposes a hybrid bias correction method by combining a commonly used bias correction method (i.e., the daily bias correction method, DBC) with the first-order two-state Markov chain. The new method takes advantage of the traditional bias correction method in correcting the distributional attributes of precipitation amounts and wet-day frequency, and it takes advantage of the Markov chain in simulating the sequence of precipitation occurrence. The performance of the proposed Markov-chain based bias correction method (MCBC) is evaluated using comparisons to the DBC over 10 stations across China.

Study Area
The performance of the proposed bias correction was demonstrated using ten meteorological stations distributed across different climatic zones in China ( Figure 1). The mean annual precipitation at these stations ranged from 185 mm to 1565 mm, with the wet day frequency varying from 15% to 44%. The detailed information, including longitude, latitude, and elevation for the ten stations, can be found in Table 1.
Atmosphere 2018, 9, x FOR PEER REVIEW 3 of 17 chain-based model in reproducing the precipitation series at five stations in Oklahoma, where they found that the sequence of the observed precipitation was well reproduced by this method for all stations. The Markov chain-based model is usually incorporated into weather generators to simulate the temporal sequence of precipitation occurrence. When using the weather generator as a downscaling tool, its parameters need to be adjusted based on climate change information between future and historical periods. However, the adjustment of weather generator parameters is not an easy task, especially for downscaling precipitation occurrence. This study proposes a hybrid bias correction method by combining a commonly used bias correction method (i.e., the daily bias correction method, DBC) with the first-order two-state Markov chain. The new method takes advantage of the traditional bias correction method in correcting the distributional attributes of precipitation amounts and wet-day frequency, and it takes advantage of the Markov chain in simulating the sequence of precipitation occurrence. The performance of the proposed Markov-chain based bias correction method (MCBC) is evaluated using comparisons to the DBC over 10 stations across China.

Study Area
The performance of the proposed bias correction was demonstrated using ten meteorological stations distributed across different climatic zones in China ( Figure 1). The mean annual precipitation at these stations ranged from 185 mm to 1565 mm, with the wet day frequency varying from 15% to 44%. The detailed information, including longitude, latitude, and elevation for the ten stations, can be found in Table 1.

Data Sources
This study used both the observed and GCM-simulated precipitation. The observed precipitation series of 10 meteorological stations were obtained from the China Meteorological Data Sharing Service System (http://data.cma.cn/data/index/0b9164954813c573.html), covering the 1961-2014 period. Precipitation simulated by the 10 GCMs was extracted from the Coupled Model Inter-comparison Project Phase 5 (CMIP5) database. Table 2 shows the detailed information of these GCMs. The grid box that overlaps the observed gauge was selected to represent the simulated precipitation for this station. All simulated precipitation time series cover the 1961-2005 period.

Methodology
The proposed bias correction method is a two-stage hybrid method combining the quantile-mapping (QM) method with the first-order two-state Markov chain. The first stage involves correcting the distribution of precipitation amounts on wet days and the wet-day frequency of the projected GCM-simulated precipitation using the QM-based method by Chen et al. [5]. This bias Atmosphere 2020, 11, 109 5 of 17 correction method is conceptually simple and widely used [10,21,29,43,44]. The second stage involves generating the temporal sequence of precipitation occurrence using the first-order two-state Markov chain based on the first stage corrected monthly precipitation amounts. More details of this method are described as follows.

Quantile Mapping Bias Correction Method
The QM method used in this study is the DBC method [10], which is a combination of the daily translation (DT) [45] and local intensity scaling (LOCI) methods [27]. The DBC method inherits the merit of the LOCI method in correcting the wet-day frequency of precipitation occurrence and that of the DT method in correcting the frequency distribution of precipitation amounts. Specifically, for each calendar month, the LOCI method first determines a threshold for the GCM-simulated precipitation time series so that the wet-day frequency of simulated precipitation series for the calibration period is equal to that of the observed precipitation. Then, different scaling factors, which are the ratios in percentiles between the observed and climate model-simulated precipitation in the calibration period, are applied to the corresponding percentiles in the validation period. One hundred percentiles are used to represent the frequency distribution of the precipitation amounts. All bias correction methods assume that the biases of climate model simulations are stationary over time, which means that the biases of climate model simulations at the calibration period do not change at the validation period. Even though the distribution of precipitation at the validation period may differ from that at the calibration period, the same biases can be removed. The formulation of DBC for correcting the distribution of precipitation is presented as follow: P cor,vali,m,p = P obs,m,p P cali,m,p * P vali,m,p , The subscripts obs, cali, and vali denote the observation, calibration period, and validation periods; and m and P denote the months and percentile, respectively, while cor denotes corrected data.

The First-Order Two-State Markov Chain
The first-order two-state Markov chain is widely used to simulate the sequence of precipitation, and its consistent performance has been demonstrated in many studies [6,34,46,47]. When using this method, the probability of precipitation on a given day is determined based on the precipitation status (dry or wet) of a previous day. The transition probability P 11 is the probability of a wet day following a wet day, and P 01 is the probability of a wet day following a dry day, as shown in the following formulas: Given the status of a previous day, for example, a dry day, if a random number derived from a standard uniform distribution is less than P 01 , a wet day is generated; and otherwise, a dry day is generated. A similar process can also be used with P 11 .

Hybrid Bias Correction Method
When using the first-order, two-state Markov chain to simulate the temporal sequence of precipitation occurrence, P 01 and P 11 must be adjusted according to the climate change signal simulated by the climate models. Previous studies showed that there are strong linear relationships between the mean monthly precipitation and the two conditional transition probabilities [16,19,30]. The linear relationship can be fitted using the observed precipitation time series. In other words, if we know the mean monthly precipitation of the validation period, two conditional transition probabilities for this period can be obtained based on the linear relationship established in the calibration period. As mentioned earlier, the QM method can provide an accurate simulation of the monthly precipitation amounts.
A linear relationship between the transition probabilities of precipitation occurrence and the mean monthly precipitation is first established using the four-point regression method of Zhang et al. [16]. The method is based on the observed precipitation time series, as in the following three steps.
(a) Taking January as an example, the total precipitation of this month over the whole period is sorted in ascending order and then separated into dry and wet groups. Then the mean monthly precipitation is calculated from each group. The daily series corresponded to the months in each group are gathered and used to calculate P 01 and P 11 . With this step, a set of P01, P11, and the mean monthly precipitation is calculated for each group. (b) The total precipitation of this month each year is also divided into two, even periods in chronological order. A set of parameters is also calculated for each period. (c) A linear equation is established for P 01 (and P 11 ) and the mean monthly precipitation based on the four samples calculated from step-(a) and step-(b).
The use of two extreme cases (the sorted dry and wet conditions) of precipitation occurrence versus the amounts is done to ensure that the fitted relationship can hold for the changing climate. In other words, including the two extreme cases in establishing the linear relationship implicitly incorporates non-stationary climate conditions when estimating the transition probabilities of precipitation occurrence for any climatic conditions within the entire range. With the bias-corrected mean monthly precipitation for the validation period, the transition probabilities P 01 and P 11 can be regressed using the established linear relationships for each calendar month. The first-order, two-state Markov chain is used to generate the sequence of precipitation occurrence. For a predicted wet day, the precipitation amount is generated by realigning the DBC-corrected precipitation amounts. Specifically, a set of random numbers is first generated for wet days simulated by the first-order, two-state Markov chain. The random numbers are then replaced by DBC-corrected precipitation amounts, according to the numerical value of the random number. In other words, the larger precipitation amount is given to the larger random number, and vice versa. If the first-order, two-state Markov chain is perfect for simulating the sequence of precipitation occurrence, then the simulated wet-and dry-day frequencies should be equal to the frequencies generated by the DBC. However, this may not be the case. Scaling the realigned daily precipitation for each month may be necessary to correct or offset the generation error of the first-order, two-state Markov chain. This is done by multiplying the ratio of the DBC-corrected mean monthly precipitation and MCBC-corrected mean monthly precipitation by the daily precipitation series of the MCBC-sampled for each specific month.

Analysis Procedure And Statistics
The performance of the proposed hybrid method was validated based on the observed and GCM-simulated precipitation time series from 1961 to 2005. Before evaluating the MCBC method, its baseline model is demonstrated using 10 stations in China, i.e., the performance of the DBC method and the linear relationships between the mean monthly precipitation and conditional probabilities of precipitation occurrence. A split-sample cross-validation method was conducted to evaluate the performance of the DBC method [8]. Specifically, the time series of the observation and GCM simulations for the 1961-2005 period is divided into five equal sub-periods (1961-1969, 1970-1978, 1979-1987, 1988-1996, and 1997-2005). Four sub-periods are used to calibrate the DBC method and the remaining sub-period is used for validation. The above processes are conducted five times to get five validation sub-periods. The five validation sub-periods are then combined to get a 45-year validation period. The data length used in the calibration of the linear relationships between the mean monthly precipitation and the conditional probabilities of precipitation occurrence was more than 90 years in prior studies by Zhang et al. [16,19]. Owing to the limited data length, we conducted this process with the total observations from 1961-2014.
The Markov chain-based bias correction (MCBC) method was applied to correct the daily precipitation series simulated by the climate models for each station. The nonparametric two-sample Kolmogorov-Smirnov (K-S) test was used to test whether the bias-corrected and observed precipitations were from the same distribution. The same method has been used in many other studies [16,19,30,35]. The statistical test D was calculated from the length of the two tested samples and the null hypothesis was that the two samples were from the same distribution. When the maximum deviation between the cumulative frequency distribution of two samples was larger than D at the 0.05 significance level, the null hypothesis was rejected, and vice versa.

Validation of the Daily Bias Correction
The performance of the DBC method for bias correcting GCM precipitation has been verified in many studies [21,29]. Since the premise of the MCBC method is to obtain the well-corrected mean monthly precipitation used to interpolate the conditional transition probabilities, only the performance of the DBC method in correcting the mean monthly precipitation is demonstrated in Figure 2.

Validation of the Linear Relationships
The data points of P01, P11, and mean monthly precipitation, which were calculated from the dry group, wet group, and two even groups of the calibration period, are plotted in Figures 3 and 4 for January and July, respectively for the ten stations. The linear curves that were regressed using the The raw GCMs considerably overestimated the mean monthly precipitation, with the medians of the relative errors (REs) ranging between 50% and 300% for most stations and months. In terms of the RE, the simulated mean monthly precipitations in July, August, and September are slightly better than in other months for most stations. The result was because most stations had higher precipitation in these three months than in other months. The biases of the mean monthly precipitation simulated by the GCMs were largely reduced by the DBC method. For example, the medians of REs across all GCMs ranged between 0.01% and 14.7%, and these errors were due to cross-validation. Although the performance of the DBC method partly depends on a specific month and location, it generally corrects the mean monthly precipitation for all GCMs accurately.

Validation of the Linear Relationships
The data points of P 01 , P 11, and mean monthly precipitation, which were calculated from the dry group, wet group, and two even groups of the calibration period, are plotted in Figures 3 and 4 for January and July, respectively for the ten stations. The linear curves that were regressed using the four data points were used to interpolate the conditional transition probabilities in the validation period. For comparison, P 01 , P 11 , and the mean monthly precipitation were also calculated for the validation period and plotted in Figures 3 and 4 as solid points. The names of the legend 'cali' and 'vali' represent calibration and validation, respectively.
Atmosphere 2018, 9, x FOR PEER REVIEW 9 of 17 stations) for P01, and for 75 out of 120 cases for P11. The reason why there was a stronger relationship for P01 than P11 may be because the occurrence of wet-following-wet events is lesser than that of wetfollowing-dry events, which contributed to greater variability and less reliability when estimating P11. As Zhang et al. (2012) presented in their study that a longer period should be used to obtain large samples of wet-following-wet events to improve the accuracy of the estimated P11.  Relationships between the conditional precipitation occurrence probabilities of wet-following-wet (P 11 ) and wet-following-dry (P 01 ) and mean monthly precipitation amounts in January for all stations.

Performance of the MCBC method
The performance of the MCBC method was first evaluated with respect to reproducing the cumulative frequencies of the observed wet-and dry-spell lengths. Figure 5 and Figure 6 present the cumulative frequencies of the wet-and dry-spell lengths generated by the DBC and MCBC methods (for one series), respectively. The cumulative frequencies of the observed wet-and dry-spell lengths were also used for comparison. Since similar results were obtained for all stations, only four stations from different climate zones were presented for illustrations. Generally, the DBC method underestimated the cumulative frequency of the wet-and dry-spells, except for the wet-spell of the first station ( Figure 5). A significant uncertainty related to the choice of a GCM simulation was observed. However, the MCBC method largely reduced the biases of the cumulative frequency for the wet and dry spells. The uncertainty related to the GCMs was considerably corrected.
The K-S test showed that about 70% of the stations and GCM combinations (n = 10 stations × 10 GCMs) could not reject the null hypothesis for the cumulative frequencies of the wet-and dry-spells at the P = 0.05 significance level when using the proposed MCBC method. In other words, about 70% of cases were observed, and the bias-corrected cumulative frequencies may be from the same distribution. However, this value was only 20% when using the DBC method. Relationships between the conditional precipitation occurrence probabilities of wet-following-wet (P 11 ) and wet-following-dry (P 01 ) and mean monthly precipitation amounts in July for all stations.
Overall, the regression lines fitted the data points very well for all months and stations, indicating a strong linear relationship between the mean monthly precipitation and conditional transition probabilities. In particular, the data points of P 01 from both the calibration and validation periods were fairly close to the regression lines. The correlation coefficient between mean monthly precipitation and transition probabilities was greater than 0.9 for 91 out of 120 cases (12 months × 10 stations) for P 01 , and for 75 out of 120 cases for P 11 . The reason why there was a stronger relationship for P 01 than P 11 may be because the occurrence of wet-following-wet events is lesser than that of wet-following-dry events, which contributed to greater variability and less reliability when estimating P 11 . As Zhang et al. (2012) presented in their study that a longer period should be used to obtain large samples of wet-following-wet events to improve the accuracy of the estimated P 11 .

Performance of the MCBC Method
The performance of the MCBC method was first evaluated with respect to reproducing the cumulative frequencies of the observed wet-and dry-spell lengths. Figures 5 and 6 present the cumulative frequencies of the wet-and dry-spell lengths generated by the DBC and MCBC methods (for one series), respectively. The cumulative frequencies of the observed wet-and dry-spell lengths were also used for comparison. Since similar results were obtained for all stations, only four stations from different climate zones were presented for illustrations. Generally, the DBC method underestimated the cumulative frequency of the wet-and dry-spells, except for the wet-spell of the first station ( Figure 5). A significant uncertainty related to the choice of a GCM simulation was observed. However, the MCBC method largely reduced the biases of the cumulative frequency for the wet and dry spells. The uncertainty related to the GCMs was considerably corrected.
The K-S test showed that about 70% of the stations and GCM combinations (n = 10 stations × 10 GCMs) could not reject the null hypothesis for the cumulative frequencies of the wet-and dry-spells at the p = 0.05 significance level when using the proposed MCBC method. In other words, about 70% of cases were observed, and the bias-corrected cumulative frequencies may be from the same distribution. However, this value was only 20% when using the DBC method.  Table 2. The X-axis is on a logarithmic scale.  Table 2. The X-axis is on a logarithmic scale.
Consecutive wet days are an essential factor for high-impact flooding [33] and consecutive dry days for an extended period are one of the causes of droughts [48]. Accurately predicting the long wetand dry-spell is vital for dealing with extreme hydrological events. Figure 7 shows the performances of the DBC (red bars) and MCBC (blue bars) methods in simulating the frequency of the spell lengths with more than or equal to 3, 5, and 7 consecutive dry days and wet days. The ratio of bias-corrected frequency to observed frequency was used as a metric to evaluate the performance of two bias correction methods. The closer the ratio was to 1, the more accurate the simulation. The results showed that the DBC method overestimated the frequencies of dry and wet spells for most stations. The result was expected, as the temporal structure was not specifically corrected by the DBC method. The MCBC method simulated the frequencies of dry and wet spells much better than the DBC method, indicated by the fact that the blue bars were shorter and closer to the baseline. Specifically, for the DBC method, the mean ratios across all stations with consecutive wet days being more than or equal to 3 days, 5 days, and 7 days ranged between 0.875 and 1.438, 0.938 and 2.471, 1.086 and 8.796, respectively. In contrast, the mean ratios for the MCBC method ranged between 0.909 and 1.074, 0.853 and 1.354, 0.772 and 1.434. Similar results were also observed for consecutive dry days. Figure 8 shows the RE of the mean (Figure 8a,b), standard deviation (Figure 8c,d), and the longest length (Figure 8e,f) of wet-and dry-spells generated by the DBC and MCBC methods. The mean spell length simulated by the DBC method was considerably biased, with the median of REs ranging between −0.45% and 21.7% for the wet spell and between 8.62% and 27.44% for the dry spell. The MCBC method largely reduced the biases of the simulated mean spell length, with the median of REs ranging from −4.36% to −1.24% for the wet spell and ranging from −4.16% to −0.64% for the dry spell. More importantly, the uncertainty related to climate models was largely reduced. Similar results were also observed for the standard deviation of spell lengths (Figure 8c,d). The MCBC method also performed slightly better than the DBC method with respect to reproducing the longest wet and dry spell lengths, with the median of RE of 10-simulation medians ranging from −33.3% to 11.1% for the wet spell and ranging from −50.8% to −2.2% for the dry spell. For the DBC method, the median of RE ranged from −23.3% to 61.1% for the wet spell and from −31.0% to 27.2% for the dry spell. Although the benefits of the MCBC method are limited, the uncertainty related to climate model is reduced, with the interquartile ranges varying from 3.3% to 17.6% for the wet spell and ranging from 3.2% to 13.6% for the dry spell. In contrast, the interquartile ranges for the RE of the longest spell length simulated by the DBC method ranged from 26.6% to 77.8% for the wet spell and ranged from 20% to 36.6% for the dry spell.  Table 2. The X-axis is on a logarithmic scale.  Table 2. The X-axis is on a logarithmic scale. DBC method, indicated by the fact that the blue bars were shorter and closer to the baseline. Specifically, for the DBC method, the mean ratios across all stations with consecutive wet days being more than or equal to 3 days, 5 days, and 7 days ranged between 0.875 and 1.438, 0.938 and 2.471, 1.086 and 8.796, respectively. In contrast, the mean ratios for the MCBC method ranged between 0.909 and 1.074, 0.853 and 1.354, 0.772 and 1.434. Similar results were also observed for consecutive dry days.  Table 2.
median of RE ranged from −23.3% to 61.1% for the wet spell and from −31.0% to 27.2% for the dry spell. Although the benefits of the MCBC method are limited, the uncertainty related to climate model is reduced, with the interquartile ranges varying from 3.3% to 17.6% for the wet spell and ranging from 3.2% to 13.6% for the dry spell. In contrast, the interquartile ranges for the RE of the longest spell length simulated by the DBC method ranged from 26.6% to 77.8% for the wet spell and ranged from 20% to 36.6% for the dry spell.

Discussions and Conclusions
Accurately predicting the temporal sequence of precipitation occurrence at field or site-specific scales is critical in assessing climate change impacts on food security [42]. Daily soil water dynamics that are influenced mainly by daily precipitation sequences significantly affect crop growth and production, especially at the critical reproductive stages. More importantly, prolonged drought or wet spells can be detrimental to crop production. Thus, adequate simulation of wet-and dry-spell lengths is essential for adequately simulating the impacts of climate change on crop production and food security. Although bias correction methods, especially QM-based methods, can correct the precipitation amount simulated by climate models, the temporal sequence of precipitation occurrence is usually not explicitly considered. This study proposed a hybrid bias correction method to address the above problem. The MCBC method can be considered as a hybrid method combining QM bias correction and the stochastic Markov chain-based model. It not only has the advantage of the commonly used bias correction method in correcting the distributional attributes of precipitation amounts but it also takes advantage of the Markov chain in simulating the temporal sequence of precipitation occurrence. An evaluation of the proposed method showed that the temporal sequence of the precipitation occurrence was well corrected, and the proposed method consistently performed better than the traditional QM bias correction method. Specifically, the distribution of MCBC corrected wet-and dry-spell lengths fitted that of the observations very well, with the K-S test showing that more than 70% of the corrected series were likely from the same distribution of the observation at the p = 0.05 significance level. The MCBC method also reasonably reproduces the observed statistics of spell length, such as the frequency of spell length, mean, standard deviation, and the longest extreme spell length.
In addition, the Markov chain-based method used in this study is conceptually simple, easy to use, and requires less computation, because its parameters can be adjusted based on QM corrected results with a linear relationship of observed data. Even though the DBC model is applied as a baseline model in this study, the Markov chain-based model can be used with any other bias correction method. Moreover, this study uses the linear relationship between the observed mean monthly precipitation and the transition probabilities to adjust Markov chain parameters. Meanwhile, other methods can also be used, such as the analytical relationship between monthly mean precipitation, probability of rainfall occurrence, and daily mean precipitation as proposed by Wilks [49], and also the delta change method by Chen et al. [36]. Moreover, higher-order Markov chains may perform better than the first-order model in simulating precipitation sequence, especially for simulating dry spell length for long-durations over the dry regions [30,36,39,50]. However, as the order increases, the number of parameters also increases exponentially, and the parameter adjustment becomes more complicated [39]. By contrast, the first-order two-state Markov chain only has two parameters, and it is easy to operate and adequate in simulating wet and dry spells for various climates [30,39,49,51].
Since the main objective of this study was to propose a new bias correction method to deal with the biased temporal sequence of climate model simulated precipitation for field and site-specific food security studies, the inter-variable correlations and the inter-site correlation of precipitation time series were not considered. While the correlation between precipitation and temperature, as well as solar radiation, is an essential input to crop and hydrological models, for some other studies, such as hydrological impact studies, the precipitation time series over multiple sites within a watershed may be needed to study the spatial variability of climate change impacts. However, the inter-variable and inter-site correlation can be easily induced using a post-processing method, such as the distribution-free shuffle approach [52][53][54] and the copula-based method [55,56]. The post-processing method can be an avenue for future studies.