Bayesian Hierarchical Model Uncertainty Quantification for Future Hydroclimate Projections in Southern Hills-Gulf Region , USA

The study investigates the hierarchical uncertainty of multi-ensemble hydroclimate projections for the Southern Hills-Gulf region, USA, considering emission pathways and a global climate model (GCM) as two main sources of uncertainty. Forty projections of downscaled daily air temperature and precipitation from 2010 to 2099 under four emission pathways and ten CMIP5 GCMs are adopted for hydroclimate modeling via the HELP3 hydrologic model. This study focuses on evapotranspiration (ET), surface runoff, and groundwater recharge projections in this century. Climate projection uncertainty is characterized by the hierarchical Bayesian model averaging (HBMA) method, which segregates emission pathway uncertainty and climate model uncertainty. HBMA is able to derive ensemble means and standard deviations, arising from individual uncertainty sources, for ET, runoff, and recharge. The model results show that future recharge in the Southern Hills-Gulf region is more sensitive to different climate projections and exhibits higher variability than ET and runoff. Overall, ET is likely to increase and runoff is likely to decrease in this century given the current emission path scenarios. Runoff are predicted to have an 18% to 20% decrease and ET is predicted to have around a 3% increase throughout the century. Groundwater recharge is likely to increase in this century with a decreasing trend. Recharge would increase about 13% in the early century and will have only a 3% increase in the late century. All hydrological projections have increasing uncertainty towards the end of the century. The HBMA result suggests that the GCM uncertainty dominates the overall hydrological projection uncertainty in the early century and the mid-century. The emission pathway uncertainty becomes important in the late century.


Introduction
Uncertainties in hydroclimate projections provide necessary information to support the planning and management of future water recourses, but those uncertainties have not been sufficiently addressed in the IPCC (Intergovernmental Panel on Climate Change) Fifth Assessment [1].It is imperative to quantify the uncertainty associated with climate projections in order for detection and attribution as well as for adaptation and mitigation strategies [2].The anthropogenic forcing impacts are investigated through future climate projections using general circulation models or global climate models (GCMs).There are different types uncertainties associated with GCM projections, which are mainly from model representations, scenarios, and internal model variability [3,4].In addition to the various sources of uncertainty inherent in GCMs, additional uncertainties are also introduced during downscaling of raw GCM output for hydrologic simulation.Uncertainties resulting from GCMs and greenhouse gas (GHG) emissions generally have been given more attention in the uncertainty assessment of hydrological impacts of climate change [5].There has been a growing interest in evaluating total uncertainty in future hydrologic projections using combinations of emission scenarios, climate models, downscaling methods, and hydrologic models [6][7][8][9][10][11][12][13][14][15].
Assessment using a single model is prone to underestimate uncertainty.There is no single model proclaimed as the best model due to the inaccuracy of each individual GCM.Hence, future climate projections should take into account multiple GCMs [16,17].Compared to a single GCM approach, a multi-model ensemble of GCMs often increases reliability of information, and its average shows more consistency with climate observations [16,18].The IPCC suggests using multi-model ensembles for evaluating climate projections [19].In addition, weighting models based on some indices of performance might further increase the reliability of projections.Model weights can be defined based on model performance in simulating historical climate [20][21][22][23][24][25][26].
Moving forward, it is of great importance to better understand the influence of individual sources of uncertainty [43] as well as their relative contribution to the full hydroclimate uncertainty space.To address this issue, we utilize a hierarchical Bayesian model averaging (HBMA) method [44] that can systematically segregate sources of uncertainty in a hierarchical structure (known as the BMA tree).
We target hydroclimate modeling uncertainties arising from downscaled precipitation and air temperature (hereafter as temperature) projections and use hydroclimate modeling on the Southern Hills-Gulf region for illustration.This study uses the downscaled precipitation and temperature projections from References [45,46], where there are 21 Coupled Model Intercomparison Project phase 5 (CMIP5) GCMs developed by 16 organizations.Some downscaled data were derived by the same GCM with different initial conditions.Some GCMs were not run for the entire four representative concentration pathways (RCP2.6,RCP4.5, RCP6.0, and RCP8.5).This study selects 10 GCMs (see Table 1) that were run for all four RCPs with one initial condition.A HELP3 model [47] was built to simulate future surface runoff, evapotranspiration, and groundwater recharge based on the selected downscaled precipitation and temperature projections.
Table 1.Modeling centers/groups and the corresponding downscaled data of 10 GCMs under the four representative concentration pathways (RCP2.6,RCP4.5, RCP6.0, and RCP8.5) [45,46] for the Southern Hills-Gulf region.This research is structured as follows: Section 2 introduces the HBMA method; Section 3 describes the case study setup including study area, climate projections and performance evaluation, and hydrologic modeling; Section 4 discusses the results and findings; and Section 5 concludes the study.

Bayesian Model Averaging (BMA) Tree
To organize the multiple sources of uncertainty embedded in hydroclimate projections, a BMA tree [44] is adopted.Figure 1 shows a hierarchical structure of CMIP5 multi-model ensembles considering the choices of emission pathways and GCMs as the two main sources of uncertainty.Each source of uncertainty creates a number of mutually exclusive propositions (i.e., unique modeling choices).The future emission scenario is placed at the first level, which contains four propositions: RCP2.6,RCP4.5, RCP6.0, and RCP8.5.Under each emission pathway, different GCMs are placed at the second level.When other sources of uncertainty are considered (e.g., different downscaling approaches, hydrologic models, or model parameter sets), the uncertainty structure in Figure 1 can be further expanded.This research is structured as follows: Section 2 introduces the HBMA method; Section 3 describes the case study setup including study area, climate projections and performance evaluation, and hydrologic modeling; Section 4 discusses the results and findings; and Section 5 concludes the study.

