**Advertisement**

ABSTRACT

Planning of mitigation and adaptation strategies to a changing climate can benefit from a good understanding of climate change impacts on human life and local society, which leads to an increasing requirement for reliable projections of future climate change at regional scales. This paper presents an ensemble of high-resolution regional climate simulations for the province of

(ProQuest: ... denotes formulae omitted.)

1. Introduction

Climate change is becoming one of the most pressing issues around the world. It has already started to affect every continent, country, community, and individual. As the largest economy in

Forecasts of future climate change with current stateof- the-art climate models can provide potential evidence for decision makers and policy makers to help determine how to adapt to or mitigate climate change. However, these forecasts are inevitably uncertain due to our incomplete understanding of the climate system in terms of its complicated physical processes and natural variability, as well as the response to rising levels of greenhouse gases (Allen et al. 2000; Murphy et al. 2004; Stainforth et al. 2005; Stott and Kettleborough 2002; Webster et al. 2003). This further results in considerable uncertainties about the rates of change that can be expected, such as changes in extremes of temperature and precipitation and sea level rise (Karl and Trenberth 2003). No single model can be powerful enough to tackle such uncertainties all at once, so it is necessary to utilize results from a range of coupled models (Houghton et al. 2001).

Previously, a number of climate research projects based on multimodel ensemble (MME) or perturbed physics ensemble (PPE) approaches have been carried out to explore techniques for quantifying uncertainties of future climate change (e.g., Barnett et al. 2006; Furrer et al. 2007; Giorgi and Francisco 2000; Giorgi and Mearns 2002, 2003; Greene et al. 2006; Harris et al. 2013; Murphy et al. 2007; Sexton et al. 2012; Tebaldi and SansÓ 2009; Tebaldi et al. 2005; Watterson and Whetton 2011). The MME approach usually consists of various GCMs developed at different modeling centers around the world to sample both structural and parametric uncertainties to a limited degree (Meehl et al. 2007; Taylor et al. 2012), but it is unable to sample the consequences of either type of uncertainty in a systematic fashion since it is assembled on an opportunistic basis from current available models (Murphy et al. 2007; Sexton et al. 2012). The PPE usually consists of variants of a single base model with perturbed parameters limited in a space of possible model configurations (Collins et al. 2006; Murphy et al. 2004; Webb et al. 2006; Yokohata et al. 2010). The main advantage of the PPE approach is that it allows greater control over the design of experiments to sample parametric uncertainties within a single model framework. Both ensemble approaches can generate a large number of projections of future climate for various scenarios, but how to combine these multiple projections and interpret them into policy-relevant information has become a major challenge in recent years, due to the lack of verification of climate projections, the problems of model dependence, bias, and tuning, and the difficulty in making sense of an ''ensemble of opportunity'' (Knutti et al. 2010; Tebaldi and Knutti 2007).

