- Email: [email protected]
A method to extend temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates
Journal Pre-proofs Research papers A method to extend temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates Yuan Li, Q.J. Wang, Hai He, Zhiyong Wu, Guihua Lu PII: DOI: Reference:
S0022-1694(19)31090-X https://doi.org/10.1016/j.jhydrol.2019.124355 HYDROL 124355
To appear in:
Journal of Hydrology
Received Date: Revised Date: Accepted Date:
17 June 2019 13 October 2019 12 November 2019
Please cite this article as: Li, Y., Wang, Q.J., He, H., Wu, Z., Lu, G., A method to extend temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates, Journal of Hydrology (2019), doi: https://doi.org/ 10.1016/j.jhydrol.2019.124355
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Published by Elsevier B.V.
Manuscript submitted to Journal of Hydrology
A method to extend temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates Yuan LI1, 2, Q. J. WANG2, Hai HE1, Zhiyong WU1, and Guihua LU1
1College
of Hydrology and Water Resources, Hohai University, Nanjing
210098, China 2Department
of Infrastructure Engineering, The University of Melbourne,
Parkville 3010, Australia Corresponding author: Zhiyong WU ([email protected])
Highlights:
A Bayesian joint probability method is presented to calibrate reanalysis estimates
The calibrated estimates are more accurate than the original reanalysis estimates
The uncertainty spread produced from the method is shown to be reliable
Key words: Bayesian joint probability method; precipitation estimates; temporal coverage; reanalysis dataset; uncertainty Declarations of interest: none 1
Manuscript submitted to Journal of Hydrology
Abstract Available high quality precipitation datasets generally cover only recent periods. To extend these datasets back in time is challenging due to historically lower rain gauge network density or poorer remote sensing. In this study, a Bayesian joint probability method is presented to extend the temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates. Relationships between precipitation estimates from a high quality dataset and precipitation estimates from a reanalysis dataset are established by using data from the overlapping period of the two datasets. The relationships are then used to calibrate reanalysis estimates for the period when only reanalysis data is available. Rain gauge observations are also used to enhance the calibration. The method brings the reanalysis dataset to the same spatial resolution as the high quality dataset, corrects bias, and makes the reanalysis data more consistent with the high quality dataset. The calibrated estimates generated are presented in the form of ensembles to represent uncertainty. The method is applied to the Han River basin in China, using the CMPA high quality dataset and the ERA-interim reanalysis dataset. The method is shown to be highly effective in improving the quality of ERA-Interim. 1 Introduction In recent times, there has been significant improvement in the monitoring of precipitation in many parts of the world. In China, the rain gauge network has grown considerably over the last few decades. The automatic rain gauge network, set up in 2006 by the China Meteorological Administration (CMA), now consists of more than 57,000 stations (http://www.observation-cma.com/). In southeastern China, the mean distance
2
Manuscript submitted to Journal of Hydrology
between two nearest automatic rain gauges is approximately 9.6 km (Shen et al., 2018). While rain gauges provide the most accurate point measurement of precipitation, remote sensing can provide better spatial coverage in precipitation estimation across the globe. The NOAA-NCEP Climate Prediction Center (CPC) morphing technique (CMORPH) uses motion vectors derived from half-hourly interval geostationary satellite IR imagery to propagate relatively high quality precipitation estimates over the globe. The Global Precipitation Measurement (GPM) Mission is also expected to provide higher quality precipitation estimates resulting from better orbital inclination, radar frequencies and channels of the passive microwave imager (Huffman et al., 2014). High quality precipitation datasets have been developed by taking advantage of both high density rain gauge network and advanced remote sensing. The Integrated Multi-satellitE Retrievals for GPM (IMERG) is intended to intercalibrate, merge, and interpolate “all” satellite microwave precipitation estimates, together with infrared satellite estimates, and rain gauge analyses at high spatial (0.1°) and temporal (30 min) resolutions over the entire globe (Gaona et al., 2016). The China Merge Precipitation Analysis Hourly V1.0 product (CMPA) used an improved probability density function-optimal interpolation (PDF-OI) method to generate 0.1° and hourly precipitation estimates over China by merging the satellite precipitation estimates from the CMORPH and more than 30,000 automatic rain gauges (Shen et al., 2014). Owing to the high density rain gauge network and rigorous quality control, the CMPA dataset for precipitation is of relatively high quality. Wang et al. (2016) evaluated the CMPA dataset against a dense rain gauge network that was not incorporated in the interpolation process of the CMPA estimates. Their results indicated that the CMPA estimates capture the
3
Manuscript submitted to Journal of Hydrology
spatial patterns of typhoon-related rain depth, and are in agreement with the spatiotemporal evolutions of hourly gauge observations by correlation coefficients from 0.93 to 0.99. A drawback of these high quality precipitation datasets is that they are available only for relatively recent periods. The GPM-IMERG product only covers the period following the launch of the GPM satellite in 2014. The CMPA product is only available after 2008 when the rain gauge network density became relatively high. For many applications, precipitation datasets of much longer temporal coverage are required. Extending the high quality datasets further back in time is challenging due to low rain gauge network density or limitations of remote sensing. The Tropical Rainfall Measuring Mission (TRMM), which was launched in 1997, provided 3-hourly precipitation analysis at a spatial resolution of 0.25° for over 17 years (Huffman et al., 2007; Joyce et al., 2004). However, the accuracy of TRMM precipitation estimates left much to be desired. Huang et al. (2014) suggested that the TRMM multisatellite precipitation analysis (TMPA) products should be used with caution for extreme precipitation events. Reanalysis
data
is
an
alternative
source
of
precipitation
estimates.
Physically-based numerical models are able to generate precipitation estimates at various spatial and temporal resolutions, and the temporal coverage could be longer than both rain gauges and remote sensing. Studies have indicated that reanalysis datasets are able to give reasonable monthly averaged surface temperature, radiation fluxes, and wind speed (Lindsay et al., 2014). However, the accuracy of the reanalysis precipitation estimates still needs to be improved due to the imperfect description of small-scale physical processes of the numerical models (Leeuw et al., 2015; Nkiaka et al., 2017). 4
Manuscript submitted to Journal of Hydrology
One way of improving the reanalysis data is to calibrate the reanalysis data against rain gauge observations. However, this requires the availability of high density rain gauge networks, which are often not available. In this study, we introduce a method to calibrate the reanalysis data by taking advantage of the availability of high quality datasets (such as CMPA) for recent times. Using the period when both datasets are available, we establish relationships between precipitation estimates of the high quality dataset and precipitation estimates of the reanalysis dataset. The relationships can then be applied to calibrate the reanalysis estimates for the period when only reanalysis data is available. The method brings the reanalysis dataset to the same spatial resolution as the high quality dataset, corrects bias, and makes the reanalysis data more consistent with the high quality dataset. The calibrated estimates generated from our calibration method are presented in the form of ensembles to represent uncertainty. The particular method we use is based on the Bayesian joint probability (BJP) approach. The BJP approach, first developed for flood record augmentation, has found wide applications (Robertson and Wang, 2012; Schepen et al., 2018; Strazzo et al., 2019; Wang, 2001; Wang et al., 2009; Zhao et al., 2019a; Zhao et al., 2019b; Zhao et al., 2016). In BJP, data normalization is used to handle skewed variables (Wang et al., 2012), missing values are permitted, and zero values are treated as censored data. In addition, multiple data sources are allowed due to their multivariate formulation. In previous studies, the BJP method is always used to establish relationships between predictors and predictand for each grid cell. However, the intersite correlations have never been considered before. In this study, the BJP method is used as a spatial model to quantifying intersite correlations between grid cells and rain gauges from different locations.
5
Manuscript submitted to Journal of Hydrology
The efficacy of this calibration method for extending the temporal coverage of the CMPA precipitation dataset by calibrating the ERA-interim reanalysis over the Han River basin in China is evaluated in this study. The Han River basin is the water source of the middle route of the South-to-North Water Diversion Project (SNWDP) in China. Therefore, availability of accurate precipitation estimates of long temporal coverage over this region is important. The remainder of the paper is structured as follows. The study area and data sources are introduced in Section 2. An introduction of the BJP modelling method is provided and the evaluation methods are introduced in Section 3. Section 4 presents the evaluation results of the ensemble mean and ensemble distribution of calibrated precipitation estimates. Conclusions are drawn in Section 5.
2 Study area and data The Han River basin is located between 106-114°E and 29-34°N with an area of nearly 170,400 km2 (Figure 1). The basin is subject to sub-tropical monsoon climate, and its topography is highly complex. The elevation of western Han River basin is over 3000 m, while the eastern Han River basin is composed of recent alluvial plains (Li et al., 2017).
6
Manuscript submitted to Journal of Hydrology
Figure 1. Location of the Han River basin, rain gauges, ERA-interim grid, and the China Merge Precipitation Analysis Hourly V1.0 product (CMPA) grid. Figure 2 illustrates the spatial characteristics of daily precipitation for each month. It is clear that the precipitation shows high spatial and temporal variability over the Han River basin, especially in the wet months from May to September. In July, the daily mean precipitation is greater than 5 mm/day with a standard deviation of greater than 10 mm/day in the upper and lower sub basins over the period of 2008 – 2016.
7
Manuscript submitted to Journal of Hydrology
Figure 2. Spatial distribution of daily mean and daily standard deviation of CMPA estimates over the period of 2008 – 2016 for each month (mm/day). Both the ERA-interim and the CMPA datasets are used in this study. A brief introduction of these two datasets has been presented in previous section. In addition, rain gauge observations of long temporal coverage are also used as inputs to bring additional information to the calibration. The rain gauge observations are extracted from the “Daily Surface Climate Variables of China” dataset (SURF_CLI_CHN_MUL_DAY V3.0). It was released by the China Meteorological Administration (CMA) with high quality control. Out of the 756 long term national rain gauges, 48 rain gauges are located in or around the Han River basin and are used in this study (Figure 1). The daily ERA-interim reanalysis precipitation (hereafter Original ERA_daily) is accumulated from two forecasts (at 00:00 and 12:00) of 12 hour precipitation each. The CMPA daily precipitation (hereafter CMPA_daily) is accumulated from hourly values. The
8
Manuscript submitted to Journal of Hydrology
rain gauge observations (hereafter GAUGE_daily) are already daily. Details of the final datasets used in this study are summarized in Table 1. Table 1. A summary of daily precipitation datasets used in this study, and the role of each dataset in the BJP method. Dataset
Spatial resolution
Temporal coverage
Role in the BJP
CMPA_daily
0.1°
2008 - 2016
Predictand
Original ERA_daily
0.75°
1979 - 2016
Predictor
GAUGE_daily
~7000 km2/gauge
1979 - 2016
Predictor
3 Methods 3.1 The BJP method 3.1.1 Model formulation In this study, the BJP method is applied to calibrate the ERA_daily using the CMPA_daily. The CMPA_daily is treated as the most accurate owing to the high rain gauge network density and advanced remote sensing technique. We denote the value of the given 0.1° CMPA_daily grid cell as the predictand Y, the corresponding value of the nearest ERA_daily as the predictor 𝑋1. To enhance the calibration, rain gauge observations are also used for calibration as well. The distances from each 0.1° CMPA_daily grid cell to the rain gauges are presented in Figure 3. The mean distance between CMPA grid cells and the corresponding nearest gauge is approximately 30 km, while the corresponding second nearest rain gauge is nearly 60 km. Thus, only the two
9
Manuscript submitted to Journal of Hydrology
nearest rain gauges are selected as predictors 𝑋2 and 𝑋3 for the given CMPA_daily grid cell here.
Figure 3. The inter-distances between CMPA grid cells and rain gauges over the Han River basin. The probability distribution of daily precipitation is always highly skewed, which makes it difficult for parameter inference. To over this problem, we transform both the predictand Y and the predictors 𝑋𝑇 = [𝑋1 𝑋2 𝑋3] to normal distributions first. The Box-cox transformation has been widely used for transformation. However, the one-parameter Box-Cox transformation is not always able to achieve variance stabilization, while the 10
Manuscript submitted to Journal of Hydrology
two-parameter Box-Cox transformation leads to unrealistically large uncertainty. The log-sinh transformation is able to stabilize the variance and normalize the variable. It has been proved to outperform both the one-parameter and two-parameter Box-Cox transformations (Wang et al., 2012). This transformation has also been tested on a wide range of hydrological data, including hydrological extremes (Bennett et al., 2014a; Bennett et al., 2014b). In this study, the predictors 𝑋𝑇 (the matrix of the three datasets) are normalized to 𝑈𝑇 = [𝑈1 𝑈2 𝑈3], and the predictand Y is normalized to 𝑉 using the log-sinh transformation, 𝑢𝑖 = 𝛽𝑋―1 ln (sinh (𝛼𝑋𝑖 + 𝛽𝑋𝑖𝑥𝑖)) 𝑖
(1)
𝑣 = 𝛽𝑌―1ln (sinh (𝛼𝑌 + 𝛽𝑌𝑦))
(2)
where 𝛼𝑋,𝑖, 𝛽𝑋,𝑖, 𝛼𝑌, and 𝛽𝑌 are unknown transformation parameters for predictors and predictand, 𝑖 = 1,⋯, 3. The matrix 𝑍𝑇 = [𝑈1 𝑈2 𝑈3 𝑉] is then assumed to follow a multivariate normal distribution. 𝑍~𝐍(𝝁, ∑)
(3)
where 𝝁 and ∑ are the mean vector and covariance matrices to be estimated. 𝝁𝑻 = [𝜇𝑈1 𝜇𝑈2 𝜇𝑈3 𝜇𝑉]
(4)
∑ = 𝝈𝑹𝝈𝑻
(5)
where 𝝈 and 𝑹 are the standard deviation vector and correlation coefficient matrix, respectively: 𝝈𝑻 = [𝜎𝑈1 𝜎𝑈2 𝜎𝑈3 𝜎𝑉] 11
(6)
Manuscript submitted to Journal of Hydrology
[
1 𝑟𝑈1 𝑈2 𝑟𝑈1 𝑈3 𝑟𝑈1 𝑉 𝑟𝑈2 𝑈1 1 𝑟𝑈2 𝑈3 𝑟𝑈2 𝑉 𝐑= 𝑟 1 𝑟𝑈3 𝑉 𝑈3 𝑈1 𝑟𝑈3 𝑈2 𝑟𝑉 𝑈1 𝑟𝑉 𝑈2 𝑟𝑉 𝑈3 1
]
(7)
Note that the correlation coefficient matrix 𝐑 is symmetric. Thus, the joint distribution of predictand and predictors has a total of 14 unknown parameters. Zero values occasionally occur in dry months for both predictors and predictand over the Han River basin. This always leads to a mixed discrete-continuous multivariate distribution, which is difficult to formulate. In this study, we follow Wang and Robertson (2011) and treat zero values as censored data with unknown values below or equal to zero. In this way, the predictors and the predictand are considered to follow continuous multivariate distributions. This treatment allows us to simplify the problem of mixed discrete-continuous multivariate distribution, and to deal with zero values. 3.1.2 Parameter inference The transformation parameters and the joint distribution parameters are estimated when the ERA_daily, GAUGE_daily, and CMPA_daily are concurrently available. The choice of training data is still debated as parameter inference is affected by seasonality effects. Pooling all data from different months is not suitable. One way to solve this problem is to pre-deseasonalize data before parameter inference. Studies have proved that deseasonalizing data could potential improve forecast skills when using neural network forecasting models (Nelson et al., 1999; Zhang and Kline, 2007). A more direct approach is to pool data with similar seasonality. Hamill (2007) suggested that including days from the same season would expand the training period, and would likely improve the performance of the model. As shown in Fig. 2, the characteristics of daily 12
Manuscript submitted to Journal of Hydrology
precipitation vary from month to month. Thus, the daily values of each month are pooled together for each variable. There are approximately 270 days in total that are selected to make parameter inferences for a given grid cell of CMPA_daily in each month. Given
a
data
series 𝐃(𝑋𝑖)
of
each
predictor 𝑋𝑖, 𝐃(𝑋𝑖) =
{𝑥𝑖(t), t = 1, ⋯, n;i = 1, 2 , 3}, the transformation parameters of 𝛼𝑋,𝑖 and 𝛽𝑋,𝑖 are estimated using a Bayesian maximum a posteriori (MAP) solution. The resulting parameter values of 𝛼𝑋,𝑖 and 𝛽𝑋,𝑖 are taken as point estimates. Parameter uncertainty is not considered for these parameters. The transformation parameters of 𝛼𝑌 and 𝛽𝑌 are estimated using the same strategy as the predictors, given a data series 𝐃(𝑌) = {𝑦(t), t = 1, ⋯, n}. The joint distribution parameters are then estimated using transformed data series
of
both
predictors
and
predictand,
𝐃(U1,𝑈2,𝑈3,𝑉) =
{(𝑢1(t), 𝑢2(t),𝑢3(t), 𝑣(t)), t = 1, ⋯, n}. Here, we denote the 14 unknown parameters of the joint distribution as 𝜽 = {𝝁, ∑}. The posterior distribution of 𝜽 is estimated using a Bayesian framework, 𝑝(𝜽|𝐃(𝑈1,𝑈2,𝑈3,𝑉)) ∝ 𝑝(𝜽)𝑝(𝐃(𝑈1,𝑈2,𝑈3,𝑉)|𝜽)
(8)
where 𝑝(𝜽) is the prior distribution of parameters, and 𝑝(𝐃(𝑈1,𝑈2,𝑈3,𝑉)|𝜽) is the likelihood. In this study, non-informative distributions are used as priors of parameters 𝜽, 𝑝(𝜽) ∝ 1. As the posterior distributions of parameters 𝜽 are not standard distributions, analytical integration is not appropriate. To overcome this problem, we use the Markov chain Monte Carlo technique to draw a sample of parameter values which represent the joint posterior distribution of parameters 𝜽. Specifically, we use the Gibbs Sampling 13
Manuscript submitted to Journal of Hydrology
algorithm to draw (1000) samples from the conditional posterior distributions. To diminish the influence of the starting values, we discard the first 5000 iterations in Markov chain simulation as burn-in period. We also set the length of the Markov chain as 30, 000, which is proportional to the number of variables to obtain sufficient number of representative samples of the posterior distribution (Zhao et al., 2016). A visual inspection of the sampled parameter chains is made to monitor the convergence of the sample (Cowles and Carlin, 1996). 3.1.3 Model use for calibration Once the parameters are sampled, the BJP model can be used to calibrate the ERA_daily when the CMPA_daily is not available. Given a new set of predictor values, 𝐃(𝑋𝑖∗ ) = {𝑥𝑖(t ∗ ), t ∗ = 1, ⋯, n}, transformed values 𝐃(𝑈𝑖∗ ) = {𝑢𝑖(t ∗ ), t ∗ = 1, ⋯, n} can be found using the estimated transformation parameters 𝛼𝑋,𝑖 and 𝛽𝑋,𝑖. The posterior predictive distribution of 𝑣(t ∗ ) is given by 𝑓𝑉(𝑣(t ∗ )) = ∫𝑝(𝑣(t ∗ )|𝑢1(t ∗ ), 𝑢2(t ∗ ),𝑢3(t ∗ ),𝜽)𝑝(𝜽)𝑝(𝐃(𝑈1,𝑈2,𝑈3,𝑉)|𝜽)d𝜽
(9)
Again, the Gibbs Sampling algorithm is used to obtain samples of 𝑓𝑉(𝑣(t ∗ )) by giving each of the 1000 sets of parameter values 𝜽. The samples are then back-transformed to produce a calibrated precipitation estimates of 𝑦(t ∗ ) using parameters 𝛼𝑌 and 𝛽𝑌. Finally, we develop the BJP model for each grid cell of the CMPA_daily and each month separately.
14
Manuscript submitted to Journal of Hydrology
3.2 Evaluation Methods 3.2.1 Cross validation To assess the performance of the BJP calibration method, a ten-fold cross validation scheme is used here. The daily values of both the predictors and the predictand in each month from 2008 to 2016 are pooled together. Both the predictors and the predictand are then partitioned into ten equal-sized subsamples. Of the ten subsamples, nine subsamples of approximately 240 daily values are used for training, while the remaining one single subsample of nearly 30 daily values are used for validation. This process is then repeated ten times, with each of the ten subsamples used exactly once as the validation data. The independently calibrated ERA daily estimates (hereafter Calibrated ERA_daily) for all the 10 sub-samples are pooled together and compared against the CMPA_daily estimates. 3.2.2 Evaluation of ensemble mean The ensemble mean of Calibrated ERA_daily for each 0.1° CMPA grid cell is calculated and evaluated for absolute value of Relative Bias (RB), root mean square error (RMSE), and Nash–Sutcliffe coefficient of efficiency (NSE): ∑𝑡 = 1|𝑦𝑡∗ ― 𝑜𝑡| 𝑇
RB =
𝑇
∑𝑡 = 1 𝑜𝑡
1 𝑇 ∗ ∑ 𝑇 𝑡 = 1(𝑦𝑡
RMSE =
2
― 𝑜𝑡)
(11)
2
∑𝑡 = 1(𝑦𝑡∗ ― 𝑜𝑡) 𝑇
NSE = 1 ―
(10)
× 100%
𝑇
∑𝑡 = 1(𝑜 ― 𝑜𝑡)2 15
(12)
Manuscript submitted to Journal of Hydrology
where 𝑦𝑡∗ is the ensemble mean of the Calibrated ERA_daily precipitation for day 𝑡, 𝑡 = 1,2,⋯,𝑇; 𝑜𝑡 is the CMPA_daily precipitation, which is treated as the most accurate precipitation estimates in this study; 𝑜 is the averaged CMPA_daily precipitation over the period of analysis. The RB is used to evaluate the agreement between the averaged value of ensemble mean of Calibrated ERA_daily and the CMPA_daily, while the RMSE is used to evaluate the average error magnitude. The NSE is selected to evaluate the relative magnitude of the residual variance of the ensemble mean of Calibrated ERA_daily compared to the CMPA_daily variance. 3.2.3 Evaluation of ensemble distribution The continuous ranked probability score (Matheson and Winkler, 1976), and the α-index (Renard et al., 2010) are used to evaluate the accuracy and the reliability of the fit of distribution in this study. The continuous ranked probability score (CRPS) is used to provide an overall evaluation of the Calibrated ERA_daily for each 0.1° CMPA grid cell: CRPS 1 = 𝑇
𝑇
∑∫[𝐹(𝑦
∗ 𝑡
) ― 𝐻(𝑦𝑡∗ ― 𝑜𝑡)]2𝑑𝑦
where 𝑦𝑡∗ is
𝑡=1
(13) the daily precipitation to be estimated for day t; 𝐹() is the cumulative distribution function of the calibrated ERA_daily ensemble, as represented by the sampled (1000) values; and 𝐻() is the Heaviside step function defined as
16
Manuscript submitted to Journal of Hydrology
𝐻(𝑦𝑡∗ ― 𝑜𝑡) =
{
0 1
𝑦𝑡∗ < 𝑜𝑡 𝑦𝑡∗ ≥ 𝑜𝑡
(14)
where 𝑜𝑡 is the CMPA_daily precipitation, which is treated as the actual precipitation. When 𝐹() is replaced by a step distribution based on a deterministic value (such as the Original ERA_daily), the CRPS is reduced to the mean absolute error (MAE). A CRPS skill score is then calculated by comparing the CRPS of the Calibrated ERA_daily with the CRPS of a reference distribution for estimating precipitation:
CRPSSS =
CRPSREF ― CRPS CRPSREF
× 100%
(15)
In this study, a cross-validated climatology of CMPA_daily is used as the reference. A skill score of 100% indicates that the calibrated ERA_daily ensemble is the same as the CMPA_daily observation, whereas a skill score of 0% suggests that the calibrated ERA_daily ensemble shows no improvement over the cross-validated climatology. A negative skill score means that the calibrated ERA_daily ensemble is inferior to the cross-validated climatology. The reliability of Calibrated ERA_daily is evaluated by the probability integral transforms (PITs) of CMPA_daily. The PIT of 𝑜𝑡 is given by 𝜋𝑡 = 𝐹(𝑜𝑡)
(16)
where F() is the cumulative distribution function of the calibrated ERA_daily ensemble. When 𝑜𝑡 is equal to zero, the precise value for 𝜋𝑡 cannot be found from equation (16). Instead, a pseudo-PIT value is sampled from a uniform distribution with a range of [0, 𝜋𝑡], and subsequently using it with other real 𝜋𝑡 values to construct the PIT plots (Schepen et
17
Manuscript submitted to Journal of Hydrology
al., 2018; Wang and Robertson, 2011). If the Calibrated ERA_daily ensemble is reliable, 𝜋𝑡 should be uniformly distributed. The 𝜋𝑡 values are then summarized into an α-index:
α = 1.0 ―
2 𝑇
∑
𝑇
|𝜋
∗ 𝑡
𝑡=1
―
|
𝑡 𝑇+1
(17) where 𝜋𝑡∗ is the
sorted 𝜋𝑡 in increasing order. The α-index ranges from 0 to 1, and represents the total deviation of 𝜋𝑡∗ from corresponding uniform quantile. A higher α-index indicates higher reliability. 3.3 Further Analysis In this study, an important assumption is made that the CMPA precipitation estimates are of high quality, as such they can be used to calibrate the ERA reanalysis estimates. It is not right to directly verify assumption as follows by comparing CMPA precipitation estimates with observations at gauges, as the observations are likely to have been used in deriving the CMPA precipitation estimates. Here we indirectly evaluate the assumption as follows. We use the 2008-2016 data to establish models to calibrate ERA reanalysis estimates against the CMPA precipitation estimates only. Unlike the models used in Section 3.1, gauged observations are not used here.
The new calibration models are then used to calibrate the ERA
reanalysis estimates for 1979-2007. The calibrated ERA reanalysis estimates are compared with gauged observations. For the calibration to lead to improvement over the raw ERA reanalysis estimates, the CMPA precipitation estimates would have to be more accurate than the ERA reanalysis estimates.
18
Manuscript submitted to Journal of Hydrology
This further analysis is to ensure that the quality of the CMPA precipitation estimates are high enough that its use adds value to the calibration and does not adversely corrupt the ERA reanalysis estimates. In practical applications, the gauged observations are also used as predictors to derive the final estimates.
4 Results 4.1 Ensemble mean The RB and RMSE of the ensemble mean of Calibrated ERA_daily for each month are plotted and compared with the Original ERA_daily in Figure 4 and Figure 5. The RB and RMSE of the Original ERA_daily are greatly reduced after the BJP calibration, especially in wet months in summer. In July, the RB of the Original ERA_daily is always greater than 150%, with RMSE values greater than 10 mm/day. In contrast, the Calibrated ERA_daily shows much higher accuracy in terms of RB values of lower than 50% and RMSE values of lower than 8 mm/day in most CMPA_daily grid cells. The reduction of RB and RMSE values are also evident in other months. The RB values reduce from nearly 100% to less than 50% while the RMSE values reduce from approximately 5 mm/day to less than 2 mm/day in months in spring and autumn.
19
Manuscript submitted to Journal of Hydrology
Figure 4. RB of Original ERA_daily and RB of ensemble mean of Calibrated ERA_daily compared to CMPA_daily for each month (%).
Figure 5. RMSE of Original ERA_daily and RMSE of ensemble mean of Calibrated ERA_daily compared to CMPA_daily for each month (mm/day). 20
Manuscript submitted to Journal of Hydrology
Figure 6 compares the NSE of the Original ERA_daily and the ensemble mean of Calibrated ERA_daily. The results also suggest that the accuracy of the Original ERA_daily has been greatly improved after the BJP calibration. In July, the NSE values of the Original ERA_daily are always less than 0.2. In contrast, the NSE values of the ensemble mean of calibrated ERA_daily are greater than 0.5 for most grid cells. The NSE
values increase from less than 0.4 to more than 0.6 after the BJP calibration in
months in spring and autumn as well. However, we should also notice that the NSE values are still low from December to Februrary owing to the relatively lower precipitation in these months shown in Fig. 2.
Figure 6. NSE of Original ERA_daily and NSE of ensemble mean of Calibrated ERA_daily compared to CMPA_daily for each month.
21
Manuscript submitted to Journal of Hydrology
4.2 Ensemble distribution The CRPS skill score of Original ERA_daily and the CRPS skill score of Calibrated ERA_daily are compared in Figure 7. The results suggest that the ensemble distribution of Calibrated ERA_daily is more accurate than the Original ERA_daily. The Calibrated ERA_daily usually exhibits positive skill in most areas in July. Whereas, negative skill is observed for Original ERA_daily. Improvement is observed in other months as well in terms of CRPS skill score.
Figure 7. CRPS skill score of Original ERA_daily and CRPS skill score of Calibrated ERA_daily compared to CMPA_daily for each month (%). The α-index of Calibrated ERA_daily for each month is presented in Figure 8. It is clear that the α-index is mostly in the range of 0.9 to 1.0 for all months, indicating that the Calibrated ERA_daily is of high reliability. The PIT diagram of a selected grid cell with lowest α-index value of 0.825 in February is also given in Figure 8. Although small 22
Manuscript submitted to Journal of Hydrology
departures from the diagonal line are observed, the PIT values for the selected grid cell lie mostly within the Kolmogorov 5% significance bands, suggesting that the uncertainty spread is reliable for all months.
Figure 8. α-index of Calibrated ERA_daily, and PIT uniform probability plot for the selected grid cell with lowest α-index value (1:1 solid line, theoretical uniform distribution; dashed lines, Kolmogorov 5% significance band; blue circles, PIT values of CMPA_daily precipitation; red circles, pseudo PIT values).
4.3 Calibrated ERA_daily precipitation for 1979 - 2007 The BJP models are subsequently established based on all available 2008 - 2016 data of 0.1° CMPA_daily precipitation, Original ERA_daily, and GAUGE_daily. The models are then used to produce Calibrated ERA_daily from the Original ERA_daily and GAUGE_daily for 1979 – 2007, for which 0.1° CMPA_daily is not available. 23
Manuscript submitted to Journal of Hydrology
Figure 9 compares the daily precipitation distribution of the Original ERA_daily, the GAUGE_daily, and the ensemble mean of Calibrated ERA_daily for an extreme precipitation event in October 1983. Beginning on 3 October 1983 and lasting to 5 October 1983, an extreme precipitation event occurred in the Han River basin. Similar to the example shown in Section 4.1, the Original ERA_daily is too low in spatial variability compared with the GAUGE_daily. The high precipitation intensity observed by the GAUGE_daily is greatly underestimated by the Original ERA_daily. The ensemble mean of Calibrated ERA_daily takes guidance from both the Original ERA_daily and GAUGE_daily. The spatial characteristics of the Original ERA_daily are broadly preserved, while high precipitation intensity of the GAUGE_daily is well captured by the ensemble mean of Calibrated ERA_daily.
24
Manuscript submitted to Journal of Hydrology
Figure 9. Daily precipitation distribution of the ensemble mean of Calibrated ERA_daily, Original ERA_daily, and GAUGE_daily for an extreme precipitation event from 0000 UTC 03 October 1983 to 0000 UTC 05 October 1983. 4.4 Further Analysis Recalling Section 3.3, the purpose of this further analysis is to ensure that the quality of the CMPA precipitation estimates is high enough that its use adds value to the calibration and does not adversely corrupt the ERA reanalysis estimates. We use the 2008-2016 data to establish models to calibrate ERA reanalysis estimates against the CMPA precipitation estimates only. Unlike the models used in Section 3.1, gauged
25
Manuscript submitted to Journal of Hydrology
observations are not used here. The new calibration models are then used to calibrate the ERA reanalysis estimates for 1979-2007. We should notice that the direct comparison between gauged observations and gridded estimates should be avoided due to the different representativeness of these two datasets. Currently, there are three indirect approaches to compare the gauged observations and gridded estimates. Zhao and Fu (2006) suggested using the kriging technique to upscale gauged observations onto the grids and then compare with the gridded estimates. However, the upscaling approach will introduce large errors and uncertainties when the gauges are sparsely distributed. Another approach is downscaling the gridded estimates at each station site using the bilinear interpolation method (Gao and Liu, 2012). Nevertheless, the downscaling approach may smooth extremes from the gridded data (Accadia et al., 2003). The third approach is to compare the gridded estimates and averaged values of gauged observations falling within the corresponding grid cell. This approach has been proved to work well with precipitation estimates comparisons (Blacutt et al., 2015; Ward et al., 2011). However, we should notice that spatial sampling error exists owing to small spatial scale of rainfall systems compared to rain gauge network density (Krajewski et al., 2000; Morrissey et al., 1995). Therefore, it is practically impossible to identify an optimal method applicable for all areas (Duan et al., 2016). In this study, the third approach is adapted here to compare the Calibrated ERA_daily and the Original ERA_daily against the gauged observations. The Calibrated ERA_daily is at a spatial resolution of 0.1° as is the CMPA_daily. Thus, there is only one gauge falling within each 0.1° grid cell of Calibrated ERA_daily. We directly compare the Calibrated ERA_daily and the corresponding gauged 26
Manuscript submitted to Journal of Hydrology
observations here. For grid cells of the Original ERA_daily containing more than one gauge, the averaged values of the gauges are calculated. The results are then compared against the Original ERA_daily. The NSE of Original ERA_daily and ensemble mean of Calibrated ERA_daily for the period of 1979-2007 are then presented in Fig. 10. The NSE of the Calibrated ERA_daily is much more consistent among gauge stations than the Original ERA_daily. The latter can be very poor in performance in some locations and months. In most months, the median NSE of the Calibrated ERA_daily is better than the Original ERA_daily. Overall, the calibration has led to considerable improvement. This analysis indicates that the use of the CMPA precipitation estimates as a high quality dataset is reasonable at least for this study case.
27
Manuscript submitted to Journal of Hydrology
Figure 10. Box plots of NSE of Original ERA_daily and ensemble mean of Calibrated ERA_daily for the period of 1979-2007. The red lines are the 50th percentiles, top and bottom of each box are the 75th and 25th percentiles, and whiskers are the 95th and 5th percentiles. 5 Conclusions This paper presents a Bayesian joint probability (BJP) method to extend the temporal coverage of high quality precipitation datasets by calibrating reanalysis
28
Manuscript submitted to Journal of Hydrology
estimates. The key concept is as follows. Relationships between precipitation estimates from a high quality dataset and precipitation estimates from a reanalysis dataset can be established by using data from the overlapping period of the two datasets. The relationships can then be used to calibrate reanalysis estimates for the period when only reanalysis data is available. Additionally, the method makes use of long records from rain gauges to enhance the calibration. The method brings the reanalysis dataset to the same spatial resolution of the high quality dataset, corrects bias, and makes the reanalysis data more consistent with the high quality dataset. The calibrated estimates generated are presented in the form of ensembles to represent uncertainty. In this study, the BJP method is used to extend the temporal coverage of a high quality precipitation dataset – CMPA (2008 - 2016) by calibrating the ERA-interim reanalysis precipitation estimates (1979 - 2016) over the Han River basin. Records from a sparsely distributed rain gauge network (1979 - 2016) are also used to enhance the calibration. To evaluate the accuracy and reliability of calibrated precipitation estimates, a ten-fold cross-validation scheme is applied to the overlapping period (2008 - 2016) of the three datasets. The original reanalysis estimates and the calibrated reanalysis estimates are evaluated against the precipitation estimates from the high quality dataset. Further analysis is performed to ensure that the quality of the CMPA precipitation estimates is high enough that its use adds value to the calibration and does not adversely corrupt the ERA reanalysis estimates. The calibrated reanalysis estimates are shown to be much more accurate than the original reanalysis estimates, in terms of RB, RMSE and CRPS skill score. The 29
Manuscript submitted to Journal of Hydrology
uncertainty spread produced from the calibration is shown to be reliable. The CMPA dataset is also proved to be effective to calibration. We conclude that the method introduced in this study is highly effective in improving the quality of ERA-Interim. The fifth generation of reanalysis dataset from ECMWF (ERA-5) will be published with a horizontal resolution of nearly 30 km from 1950 to the present by early 2019 (Brönnimann et al., 2018). The application of the calibration method developed in this study to the new reanalysis products will be straightforward.
Acknowledgments We would like to thank Geoff Pegram and another two anonymous reviewers for their insightful comments. This work was funded by the National Key R&D Program of China (Grant No. 2018 YFC0407701, 2017YFC1502403), the Fundamental Research Funds for the Central Universities (Grants No. 2019B10214, 2017B681X14), the National Natural Science Foundation of China (Grants No. 51579065, 51779071), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant No. KYCX17_0413), and the China Scholarship Council scholarship. We also thank Dr Tongtiegang Zhao for providing technical support to the work. References Accadia, C., Mariani, S., Casaioli, M., Lavagnini, A., Speranza, A., 2003. Sensitivity of precipitation forecast skill scores to bilinear interpolation and a simple nearest-neighbor average method on high-resolution verification grids. Weather
30
Manuscript submitted to Journal of Hydrology
and
forecasting,
18(5):
918-932.
https://doi.org/10.1175/1520-0434(2003)018<0918:SOPFSS>2.0.CO;2. Bennett, J.C. et al., 2014a. A System for Continuous Hydrological Ensemble Forecasting (SCHEF)
to
lead
times
of
9
days.
J
Hydrol,
519:
2832-2846.
https://doi.org/10.5194/nhess-14-219-2014. Bennett, J.C., Wang, Q., Pokhrel, P., Robertson, D.E., 2014b. The challenge of forecasting high streamflows 1–3 months in advance with lagged climate indices in southeast Australia. Natural Hazards and Earth System Sciences, 14(2): 219-233. https://doi.org/10.1016/j.jhydrol.2014.08.010. Brönnimann, S. et al., 2018. Observations for reanalyses. Bulletin of the American Meteorological
Society,
99(9):
1851-1866.
https://doi.org/10.1175/BAMS-D-17-0229.1. Blacutt, L.A., Herdies, D.L., de Gonçalves, L.G.G., Vila, D.A., Andrade, M., 2015. Precipitation comparison for the CFSR, MERRA, TRMM3B42 and Combined Scheme
datasets
in
Bolivia.
Atmos
Res,
163:
117-131.
https://doi.org/10.1016/j.atmosres.2015.02.002. China Merge Precipitation Analysis Hourly V1.0 product (2018). China Meteorological Administration,
retrieved
from
http://data.cma.cn/data/detail/dataCode/SEVP_CLI_CHN_MERGE_CMP_PRE_ HOUR_GRID_0.10/. Cowles, M., and Carlin B. P., 1996. Markov chain Monte Carlo convergence diagnostics: A
comparative
review,
J.
Am.
https://doi.org/10.2307/2291683.
31
Stat.
Assoc.,
91,
883–904.
Manuscript submitted to Journal of Hydrology
Daily Surface Climate Variables of China (2018). China Meteorological Administration, retrieved
from
http://data.cma.cn/data/detail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0/. Duan, Z., Liu, J., Tuo, Y., Chiogna, G., Disse, M., 2016. Evaluation of eight high spatial resolution gridded precipitation products in Adige Basin (Italy) at multiple temporal and spatial scales. Science of the Total Environment, 573: 1536-1553. https://doi.org/10.1016/j.scitotenv.2016.08.213. Gao, Y., Liu, M., 2012. Evaluation of high-resolution satellite precipitation products using rain gauge observations over the Tibetan Plateau. Hydrology & Earth System Sciences Discussions, 9(8). https://doi.org/10.5194/hess-17-837-2013. Gaona, M.F.R., Overeem, A., Leijnse, H., Uijlenhoet, R., 2016. First-Year Evaluation of GPM Rainfall over the Netherlands: IMERG Day 1 Final Run (VO3D). J Hydrometeorol, 17(11): 2799-2814. https://doi.org/10.1175/JHM-D-16-0087.1. Hamill, T. M. (2007), Comments on “Calibrated Surface Temperature Forecasts from the Canadian Ensemble Prediction System Using Bayesian Model Averaging”, Monthly
Weather
Review,
135(12),
4226-4230,
https://doi.org/10.1175/2007mwr1963.1. Huang, Y. et al., 2014. Evaluation of Version-7 TRMM Multi-Satellite Precipitation Analysis Product during the Beijing Extreme Heavy Rainfall Event of 21 July 2012. Water, 6(1): 32-44. https://doi.org/10.3390/w6010032. Huffman, G.H. et al., 2014. NASA Global Precipitation Measurement (GPM) Integrated Multi-satellitE Retrievals for GPM (IMERG). Algorithm theoretical basis document, version, 4: 30.
32
Manuscript submitted to Journal of Hydrology
Huffman, G.J. et al., 2007. The TRMM multisatellite precipitation analysis (TMPA): Quasi-global, multiyear, combined-sensor precipitation estimates at fine scales. J Hydrometeorol, 8(1): 38-55. https://doi.org/10.1175/JHM560.1. Joyce, R.J., Janowiak, J.E., Arkin, P.A., Xie, P.P., 2004. CMORPH: A method that produces global precipitation estimates from passive microwave and infrared data at high spatial and temporal resolution. J Hydrometeorol, 5(3): 487-503. https://doi.org/10.1175/1525-7541(2004)005<0487:CAMTPG>2.0.CO;2. Krajewski, W.F., Ciach, G.J., McCollum, J.R., Bacotiu, C., 2000. Initial validation of the Global Precipitation Climatology Project monthly rainfall over the United States. Journal
of
Applied
Meteorology,
39(7):
1071-1086.
https://doi.org/10.1175/1520-0450(2000)039<1071:IVOTGP>2.0.CO;2. Leeuw, J.D., Methven, J., Blackburn, M., 2015. Evaluation of ERA-Interim reanalysis precipitation products using England and Wales observations. Q J Roy Meteor Soc, 141(688): 798-806. https://doi.org/10.1002/qj.2395. Li, Y., Lu, G.H., Wu, Z.Y., He, H., He, J., 2017. High-Resolution Dynamical Downscaling of Seasonal Precipitation Forecasts for the Hanjiang Basin in China Using the Weather Research and Forecasting Model. J Appl Meteorol Clim, 56(5): 1515-1536. https://doi.org/10.1175/Jamc-D-16-0268.1. Lindsay, R., Wensnahan, M., Schweiger, A., Zhang, J., 2014. Evaluation of Seven Different Atmospheric Reanalysis Products in the Arctic. J Climate, 27(7): 2588-2606. https://doi.org/10.1175/JCLI-D-13-00014.1.
33
Manuscript submitted to Journal of Hydrology
Matheson, J.E., Winkler, R.L., 1976. Scoring rules for continuous probability distributions.
Management
science,
22(10):
1087-1096.
https://doi.org/10.1287/mnsc.22.10.1087. Morrissey, M.L., Maliekal, J.A., Greene, J.S., Wang, J., 1995. The uncertainty of simple spatial averages using rain gauge networks. Water Resources Research, 31(8): 2011-2017. https://doi.org/10.1029/95WR01232. Nelson, M., Hill, T., Remus, W., O'Connor, M., 1999. Time series forecasting using neural networks: Should the data be deseasonalized first? Journal of forecasting, 18(5):
359-367.
https://doi.org/10.1002/(SICI)1099-131X(199909)18:5<359::AID-FOR746>3.0.C O;2-P. Nkiaka, E., Nawaz, N.R., Lovett, J.C., 2017. Evaluating global reanalysis precipitation datasets with rain gauge measurements in the Sudano‐Sahel region: case study of the Logone catchment, Lake Chad Basin. Meteorol Appl, 24(1): 9-18. https://doi.org/10.1002/met.1600. Poli, P. et al., 2016. ERA-20C: An Atmospheric Reanalysis of the Twentieth Century. J Climate, 29(11): 4083-4097. https://doi.org/10.1175/Jcli-D-15-0556.1. Renard, B., Kavetski, D., Kuczera, G., Thyer, M., Franks, S.W., 2010. Understanding predictive uncertainty in hydrologic modeling: The challenge of identifying input and
structural
errors.
Water
https://doi.org/10.1029/2009wr008328.
34
Resour
Res,
46.
Manuscript submitted to Journal of Hydrology
Robertson, D.E., Wang, Q.J., 2012. Bayesian Approach to Predictor Selection for Seasonal
Strearnflow
Forecasting.
J
Hydrometeorol,
13(1):
155-171.
https://doi.org/10.1175/Jhm-D-10-05009.1. Schepen, A., Zhao, T.T.G., Wang, Q.J., Robertson, D.E., 2018. A Bayesian modelling method for post-processing daily sub-seasonal to seasonal rainfall forecasts from global climate models and evaluation for 12 Australian catchments. Hydrol Earth Syst Sc, 22(2): 1615-1628. https://doi.org/10.5194/hess-22-1615-2018. Shen, Y., Hong, Z., Pan, Y., Yu, J., Maguire, L., 2018. China’s 1 km Merged Gauge, Radar and Satellite Experimental Precipitation Dataset. Remote Sensing, 10(2): 264. https://doi.org/10.3390/rs10020264. Shen, Y., Zhao, P., Pan, Y., Yu, J., 2014. A high spatiotemporal gauge-satellite merged precipitation analysis over China. J Geophys Res-Atmos, 119(6): 3063-3075. https://doi.org/10.1002/2013JD020686. Strazzo, S. et al., 2019. Application of a Hybrid Statistical–Dynamical System to Seasonal Prediction of North American Temperature and Precipitation. Mon Weather Rev, 147(2): 607-625. https://doi.org/10.1175/MWR-D-18-0156.1. Wang, D.S. et al., 2016. Evaluation of CMPA precipitation estimate in the evolution of typhoon-related storm rainfall in Guangdong, China. J Hydroinform, 18(6): 1055-1068. https://doi.org/10.2166/hydro.2016.241. Wang, Q.J., 2001. A Bayesian joint probability approach for flood record augmentation. Water Resour Res, 37(6): 1707-1712. https://doi.org/10.1029/2000WR900401.
35
Manuscript submitted to Journal of Hydrology
Wang, Q.J., Robertson, D.E., 2011. Multisite probabilistic forecasting of seasonal flows for
streams
with
zero
value
occurrences.
Water
Resour
Res,
47.
https://doi.org/10.1029/2010wr009333. Wang, Q.J., Robertson, D.E., Chiew, F.H.S., 2009. A Bayesian joint probability modeling approach for seasonal forecasting of streamflows at multiple sites. Water Resour Res, 45. https://doi.org/10.1029/2008wr007355. Wang, Q.J., Shrestha, D.L., Robertson, D.E., Pokhrel, P., 2012. A log-sinh transformation for data normalization and variance stabilization. Water Resour Res, 48(5): 1-7. https://doi.org/10.1029/2011wr010973. Ward, E., Buytaert, W., Peaver, L., Wheater, H., 2011. Evaluation of precipitation products over complex mountainous terrain: A water resources perspective. Advances
in
Water
Resources,
34(10):
1222-1231.
https://doi.org/10.1016/j.advwatres.2011.05.007. Zhang, G.P., Kline, D.M., 2007. Quarterly time-series forecasting with neural networks. IEEE
transactions
on
neural
networks,
18(6):
1800-1814.
https://doi.org/10.1109/TNN.2007.896859. Zhao, T., Fu, C., 2006. Comparison of products from ERA-40, NCEP-2, and CRU with station data for summer precipitation over China. Adv Atmos Sci, 23(4): 593-604. https://doi.org/10.1007/s00376-006-0593-1. Zhao, T., Wang, Q.J., Schepen, A., 2019a. A Bayesian modelling approach to forecasting short-term reference crop evapotranspiration from GCM outputs. Agricultural
and
Forest
Meteorology,
https://doi.org/10.1016/j.agrformet.2019.02.003.
36
269:
88-101.
Manuscript submitted to Journal of Hydrology
Zhao, T., Wang, Q.J., Schepen, A., Griffiths, M., 2019b. Ensemble forecasting of monthly and seasonal reference crop evapotranspiration based on global climate model
outputs.
Agricultural
and
forest
meteorology,
264:
114-124.
https://doi.org/10.1016/j.agrformet.2018.10.001. Zhao, T.T.G., Schepen, A., Wang, Q.J., 2016. Ensemble forecasting of sub-seasonal to seasonal streamflow by a Bayesian joint probability modelling approach. J Hydrol, 541: 839-849. https://doi.org/10.1016/j.jhydrol.2016.07.040.
37
S0022-1694(19)31090-X https://doi.org/10.1016/j.jhydrol.2019.124355 HYDROL 124355
To appear in:
Journal of Hydrology
Received Date: Revised Date: Accepted Date:
17 June 2019 13 October 2019 12 November 2019
Please cite this article as: Li, Y., Wang, Q.J., He, H., Wu, Z., Lu, G., A method to extend temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates, Journal of Hydrology (2019), doi: https://doi.org/ 10.1016/j.jhydrol.2019.124355
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 Published by Elsevier B.V.
Manuscript submitted to Journal of Hydrology
A method to extend temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates Yuan LI1, 2, Q. J. WANG2, Hai HE1, Zhiyong WU1, and Guihua LU1
1College
of Hydrology and Water Resources, Hohai University, Nanjing
210098, China 2Department
of Infrastructure Engineering, The University of Melbourne,
Parkville 3010, Australia Corresponding author: Zhiyong WU ([email protected])
Highlights:
A Bayesian joint probability method is presented to calibrate reanalysis estimates
The calibrated estimates are more accurate than the original reanalysis estimates
The uncertainty spread produced from the method is shown to be reliable
Key words: Bayesian joint probability method; precipitation estimates; temporal coverage; reanalysis dataset; uncertainty Declarations of interest: none 1
Manuscript submitted to Journal of Hydrology
Abstract Available high quality precipitation datasets generally cover only recent periods. To extend these datasets back in time is challenging due to historically lower rain gauge network density or poorer remote sensing. In this study, a Bayesian joint probability method is presented to extend the temporal coverage of high quality precipitation datasets by calibrating reanalysis estimates. Relationships between precipitation estimates from a high quality dataset and precipitation estimates from a reanalysis dataset are established by using data from the overlapping period of the two datasets. The relationships are then used to calibrate reanalysis estimates for the period when only reanalysis data is available. Rain gauge observations are also used to enhance the calibration. The method brings the reanalysis dataset to the same spatial resolution as the high quality dataset, corrects bias, and makes the reanalysis data more consistent with the high quality dataset. The calibrated estimates generated are presented in the form of ensembles to represent uncertainty. The method is applied to the Han River basin in China, using the CMPA high quality dataset and the ERA-interim reanalysis dataset. The method is shown to be highly effective in improving the quality of ERA-Interim. 1 Introduction In recent times, there has been significant improvement in the monitoring of precipitation in many parts of the world. In China, the rain gauge network has grown considerably over the last few decades. The automatic rain gauge network, set up in 2006 by the China Meteorological Administration (CMA), now consists of more than 57,000 stations (http://www.observation-cma.com/). In southeastern China, the mean distance
2
Manuscript submitted to Journal of Hydrology
between two nearest automatic rain gauges is approximately 9.6 km (Shen et al., 2018). While rain gauges provide the most accurate point measurement of precipitation, remote sensing can provide better spatial coverage in precipitation estimation across the globe. The NOAA-NCEP Climate Prediction Center (CPC) morphing technique (CMORPH) uses motion vectors derived from half-hourly interval geostationary satellite IR imagery to propagate relatively high quality precipitation estimates over the globe. The Global Precipitation Measurement (GPM) Mission is also expected to provide higher quality precipitation estimates resulting from better orbital inclination, radar frequencies and channels of the passive microwave imager (Huffman et al., 2014). High quality precipitation datasets have been developed by taking advantage of both high density rain gauge network and advanced remote sensing. The Integrated Multi-satellitE Retrievals for GPM (IMERG) is intended to intercalibrate, merge, and interpolate “all” satellite microwave precipitation estimates, together with infrared satellite estimates, and rain gauge analyses at high spatial (0.1°) and temporal (30 min) resolutions over the entire globe (Gaona et al., 2016). The China Merge Precipitation Analysis Hourly V1.0 product (CMPA) used an improved probability density function-optimal interpolation (PDF-OI) method to generate 0.1° and hourly precipitation estimates over China by merging the satellite precipitation estimates from the CMORPH and more than 30,000 automatic rain gauges (Shen et al., 2014). Owing to the high density rain gauge network and rigorous quality control, the CMPA dataset for precipitation is of relatively high quality. Wang et al. (2016) evaluated the CMPA dataset against a dense rain gauge network that was not incorporated in the interpolation process of the CMPA estimates. Their results indicated that the CMPA estimates capture the
3
Manuscript submitted to Journal of Hydrology
spatial patterns of typhoon-related rain depth, and are in agreement with the spatiotemporal evolutions of hourly gauge observations by correlation coefficients from 0.93 to 0.99. A drawback of these high quality precipitation datasets is that they are available only for relatively recent periods. The GPM-IMERG product only covers the period following the launch of the GPM satellite in 2014. The CMPA product is only available after 2008 when the rain gauge network density became relatively high. For many applications, precipitation datasets of much longer temporal coverage are required. Extending the high quality datasets further back in time is challenging due to low rain gauge network density or limitations of remote sensing. The Tropical Rainfall Measuring Mission (TRMM), which was launched in 1997, provided 3-hourly precipitation analysis at a spatial resolution of 0.25° for over 17 years (Huffman et al., 2007; Joyce et al., 2004). However, the accuracy of TRMM precipitation estimates left much to be desired. Huang et al. (2014) suggested that the TRMM multisatellite precipitation analysis (TMPA) products should be used with caution for extreme precipitation events. Reanalysis
data
is
an
alternative
source
of
precipitation
estimates.
Physically-based numerical models are able to generate precipitation estimates at various spatial and temporal resolutions, and the temporal coverage could be longer than both rain gauges and remote sensing. Studies have indicated that reanalysis datasets are able to give reasonable monthly averaged surface temperature, radiation fluxes, and wind speed (Lindsay et al., 2014). However, the accuracy of the reanalysis precipitation estimates still needs to be improved due to the imperfect description of small-scale physical processes of the numerical models (Leeuw et al., 2015; Nkiaka et al., 2017). 4
Manuscript submitted to Journal of Hydrology
One way of improving the reanalysis data is to calibrate the reanalysis data against rain gauge observations. However, this requires the availability of high density rain gauge networks, which are often not available. In this study, we introduce a method to calibrate the reanalysis data by taking advantage of the availability of high quality datasets (such as CMPA) for recent times. Using the period when both datasets are available, we establish relationships between precipitation estimates of the high quality dataset and precipitation estimates of the reanalysis dataset. The relationships can then be applied to calibrate the reanalysis estimates for the period when only reanalysis data is available. The method brings the reanalysis dataset to the same spatial resolution as the high quality dataset, corrects bias, and makes the reanalysis data more consistent with the high quality dataset. The calibrated estimates generated from our calibration method are presented in the form of ensembles to represent uncertainty. The particular method we use is based on the Bayesian joint probability (BJP) approach. The BJP approach, first developed for flood record augmentation, has found wide applications (Robertson and Wang, 2012; Schepen et al., 2018; Strazzo et al., 2019; Wang, 2001; Wang et al., 2009; Zhao et al., 2019a; Zhao et al., 2019b; Zhao et al., 2016). In BJP, data normalization is used to handle skewed variables (Wang et al., 2012), missing values are permitted, and zero values are treated as censored data. In addition, multiple data sources are allowed due to their multivariate formulation. In previous studies, the BJP method is always used to establish relationships between predictors and predictand for each grid cell. However, the intersite correlations have never been considered before. In this study, the BJP method is used as a spatial model to quantifying intersite correlations between grid cells and rain gauges from different locations.
5
Manuscript submitted to Journal of Hydrology
The efficacy of this calibration method for extending the temporal coverage of the CMPA precipitation dataset by calibrating the ERA-interim reanalysis over the Han River basin in China is evaluated in this study. The Han River basin is the water source of the middle route of the South-to-North Water Diversion Project (SNWDP) in China. Therefore, availability of accurate precipitation estimates of long temporal coverage over this region is important. The remainder of the paper is structured as follows. The study area and data sources are introduced in Section 2. An introduction of the BJP modelling method is provided and the evaluation methods are introduced in Section 3. Section 4 presents the evaluation results of the ensemble mean and ensemble distribution of calibrated precipitation estimates. Conclusions are drawn in Section 5.
2 Study area and data The Han River basin is located between 106-114°E and 29-34°N with an area of nearly 170,400 km2 (Figure 1). The basin is subject to sub-tropical monsoon climate, and its topography is highly complex. The elevation of western Han River basin is over 3000 m, while the eastern Han River basin is composed of recent alluvial plains (Li et al., 2017).
6
Manuscript submitted to Journal of Hydrology
Figure 1. Location of the Han River basin, rain gauges, ERA-interim grid, and the China Merge Precipitation Analysis Hourly V1.0 product (CMPA) grid. Figure 2 illustrates the spatial characteristics of daily precipitation for each month. It is clear that the precipitation shows high spatial and temporal variability over the Han River basin, especially in the wet months from May to September. In July, the daily mean precipitation is greater than 5 mm/day with a standard deviation of greater than 10 mm/day in the upper and lower sub basins over the period of 2008 – 2016.
7
Manuscript submitted to Journal of Hydrology
Figure 2. Spatial distribution of daily mean and daily standard deviation of CMPA estimates over the period of 2008 – 2016 for each month (mm/day). Both the ERA-interim and the CMPA datasets are used in this study. A brief introduction of these two datasets has been presented in previous section. In addition, rain gauge observations of long temporal coverage are also used as inputs to bring additional information to the calibration. The rain gauge observations are extracted from the “Daily Surface Climate Variables of China” dataset (SURF_CLI_CHN_MUL_DAY V3.0). It was released by the China Meteorological Administration (CMA) with high quality control. Out of the 756 long term national rain gauges, 48 rain gauges are located in or around the Han River basin and are used in this study (Figure 1). The daily ERA-interim reanalysis precipitation (hereafter Original ERA_daily) is accumulated from two forecasts (at 00:00 and 12:00) of 12 hour precipitation each. The CMPA daily precipitation (hereafter CMPA_daily) is accumulated from hourly values. The
8
Manuscript submitted to Journal of Hydrology
rain gauge observations (hereafter GAUGE_daily) are already daily. Details of the final datasets used in this study are summarized in Table 1. Table 1. A summary of daily precipitation datasets used in this study, and the role of each dataset in the BJP method. Dataset
Spatial resolution
Temporal coverage
Role in the BJP
CMPA_daily
0.1°
2008 - 2016
Predictand
Original ERA_daily
0.75°
1979 - 2016
Predictor
GAUGE_daily
~7000 km2/gauge
1979 - 2016
Predictor
3 Methods 3.1 The BJP method 3.1.1 Model formulation In this study, the BJP method is applied to calibrate the ERA_daily using the CMPA_daily. The CMPA_daily is treated as the most accurate owing to the high rain gauge network density and advanced remote sensing technique. We denote the value of the given 0.1° CMPA_daily grid cell as the predictand Y, the corresponding value of the nearest ERA_daily as the predictor 𝑋1. To enhance the calibration, rain gauge observations are also used for calibration as well. The distances from each 0.1° CMPA_daily grid cell to the rain gauges are presented in Figure 3. The mean distance between CMPA grid cells and the corresponding nearest gauge is approximately 30 km, while the corresponding second nearest rain gauge is nearly 60 km. Thus, only the two
9
Manuscript submitted to Journal of Hydrology
nearest rain gauges are selected as predictors 𝑋2 and 𝑋3 for the given CMPA_daily grid cell here.
Figure 3. The inter-distances between CMPA grid cells and rain gauges over the Han River basin. The probability distribution of daily precipitation is always highly skewed, which makes it difficult for parameter inference. To over this problem, we transform both the predictand Y and the predictors 𝑋𝑇 = [𝑋1 𝑋2 𝑋3] to normal distributions first. The Box-cox transformation has been widely used for transformation. However, the one-parameter Box-Cox transformation is not always able to achieve variance stabilization, while the 10
Manuscript submitted to Journal of Hydrology
two-parameter Box-Cox transformation leads to unrealistically large uncertainty. The log-sinh transformation is able to stabilize the variance and normalize the variable. It has been proved to outperform both the one-parameter and two-parameter Box-Cox transformations (Wang et al., 2012). This transformation has also been tested on a wide range of hydrological data, including hydrological extremes (Bennett et al., 2014a; Bennett et al., 2014b). In this study, the predictors 𝑋𝑇 (the matrix of the three datasets) are normalized to 𝑈𝑇 = [𝑈1 𝑈2 𝑈3], and the predictand Y is normalized to 𝑉 using the log-sinh transformation, 𝑢𝑖 = 𝛽𝑋―1 ln (sinh (𝛼𝑋𝑖 + 𝛽𝑋𝑖𝑥𝑖)) 𝑖
(1)
𝑣 = 𝛽𝑌―1ln (sinh (𝛼𝑌 + 𝛽𝑌𝑦))
(2)
where 𝛼𝑋,𝑖, 𝛽𝑋,𝑖, 𝛼𝑌, and 𝛽𝑌 are unknown transformation parameters for predictors and predictand, 𝑖 = 1,⋯, 3. The matrix 𝑍𝑇 = [𝑈1 𝑈2 𝑈3 𝑉] is then assumed to follow a multivariate normal distribution. 𝑍~𝐍(𝝁, ∑)
(3)
where 𝝁 and ∑ are the mean vector and covariance matrices to be estimated. 𝝁𝑻 = [𝜇𝑈1 𝜇𝑈2 𝜇𝑈3 𝜇𝑉]
(4)
∑ = 𝝈𝑹𝝈𝑻
(5)
where 𝝈 and 𝑹 are the standard deviation vector and correlation coefficient matrix, respectively: 𝝈𝑻 = [𝜎𝑈1 𝜎𝑈2 𝜎𝑈3 𝜎𝑉] 11
(6)
Manuscript submitted to Journal of Hydrology
[
1 𝑟𝑈1 𝑈2 𝑟𝑈1 𝑈3 𝑟𝑈1 𝑉 𝑟𝑈2 𝑈1 1 𝑟𝑈2 𝑈3 𝑟𝑈2 𝑉 𝐑= 𝑟 1 𝑟𝑈3 𝑉 𝑈3 𝑈1 𝑟𝑈3 𝑈2 𝑟𝑉 𝑈1 𝑟𝑉 𝑈2 𝑟𝑉 𝑈3 1
]
(7)
Note that the correlation coefficient matrix 𝐑 is symmetric. Thus, the joint distribution of predictand and predictors has a total of 14 unknown parameters. Zero values occasionally occur in dry months for both predictors and predictand over the Han River basin. This always leads to a mixed discrete-continuous multivariate distribution, which is difficult to formulate. In this study, we follow Wang and Robertson (2011) and treat zero values as censored data with unknown values below or equal to zero. In this way, the predictors and the predictand are considered to follow continuous multivariate distributions. This treatment allows us to simplify the problem of mixed discrete-continuous multivariate distribution, and to deal with zero values. 3.1.2 Parameter inference The transformation parameters and the joint distribution parameters are estimated when the ERA_daily, GAUGE_daily, and CMPA_daily are concurrently available. The choice of training data is still debated as parameter inference is affected by seasonality effects. Pooling all data from different months is not suitable. One way to solve this problem is to pre-deseasonalize data before parameter inference. Studies have proved that deseasonalizing data could potential improve forecast skills when using neural network forecasting models (Nelson et al., 1999; Zhang and Kline, 2007). A more direct approach is to pool data with similar seasonality. Hamill (2007) suggested that including days from the same season would expand the training period, and would likely improve the performance of the model. As shown in Fig. 2, the characteristics of daily 12
Manuscript submitted to Journal of Hydrology
precipitation vary from month to month. Thus, the daily values of each month are pooled together for each variable. There are approximately 270 days in total that are selected to make parameter inferences for a given grid cell of CMPA_daily in each month. Given
a
data
series 𝐃(𝑋𝑖)
of
each
predictor 𝑋𝑖, 𝐃(𝑋𝑖) =
{𝑥𝑖(t), t = 1, ⋯, n;i = 1, 2 , 3}, the transformation parameters of 𝛼𝑋,𝑖 and 𝛽𝑋,𝑖 are estimated using a Bayesian maximum a posteriori (MAP) solution. The resulting parameter values of 𝛼𝑋,𝑖 and 𝛽𝑋,𝑖 are taken as point estimates. Parameter uncertainty is not considered for these parameters. The transformation parameters of 𝛼𝑌 and 𝛽𝑌 are estimated using the same strategy as the predictors, given a data series 𝐃(𝑌) = {𝑦(t), t = 1, ⋯, n}. The joint distribution parameters are then estimated using transformed data series
of
both
predictors
and
predictand,
𝐃(U1,𝑈2,𝑈3,𝑉) =
{(𝑢1(t), 𝑢2(t),𝑢3(t), 𝑣(t)), t = 1, ⋯, n}. Here, we denote the 14 unknown parameters of the joint distribution as 𝜽 = {𝝁, ∑}. The posterior distribution of 𝜽 is estimated using a Bayesian framework, 𝑝(𝜽|𝐃(𝑈1,𝑈2,𝑈3,𝑉)) ∝ 𝑝(𝜽)𝑝(𝐃(𝑈1,𝑈2,𝑈3,𝑉)|𝜽)
(8)
where 𝑝(𝜽) is the prior distribution of parameters, and 𝑝(𝐃(𝑈1,𝑈2,𝑈3,𝑉)|𝜽) is the likelihood. In this study, non-informative distributions are used as priors of parameters 𝜽, 𝑝(𝜽) ∝ 1. As the posterior distributions of parameters 𝜽 are not standard distributions, analytical integration is not appropriate. To overcome this problem, we use the Markov chain Monte Carlo technique to draw a sample of parameter values which represent the joint posterior distribution of parameters 𝜽. Specifically, we use the Gibbs Sampling 13
Manuscript submitted to Journal of Hydrology
algorithm to draw (1000) samples from the conditional posterior distributions. To diminish the influence of the starting values, we discard the first 5000 iterations in Markov chain simulation as burn-in period. We also set the length of the Markov chain as 30, 000, which is proportional to the number of variables to obtain sufficient number of representative samples of the posterior distribution (Zhao et al., 2016). A visual inspection of the sampled parameter chains is made to monitor the convergence of the sample (Cowles and Carlin, 1996). 3.1.3 Model use for calibration Once the parameters are sampled, the BJP model can be used to calibrate the ERA_daily when the CMPA_daily is not available. Given a new set of predictor values, 𝐃(𝑋𝑖∗ ) = {𝑥𝑖(t ∗ ), t ∗ = 1, ⋯, n}, transformed values 𝐃(𝑈𝑖∗ ) = {𝑢𝑖(t ∗ ), t ∗ = 1, ⋯, n} can be found using the estimated transformation parameters 𝛼𝑋,𝑖 and 𝛽𝑋,𝑖. The posterior predictive distribution of 𝑣(t ∗ ) is given by 𝑓𝑉(𝑣(t ∗ )) = ∫𝑝(𝑣(t ∗ )|𝑢1(t ∗ ), 𝑢2(t ∗ ),𝑢3(t ∗ ),𝜽)𝑝(𝜽)𝑝(𝐃(𝑈1,𝑈2,𝑈3,𝑉)|𝜽)d𝜽
(9)
Again, the Gibbs Sampling algorithm is used to obtain samples of 𝑓𝑉(𝑣(t ∗ )) by giving each of the 1000 sets of parameter values 𝜽. The samples are then back-transformed to produce a calibrated precipitation estimates of 𝑦(t ∗ ) using parameters 𝛼𝑌 and 𝛽𝑌. Finally, we develop the BJP model for each grid cell of the CMPA_daily and each month separately.
14
Manuscript submitted to Journal of Hydrology
3.2 Evaluation Methods 3.2.1 Cross validation To assess the performance of the BJP calibration method, a ten-fold cross validation scheme is used here. The daily values of both the predictors and the predictand in each month from 2008 to 2016 are pooled together. Both the predictors and the predictand are then partitioned into ten equal-sized subsamples. Of the ten subsamples, nine subsamples of approximately 240 daily values are used for training, while the remaining one single subsample of nearly 30 daily values are used for validation. This process is then repeated ten times, with each of the ten subsamples used exactly once as the validation data. The independently calibrated ERA daily estimates (hereafter Calibrated ERA_daily) for all the 10 sub-samples are pooled together and compared against the CMPA_daily estimates. 3.2.2 Evaluation of ensemble mean The ensemble mean of Calibrated ERA_daily for each 0.1° CMPA grid cell is calculated and evaluated for absolute value of Relative Bias (RB), root mean square error (RMSE), and Nash–Sutcliffe coefficient of efficiency (NSE): ∑𝑡 = 1|𝑦𝑡∗ ― 𝑜𝑡| 𝑇
RB =
𝑇
∑𝑡 = 1 𝑜𝑡
1 𝑇 ∗ ∑ 𝑇 𝑡 = 1(𝑦𝑡
RMSE =
2
― 𝑜𝑡)
(11)
2
∑𝑡 = 1(𝑦𝑡∗ ― 𝑜𝑡) 𝑇
NSE = 1 ―
(10)
× 100%
𝑇
∑𝑡 = 1(𝑜 ― 𝑜𝑡)2 15
(12)
Manuscript submitted to Journal of Hydrology
where 𝑦𝑡∗ is the ensemble mean of the Calibrated ERA_daily precipitation for day 𝑡, 𝑡 = 1,2,⋯,𝑇; 𝑜𝑡 is the CMPA_daily precipitation, which is treated as the most accurate precipitation estimates in this study; 𝑜 is the averaged CMPA_daily precipitation over the period of analysis. The RB is used to evaluate the agreement between the averaged value of ensemble mean of Calibrated ERA_daily and the CMPA_daily, while the RMSE is used to evaluate the average error magnitude. The NSE is selected to evaluate the relative magnitude of the residual variance of the ensemble mean of Calibrated ERA_daily compared to the CMPA_daily variance. 3.2.3 Evaluation of ensemble distribution The continuous ranked probability score (Matheson and Winkler, 1976), and the α-index (Renard et al., 2010) are used to evaluate the accuracy and the reliability of the fit of distribution in this study. The continuous ranked probability score (CRPS) is used to provide an overall evaluation of the Calibrated ERA_daily for each 0.1° CMPA grid cell: CRPS 1 = 𝑇
𝑇
∑∫[𝐹(𝑦
∗ 𝑡
) ― 𝐻(𝑦𝑡∗ ― 𝑜𝑡)]2𝑑𝑦
where 𝑦𝑡∗ is
𝑡=1
(13) the daily precipitation to be estimated for day t; 𝐹() is the cumulative distribution function of the calibrated ERA_daily ensemble, as represented by the sampled (1000) values; and 𝐻() is the Heaviside step function defined as
16
Manuscript submitted to Journal of Hydrology
𝐻(𝑦𝑡∗ ― 𝑜𝑡) =
{
0 1
𝑦𝑡∗ < 𝑜𝑡 𝑦𝑡∗ ≥ 𝑜𝑡
(14)
where 𝑜𝑡 is the CMPA_daily precipitation, which is treated as the actual precipitation. When 𝐹() is replaced by a step distribution based on a deterministic value (such as the Original ERA_daily), the CRPS is reduced to the mean absolute error (MAE). A CRPS skill score is then calculated by comparing the CRPS of the Calibrated ERA_daily with the CRPS of a reference distribution for estimating precipitation:
CRPSSS =
CRPSREF ― CRPS CRPSREF
× 100%
(15)
In this study, a cross-validated climatology of CMPA_daily is used as the reference. A skill score of 100% indicates that the calibrated ERA_daily ensemble is the same as the CMPA_daily observation, whereas a skill score of 0% suggests that the calibrated ERA_daily ensemble shows no improvement over the cross-validated climatology. A negative skill score means that the calibrated ERA_daily ensemble is inferior to the cross-validated climatology. The reliability of Calibrated ERA_daily is evaluated by the probability integral transforms (PITs) of CMPA_daily. The PIT of 𝑜𝑡 is given by 𝜋𝑡 = 𝐹(𝑜𝑡)
(16)
where F() is the cumulative distribution function of the calibrated ERA_daily ensemble. When 𝑜𝑡 is equal to zero, the precise value for 𝜋𝑡 cannot be found from equation (16). Instead, a pseudo-PIT value is sampled from a uniform distribution with a range of [0, 𝜋𝑡], and subsequently using it with other real 𝜋𝑡 values to construct the PIT plots (Schepen et
17
Manuscript submitted to Journal of Hydrology
al., 2018; Wang and Robertson, 2011). If the Calibrated ERA_daily ensemble is reliable, 𝜋𝑡 should be uniformly distributed. The 𝜋𝑡 values are then summarized into an α-index:
α = 1.0 ―
2 𝑇
∑
𝑇
|𝜋
∗ 𝑡
𝑡=1
―
|
𝑡 𝑇+1
(17) where 𝜋𝑡∗ is the
sorted 𝜋𝑡 in increasing order. The α-index ranges from 0 to 1, and represents the total deviation of 𝜋𝑡∗ from corresponding uniform quantile. A higher α-index indicates higher reliability. 3.3 Further Analysis In this study, an important assumption is made that the CMPA precipitation estimates are of high quality, as such they can be used to calibrate the ERA reanalysis estimates. It is not right to directly verify assumption as follows by comparing CMPA precipitation estimates with observations at gauges, as the observations are likely to have been used in deriving the CMPA precipitation estimates. Here we indirectly evaluate the assumption as follows. We use the 2008-2016 data to establish models to calibrate ERA reanalysis estimates against the CMPA precipitation estimates only. Unlike the models used in Section 3.1, gauged observations are not used here.
The new calibration models are then used to calibrate the ERA
reanalysis estimates for 1979-2007. The calibrated ERA reanalysis estimates are compared with gauged observations. For the calibration to lead to improvement over the raw ERA reanalysis estimates, the CMPA precipitation estimates would have to be more accurate than the ERA reanalysis estimates.
18
Manuscript submitted to Journal of Hydrology
This further analysis is to ensure that the quality of the CMPA precipitation estimates are high enough that its use adds value to the calibration and does not adversely corrupt the ERA reanalysis estimates. In practical applications, the gauged observations are also used as predictors to derive the final estimates.
4 Results 4.1 Ensemble mean The RB and RMSE of the ensemble mean of Calibrated ERA_daily for each month are plotted and compared with the Original ERA_daily in Figure 4 and Figure 5. The RB and RMSE of the Original ERA_daily are greatly reduced after the BJP calibration, especially in wet months in summer. In July, the RB of the Original ERA_daily is always greater than 150%, with RMSE values greater than 10 mm/day. In contrast, the Calibrated ERA_daily shows much higher accuracy in terms of RB values of lower than 50% and RMSE values of lower than 8 mm/day in most CMPA_daily grid cells. The reduction of RB and RMSE values are also evident in other months. The RB values reduce from nearly 100% to less than 50% while the RMSE values reduce from approximately 5 mm/day to less than 2 mm/day in months in spring and autumn.
19
Manuscript submitted to Journal of Hydrology
Figure 4. RB of Original ERA_daily and RB of ensemble mean of Calibrated ERA_daily compared to CMPA_daily for each month (%).
Figure 5. RMSE of Original ERA_daily and RMSE of ensemble mean of Calibrated ERA_daily compared to CMPA_daily for each month (mm/day). 20
Manuscript submitted to Journal of Hydrology
Figure 6 compares the NSE of the Original ERA_daily and the ensemble mean of Calibrated ERA_daily. The results also suggest that the accuracy of the Original ERA_daily has been greatly improved after the BJP calibration. In July, the NSE values of the Original ERA_daily are always less than 0.2. In contrast, the NSE values of the ensemble mean of calibrated ERA_daily are greater than 0.5 for most grid cells. The NSE
values increase from less than 0.4 to more than 0.6 after the BJP calibration in
months in spring and autumn as well. However, we should also notice that the NSE values are still low from December to Februrary owing to the relatively lower precipitation in these months shown in Fig. 2.
Figure 6. NSE of Original ERA_daily and NSE of ensemble mean of Calibrated ERA_daily compared to CMPA_daily for each month.
21
Manuscript submitted to Journal of Hydrology
4.2 Ensemble distribution The CRPS skill score of Original ERA_daily and the CRPS skill score of Calibrated ERA_daily are compared in Figure 7. The results suggest that the ensemble distribution of Calibrated ERA_daily is more accurate than the Original ERA_daily. The Calibrated ERA_daily usually exhibits positive skill in most areas in July. Whereas, negative skill is observed for Original ERA_daily. Improvement is observed in other months as well in terms of CRPS skill score.
Figure 7. CRPS skill score of Original ERA_daily and CRPS skill score of Calibrated ERA_daily compared to CMPA_daily for each month (%). The α-index of Calibrated ERA_daily for each month is presented in Figure 8. It is clear that the α-index is mostly in the range of 0.9 to 1.0 for all months, indicating that the Calibrated ERA_daily is of high reliability. The PIT diagram of a selected grid cell with lowest α-index value of 0.825 in February is also given in Figure 8. Although small 22
Manuscript submitted to Journal of Hydrology
departures from the diagonal line are observed, the PIT values for the selected grid cell lie mostly within the Kolmogorov 5% significance bands, suggesting that the uncertainty spread is reliable for all months.
Figure 8. α-index of Calibrated ERA_daily, and PIT uniform probability plot for the selected grid cell with lowest α-index value (1:1 solid line, theoretical uniform distribution; dashed lines, Kolmogorov 5% significance band; blue circles, PIT values of CMPA_daily precipitation; red circles, pseudo PIT values).
4.3 Calibrated ERA_daily precipitation for 1979 - 2007 The BJP models are subsequently established based on all available 2008 - 2016 data of 0.1° CMPA_daily precipitation, Original ERA_daily, and GAUGE_daily. The models are then used to produce Calibrated ERA_daily from the Original ERA_daily and GAUGE_daily for 1979 – 2007, for which 0.1° CMPA_daily is not available. 23
Manuscript submitted to Journal of Hydrology
Figure 9 compares the daily precipitation distribution of the Original ERA_daily, the GAUGE_daily, and the ensemble mean of Calibrated ERA_daily for an extreme precipitation event in October 1983. Beginning on 3 October 1983 and lasting to 5 October 1983, an extreme precipitation event occurred in the Han River basin. Similar to the example shown in Section 4.1, the Original ERA_daily is too low in spatial variability compared with the GAUGE_daily. The high precipitation intensity observed by the GAUGE_daily is greatly underestimated by the Original ERA_daily. The ensemble mean of Calibrated ERA_daily takes guidance from both the Original ERA_daily and GAUGE_daily. The spatial characteristics of the Original ERA_daily are broadly preserved, while high precipitation intensity of the GAUGE_daily is well captured by the ensemble mean of Calibrated ERA_daily.
24
Manuscript submitted to Journal of Hydrology
Figure 9. Daily precipitation distribution of the ensemble mean of Calibrated ERA_daily, Original ERA_daily, and GAUGE_daily for an extreme precipitation event from 0000 UTC 03 October 1983 to 0000 UTC 05 October 1983. 4.4 Further Analysis Recalling Section 3.3, the purpose of this further analysis is to ensure that the quality of the CMPA precipitation estimates is high enough that its use adds value to the calibration and does not adversely corrupt the ERA reanalysis estimates. We use the 2008-2016 data to establish models to calibrate ERA reanalysis estimates against the CMPA precipitation estimates only. Unlike the models used in Section 3.1, gauged
25
Manuscript submitted to Journal of Hydrology
observations are not used here. The new calibration models are then used to calibrate the ERA reanalysis estimates for 1979-2007. We should notice that the direct comparison between gauged observations and gridded estimates should be avoided due to the different representativeness of these two datasets. Currently, there are three indirect approaches to compare the gauged observations and gridded estimates. Zhao and Fu (2006) suggested using the kriging technique to upscale gauged observations onto the grids and then compare with the gridded estimates. However, the upscaling approach will introduce large errors and uncertainties when the gauges are sparsely distributed. Another approach is downscaling the gridded estimates at each station site using the bilinear interpolation method (Gao and Liu, 2012). Nevertheless, the downscaling approach may smooth extremes from the gridded data (Accadia et al., 2003). The third approach is to compare the gridded estimates and averaged values of gauged observations falling within the corresponding grid cell. This approach has been proved to work well with precipitation estimates comparisons (Blacutt et al., 2015; Ward et al., 2011). However, we should notice that spatial sampling error exists owing to small spatial scale of rainfall systems compared to rain gauge network density (Krajewski et al., 2000; Morrissey et al., 1995). Therefore, it is practically impossible to identify an optimal method applicable for all areas (Duan et al., 2016). In this study, the third approach is adapted here to compare the Calibrated ERA_daily and the Original ERA_daily against the gauged observations. The Calibrated ERA_daily is at a spatial resolution of 0.1° as is the CMPA_daily. Thus, there is only one gauge falling within each 0.1° grid cell of Calibrated ERA_daily. We directly compare the Calibrated ERA_daily and the corresponding gauged 26
Manuscript submitted to Journal of Hydrology
observations here. For grid cells of the Original ERA_daily containing more than one gauge, the averaged values of the gauges are calculated. The results are then compared against the Original ERA_daily. The NSE of Original ERA_daily and ensemble mean of Calibrated ERA_daily for the period of 1979-2007 are then presented in Fig. 10. The NSE of the Calibrated ERA_daily is much more consistent among gauge stations than the Original ERA_daily. The latter can be very poor in performance in some locations and months. In most months, the median NSE of the Calibrated ERA_daily is better than the Original ERA_daily. Overall, the calibration has led to considerable improvement. This analysis indicates that the use of the CMPA precipitation estimates as a high quality dataset is reasonable at least for this study case.
27
Manuscript submitted to Journal of Hydrology
Figure 10. Box plots of NSE of Original ERA_daily and ensemble mean of Calibrated ERA_daily for the period of 1979-2007. The red lines are the 50th percentiles, top and bottom of each box are the 75th and 25th percentiles, and whiskers are the 95th and 5th percentiles. 5 Conclusions This paper presents a Bayesian joint probability (BJP) method to extend the temporal coverage of high quality precipitation datasets by calibrating reanalysis
28
Manuscript submitted to Journal of Hydrology
estimates. The key concept is as follows. Relationships between precipitation estimates from a high quality dataset and precipitation estimates from a reanalysis dataset can be established by using data from the overlapping period of the two datasets. The relationships can then be used to calibrate reanalysis estimates for the period when only reanalysis data is available. Additionally, the method makes use of long records from rain gauges to enhance the calibration. The method brings the reanalysis dataset to the same spatial resolution of the high quality dataset, corrects bias, and makes the reanalysis data more consistent with the high quality dataset. The calibrated estimates generated are presented in the form of ensembles to represent uncertainty. In this study, the BJP method is used to extend the temporal coverage of a high quality precipitation dataset – CMPA (2008 - 2016) by calibrating the ERA-interim reanalysis precipitation estimates (1979 - 2016) over the Han River basin. Records from a sparsely distributed rain gauge network (1979 - 2016) are also used to enhance the calibration. To evaluate the accuracy and reliability of calibrated precipitation estimates, a ten-fold cross-validation scheme is applied to the overlapping period (2008 - 2016) of the three datasets. The original reanalysis estimates and the calibrated reanalysis estimates are evaluated against the precipitation estimates from the high quality dataset. Further analysis is performed to ensure that the quality of the CMPA precipitation estimates is high enough that its use adds value to the calibration and does not adversely corrupt the ERA reanalysis estimates. The calibrated reanalysis estimates are shown to be much more accurate than the original reanalysis estimates, in terms of RB, RMSE and CRPS skill score. The 29
Manuscript submitted to Journal of Hydrology
uncertainty spread produced from the calibration is shown to be reliable. The CMPA dataset is also proved to be effective to calibration. We conclude that the method introduced in this study is highly effective in improving the quality of ERA-Interim. The fifth generation of reanalysis dataset from ECMWF (ERA-5) will be published with a horizontal resolution of nearly 30 km from 1950 to the present by early 2019 (Brönnimann et al., 2018). The application of the calibration method developed in this study to the new reanalysis products will be straightforward.
Acknowledgments We would like to thank Geoff Pegram and another two anonymous reviewers for their insightful comments. This work was funded by the National Key R&D Program of China (Grant No. 2018 YFC0407701, 2017YFC1502403), the Fundamental Research Funds for the Central Universities (Grants No. 2019B10214, 2017B681X14), the National Natural Science Foundation of China (Grants No. 51579065, 51779071), the Postgraduate Research & Practice Innovation Program of Jiangsu Province (Grant No. KYCX17_0413), and the China Scholarship Council scholarship. We also thank Dr Tongtiegang Zhao for providing technical support to the work. References Accadia, C., Mariani, S., Casaioli, M., Lavagnini, A., Speranza, A., 2003. Sensitivity of precipitation forecast skill scores to bilinear interpolation and a simple nearest-neighbor average method on high-resolution verification grids. Weather
30
Manuscript submitted to Journal of Hydrology
and
forecasting,
18(5):
918-932.
https://doi.org/10.1175/1520-0434(2003)018<0918:SOPFSS>2.0.CO;2. Bennett, J.C. et al., 2014a. A System for Continuous Hydrological Ensemble Forecasting (SCHEF)
to
lead
times
of
9
days.
J
Hydrol,
519:
2832-2846.
https://doi.org/10.5194/nhess-14-219-2014. Bennett, J.C., Wang, Q., Pokhrel, P., Robertson, D.E., 2014b. The challenge of forecasting high streamflows 1–3 months in advance with lagged climate indices in southeast Australia. Natural Hazards and Earth System Sciences, 14(2): 219-233. https://doi.org/10.1016/j.jhydrol.2014.08.010. Brönnimann, S. et al., 2018. Observations for reanalyses. Bulletin of the American Meteorological
Society,
99(9):
1851-1866.
https://doi.org/10.1175/BAMS-D-17-0229.1. Blacutt, L.A., Herdies, D.L., de Gonçalves, L.G.G., Vila, D.A., Andrade, M., 2015. Precipitation comparison for the CFSR, MERRA, TRMM3B42 and Combined Scheme
datasets
in
Bolivia.
Atmos
Res,
163:
117-131.
https://doi.org/10.1016/j.atmosres.2015.02.002. China Merge Precipitation Analysis Hourly V1.0 product (2018). China Meteorological Administration,
retrieved
from
http://data.cma.cn/data/detail/dataCode/SEVP_CLI_CHN_MERGE_CMP_PRE_ HOUR_GRID_0.10/. Cowles, M., and Carlin B. P., 1996. Markov chain Monte Carlo convergence diagnostics: A
comparative
review,
J.
Am.
https://doi.org/10.2307/2291683.
31
Stat.
Assoc.,
91,
883–904.
Manuscript submitted to Journal of Hydrology
Daily Surface Climate Variables of China (2018). China Meteorological Administration, retrieved
from
http://data.cma.cn/data/detail/dataCode/SURF_CLI_CHN_MUL_DAY_V3.0/. Duan, Z., Liu, J., Tuo, Y., Chiogna, G., Disse, M., 2016. Evaluation of eight high spatial resolution gridded precipitation products in Adige Basin (Italy) at multiple temporal and spatial scales. Science of the Total Environment, 573: 1536-1553. https://doi.org/10.1016/j.scitotenv.2016.08.213. Gao, Y., Liu, M., 2012. Evaluation of high-resolution satellite precipitation products using rain gauge observations over the Tibetan Plateau. Hydrology & Earth System Sciences Discussions, 9(8). https://doi.org/10.5194/hess-17-837-2013. Gaona, M.F.R., Overeem, A., Leijnse, H., Uijlenhoet, R., 2016. First-Year Evaluation of GPM Rainfall over the Netherlands: IMERG Day 1 Final Run (VO3D). J Hydrometeorol, 17(11): 2799-2814. https://doi.org/10.1175/JHM-D-16-0087.1. Hamill, T. M. (2007), Comments on “Calibrated Surface Temperature Forecasts from the Canadian Ensemble Prediction System Using Bayesian Model Averaging”, Monthly
Weather
Review,
135(12),
4226-4230,
https://doi.org/10.1175/2007mwr1963.1. Huang, Y. et al., 2014. Evaluation of Version-7 TRMM Multi-Satellite Precipitation Analysis Product during the Beijing Extreme Heavy Rainfall Event of 21 July 2012. Water, 6(1): 32-44. https://doi.org/10.3390/w6010032. Huffman, G.H. et al., 2014. NASA Global Precipitation Measurement (GPM) Integrated Multi-satellitE Retrievals for GPM (IMERG). Algorithm theoretical basis document, version, 4: 30.
32
Manuscript submitted to Journal of Hydrology
Huffman, G.J. et al., 2007. The TRMM multisatellite precipitation analysis (TMPA): Quasi-global, multiyear, combined-sensor precipitation estimates at fine scales. J Hydrometeorol, 8(1): 38-55. https://doi.org/10.1175/JHM560.1. Joyce, R.J., Janowiak, J.E., Arkin, P.A., Xie, P.P., 2004. CMORPH: A method that produces global precipitation estimates from passive microwave and infrared data at high spatial and temporal resolution. J Hydrometeorol, 5(3): 487-503. https://doi.org/10.1175/1525-7541(2004)005<0487:CAMTPG>2.0.CO;2. Krajewski, W.F., Ciach, G.J., McCollum, J.R., Bacotiu, C., 2000. Initial validation of the Global Precipitation Climatology Project monthly rainfall over the United States. Journal
of
Applied
Meteorology,
39(7):
1071-1086.
https://doi.org/10.1175/1520-0450(2000)039<1071:IVOTGP>2.0.CO;2. Leeuw, J.D., Methven, J., Blackburn, M., 2015. Evaluation of ERA-Interim reanalysis precipitation products using England and Wales observations. Q J Roy Meteor Soc, 141(688): 798-806. https://doi.org/10.1002/qj.2395. Li, Y., Lu, G.H., Wu, Z.Y., He, H., He, J., 2017. High-Resolution Dynamical Downscaling of Seasonal Precipitation Forecasts for the Hanjiang Basin in China Using the Weather Research and Forecasting Model. J Appl Meteorol Clim, 56(5): 1515-1536. https://doi.org/10.1175/Jamc-D-16-0268.1. Lindsay, R., Wensnahan, M., Schweiger, A., Zhang, J., 2014. Evaluation of Seven Different Atmospheric Reanalysis Products in the Arctic. J Climate, 27(7): 2588-2606. https://doi.org/10.1175/JCLI-D-13-00014.1.
33
Manuscript submitted to Journal of Hydrology
Matheson, J.E., Winkler, R.L., 1976. Scoring rules for continuous probability distributions.
Management
science,
22(10):
1087-1096.
https://doi.org/10.1287/mnsc.22.10.1087. Morrissey, M.L., Maliekal, J.A., Greene, J.S., Wang, J., 1995. The uncertainty of simple spatial averages using rain gauge networks. Water Resources Research, 31(8): 2011-2017. https://doi.org/10.1029/95WR01232. Nelson, M., Hill, T., Remus, W., O'Connor, M., 1999. Time series forecasting using neural networks: Should the data be deseasonalized first? Journal of forecasting, 18(5):
359-367.
https://doi.org/10.1002/(SICI)1099-131X(199909)18:5<359::AID-FOR746>3.0.C O;2-P. Nkiaka, E., Nawaz, N.R., Lovett, J.C., 2017. Evaluating global reanalysis precipitation datasets with rain gauge measurements in the Sudano‐Sahel region: case study of the Logone catchment, Lake Chad Basin. Meteorol Appl, 24(1): 9-18. https://doi.org/10.1002/met.1600. Poli, P. et al., 2016. ERA-20C: An Atmospheric Reanalysis of the Twentieth Century. J Climate, 29(11): 4083-4097. https://doi.org/10.1175/Jcli-D-15-0556.1. Renard, B., Kavetski, D., Kuczera, G., Thyer, M., Franks, S.W., 2010. Understanding predictive uncertainty in hydrologic modeling: The challenge of identifying input and
structural
errors.
Water
https://doi.org/10.1029/2009wr008328.
34
Resour
Res,
46.
Manuscript submitted to Journal of Hydrology
Robertson, D.E., Wang, Q.J., 2012. Bayesian Approach to Predictor Selection for Seasonal
Strearnflow
Forecasting.
J
Hydrometeorol,
13(1):
155-171.
https://doi.org/10.1175/Jhm-D-10-05009.1. Schepen, A., Zhao, T.T.G., Wang, Q.J., Robertson, D.E., 2018. A Bayesian modelling method for post-processing daily sub-seasonal to seasonal rainfall forecasts from global climate models and evaluation for 12 Australian catchments. Hydrol Earth Syst Sc, 22(2): 1615-1628. https://doi.org/10.5194/hess-22-1615-2018. Shen, Y., Hong, Z., Pan, Y., Yu, J., Maguire, L., 2018. China’s 1 km Merged Gauge, Radar and Satellite Experimental Precipitation Dataset. Remote Sensing, 10(2): 264. https://doi.org/10.3390/rs10020264. Shen, Y., Zhao, P., Pan, Y., Yu, J., 2014. A high spatiotemporal gauge-satellite merged precipitation analysis over China. J Geophys Res-Atmos, 119(6): 3063-3075. https://doi.org/10.1002/2013JD020686. Strazzo, S. et al., 2019. Application of a Hybrid Statistical–Dynamical System to Seasonal Prediction of North American Temperature and Precipitation. Mon Weather Rev, 147(2): 607-625. https://doi.org/10.1175/MWR-D-18-0156.1. Wang, D.S. et al., 2016. Evaluation of CMPA precipitation estimate in the evolution of typhoon-related storm rainfall in Guangdong, China. J Hydroinform, 18(6): 1055-1068. https://doi.org/10.2166/hydro.2016.241. Wang, Q.J., 2001. A Bayesian joint probability approach for flood record augmentation. Water Resour Res, 37(6): 1707-1712. https://doi.org/10.1029/2000WR900401.
35
Manuscript submitted to Journal of Hydrology
Wang, Q.J., Robertson, D.E., 2011. Multisite probabilistic forecasting of seasonal flows for
streams
with
zero
value
occurrences.
Water
Resour
Res,
47.
https://doi.org/10.1029/2010wr009333. Wang, Q.J., Robertson, D.E., Chiew, F.H.S., 2009. A Bayesian joint probability modeling approach for seasonal forecasting of streamflows at multiple sites. Water Resour Res, 45. https://doi.org/10.1029/2008wr007355. Wang, Q.J., Shrestha, D.L., Robertson, D.E., Pokhrel, P., 2012. A log-sinh transformation for data normalization and variance stabilization. Water Resour Res, 48(5): 1-7. https://doi.org/10.1029/2011wr010973. Ward, E., Buytaert, W., Peaver, L., Wheater, H., 2011. Evaluation of precipitation products over complex mountainous terrain: A water resources perspective. Advances
in
Water
Resources,
34(10):
1222-1231.
https://doi.org/10.1016/j.advwatres.2011.05.007. Zhang, G.P., Kline, D.M., 2007. Quarterly time-series forecasting with neural networks. IEEE
transactions
on
neural
networks,
18(6):
1800-1814.
https://doi.org/10.1109/TNN.2007.896859. Zhao, T., Fu, C., 2006. Comparison of products from ERA-40, NCEP-2, and CRU with station data for summer precipitation over China. Adv Atmos Sci, 23(4): 593-604. https://doi.org/10.1007/s00376-006-0593-1. Zhao, T., Wang, Q.J., Schepen, A., 2019a. A Bayesian modelling approach to forecasting short-term reference crop evapotranspiration from GCM outputs. Agricultural
and
Forest
Meteorology,
https://doi.org/10.1016/j.agrformet.2019.02.003.
36
269:
88-101.
Manuscript submitted to Journal of Hydrology
Zhao, T., Wang, Q.J., Schepen, A., Griffiths, M., 2019b. Ensemble forecasting of monthly and seasonal reference crop evapotranspiration based on global climate model
outputs.
Agricultural
and
forest
meteorology,
264:
114-124.
https://doi.org/10.1016/j.agrformet.2018.10.001. Zhao, T.T.G., Schepen, A., Wang, Q.J., 2016. Ensemble forecasting of sub-seasonal to seasonal streamflow by a Bayesian joint probability modelling approach. J Hydrol, 541: 839-849. https://doi.org/10.1016/j.jhydrol.2016.07.040.
37