Bayesian Model Averaging (BMA) Tree
To organize the multiple sources of uncertainty embedded in hydroclimate projections, a BMA tree [44] is adopted.Figure 1 shows a hierarchical structure of CMIP5 multi-model ensembles considering the choices of emission pathways and GCMs as the two main sources of uncertainty.Each source of uncertainty creates a number of mutually exclusive propositions (i.e., unique modeling choices).The future emission scenario is placed at the first level, which contains four propositions: RCP2.6,RCP4.5, RCP6.0, and RCP8.5.Under each emission pathway, different GCMs are placed at the second level.When other sources of uncertainty are considered (e.g., different downscaling approaches, hydrologic models, or model parameter sets), the uncertainty structure in Figure 1 can be further expanded.In the BMA tree, a parent model is a model at a vertex of a level, and a model one level below the parent model is a child model.The model at the apex of the BMA tree is called the hierarch model, which is a union of all sub-models.Each source of uncertainty corresponds to each level in the BMA tree and builds on top of other sources of uncertainty lower on the tree.The base level of the BMA tree depicts base models resulting from combinations of propositions at all different levels.Increasing the number of sources of uncertainty increases the number of levels of the BMA tree and the number of models.Through a BMA tree, the collective uncertainty from all child models can be hierarchically summarized into their parent models.This approach may enable the illustration of different sources of uncertainty in a clear and organized fashion.

Hierarchical Bayesian Model Averaging (HBMA)
Given p sources of uncertainty, base models, at level p are created, where the subscript ( ... ) p ij lm    locates a base model top down in the hierarchy (see Figure 1 for an example) from the first level to the second level and so forth to reach to base level p.For instance, the ith model at the first level of a BMA tree is denoted by

M at
level p-1, and so forth, until the hierarch level is reached [44].At the hierarch level (level zero), it is called "the hierarch BMA model," where predictions are based on all base models and thus ensemble averaging results are the same as using the BMA model in Hoeting et al. [27].
The posterior probability for a predicted quantity Δ at level n given observational data D and models

M
at level n + 1 is: In the BMA tree, a parent model is a model at a vertex of a level, and a model one level below the parent model is a child model.The model at the apex of the BMA tree is called the hierarch model, which is a union of all sub-models.Each source of uncertainty corresponds to each level in the BMA tree and builds on top of other sources of uncertainty lower on the tree.The base level of the BMA tree depicts base models resulting from combinations of propositions at all different levels.Increasing the number of sources of uncertainty increases the number of levels of the BMA tree and the number of models.Through a BMA tree, the collective uncertainty from all child models can be hierarchically summarized into their parent models.This approach may enable the illustration of different sources of uncertainty in a clear and organized fashion.

Hierarchical Bayesian Model Averaging (HBMA)
Given p sources of uncertainty, base models, M (ij...lm) ∈ M p at level p are created, where the subscript (ij...lm) p locates a base model top down in the hierarchy (see Figure 1 for an example) from the first level to the second level and so forth to reach to base level p.For instance, the ith model at the first level of a BMA tree is denoted by 2 ∈ M 2 , and so forth, until the base level p is reached.Furthermore, models M p−1 at level p − 1 are BMA models of their corresponding child base models M p at level p, models M p−2 at level p − 2 are BMA models of their corresponding child BMA models M p−1 at level p − 1, and so forth, until the hierarch level is reached [44].At the hierarch level (level zero), it is called "the hierarch BMA model", where predictions are based on all base models and thus ensemble averaging results are the same as using the BMA model in Hoeting et al. [27].The posterior probability for a predicted quantity ∆ at level n given observational data D and models M n+1 at level n + 1 is: where E M n+1 is the expectation operator, which calculates the mean over models M n+1 as follows: where Pr M (m) n+1 |D, M n are the conditional posterior model probabilities for models at level n + 1 under their parent models at level n, which are calculated by: where Pr M (m) n+1 |D are the posterior model probabilities for models at level n + 1, and posterior model probabilities at level n are given by: Equation ( 4) shows that posterior model probabilities at level n can be calculated as long as the posterior model probabilities of base models are known.
Note that the posterior model probabilities in the BMA paradigm do not represent absolute importance, but they are relative probabilities for the models selected by the analyst.Posterior model probabilities change when the pool of models changes.Moreover, Using BMA or HBMA takes place under the mutually exclusive and collectively exhaustive (MECE) principle, where models are considered to be mutually independent of each other and the combination of all models covers all possible outcomes.Analysts should exhaust all possible models for ensemble analysis.
Since Equation ( 1) is a recursive equation, one can expand the right-hand side of Equation (1) up to the base models: Using the law of total expectation and law of total variance, the expectation and variance of predicted quantity ∆ at level n are: where E ∆|D, M p is the prediction mean of ∆ given observational data D and models M p at level p. Var(∆|D, M n+1 ) is the prediction variance using models at level n + 1.It includes the within-model variance E M n+1 [Var(∆|D, M n+1 )] and the between-model variance Var is the variance operator, which calculates the between-model variance using models at level n + 1.For example, given an RCP at level 1 (see Figure 1), Var(∆|D, M 1 ) is the projection variance for ∆ using various GCMs (M 2 ) at level 2 under the same RCP.The within-model variance is the expectation of projection variances of individual GCMs, Var(∆|D, M 2 ), under the same RCP.The between-model variance Var M 1 [E(∆|D, M 1 )] is due to the variation of the mean predicted ∆ by different GCMs at level 1.Since Equation ( 7) is a recursive equation, the within-model variance is obtained using: Water 2019, 11, 268 which contains the within-model variance and the between-model variance for ∆ using models at level n + 2. The between-model variance is:

Mean and Variance at the Hierarch Level
The prediction mean of ∆ at the hierarch level can be obtained by averaging the prediction means of models at level 1 or using Equation ( 2) recursively from the base level to level 1: The total prediction variance of ∆ at the hierarch level is: Equation ( 11) is used to evaluate the contribution of individual sources of uncertainty to the total uncertainty.Using Figure 1 as an example, the total variance of multi-model hydroclimate projections can be decomposed into three components: The first term on the right-hand side of Equation ( 12) is the projection variance due to emission pathway.The second term is the projection variance caused by different GCMs.The third term is the projection variance caused by parameter uncertainty in individual climate models, which can be assessed by the perturbations of parameters over the range of plausible values [48].However, given that such information is unavailable for most of the CMIP5 models, the third term is not considered in this study.
The relative contributions of these sources of uncertainty are defined as the proportions of their corresponding variance terms to the total variance.Building upon the conventional BMA method [27], this HBMA method provides a theoretical sound foundation for the quantification of relative contributions from different sources of uncertainty in multi-ensemble hydroclimate projections.

Study Area and Model Data
Figure 2 shows the study area, the Southern Hills-Gulf region.The study area covers the Southern Hills regional aquifer system [49] and the Gulf region of southeastern Louisiana.We analyze hydroclimate projection uncertainty due to emission scenarios and GCM uncertainties.The study area lies between latitudes 28.93 • N and 33.06 • N and between longitudes 88.81 • W and 91.98 • W and covers 79,126 km 2 .It includes 28 parishes of Louisiana and 20 counties of Mississippi.The area contains 26 U.S. Geological Survey (USGS) 8-digit hydrologic units (HUC8 or subbasin) and 728 12-digit hydrologic units (HUC12 or subwatershed) with average sizes of 588 and 33 km 2 , respectively.Precipitation at the outcrops in southwestern Mississippi is the main source of water to the deep aquifers in southeastern Louisiana.
To evaluate the posterior model probability for each GCM, the gridded 1950-2006 monthly precipitation and temperature observation over the study area were obtained from Maurer et al. [50] and Maurer [51].The horizontal grid resolution of precipitation and temperature is 1/8 • (≈12 km), which covers 439 cells in the study area.To evaluate the posterior model probability for each GCM, the gridded 1950-2006 monthly precipitation and temperature observation over the study area were obtained from Maurer et al. [50] and Maurer [51].The horizontal grid resolution of precipitation and temperature is 1/8° (≈12 km), which covers 439 cells in the study area.
To enable high spatial-resolution hydrologic projections, this study created vector-based data maps that included land use/land cover, topographic slope, base flow index, HUC8 boundaries, surficial soil type, leaf area index (LAI), lithology, precipitation, and temperature, and intersected the data maps using Esri ArcGIS 10.2.As a result, the study area was divided into 2,669,533 subdivisions, each of which had homogeneous model parameters for high-performance hydrologic simulation.Detailed information on the various data and maps for the study area can be found in Beigi and Tsai [52,53].

Downscaled Climate Projection Data Set
The downscaled daily precipitation and temperature projections from 1950 to 2099 were obtained from the Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections data archive [45,46], where the bias corrected constructed analogs (BCCA) method was used to downscale the raw GCM outputs to 1/8° (≈12 km) horizontal resolution, which is consistent with the gridded To enable high spatial-resolution hydrologic projections, this study created vector-based data maps that included land use/land cover, topographic slope, base flow index, HUC8 boundaries, surficial soil type, leaf area index (LAI), lithology, precipitation, and temperature, and intersected the data maps using Esri ArcGIS 10.2.As a result, the study area was divided into 2,669,533 subdivisions, each of which had homogeneous model parameters for high-performance hydrologic simulation.Detailed information on the various data and maps for the study area can be found in Beigi and Tsai [52,53].

Downscaled Climate Projection Data Set
The downscaled daily precipitation and temperature projections from 1950 to 2099 were obtained from the Downscaled CMIP3 and CMIP5 Climate and Hydrology Projections data archive [45,46], where the bias corrected constructed analogs (BCCA) method was used to downscale the raw GCM outputs to 1/8 • (≈12 km) horizontal resolution, which is consistent with the gridded observations [50,51].The archive integrates the state-of-the-art climate models and the most recent GHG emission scenarios [54].The historical GHG forcing was used during the baseline control simulation (before 2006), and followed by four RCPs-RCP2.6,RCP4.5, RCP6.0, and RCP8.5-representing the lowest, lower, midrange, and highest emission pathways, respectively.RCP2.6 is the most optimistic scenario that features the least warming (0.3-1.7 • C) and the lowest GHG concentration.On the contrary, RCP8.5 portrays the most pessimistic scenario that features the highest GHG concentration and the greatest warming (2.6-4.8 • C) across the 21st century [55].
A total of 40 sets of climate projections that are combinations of four RCPs and 10 CMIP5 GCMs were comprehensively collected (see Table 1).Referring to the hierarchical structure illustrated in Figure 1, there are four propositions (RCPs) in the first level.Under each RCP, there are 10 propositions (GCMs) in the second level.Note also that because downscaled climate projections were available only under the BCCA downscaling method, this case study did not consider uncertainty from other downscaling techniques.