There are obviously different ways to synthesize modeling results of MMEs or PPEs. A straightforward way is to calculate the multimodel averages for given diagnostics or variables where each model is weighted equally (Reuroaiseuroanen and Palmer 2001). In many cases, combining ensemble results through Bayesian methods or weighted averages, where weights are determined by comparing model predictions to observations, shows improved performance better than simple averages (e.g.,

However, previous probabilistic studies have mainly focused on quantifying climate change uncertainties at global scales based on various GCM ensembles (e.g., Furrer et al. 2007; Greene et al. 2006; Tebaldi and Knutti 2007; Tebaldi and SansÓ 2009). Additional downscaling work needs to be done before feeding the probabilistic projections into impact models (e.g., hydrological modes or crop models) because of the coarse resolution of GCM outputs, by either statistical approaches (e.g., Wang et al. 2013; Wilby et al. 2004) or regional climate models (RCMs) . A small number of papers published in recent years (e.g., Benestad 2004; Luo et al. 2005; Manning et al. 2009; New et al. 2007) have approached the problem of deriving regional probabilistic projections at a finer resolution to facilitate the assessment of regional and local climate change impacts, but they are characterized by a less general approach and are tailored to specific regions or impacts studies.

The objective of this research is to develop highresolution probabilistic projections of temperature changes for the province of

2. Development of RCM ensemble

a. Study area

As the second largest province in

In 2007,

b. Regional climate modeling using PRECIS

Regional climate modeling is an essential way for conducting thorough assessments of climate change impacts by providing impact researchers with regional detail of how future climate might change. An RCM is a powerful tool to add small-scale detailed information of future climate change to the large-scale projections of a GCM. RCMs are physics-based climate models and as such represent most or all of the processes, interactions, and feedbacks between the climate system components that are represented in GCMs (Jones et al. 2004). By taking coarse-resolution information from GCMs, RCMs can develop temporally and spatially finescale information using their higher-resolution representation of the climate system. In this study, we apply a widely used RCM developed at the Hadley Centre, PRECIS, to facilitate regional climate modeling over

Using a single GCM projection to drive RCMs can provide us with some information about the expected changes under a given emission scenario for the region of interest, but it does not provide us with more information about how confident we should be in those changes (Bellprat et al. 2012). The ensemble approach through either MMEs or PPEs is widely accepted as an effective way to explore the range or spread of projections from multiple members, which enables us to gain a better understanding of the uncertainties in climate modeling. The Hadley Centre has published 17 sets of boundary data from a perturbed physics ensemble (i.e., HadCM3Q0-Q16, known as QUMP), which is based on HadCM3 under the Special Report on Emissions Scenarios (SRES) A1B emission scenario, for use with PRECIS in order to allow users to generate an ensemble of high-resolution regional simulations (McSweeney and Jones 2010). Downscaling the 17- member PPE ensemble with PRECIS would require very large inputs of computing resources, data storage, and data analyses. To explore the range of uncertainties while minimizing these requirements, we select a subset of five members (i.e., HadCM3Q0, Q3, Q10, Q13, and Q15) from the QUMP ensemble according to the Hadley Centre's recommendation (see http://www.metoffice.gov.uk/precis/qump). HadCM3Q0 is first selected as it is the standard, unperturbed model using the original parameter settings as applied in the atmospheric component of HadCM3. Selection of the remaining four members is based on 1) their performances in simulating key features of the climate over

We run five PRECIS experiments driven by boundary conditions from the selected GCM members from 1950 to 2099 at its highest horizontal resolution (i.e., 25 km). This allows us to carry out undermentioned probabilistic analyses by providing full simulation coverage from the present day to the future. The PRECIS model outputs are extracted and divided into four 30-yr periods: one baseline period (1961-90), and three future periods (2020-49, 2040-69, and 2070-99), representing its simulations for the province of

c. Data

We focus on the temperature responses over the province of

3. Bayesian hierarchical model

Suppose we have observation data for the present climate, denoted as x0, and the PRECIS simulations are expressed as xi (for current climate) and yi (for future climate), for i51, 2, . . . , N; here N 55, indicating the total number of PRECIS runs. We further assume that the PRECIS outputs depend on some parameters that are unknown due to uncertainties in climate models. In addition, the quality of observation data may be affected by some unexpected factors such as measurement error, equipment failure, and so on. Let Q be the vector of all unknown parameters involved in both observations and model simulations. The Bayesian viewpoint allows us to treat these parameters as random variables in order to quantify the uncertainties of interest in a statistical way. Specifically, we can construct a probabilistic model for random parameters Q, which are conditional to existing data D consisting of observations x0 and model simulations xi and yi, as follows:

... (1)

where p(QjD) is the posterior distribution of Q given our best understanding of the climate system based on existing observations and model simulations (it is a probabilistic representation of what we can conclude about the unknown parameters after we observe and model the climate system); p(Q) is the prior distribution of Q indicating what we know about the unknown parameters before we obtain data D; p(DjQ) is the likelihood specifying the conditional distribution of the data given all involved parameters, which is formulated under some statistical assumptions; and the symbol } means a proportional relationship up to a normalizing constant (namely marginal distribution), which is usually intractable integral and not necessary to compute in Markov chain

a. Likelihoods

We assume Gaussian distribution for observations x0:

... (2)