Climate Model Evaluation
Although there are various approaches to evaluating the reasonableness of a GCM from different perspectives (e.g., using observed sea surface temperature [56,57]), given our specific hydrologic simulation needs, we focused on examining the similarity of simulated 1950-2006 monthly temperature and precipitation observations to historical observations.Based on the CMIP5 experimental design, the climate projections across different RCPs from a particular GCM will have a common baseline period (before 2006).Therefore, we evaluated the 10 GCMs listed in Table 1 as a basis to estimate the posterior model probabilities in Equation ( 4).Because statistical bias-correction has been a part of BCCA (i.e., the long-term mean has been corrected by observation), the performance evaluated herein was mainly governed by the similarity of interannual variability of temperature and precipitation to historical observations.The study also assumed equal prior model probability for the 10 climate models before their outputs were compared to the historical precipitation and temperature data.
By adopting the Bayesian information criterion (BIC) and the variance window [58], the posterior model probabilities of the climate models were approximated as follows: where p − BIC min and α is the scaling factor that determines the size of the variance window.BIC (m) p is the BIC value of climate model M (m) p , and BIC min is the minimum BIC value among the climate models.Focusing on monthly precipitation and temperature, the BIC can be written as: where P p values or posterior model probabilities.

Hydrologic Modeling
In this study, HELP3 [47] was adopted to estimate evapotranspiration (ET), surface runoff, and groundwater recharge in the Southern Hills-Gulf region given the aforementioned climate projections.
HELP3 is a quasi-two-dimensional hydrologic model that simulates processes including ET, surface runoff, soil moisture storage, and vertical unsaturated drainage for individual layered soil columns.The HELP3 model has been used in many hydrological studies to estimate groundwater recharge, runoff, and ET (e.g., References [59][60][61][62][63]) and has been extensively verified to be a suitable model for hydrologic studies [47,64,65].

HELP3 Input Data
Daily precipitation, temperature, solar radiation, annual wind speed, and quarterly relative humidity are required inputs to the HELP3 model.For this study, we compiled historical daily precipitation and temperature data from 1950 to 2009 at a 12-km resolution from Maurer et al. [50] and Maurer [51].We also compiled averaged quarterly relative humidity and annual wind speed that were obtained from the Southern Regional Climate Center [66].The downscaled solar radiation data were generated synthetically using the weather generator (WGEN) model of Richardson and Wright [67].
Detailed surficial soil texture classes were obtained from the Natural Resources Conservation Service [68] for each soil layer in the HELP3 model.Land use and land cover data at high spatial resolution (30 m) were obtained from the National Land Cover Dataset 2011 [69], which is the latest national land cover product produced by the Multi-Resolution Land Characteristics (MRLC) Consortium [70].The status quo of land use/land cover and soil type is assumed as inputs to HELP3 for 2010-2099.

Parallel Computation for High-Resolution Hydrologic Prediction
Forty sets of downscaled daily precipitation and temperature projections were used as input to drive the HELP3 model that projects future evapotranspiration, surface runoff and groundwater recharge.For each of the 40 climate projections, a large number (2,669,533) of HELP3 model runs were needed for hydrologic prediction for the study area.It was computationally impractical to execute millions of HELP3 model runs in a single-core processor.The Python programming language was used in this study to feed 2,669,533 HELP3 model runs to SuperMike-II, a supercomputer with multiple cores at Louisiana State University for parallel computation.SuperMike-II is a 146-TFlop Peak Performance cluster running on the Red Hat Enterprise Linux 6 operating system.SuperMike-II contains 440 compute nodes, each of which consists of two 8-Core Sandy Bridge Xeon 64-bit processors with 2.6 GHz core frequency.
We used the USGS WaterWatch runoff data of the hydrologic units (see Figure 2) to adjust the curve number (CN) for each subdivision in the HELP3 model.Yearly USGS computed runoffs from 1950 to 2009 (60 years) were used for model calibration.Based on the regression Equation (34) in HELP Engineering Documentation for Version 3 [47], an equation was proposed to adjust user-specified antecedent moisture condition II (AMC-II) curve number for different topographic slope and slope length conditions for each subdivision.The optimal curve numbers (CNs) for each hydrologic unit were obtained by minimizing the sum of squared errors of the HELP3 calculated to USGS computed direct runoffs as follows: where Q HELP3 u is the calculated direct runoff at year t for hydrologic unit u from the HELP3 model and is the computed direct runoff at year t for hydrologic unit u from the USGS WaterWatch database.Readers are referred to Beigi and Tsai [52] for detailed model calibration techniques.The root mean square error (RMSE) on yearly runoff ranges from 90.32 to 137.45 mm, which are much smaller than the range of the WaterWatch yearly runoff data.The Nash-Sutcliffe model efficiency coefficient (NSE) ranged from −0.10 to 0.50.The lowest three NSE values are −0.10,0.10, and 0.11 for Coles Creek, Bayou Cocodrie, and Lake Maurepas HUC8s, respectively.Other HUC8s have an NSE higher than 0.15.Watersheds with low NSE show high runoff variability in the WaterWatch data.Low NSE indicates that the HELP3 predictions are as accurate as the mean of the WaterWatch data but cannot capture the data variability.The highest three NSE values are 0.50, 0.47, and 0.44 for Tangipahoa, Middle Pearl-Silver, and Middle Pearl-Strong HUC8s, respectively.We used the estimated ET from the MOD16 evapotranspiration dataset [71,72] to verify the calibrated HELP3 model.The RMSE of the HELP3 yearly estimated ET from 2000 to 2010 versus the MOD16 yearly computed ET is 93.52 mm, which is much smaller than the range of the MOD16 ET estimates.The NSE is 0.59.

Posterior Model Probabilities
Based on Equations ( 13) and ( 14), we evaluated the posterior model probabilities for the 10 GCMs listed in Table 1 using simulated and observed 1950-2006 monthly temperature and precipitation observations.Table 2 summarizes the weighted fitting residuals and ∆BIC values for the 10 climate models with individual posterior model probabilities.The GCM model, gfdl-esm2g.1, had the highest posterior model probability (23.0%), followed by miroc-esm-chem.1 (19.4%) and miroc5.1 (14.7%).The posterior model probabilities of GCMs in Table 2 represent conditional posterior model probabilities at level 2 under each RCP shown in Figure 3.The GCM probability distributions were the same under each RCP; therefore, each RCP had a posterior model probability 25% at level 1.The posterior model probabilities of GCMs at level 2, according to Equation (3), are shown in Figure 3.

Temporal Analysis with HBMA
Recharge, runoff, and ET projections in the Southern Hills-Gulf region for 2010-2099 using the first and the second best GCM members at level 2 for four emission pathways are shown in Figure 4. Projected recharge showed much larger variability compared to runoff and ET projections.Recharge

Temporal Analysis with HBMA
Recharge, runoff, and ET projections in the Southern Hills-Gulf region for 2010-2099 using the first and the second best GCM members at level 2 for four emission pathways are shown in Figure 4. Projected recharge showed much larger variability compared to runoff and ET projections.Recharge quantities projected by different climate models were quite different.The ET projection for the 21st century was consistently high regardless of climate models and emission pathways.

Temporal Analysis with HBMA
Recharge, runoff, and ET projections in the Southern Hills-Gulf region for 2010-2099 using the first and the second best GCM members at level 2 for four emission pathways are shown in Figure 4. Projected recharge showed much larger variability compared to runoff and ET projections.Recharge quantities projected by different climate models were quite different.The ET projection for the 21st century was consistently high regardless of climate models and emission pathways.Using BMA at level 1, the mean projection and standard deviation of recharge, runoff, and ET due to using the 10 GCMs are shown in Figure 5. Recharge projection showed the highest uncertainty, and ET projection showed the lowest uncertainty.Using linear trend analysis, RCP2.6 projected an increasing trend (the slope of a linear model = 0.31 mm/year) in recharge while RCP4.5 (slope = −0.46Using BMA at level 1, the mean projection standard deviation of recharge, runoff, and ET due to using the 10 GCMs are shown in Figure 5. Recharge projection showed the highest uncertainty, and ET projection showed the lowest uncertainty.Using linear trend analysis, RCP2.6 projected an increasing trend (the slope of a linear model = 0.31 mm/year) in recharge while RCP4.5 (slope = −0.46mm/year), RCP6.0 (slope = −1.01mm/year), and RCP8.5 (slope = −1.16mm/year) each showed a decreasing trend in recharge.For runoff projection, RCP2.6 projected an increasing trend (slope = 0.25 mm/year), and RCP4.5 (slope = −0.06mm/year), RCP6.0 (slope = −0.31mm/year), and RCP8.5 (slope = −0.24mm/year) each showed a decreasing trend.None of the emission pathways showed a noticeable trend (slope was positive, but less than 0.1 mm/year) in the ET projections.Generally speaking, ET tended to increase as recharge and runoff were likely to decrease as RCP tended toward more pessimistic scenarios.
Figure 6 shows the annual recharge, runoff, and ET predictions from 2010 to 2099.The BMA mean and one standard deviation interval resulted from the 10 GCMs and 4 emission pathways in a hierarchical order (Figure 3).The results of all 40 climate projections indicate that mean recharge would slightly decrease while mean runoff and ET would remain the same throughout 2099.The temporal averages (and ranges) of the coefficient of variation (CV) for yearly recharge, runoff, and ET projections were 40.1% (31.5-54.1%),28.6% (23.3-35.6%),and 5.4% (3.6-8.9%),respectively.Projected ET and runoff have smaller CVs than projected recharge, indicating that recharge estimation was more sensitive to climate projections.Moreover, as shown in Figures 4-6, the loss in the variability of mean hydrologic projections compared to individual hydrologic projections was compensated by the increase in projection variances.Therefore, the variability of future projections was implicitly preserved in the HBMA method.
Water 2019, 11, x FOR PEER REVIEW 11 of 20 mm/year), RCP6.0 (slope = −1.01mm/year), and RCP8.5 (slope = −1.16mm/year) each showed a decreasing trend in recharge.For runoff projection, RCP2.6 projected an increasing trend (slope = 0.25 mm/year), and RCP4.5 (slope = −0.06mm/year), RCP6.0 (slope = −0.31mm/year), and RCP8.5 (slope = −0.24mm/year) each showed a decreasing trend.None of the emission pathways showed a noticeable trend (slope was positive, but less than 0.1 mm/year) in the ET projections.Generally speaking, ET tended to increase as recharge and runoff were likely to decrease as RCP tended toward more pessimistic scenarios.Figure 6 shows the annual recharge, runoff, and ET predictions from 2010 to 2099.The BMA mean and one standard deviation interval resulted from the 10 GCMs and 4 emission pathways in a hierarchical order (Figure 3).The results of all 40 climate projections indicate that mean recharge would slightly decrease while mean runoff and ET would remain the same throughout 2099.The temporal averages (and ranges) of the coefficient of variation (CV) for yearly recharge, runoff, and ET projections were 40.1% (31.5-54.1%),28.6% (23.3-35.6%),and 5.4% (3.6-8.9%),respectively.Projected ET and runoff have smaller CVs than projected recharge, indicating that recharge estimation was more sensitive to climate projections.Moreover, as shown in Figures 4 through 6, the loss in the variability of mean hydrologic projections compared to individual hydrologic projections was compensated by the increase in projection variances.Therefore, the variability of future projections was implicitly preserved in the HBMA method.