where the notation N(m, l21 0 ) indicates a Gaussian distribution with mean m and variance l21 0 . Here m represents the true value of current climate mean, and l21 0 is treated as a random variable to indicate that the observations are centered on the true value of current climate with a random error. Observations may suffer from both random errors (i.e., measurement and sampling) and systematic errors due to different measurement platforms and practices (

... (3)

where x;N(0, l21 0 ).

Similarly, we assume Gaussian distribution for xi:

... (4)

where li is referred to as the precision of distributions xi in estimating current climate, following the definition proposed by Tebaldi et al. (2005). The assumption underlying Eq. (4) is that the model simulation is a symmetric distribution, whose center is the true value of current climate, but with an individual variability.As such, li can be treated as a quantity for assessing the performance of the ith PRECIS run in reproducing current climate, while driven by the corresponding boundary conditions from the QUMP ensemble. Accordingly, the statistical assumption for xi can be expressed as follows:

... (5)

where hi ;N(0, l21 i ). Projecting future climate with a climate model is to some extent correlated to its capability in hindcasting observed climate, so we therefore treat yi and xi as dependent distributions through a linear regression equation. Thus yi can be formulated as follows:

... (6)

where n represents the true value of future climate mean; ji ;N[0, (uli)21], the product of uli is referred to as the precision of distribution yi in terms of simulating future climate, while u is introduced as an additional parameter to allow for different precision for yi from xi in all PRECIS runs; and b is an unknown regression coefficient. A value of b equal to 0 indicates independence between yi and xi; otherwise, positive values imply direct relations and negative for inverse ones between these two quantities. Likewise, we assume the likelihood of yi as a Gaussian distribution:

... (7)

where we use the term of b(xi 2m) to imply a linear adjustment to the projection of future climate based on the model bias for current climate simulation.

b. Prior distributions

The statistical model described by Eqs. (2), (4), and (7) are formulated using a set of parameters: fm, n, b, u, l0, l1, . . . , lNg. We assume uninformative prior distributions for these parameters as follows:

(i) The true values of current and future climate means, m and n, are assumed to have uniform priors on the real line (i.e., [2', 1']).

(ii) The regression coefficient (b) is presumed to be freely varying between 21 and 11, thus a uniform distribution on [21, 1].

(iii) Based on the estimation of Giorgi and Mearns (2002) in terms of natural variability of observed temperature at different regions for winter and summer seasons, we assume gamma prior density for l0 and the first guesses for its mean and variance are 4.5 and 19.3 respectively. Thus we formulate the prior distribution of l0 as follows:

... (8)

where m51:05 and n50:23.

(iv) We assume gamma distributions for l1, . . . , lN:

... (9)

Similarly, gamma distribution for u:

... (10)

Here we set a5b5c5d50:001; this is to translate these assumed priors into gamma distributions with mean 1 and variance 1000 (Tebaldi et al. 2005). By doing so, we can get extremely diffused priors to reflect our poor understanding about these unknown parameters.

c. Posterior distributions

Inferences for the statistical model defined in Eq. (1) can be achieved by applying Bayesian theorem to the likelihoods and priors described above. The joint posterior distribution is obtained up to a constant by taking the product of all conditional distributions. Thus, we have

... (11)

From Eq. (11) we can obtain full conditional distribution for each parameter by ignoring all other parameters that are constantwith respect to the parameter of interest. In our case, full conditional distributions for all parameters are well-known distributions such as the gamma or Gaussian ones. We can therefore perform

Here we show how to obtain the probabilistic climate change projections based on the PRECIS experiments driven by different boundary conditions. First, the full conditional distribution of m is derived from Eq. (11) by fixing all other parameters, as a Gaussian distribution:

... (12)

In a similar way, the full conditional distribution of n can be obtained as follows:

... (13)

Likewise, we can derive the full conditional posterior distributions for the remaining parameters (as presented in appendix A). By doing a series of random drawings using Gibbs sampler, we have a large number of samples for m and n. The densities of these

... (14)

Thus, the density of D can be estimated using the differences between two samples of n and m. Given the limited capability of climate models in representing the real climate system, we only can give plausible distribution for future climate change. We cannot say the absolute probability of climate changing by some exact values. Instead we talk about the probability of climate change being less than or greater than a certain value. Following the approach employed in the fifth generation of climate change information for the

4. Results

a. Capability of hindcasting current climate

As mentioned above, a five-member subset was screened out from the QUMP ensemble as boundary conditions to drive five individual PRECIS runs. To validate the capability of these five runs in simulating key features of current climate over

Figures 2-4 show the boxplots for Tmax, Tmean, and Tmin respectively. Observed monthly means of Tmax and Tmean at 12 weather stations are well captured by the aggregated simulations from the ensemble runs. However, the ensemble simulations slightly overestimate Tmin in some weather stations located in the middle and northern areas. In particular, the observed means of Tmin at stations SLA,

b. Posterior distributions at selected weather stations

By applying the Bayesian model, we derive posterior distributions for all random parameters based on the assumptions for the likelihoods and priors. Here we focus on analyzing the

1) TEMPERATURE CHANGES