Future Hydrologic Projection Anomalies
Future hydrologic projection anomalies with respect to historical means were assessed by looking into the future 30-year mean projected changes with respect to the mean annuals between 1950 and 2009.The future recharge, runoff, and ET projection anomalies at different BMA levels is shown in Table 3.The estimated mean annual ET, runoff, and recharge between 1950 and 2009 were 832.9, 352.8, and 337.4 mm, respectively, which accounted for 54.4%, 23.1%, and 22.1% of the mean annual precipitation.The best and the second-best climate models at level 2 in Table 3 under each emission pathway gave more positive recharge anomalies than negative recharge anomalies in 2010-

Future Hydrologic Projection Anomalies
Future hydrologic projection anomalies with respect to historical means were assessed by looking into the future 30-year mean projected changes with respect to the mean annuals between 1950 and 2009.The future recharge, runoff, and ET projection anomalies at different BMA levels is shown in Table 3.The estimated mean annual ET, runoff, and recharge between 1950 and 2009 were 832.9, 352.8, and 337.4 mm, respectively, which accounted for 54.4%, 23.1%, and 22.1% of the mean annual precipitation.The best and the second-best climate models at level 2 in Table 3 under each emission pathway gave more positive recharge anomalies than negative recharge anomalies in 2010-2069, indicating that recharge is likely to increase in this 60-year period.From 2070 to 2099, recharge is likely to increase under RCP2.6 and is likely to decrease under RCPs 6.0 and 8.5.The rcp26.gfdl-esm2g.1 model projected the highest recharge increase (13.49%), while the rcp85.miroc-esm-chem.1 model projected the largest recharge decrease (−33.0%) between 2070 and 2099.With respect to the historical mean, the runoff is predicted to decrease and ET is predicted to increase throughout 2099.Consistent results were found between the projection anomalies at level 1 (Table 3) and the established future emissions.Under the assumption of lesser warming and lower GHG concentration, the RCP2.6 and RCP4.5 projected recharge increase was between 4.5% and 19.6% in the 21st century.RCP6.0 and RCP8.5 also projected more than a 10% recharge increase in 2010-2039.However, the recharge could be significantly reduced by as much as 16% toward the end of the century.All emission pathways (level 1) showed consistent projections that runoff is predicted to decrease between 15% and 25%, while ET is predicted to increase between 2% and 4% in this century.
Using all the 40 climate projections (the hierarch model), the modeling results show that recharge in the Southern Hills-Gulf region is likely to increase 13% in 2010-2039.This is due to increase in precipitation and decrease in surface runoff in 2010-2039.However, the recharge increase would reduce to only 3% in 2070-2099 because of the decrease in precipitation.Runoff is likely to decrease between 18% and 20% and ET is likely to increase around 3% in this century.

Spatial Analysis
Figure 7 shows the mean annual ET, recharge, and runoff distributions over the 728 HUC12s using all 40 climate projections.The maximum changes of ET, runoff, and recharge projections over the three 30-year periods were only 38.6 mm, 56.2 mm, and 44.7 mm, respectively.Outcrops of Miocene deposits in southwestern Mississippi showed a high recharge rate, which indicated good future groundwater recharge to deep sands (e.g., the "2000-foot" sand and the "2400-foot" sand) in southeastern Louisiana.Figure 8 shows the spatial distributions of BMA standard deviations.Projection uncertainty increased toward the end of the century for all hydrologic projections.Annual recharge projection had the greatest uncertainty.In 2070-2099, the standard deviation of projected annual recharge had a large range from 4.1 to 123.2 mm, followed by runoff (1.7-96.7 mm) and ET (7.8-45.6 mm).
future groundwater recharge to deep sands (e.g., the "2000-foot" sand and the "2400-foot" sand) in southeastern Louisiana.Figure 8 shows the spatial distributions of BMA standard deviations.Projection uncertainty increased toward the end of the century for all hydrologic projections.Annual recharge projection had the greatest uncertainty.In 2070-2099, the standard deviation of projected annual recharge had a large range from 4.1 to 123.2 mm, followed by runoff (1.7-96.7 mm) and ET (7.8-45.6 mm).

Contributions of Sources of Uncertainty
Contributions of individual sources of uncertainty from the emission paths and climate models to the total recharge projection uncertainty are shown in Figure 9.The contributions are the proportions of variances in the first and the second terms in Equation ( 12) to the sum of the three variances.Figure 9 shows that the dominant source of uncertainty in hydrologic projections came

Contributions of Sources of Uncertainty
Contributions of individual sources of uncertainty from the emission paths and climate models to the total recharge projection uncertainty are shown in Figure 9.The contributions are the proportions of variances in the first and the second terms in Equation ( 12) to the sum of the three variances.Figure 9 shows that the dominant source of uncertainty in hydrologic projections came from GCM uncertainty.Moreover, emission pathway uncertainty increased toward the end of the century because under different emission paths, hydrologic projections started to diverge significantly after the mid-century.Table 4 shows the uncertainty contributions from GCMs and emission paths to the projected recharge.In 2010-2039, almost all uncertainty in recharge projection came from GCM uncertainty, which reduced to 66.0% in 2070-2099.Emission path uncertainty grew from 5.4% in 2010-2039 to 34% in 2070-2099.The standard deviations were small, indicating that variations of individual source uncertainty contributions across the 728 HUC12s were small.Although not presented, the GCM and emission path uncertainties had similar contributions to the projected runoff and ET.
Water 2019, 11, x FOR PEER REVIEW 15 of 20 downscaling method, hydrological model structure, hydrological model parameters, and internal variability of climate system.Their results suggested that uncertainty arising from GCM structure is the largest source of uncertainty in studying impacts of climate change on flood frequency in England.In summary, regardless of the number of sources contributing to total uncertainty, there is a consensus among previous studies that GCM uncertainty is the major source of uncertainty when assessing the impact of climate change on hydrologic projections (e.g., References [8,79]).The other finding is also consistent with Northrop and Chandler [80] where uncertainty resulting from the choice of emission pathways becomes more important toward the end of the twenty-first century.This indicates that climate model parameterizations, sensitivities, and responses to greenhouse gases and other forcings grow over time and become evident.Although the contribution of GCM uncertainty to the total uncertainty decreases every 30-year period, this does not mean that GCM uncertainty decreases over time, but reflects the significant growth of emission pathway uncertainty.The total uncertainty actually increased towards the end of the century as shown in Figure 8.