Figure 5 shows posterior distributions of m and n for Tmax, Tmean, and Tmin at 12 weather stations. In particular, we draw the density curves based on the

To further analyze the temperature changes (i.e., DT 5n2m) under future climate, we extract the estimated values at three probability levels (i.e., 10%, 50%, and 90%) for Tmax, Tmean, and Tmin at 12 stations. Meanwhile, a 95% confidence interval (CI) is computed to give an estimation for the range of values that is likely to cover the possible changes under presumed uncertainties (Goldstein and Healy 1995). The results are listed in Table 2. Thus, we can interpret the temperature changes at each station with the corresponding confidence levels. For example, the Tmax change at stationWDA for the period of 2020-49 is very likely to be greater than 2.648C and is very likely to be less than 3.198C, while the median change tends to be 2.928C and the likely range for Tmax change would be [2.91, 2.93]8C for the same period; the Tmean change at station TCC under future climate of 2070-99 would be around 5.628C and it is very likely to be greater than 5.228C and less than 6.028C, while the credible estimate of changing interval would be [5.60, 5.64]8C. The results also report a gradient increase with time to the change of each temperature index at each weather station. For example, the median values of Tmean change at all stations range from 2.668C (at SLA) to 2.988C (at MUA) for the period of 2020-49; the range of these values for the period of 2040-69 would rise up to [3.72, 4.76]8C, with the lower bound at BTL and the upper one at MUA; to the end of the twenty-first century, the change of Tmean would jump to as high as 7.218C at MUA, with the least change being 5.448C at SLA.

2) MODEL RELIABILITY AND INFLATION-DEFLATION EFFECTS

We use li to reflect the precision of each PRECIS model run in terms of reproducing the current observed climate, while an inflation/deflation parameter (u) is introduced to represent its performance in simulating climate under future forcing, with the form of uli. Figures 6 and 7 present the posterior distributions for li and u in the form of box plots, with logarithmically scaled horizontal axes because of their high degree of skewness (Tebaldi et al. 2005). In our analysis, li is used as a scoring criterion to evaluate the overall performance of each PRECIS run in terms of simulating current climate forcing. It will be further treated as a weighting coefficient to reflect the contribution of each model run to the final combined estimate. However, li here is a random variable with gamma posterior distribution; the scoring of five PRECIS runs should be assessed through the relative positions of the box plots at 12 stations, rather than by comparing point estimates. To facilitate our analysis, we draw a vertical line at 1 for each box plot in Fig. 6, which is helpful to compare the overall performance of all model runs. Specifically, we use the location of the median line within a box plot relative to line 1 to represent its relative position. Evaluating the overall performance of a model run may suggest a complete consideration of its individual precisions in reproducing the current climate conditions of 12 weather stations. For example, a model with box plots mostly distributed to the right side of 1 receives a higher score than the one with left-distributed box plots gains. Thus, we can rank the precisions of five PRECIS runs by computing the proportion of stations with box plots shifted to the right (shown in Table 3). For the run driven by HadCM3Q0, the proportions for Tmax, Tmean, and Tmin are all as high as 100%, which leads to the highest rank in its overall performance. By contrast, the run driven by HadCM3Q3 gains the lowest rank because all box plots are distributed in the leftside of 1. Among the remaining three runs, the ones driven by HadCM3Q15 and Q13 are both well rated owing to their good performance in more than half the stations (ranking second and third respectively), whereas the one driven by HadCM3Q10 performs relatively poorly (ranking the fourth) in the majority of the 12 weather stations.

The inflation/deflation parameter (u) is useful to help understand if a well-performed model run in the control simulation period would behave similarly or inversely in the future climate simulation. Figure 7 indicates that almost all of the weather stations report varying levels of inflation in simulating the three temperature variables, with an average inflation factor on the order of 10. It is very interesting to realize that station MUA is the only exception, oppositely expressing a deflation behavior in the periods of 2040-69 and 2070-99. However, the deflation of this individual station is not likely to squash or offset the overall magnifying trend. The consistence of inflated performance in eleven-twelfths of weather stations may imply that the PRECIS model would project future climate with higher precision than it does in the current climate simulation. But this implication is conditional on our single-model ensemble runs driven by the boundary conditions from HadCM3; further exploration in terms of such inflation or deflation effects can be done by combining different RCMs and GCMs (Hewitt 2004; Kendon et al. 2010; Mearns et al. 2009, 2012).