Conclusions
The hierarchical Bayesian model averaging (HBMA) method is able to characterize the hierarchical natural of sources of uncertainties in climate projections.The HBMA was shown be to a good approach to analyze 40 sets of CMIP5 BCCA-downscaled daily precipitation and temperature projections that were generated using 10 GCMs and 4 emission paths.Through the HELP3 model,   The finding that GCM uncertainty dominated in hydrological projections is consistent with many studies [6,8,[73][74][75], although numerical quantification of contributions in the literature is not as available as in this study.Hosseinzadehtalaei et al. [76] found that the choice of GCMs is the major contributor (up to 65% for some cases) toward intensifying precipitation change uncertainty for all return periods.Also, Mandal et al. [77] concluded that the choice of GCMs is the largest sources of uncertainty of the climate change impacts on total monthly precipitation when uncertainty arising from downscaling methods are excluded.Taylor and Kingston [78] found that GCM-related uncertainty substantially dominates uncertainty associated with climate sensitivity and hydrological model parametrization.Also, Teng et al. [75] indicated that uncertainty due to GCMs is significantly greater than uncertainty sources from rainfall-runoff models in assessing climate change impacts on runoffs.Kay et al. [8] investigated six different sources of uncertainty including: GCM structure, downscaling method, hydrological model structure, hydrological model parameters, and internal variability of climate system.Their results suggested that uncertainty arising from GCM structure is the largest source of uncertainty in studying impacts of climate change on flood frequency in England.In summary, regardless of the number of sources contributing to total uncertainty, there is a consensus among previous studies that GCM uncertainty is the major source of uncertainty when assessing the impact of climate change on hydrologic projections (e.g., References [8,79]).
The other finding is also consistent with Northrop and Chandler [80] where uncertainty resulting from the choice of emission pathways becomes more important toward the end of the twenty-first century.This indicates that climate model parameterizations, sensitivities, and responses to greenhouse gases and other forcings grow over time and become evident.Although the contribution of GCM uncertainty to the total uncertainty decreases every 30-year period, this does not mean that GCM uncertainty decreases over time, but reflects the significant growth of emission pathway uncertainty.The total uncertainty actually increased towards the end of the century as shown in Figure 8.

Conclusions
The hierarchical Bayesian model averaging (HBMA) method is able to characterize the hierarchical natural of sources of uncertainties in climate projections.The HBMA was shown be to a good approach to analyze 40 sets of CMIP5 BCCA-downscaled daily precipitation and temperature projections that were generated using 10 GCMs and 4 emission paths.Through the HELP3 model, the HBMA was able to derive means and variances of ET, runoff, and groundwater recharge projections under climate model uncertainty and emission path scenario uncertainty.
This study conducted hydroclimate modeling for the Southern Hills-Gulf region, which is computationally demanding because the HELP3 model is highly parameterized, where 40 sets of daily precipitation and temperature projections are used, and the simulation period is a century long.By using the CPU-based multi-core supercomputer, this study was able to efficiently execute several million runs of the HELP3 model in order to quantify climate-related hydrologic projection uncertainty.
Projected groundwater recharge shows higher spatial and temporal variability and variance than ET and runoff in the Southern Hills-Gulf region, indicating that recharge estimation is more sensitive to climate projections.Evapotranspiration projection shows the least variability and variance in space and time.The uncertainty (in terms of standard deviation) in projected ET, runoff, and recharge increased toward the end of the century, indicating that the 40 climate projections diverge chronologically.
The anomalies with respect to the past record  show that groundwater recharge in the Southern Hills-Gulf region is predicted to increase 13% in the early century, but is predicted to reduce to 3% increase in the late century.Runoff is predicted to decrease between 18% and 20% and ET is predicted to increase around 3% throughout the century.
The study found that the GCM uncertainty prevailed the total uncertainty in future hydrologic projections in the Southern Hills-Gulf region.For recharge projection in the early century, 95% of total uncertainty came from the GCM uncertainty.The recharge project uncertainty in the late century had 66% from the GCM uncertainty.The emission pathway uncertainty was negligible in the early century, but become significant towards the end of the century.

Figure 1 .
Figure 1.Hierarchical structure of sources of uncertainty in the CMIP5 multi-model ensembles.

M
and so forth, until the base level p is reached.Furthermore, models 1 p − M at level p − 1 are BMA models of their corresponding child base models p at level p − 2 are BMA models of their corresponding child BMA models 1 p −

Figure 1 .
Figure 1.Hierarchical structure of sources of uncertainty in the CMIP5 multi-model ensembles.

Figure 2 .
Figure 2. The study area of the Southern Hills-Gulf region, which includes 26 HUC8 (dark gray lines) and 728 HUC12 (light gray lines).

Figure 2 .
Figure 2. The study area of the Southern Hills-Gulf region, which includes 26 HUC8 (dark gray lines) and 728 HUC12 (light gray lines).
are the ith simulated monthly precipitation and temperature by model M ith observed monthly precipitation and temperature, σ 2 P,i and σ 2 T,i are the ith error variance of precipitation and temperature, N is the total number of monthly data from the 439 cells and from 1950 to 2006 (N = 439 × 57 × 12), and k (m) p is the number of model (GCM) parameters.The first two terms on the right-hand side of Equation (14) are the weighted fitting residuals.The term N ln 2π + k (m) p ln(N) represents the complexity of climate model M (m)p .Due to the lack of information on the complexity term in Equation (14), this study assumed equal model complexity for all climate models such that the complexity term had no effect on ∆BIC (m)

Water 2019 , 20 Figure 3 .
Figure 3.A BMA tree of posterior model probabilities for emission pathways and GCMs at levels 1 and 2.

Figure 3 .
Figure 3.A BMA tree of posterior model probabilities for emission pathways and GCMs at levels 1 and 2.

Figure 3 .
Figure 3.A BMA tree of posterior model probabilities for emission pathways and GCMs at levels 1 and 2.

Figure 4 .
Figure 4. Annual recharge, ET and runoff projections (2010-2099) in Southern Hills-Gulf region using the best and the second-best climate models at level 2.

Figure 4 .
Figure 4. Annual recharge, ET and runoff projections (2010-2099) in Southern Hills-Gulf region using the best and the second-best climate models at level 2.

Figure 5 .
Figure 5. Means and one standard deviation (SD) intervals of annual recharge, ET and runoff projections (2010−2099) in the Southern Hills-Gulf region using all climate models under corresponding emission pathways at level 1.

Figure 5 .
Figure 5. Means and one standard deviation (SD) intervals of annual recharge, ET and runoff projections (2010−2099) in the Southern Hills-Gulf region using all climate models under corresponding emission pathways at level 1. Water 2019, 11, x FOR PEER REVIEW 12 of 20

Figure 6 .
Figure 6.Means and one standard deviation (SD) intervals of (a) annual recharge, and (b) ET and runoff projections (2010-2099) in the Southern Hills-Gulf region at the hierarch level.

Figure 6 .
Figure 6.Means and one standard deviation (SD) intervals of (a) annual recharge, and (b) ET and runoff projections (2010-2099) in the Southern Hills-Gulf region at the hierarch level.

Table 4 .
Means and standard deviations (SD) of uncertainty contributions in 2010-2039, 2040-2069, and 2070-2099 from emission paths and GCMs to the projected recharge (including all 728 HUC12) in the Southern Hills-Gulf region.

Table 4 .
Means and standard deviations (SD) of uncertainty contributions in 2010-2039, 2040-2069, and 2070-2099 from emission paths and GCMs to the projected recharge (including all 728 HUC12) in the Southern Hills-Gulf region.

Table 2 .
Posterior model probabilities of 10 climate models based on monthly precipitation and temperature data 1950-2006 in the Southern Hills-Gulf region.

Table 3 .
Anomalies of mean annuals in 2010-2039, 2040-2069, and 2070-2099 to mean annuals of 1950-2009 in recharge, runoff, and evapotranspiration (ET) in the Southern Hills-Gulf region.Only the best and the second best GCMs under each emission path are listed in level 2.