3) TEMPORAL CORRELATION

By introducing parameter b in the Bayesian model, we are able to investigate the degree of correlation of the model simulations between future climate (yi) and current climate (xi). Figure 8 shows the distributions of b at 12 weather stations for three temperature variables and three future 30-yr periods, in the form of box plots. It is noticeable that most of the mass of the posterior densities is shifted in the right side of 0, with the mean of median values concentrated around 1.4. This implies a strong positive correlation between future climate projections and controlled climate simulation in each PRECIS experiment.

4) OBSERVATIONAL ERRORS

Figure 9 shows the posterior densities of l0, which is introduced here to quantify the uncertainties in observational data owing to system errors or natural variability, at 12 stations for three temperature variables. Similarly, we apply a logarithmically scaled axis in the box plots to help understand how the mass of the posterior densities is distributed. The plots for Tmax, Tmean, and Tmin show different patterns at 12 stations, but there is a common left-shifting tendency from Tmax to Tmean, then to Tmin while we compare the mass distributions for these three variables one station by one. Since l21 0 is referred to as the variance of observational data, bigger l0 suggests smaller variance and better quality of the observed data used in our analysis. Thus, these posterior distributions of l0 may give some insights into the quality of observational data. For example, the observation for Tmax used in this study is more likely to be concentrated on the true value of current climate with smaller deviation than the observed Tmin.

c. Maps of projected temperature changes

The purpose of this section is to display a number of maps to help understand the geographical patterns of temperature changes across the domain of

Figure 10 shows the projected changes of Tmax at different probability levels. The central estimates of change are projected to increase significantly with time periods, and so are the changes at 10% and 90% probability levels. Specifically, in the period of 2020-49, the central estimates of change are projected to be generally between 28 and 38C acrossmost of the province, with slightly larger changes in the northeast next to the

The changes of Tmin, as shown in Fig. 11, reveal similar patterns as those of Tmax, but with relatively smaller magnitudes compared to the changes of Tmax. The central estimates of changes in mean daily minimum temperature are projected to be [2, 3]8C in 2020- 49, [4, 5]8C in 2040-69, and [5, 7]8C in 2070-99, respectively. The changes in Tmin are very unlikely to be less than 18C (in 2020-49) and very unlikely to be greater than 108C (in 2070-99).

Figure 12 shows that, as with Tmax and Tmin, there is a temporally gradient increase in terms of changing magnitudes of daily mean temperature from 2020-49 to 2070-99. In addition, each map shows a gradient decreasing pattern for Tmean changes from the northeast, where the central estimates of changes in Tmean can be as high as 88C in 2070-99, to the southwest, where the projected changes in the same period at 50% probability level can be 68C or less. Such a spatially distributed pattern may be somewhat related to the large water bodies (i.e.,

5. Summary and conclusions

In this study, we developed a high-resolution RCM ensemble for the province of

The results show that there is likely to be a significant warming trend throughout the twenty-first century over

The Bayesian hierarchical model presented here is not restricted to the high-resolution PRECIS ensemble developed for the province of

Acknowledgments. This research was supported by the

REFERENCES

Allen, M. R.,

Barnett, D. N.,

Bellprat, O., S. Kotlarski, D. Leurouthi, and C. Scheuroar, 2012: Exploring perturbed physics ensembles in a regional climate model. J. Climate, 25, 4582-4599, doi:10.1175/JCLI-D-11-00275.1.

Benestad, R. E., 2004: Tentative probabilistic temperature scenarios for northern

Brooks, S., 1998: Markov chain

Collins, M.,

Furrer,

Giorgi, F., and

-, and

-, and -, 2003: Probability of regional climate change based on the Reliability Ensemble Averaging (REA) method. Geophys. Res. Lett., 30, 1629, doi:10.1029/2003GL017130.

Goldstein, H., and

Greene, A. M.,

Harris, G.

Hewitt, C. D., 2004: ENSEMBLES-based predictions of climate changes and their impacts. Eos, Trans.

Houghton, J. T., Y. Ding,

Jones, R. G.,

Karl, T. R., and

Kendon, E. J.,

Knutti, R.,

Krishnamurti, T.N., C.M. Kishtawal, T.E. LaRow,D.

-,-, Z. Zhang, T.E. LaRow,D.

Luo, Q.,

Manning, L. J., J.W.Hall,H.

McSweeney, C. F., and

-, -, and

Mearns, L. O., W. Gutowski,

-, and Coauthors, 2012: The North American Regional Climate Change Assessment Program: Overview of phase I results. Bull.

Meehl, G. A.,

Murphy, J.

-,

-, and Coauthors, 2009:

Nakicenovic, N., and

New, M., A. Lopez, S. Dessai, and

NLWIS, cited 2007: Daily 10 km gridded climate dataset: 1961- 2003 version 1.0. National Land and

-, 2011b: Climate ready:

-, 2011c: Climate progress:

Reuroaiseuroanen, J., and

Rayner,N.A., P. Brohan,D.

Robertson, A. W., U. Lall,

Sexton, D. M. H.,

Stainforth, D. A., and Coauthors, 2005: Uncertainty in predictions of the climate response to rising levels of greenhouse gases. Nature, 433, 403-406, doi:10.1038/nature03301.

Stott, P. A., and

Taylor, K. E.,

Tebaldi, C., and

-, and B. SansÓ, 2009: Joint projections of temperature and precipitation change from multiple climate models: A hierarchical Bayesian approach.

-,

Wang, X., and Coauthors, 2013: A stepwise cluster analysis approach for downscaled climate projection-A Canadian case study. Environ. Modell. Software, 49, 141-151, doi:10.1016/ j.envsoft.2013.08.006.

Watterson, I. G., and

Webb, M. J., and Coauthors, 2006: On the contribution of local feedback mechanisms to the range of climate sensitivity in two GCM ensembles. Climate Dyn., 27, 17-38, doi:10.1007/ s00382-006-0111-2.

Webster, M., and Coauthors, 2003: Uncertainty analysis of climate change and policy response. Climatic Change, 61, 295-320, doi:10.1023/B:CLIM.0000004564.09961.9f.

Wilby, R. L.,

Wilson, S.,

Yokohata, T.,

QIANGUO LIN

(Manuscript received

Corresponding author address:

E-mail: huang@iseis.org.

DOI: 10.1175/JCLI-D-13-00717.1

APPENDIX A

Bayesian Model Formulation

1) Likelihoods for existing data x0, xi, yi:

...

2) Prior distributions for parameters m, n, b, u, l0, l1, . . . , lN:

Assume uniform prior densities on the real line for m and n, uniform prior distribution on [21, 1] for b, respectively. For the remaining parameters, we assume gamma distributions as follows:

...

3) Posterior is given by, up to a normalizing constant,

...

Thus, we can obtain full conditional distributions for each parameter by ignoring all terms that are constant with respect to the parameter:

Full conditional for u:

...

Full conditional for b:

...

Full conditional for l0:

...

Full conditional for li, i51, . . . , N:

...

Full conditional for m:

...

Full conditional for n:

...

APPENDIX B

Gibbs Sampling Algorithm

Step 1. Pick a starting value for the Markov chain: (m, n, b, u, l0, l1, . . . , lN), say

...

Step 2. Update each parameter in turn:

(i) Sample a value of u from its full conditional distribution: p(u j m, n, b, l0, l1, . . . , lN, x0, x1, . . . , xN,

(ii) Sample a value of b from its full conditional distribution: p(b j m, n, u, l0, l1, . . . , lN, x0, x1, . . . , xN,

(iii) Sample a value of l0 from its full conditional distribution: p(l0 j m, n, b, u, l1, . . . , lN, x0, x1, . . . , xN,

(iv) Sample a value of li from its full conditional distribution: p(li j m, n, b, u, l0, l1, . . . , li21, li11, . . . , lN, x0, x1, . . . , xN,

(v) Sample a value of m from its full conditional distribution: p(m j n, b, u, l0, l1, . . . , lN, x0, x1, . . . , xN,

(vi) Sample a value of n from its full conditional distribution: p(n j m, b, u, l0, l1, . . . , lN, x0, x1, . . . , xN,

Step 3.Repeat step 2M21 times to produce a Markov chain of length M.

APPENDIX C

MCMC Simulation

We run the Gibbs sampling for a total of 260 000 iterations for all parameters at each 25-km grid specified by the PRECIS model. The first 10 000 iterations are treated as random drawings in the burn-in period during which the