Regionalised spatiotemporal rainfall and temperature models for flood studies in the Basque Country , Spain

A spatiotemporal point process model of rainfall is fitted to data taken from three homogeneous regions in the Basque Country, Spain. The model is the superposition of two spatiotemporal Neyman–Scott processes, in which rain cells are modelled as discs with radii that follow exponential distributions. In addition, the model includes a parameter for the radius of storm discs, so that rain only occurs when both a cell and a storm disc overlap a point. The model is fitted to data for each month, taken from each of the three homogeneous regions, using a modified method of moments procedure that ensures a smooth seasonal variation in the parameter estimates. Daily temperature data from 23 sites are used to fit a stochastic temperature model. A principal component analysis of the maximum daily temperatures across the sites indicates that 92 % of the variance is explained by the first component, implying that this component can be used to account for spatial variation. A harmonic equation with autoregressive error terms is fitted to the first principal component. The temperature model is obtained by regressing the maximum daily temperature on the first principal component, an indicator variable for the region, and altitude. This, together with scaling and a regression model of temperature range, enables hourly temperatures to be predicted. Rainfall is included as an explanatory variable but has only a marginal influence when predicting temperatures. A distributed model (TETIS; Francés et al., 2007) is calibrated for a selected catchment. Five hundred years of data are simulated using the rainfall and temperature models and used as input to the calibrated TETIS model to obtain simulated discharges to compare with observed discharges. Kolmogorov–Smirnov tests indicate that there is no significant difference in the distributions of observed and simulated maximum flows at the same sites, thus supporting the use of the spatiotemporal models for the intended application.


Introduction
Rainfall and temperature data are required in the study of hydrological systems -for example, in flood studies or in the analysis of urban drainage networks. However, historical records of data are always limited; for example, record lengths may be too short to predict high return period events, or data may be unavailable at sites of interest or only available at time scales that are too coarse for the intended application. Hence, stochastic models are used to simulate data to supplement or extend existing historical records; see, for example, Gyasi-Agyei (2005), Cowpertwait (2006), or Burton et al. (2008).
There is extensive literature on stochastic rainfall models that includes models based on spatiotemporal point processes similar to the model used here. One of the earlier models was based on a spatiotemporal Poisson process; see Cox and Isham (1988). This model does not explicitly allow for the clustering of rainfall events. Hence, most subsequent developments have allowed for the temporal clustering of rain cells using a Neyman-Scott or Bartlett-Lewis point process of cell arrival times. These build on the work of Rodriguez-Iturbe et al. (1987), who developed a methodology for the derivation of statistical properties that can be used in model fitting. For example, Northrop (1998)

developed a spatial
Published by Copernicus Publications on behalf of the European Geosciences Union. extension in which elliptical cells have occurrence times that follow a Bartlett-Lewis process, whilst Cowpertwait (1995) developed a spatiotemporal model with cell arrival times occurring in a Neyman-Scott process. Burton et al. (2010) extended the latter model to include a non-stationary mean number of rain cells, which can approximate orographic effects. Cowpertwait (2010) developed a spatial generalisation that allows a continuous distribution of storm types; this model also allows storms to have a defined spatial extent, to provide an improved fit to sample cross-correlations, and subsequently also extends the model in Leonard et al. (2008). The reader can find further information on the development of spatiotemporal point process rainfall models in Wheater et al. (2005), Cowpertwait (2010), Burton et al. (2010), and the references therein.
The rainfall model described here is the superposed discrete storm-type analogy to the continuous storm-type model used in Cowpertwait (2010), and it is more mathematically tractable and easier to implement. The model is a spatial extension to the discrete superposed temporal point process models discussed in Cowpertwait (2004) and Morrissey (2009). In the version used here, two independent superposed point processes are used to allow for different storm types (e.g. convective and stratiform rain). The methodology extends that used for the Thames study described in Cowpertwait (2006): to account for larger geographical regions through the inclusion of a storm centre and radius for each of the two superposed processes; a regionalisation procedure; and a temperature model to allow for further conditions that may give rise to flooding (e.g. snowmelt). A modified method of moments fitting procedure ensures that the parameter estimates vary smoothly across seasons. The models are intended for use in flood studies across the Basque Country, including studies for large catchment areas that require spatially representative series. In the last part of the paper, simulated rainfall and temperature data are input into the TETIS distributed catchment model (Francés et al., 2007) and observed and simulated annual maximum flows compared.

Superposed point processes
Independent stochastic point processes may be superposed (e.g. see Cox and Isham, 1988) to give additional model parameters and, therefore, increased flexibility in model fitting. For example, n independent Poisson point processes with rates λ 1 , λ 2 , ..., λ n may be superposed to give a Poisson process with overall rate λ = λ 1 + λ 2 + · · · + λ n . Equivalently, a Poisson process with rate λ may be the composition of n independent Poisson processes with rates λ i (i = 1, ..., n), such that a Poisson process selected at random has rate λ i with probability α i = λ i /λ, where n P i=1 α i = 1. When the Poisson processes each form the basis of an independent cluster process, such as a Neyman-Scott or Bartlett-Lewis point process, the additional parameters due to the clustering can be indexed to give different parameter sets for the different superposed cluster point processes (Cowpertwait, 2004;Morrissey, 2009). This can be generalised to a continuous probability density function α(x) with corresponding Poisson rates λ(x) (and other functions for any parameters used for clustering). The first approach ("discrete superposition") is the superposition of a countable set of point processes, whilst the second approach ("continuous superposition") is the superposition of an uncountable set of point processes. Continuous superposition lends itself to an analysis of the functional relationships between different model parameters, and it has been used to generalise a spatiotemporal point process model of rainfall (Cowpertwait, 2010). Discrete superposition has been used to model the temporal rainfall process (Cowpertwait, 2004;Morrissey, 2009) and is easier to implement in practice because most of the model properties of the independent processes can be summed to obtain those of the superposed process. Discrete superposition is suitable if storms are known to be of discrete distinct types. In this paper, discrete superposition is applied to the spatiotemporal Neyman-Scott model developed by Cowpertwait (1995) and Leonard et al. (2008). The superposed model is defined as follows.

Model definition
This is a straightforward extension of the previous model (Cowpertwait, 1995;Leonard et al., 2008) that essentially includes indices to allow for different storm types.
Let storm origins occur in a spatiotemporal Poisson process with rate ζ s per unit time per unit area. Suppose storms can be of n independent types and that each storm origin is of type i with probability α i (i = 1, ..., n).
Associated with each type i storm origin is a disc of random radius R s,i , where R s,i is an independent exponential random variable with parameter φ s,i . The arrival times {T ij } of type i storm origins at an arbitrary point in the plane occur in a temporal Poisson process with rate λ i (per hour), where λ i = 2 π ζ s,i /φ 2 s,i and ζ s,i is the spatiotemporal rate for a type i storm (e.g. see Cowpertwait, 2010). As outlined above, this process is equivalent to the superposition of n independent Poisson processes, and so the overall process is Poisson with Each type i storm consists of a marked point process of rain cells, denoted as {(U ij k , V ij k ), S ij k , L ij k , X ij k , R C ij k } for the j -th occurrence of a type i storm, where the marked process satisfies the following: a. {(U ij k , V ij k )} forms a 2-D Poisson process with rate ζ c,i (per unit area per storm); Hydrol. Earth Syst. Sci., 17, 479-494, 2013 www.hydrol-earth-syst-sci.net/17/479/2013/ b. (U ij k , V ij k ) and R C ij k form discs in 2-D space, where (U ij k , V ij k ) is the disc centre and R C ij k is the disc radius that is taken to be an independent exponential random variable with parameter φ c,i ; c. S ij k is the arrival time of the k-th cell in the j -th occurrence of a type i storm, where S ij k − T ij are independent exponential random variables with parameter β i ; d. L ij k is the lifetime of the k-th cell, which is taken to be an independent exponential random variable with parameter η i , so that the k-th cell in the j -th occurrence of a type i storm terminates at a time S ij k + L ij k ; and e. X ij k is a random variable representing the rain intensity (depth per unit time) of the k-th cell in the j -th occurrence of a type i storm, where X ij k has mean θ i and remains constant throughout the cell lifetime and over the area of the cell disc defined by Rain occurs at a spatiotemporal point (t, x, y) if, and only if, both a cell and storm disc overlap the point. The total intensity at spatiotemporal point (t, x, y) is then the sum of the intensities of all cells alive at time t that have discs overlapping (x, y) ∈ R 2 . Note that (c) above implies the cell arrival time point process {S ij k } is the superposition of n Neyman-Scott point processes.
For the purpose of model fitting and simulation, the cell intensities X ij k are taken to be independent exponential random variables each with survivor function P (X ij k > x) = e −x/θ i and moments given by E[X r ij k ] = r!θ r i (r = 0, 1, 2, ...). For each type i storm the number of cells C i that overlap a point in 2-D space R 2 is a Poisson random variable with mean ν i = 2 πζ c,i /φ 2 c,i (Cowpertwait, 1995). Furthermore, since the rate of storm arrivals is also related to the spatiotemporal rate of storms ζ s,i and the mean storm radius 1/φ s,i by λ i = 2 π ζ s,i /φ 2 s,i , it follows that both spatiotemporal rates, ζ s,i and ζ c,i , are functions of other model parameters and do not need to be fitted separately. Hence, the superposed spatiotemporal model is summarised by the following set of independent parameters: Although there are many possible storm types that could be envisaged, usually it is sufficient in practice to have just two storm types (n = 2 in the above) to broadly correspond to two distinct types of storms: convective and frontal systems. However, a further storm type can be associated with a spatial region -for example, to account for a third type of precipitation event caused by orography.

Properties used in model fitting
The statistical properties of the spatiotemporal rainfall model are given in the Appendix (Eqs. A2-A6). These properties are stationary, whilst the physical rainfall process is nonstationary in space and time. Non-stationarity in the physical process can be accounted for by fitting the model to discrete spatiotemporal intervals that are sufficiently small to be approximately stationary. For example, the model can be fitted to data taken over the period of a calendar month to account for seasonal changes in rainfall. Analogously, the model can be fitted to spatially homogeneous regions that are approximately stationary. We adopt this approach, using the homogeneous regions found in the earlier study of daily rainfall data from the Basque Country (Cowpertwait, 2011).
Based on the available model functions -see the Appendix for details -the following properties are used to fit the model: the mean rainfall, µ(h); the coefficient of variation, υ(h) = σ (h)/µ(h); the coefficient of skewness, κ(h) = ξ(h)/σ 3 (h); the autocorrelation, ρ(l, h) = γ (l, h)/σ 2 (h); the cross-correlation, ρ(d, l, h) = γ (d, l, h)/σ 2 (h); and the proportion of dry intervals of width h, ϕ(h). These functions are used at the hourly, six-hourly, and daily aggregation levels (h = 1, 6, and 24) in Eqs. (A2)-(A5) and at the daily level for the proportion dry Eq. (A6). The autocorrelation is used at lag l = 1 and the cross-correlation at lag l = 0. In summary, the set of properties for model fitting is

Historical data and sample estimates
The data for the project came from 357 sites: 123 records of hourly data for the period 1985-2010, and 234 records of daily data for the period 1914-2010 ( Fig. 1). About 73 % of the data were missing over this period, with most missing in the earlier part . (The fitted model can be used to fill in the missing data and disaggregate the daily data to hourly values; Cowpertwait, 2006.) Sample estimates of the properties for model fitting were found for each homogeneous region and for each month by pooling all available data for the month and region. A weighted average, based on the number of available observations, of the sample properties at the daily level of aggregation was calculated by combining the estimates from the daily data with those from the hourly series. For sample properties at aggregation levels smaller than 24 h (i.e. at the 1-and 6-h levels), estimates from the hourly series were used. In summary, the following set of sample estimates was found for each month and region:

Modified method of moments
For each homogeneous region, the model parameters are estimated using a modified method of moments procedure that minimises the sum of squares:  (Cowpertwait, 2011)): hourly data (), daily data (lightly shaded), temperature data (•) and gauges used to calibrate the TETIS catchment model ().
where g is a model function andĝ is the equivalent sample estimate taken from the historical data. The model parameters are estimated for each homogeneous region in three stages.

Stage 1: estimation of temporal parameters based on dimensionless properties
Initially, G is F minus the mean and cross-correlation ρ(d, 0, h), so that dimensionless temporal properties are fitted to estimate the parameters {λ i , ν i , β i , η i : i = 1, 2} for 2 storm types. When minimising Eq. (1), it is likely that there are a number of local minima so that two sets of parameter estimates may give similar fits to the sample properties (e.g. see Vanhaute et al., 2012). This can result in the estimates not following a smooth variation over seasons, which is particularly undesirable in a regionalisation where comparisons between the estimates from different regions (and months) are required. Hence, the procedure below is a modification of the minimisation procedure used in previous work (e.g. Cowpertwait, 2006), to ensure a smoother variation of the parameter estimates over seasons whilst retaining a good fit to the sample properties. Fitting to the dimensionless temporal properties is carried out in a series of steps as follows.
a. At the very least the parameter estimates must all be greater than zero, and so a constrained minimisation of Eq.
(1) is necessary. So the first step in fitting is to place wide arbitrary bounds on the parameters and then minimise Eq. (1). b.
Step (a) is repeated for each month to gives twelve estimates for each of λ i , ν i , β i , and η i (i = 1, 2). c. For each parameter, the estimate for the j -th month is adjusted to be the mean of the equivalent estimate of adjacent months (e.g. the January estimate of λ 1 is the mean of the February and December estimates of λ 1 ).
d. The parameters are re-estimated using the adjusted estimates from (b) as starting values when minimising Eq.
e. Steps (b)-(c) are iterated n times to give a sample of estimates for each parameter.
f. For each parameter, the median of the n iterated estimates from step (d) is used as the final estimate.

Stage 2: estimation of spatial parameters
Using the estimates of {λ i , ν i , β i , η i : i = 1, ..., 2} obtained in Stage 1 above, the cell and storm radii parameters (φ c,i and φ s,i ) are estimated for each month by minimising Eq. (1) with G = {ρ(d, 0, h) : h = 1, 6, 24; d ∈ D}, by pooling data from all available pairs of sites in the region when calculating the sample cross-correlations and set D of corresponding distances. The estimates for the spatial parameters can be taken to be the same in each region (so the estimates only vary with season), in which case : h = 1, 6, 24; d ∈ D j } and D j is the set of distances for all pairs of sites in the j -th region (j = 1, 2, 3). The estimates of {λ i , ν i , β i , η i : i = 1, ..., 2} for the j -th region are used in the calculation of G j = {ρ(d, 0, h) : h = 1, 6, 24; d ∈ D j }.

Stage 3: estimation of the scale parameter
Finally, the scale parameter θ is taken to be the same for each storm type and is estimated for each month from the sample mean hourly rainfallμ(1) using the relation: This can vary for different sites or, when the regions are clearly homogeneous, can be kept the same over a whole region. Sinceθ is a function of the mean rainfall, it can be adjusted on a site-by-site basis when good estimates of the site means are available.

Fitted rainfall model
Using the fitting procedure in Stage 1 above, a range of different parameter sets were used to fit two storm types to the dimensionless sample properties from each region. It was found that the mean number of cells per storm per site could be kept the same for the different storm types (i.e. ν 1 = ν 2 ) Hydrol without a loss in overall goodness-of-fit to the sample properties. The fit to the dimensionless properties gave estimates of λ 1 , β 1 , η 1 , λ 2 , β 2 , η 2 and ν (= ν 1 = ν 2 ), for each month and each region, which are shown in Fig. 3 and tabulated for the central region in Table 2. To illustrate the iterations from Step (d), plots are given for n = 50 iterations for the January estimates in the central region, where correlation between pairs of estimates, probably due to local minima, is evident as the peaks tend to occur at the same iteration (Fig. 2). Following Stage 2 in Sect. 3.3, φ c and φ s were estimated and are given in Table 2. The fitted and sample properties are given in Figs. 4-6.
Initially, it seemed appropriate that the scale parameter θ, which is a function of the mean rainfall (Eq. 2), should depend on altitude as found in previous studies (e.g. Cowpertwait, 2006). So the mean rainfall was found for each month for each of the 357 sites ( Fig. 1) and regressed on altitude and indicator variables for the month and region. The regression models were fitted using weighted least squares with weights corresponding to the number of observations used to calculate the mean rainfall. The results indicated that altitude was not needed as an explanatory variable when the region and month were included in the model (Table 1). (Even a regression model with 72 interacting terms only improved R 2 by 2 % when altitude was included.) This result provided further support for the regionalisation procedure (Cowpertwait, 2011). The scale parameter was therefore treated as the same within a homogeneous region and is shown plotted in Fig. 3. However, as noted in Stage 3 of the fitting procedure, this parameter can be estimated directly from the site mean  rainfall when this is available, which is the approach used in Sect. 5. The resulting parameter estimates reflect some expected characteristics of the rainfall process. For example, over the summer months, there is a general increase in both the cell intensity parameter estimate (θ) and the cell duration parameters (η), whilst there is a general decrease in the number of cells (ν) per storm (Fig. 3). This corresponds with the expected increase in summer convective storms that tend to have fewer, shorter duration (1/η) raincells, but of higher intensity (Fig. 3).
Differences between the two fitted storm types can also be discerned. For example, type 1 storms are less frequent than type 2 storms because, in general,λ 1 <λ 2 (Fig. 3). In addition, type 1 storms have cells that are less clustered (β 1 <β 2 ) and of longer duration (η 1 <η 2 ). Type 1 storms are therefore more representative of stratiform rain, which would tend to be more persistent. It should, however, be mentioned that these classifications of storm types are partially for additional flexibility in the model parameterisation, to obtain a good fit to the data, and that properties similar to those observed for convective storms could result, by chance, in simulations for any storm type in the model (because all the variables in the model are random). However, a general tendency for storm characteristics of a particular type is ensured through this classification and is supported by observing the clear distinction in the resulting estimates for the different storm types (Fig. 3).
There are also regional differences in the parameter estimates. For example, the north-eastern (green) region has a consistently higher rate of storm arrivals (λ), whilst the central (red) region has a higher number of cells per storm during the winter months (lower left-hand plot in Fig. 3). Also, whilst the north-eastern region has a tendency to experience more storms, over the summer the cells are generally less intense, as indicated by a lower value ofθ (lower right-hand plot; Fig. 3). If λ 1 is interpreted as a rate for stratiform storm occurrence, then the north-eastern region receives the highest rainfall of this type. The southern region has notably fewer storms over the summer months, having the lowest value of λ 2 (upper right-hand plot; Fig. 3). These regional differences are most likely to be due to the presence of both the ocean to the north and the Cantabrian Mountains separating the northern and southern regions.
In general, the fit to the sample properties is very good (Figs. 4 and 5). Again, seasonal and regional difference are evident in the figures. For example, there is a notable difference between the southern (blue) and north-eastern (green) regions: the southern region has lower values of autocorrelation but notably higher values of skewness, the proportion dry, and the coefficient of variation, especially over the summer months (Figs. 4 and 5). This corresponds to fewer more intense events (lower λ 2 , ν, higher θ; Fig. 3) characteristic of summer convective storms. However, a clear correspondence between any selected parameter estimate and the sample estimates is generally lacking, but this is essentially due    of constructing a full spatiotemporal temperature model that can be used to simulate multisite hourly temperature series at any location in a region.

470
There were 23 sites with (near) complete records of daily maximum temperatures, and total rainfall, for the period 1985-2010 (Fig. 1). This was used to form a matrix of 9490 rows and 23 columns of maximum daily temperatures, from which the principal components of the daily temperature data 475 were extracted (based on the correlations between the sites). In addition, the principal components of the daily rainfall totals were also found and the percentage of variance associated with each component extracted (Table 3).
The first principal component accounts for 92% of the 480 variance in the temperature data, whilst the second component only accounts for 2.4% of the variability (Table 3). This high first value, followed by a subsequent low value, indicates that the spatial variation of the temperatures is mainly accounted for in a single component. In contrast, daily rain-485 fall has a lower first component of 71% followed by values that are higher than the corresponding values for temperature (Table 3). This is due to the rainfall process hav- to the model properties (Eqs. A2-A6) being functions of all the model parameters.
The fit to the sample cross-correlations is given in Fig. 6, where it can be seen that the curves decay with distance, as expected. The spatial estimates (φ c andφ s ) are the same for all three regions, so the slight difference in the curves for the different regions is due to the differences in the temporal parameter estimates obtained in Stage 1 of the fitting procedure (and given in Fig. 3).
Some exceptions to the general goodness-of-fit may be observed -the most notable being an over-estimation in the proportion of dry days in the summer months for the central and southern regions (lower right-hand plot in Fig. 5). A dry period requires an arbitrary definition, since the model may generate very small values (and gauges usually require a small accumulation of rainfall before tipping). In this analysis, an observed daily rainfall of less than 0.1 mm is taken to be dry -an amount that could easily be lost due to evaporation. The distribution of daily rainfall is largely determined by the first three sample moments. These are closely matched by the model in the central region; hence, the lack-of-fit to the proportion dry may be of small practical significance.
It is useful to consider the fit to observed extreme values, since these are not used in the fitting procedure and are important when simulating high flows. Using the estimates in Table 2, one thousand years of hourly rainfall were simulated for the daily sites in the central region. These were aggregated to daily values and the annual maxima found for each site. The median of the annual maxima across the sites was then calculated; these were ordered in 20-year blocks for simulated and observed   ing more between-site spatial variability than temperatures. Consequently, a spatial temperature model need not contain 490 the same level of complexity as the spatial-temporal rainfall model of §2.2. The first principal component was therefore used as the basis for a spatial-temporal temperature model. Harmonic curves were fitted to the first principal component to account for seasonality. The i-th harmonic is given by: z t = sin(2iπt/365 + φ) = s i sin(2iπt/365) + c i cos(2iπt/365), where t is time measured in days and z t is the first principal component score of the maximum daily temperatures. The first harmonic was fitted by least squares and is given by: where a t is the residual error series. (Two harmonic models were fitted to the first component score of maximum daily 495 temperature by least squares regression. The first, given above, contained only the first harmonic whilst the second contained the first 20 harmonics and was thus notably more  against the standardised Gumbel variate. (Fig. 7). Although some lack-of-fit can be seen in the very largest value, where the observed value exceeds the simulated values, in general the results indicate that the model performs satisfactorily with respect to the extremes, because the observed values fall within the range of simulated values (Fig. 7). This, together with the goodness-of-fit to properties up to third order (Figs. 4 and 5), indicates that the fitted model may be suitable for flood studies. In Sect. 5, the model is further validated against properties, such as peak flow discharges, that are important in the intended application.

Models for the first principal component
In weather generators, stochastic temperature models have been related to rainfall (e.g. see Kilsby et al., 2007). Although a rainfall variable is used in the following, our results indicate that it is a poor predictor of temperature for the Basque Country, and so we consider an alternative approach, based on principal components, that can also be used to generate multisite temperature series. Principal components have been used in other studies of the spatial variation in temperature series, e.g. see Benzi et al. (1997). The approach described herein differs in that it provides a method of constructing a full spatiotemporal temperature model that can be used to simulate multisite hourly temperature series at any location in a region. There were 23 sites with (near-) complete records of daily maximum temperatures, and total rainfall, for the period 1985-2010 (Fig. 1). This was used to form a matrix of 9490 rows and 23 columns of maximum daily temperatures, from which the principal components of the daily temperature data were extracted (based on the correlations between the sites). In addition, the principal components of the daily rainfall totals were also found and the percentage of variance associated with each component extracted (Table 3).
The first principal component accounts for 92 % of the variance in the temperature data, whilst the second component only accounts for 2.4 % of the variability (Table 3). This high first value, followed by a subsequent low value, indicates that the spatial variation of the temperatures is mainly accounted for in a single component. In contrast, daily rainfall has a lower first component of 71 % followed by values that are higher than the corresponding values for temperature (Table 3). This is due to the rainfall process having more between-site spatial variability than temperatures. Consequently, a spatial temperature model need not contain the same level of complexity as the spatiotemporal rainfall model of Sect. 2.2. The first principal component was therefore used as the basis for a spatiotemporal temperature model.
Harmonic curves were fitted to the first principal component to account for seasonality. The i-th harmonic is given by z t = sin (2 i πt/365 + φ) = s i sin (2 i π t/365) + c i cos (2 i π t/365), where t is time measured in days and z t is the first principal component score of the maximum daily temperatures. The first harmonic was fitted by least squares and is given by z t = 2.1 sin (2 π t/365) + 4.7 cos (2 π t/365) + a t , where a t is the residual error series. (Two harmonic models were fitted to the first component score of maximum daily temperature by least squares regression. The first, given above, contained only the first harmonic, whilst the second contained the first 20 harmonics and was thus notably more complex. When compared to the simpler model, the more complex harmonic model only improved the adjusted R 2 by 2 %, and so the simpler model above was selected in preference.) A plot of the first ten years of the fitted values against the observed scores is shown in Fig. 8, where it is evident that the harmonic model successfully follows the seasonal variation in the first principal component. There is evidence of random variation about the curve, and so a stochastic component is required for a t . A best-fitting autoregressive model, based on maximum likelihood, was therefore found and is given by a t = 0.77 a t−1 − 0.13 a t−2 + 0.074 a t−3 + 0.039 a t−4 + w t , (4) where w t is the residual error series.
The correlogram of the residuals of the fitted AR(4) model indicates that the model successfully accounts for the autocorrelation in the residual series of the fitted harmonic model ( Fig. 9). Furthermore, there is no evidence of any persistence or seasonal variation, indicating that the harmonic model combined with the AR(4) model provides a good fit to the first principal component score of the data (Fig. 9). The sample variance of the AR(4) residuals, which is needed in simulations, is 3.7.

Spatial daily temperature model
The maximum daily temperature at each site was regressed on the following variables: altitude, region (as an indicator variable), the first principal component score (PC1; z t ), and daily rainfall. The fitted regression model is shown in Table 4. All variables are statistically significant (Table 4). However, some of the variables are statistically significant because of the large numbers of observations, but they are not of practical significance. For example, if the rainfall variable is removed, there is no change in R 2 to four significant figures, implying that the rainfall variable can be left out when predicting maximum daily temperatures. (This is also reflected in the very low coefficient for the rainfall variable; Table 4.) Overall, the fitted model explains more than 90 % of the variance in the daily temperatures, which implies the model provides a good fit to the data. In addition, as expected, the model predicts lower temperatures at higher altitudes and higher temperatures in the south (and marginally lower temperatures for higher rainfall).
Maximum daily temperatures can therefore be simulated at any site by first simulating a series of autoregressive terms (Eq. 4), adding these to the harmonic equation for the first principal component (Eq. 3) and then using the simulated first component (PC1) as a predictor in the regression model (Table 4). Rainfall can probably be left out of this procedure without any practical effect on the results.

Hourly temperature model
For this part of the study hourly temperatures were required. From the available temperature data, most of the data were at the daily level. Nevertheless, a total of 45 607 values of hourly temperatures were available for model fitting.
Let T max be the maximum daily temperature. Hourly temperatures can be related to this using the following equation: where T t is the temperature at time t (in hours; t = 1, ..., 24), A is the temperature range over the day and h t is the mean scale factor at time t (0 ≤ h t ≤ 1). Note that the scale factor h t ensures the maximum daily temperature is retained. The scale factor h t is taken to be constant for different months and sites and is estimated by taking the mean of (T max − T t )/A for each hour over all the records. This does result in the maximum temperature occurring at the same time each day, but this is unlikely to be of practical importance in the simulations. The estimates for h t are as follows: 0.906, 0.930, 0.953, 0.974, 0.993, 1.000, 0.966, 0.879, 0.724, 0.520, 0.318, 0.154, 0.0457, 0, 0.0195, 0.0932, 0.208, 0.352, 0.502, 0.633, 0.731, 0.797, 0.843, 0.878 (in order from the first hour, 00:00-01:00 UTC, after midnight).
The range A was regressed on the maximum daily temperature, month (as an indicator variable), and the daily rainfall (Table 5). The regression model for predicting the temperature range A has a lower R 2 compared to the model for predicting maximum daily temperature (Table 5). Nevertheless, it provides a method for distributing the hourly temperatures over a day, conditional on the predicted maximum temperatures, and may be adequate in practice. As expected, the model predicts a greater range of temperatures when the maximum temperature is high, which would typically occur during the summer months (see the coefficient for the maximum temperature in Table 5). Also, the predicted range (when the explanatory variable for maximum temperature takes the same value in both winter and summer) is greater in winter, as would be expected for a maritime climate with higher temperature drops at night during winter months. For example, the coefficient for February is 0.62, compared to −4.45 for August, indicating a lower range for August and, hence, a higher temperature drop in February (Table 5).

Urola catchment
In forthcoming projects, the fitted rainfall and temperature models will be coupled with a hydrological catchment model to simulate long flow series for the purpose of evaluating design discharges in flood studies. This approach removes the need to explicitly assign antecedent conditions, such as soil moisture, which are a significant source of uncertainty (Michele and Salvadore, 2002;Boughton and Droop, 2003;Camici et al., 2011). In the following case study, this approach is adopted using the spatiotemporal rainfall and temperature models to simulate data for input to a distributed catchment model. A flood risk assessment is needed for the Urola basin, which requires the establishment of different design discharges. The catchment is located in the north of the Iberian Peninsula. The region is characterised by an oceanic climate, which is humid and temperate without a dry season. The mean annual rainfall in the catchment ranges from 1200 to 1600 mm, while the mean annual temperature varies between 11.5 • C in the upper part of the valley and 13.5 • C in the lower part. Usually rainfall occurs due to the advection of North Atlantic fronts coming from the north-west and hitting the slopes of the Cantabrian Mountains, located only 30-60 km from the coast. This results in uniform and moderate precipitation. However, persistent and very intensive convective phenomena can also take place due to the combination of polar air and high sea surface temperatures, leading to heavy rainfall and flooding. The catchment covers a surface of 342 km 2 running from south to north along 65 km of main stream. Upstream, close to the Aitzkorri Sierra, the valley is narrow and steep (average slope of 1.5 %). Downwards and after receiving its main tributary, the Ibai-Eder River, the valley widens and flattens before reaching the Atlantic (Fig. 10). There are three gauges in the catchment that are known to provide reliable data (Fig. 10). However, the available record lengths, which range from 13 to 22 yr, are too short to confidently extrapolate high return period events that are needed to design protection measures against flood damage. Consequently, the Urola basin is a good example for a case study.

Calibration of the distributed catchment model
The TETIS model was selected to simulate the hydrological processes in the River Urola basin (Francés et al., 2007). The model is a conceptual distributed catchment model that divides the catchment into square cells, each characterised by six tanks linked vertically. Each tank represents the different water storages in the terrestrial phase of the hydrological cycle. Water flows downstream from each cell until reaching the river channel. A more detailed description of the model is available in Francés et al. (2007) or Vélez et al. (2009).
Three main parameters are responsible for the model output: (1) the static storage capacity, which controls the amount of water lost due to evapotranspiration in the long-term and the initial losses in the case of a flood; (2) the soil hydraulic conductivity, which affects the amount of infiltration and the velocity of interflow discharge; and (3) the subsoil hydraulic conductivity, which determines the amount of percolation and the velocity of the base flow discharge. These parameters have been estimated for the region and the model shown to successfully predict surface water runoff (Vélez et al., 2009).
There are five variables that define the initial conditions for each simulation, corresponding to the initial level in each tank. To avoid over-parameterising the model, following Beven and Binley (1992), values are obtained for the initial conditions by simulating a sufficiently long antecedent period; soil moisture conditions for long daily simulations were used (in part to reduce computational times) as initial states in the hourly simulations. In addition to these variables, the model has nine further correction factors that are used to adjust the way the three main parameters affect the hydrological processes and the value of other variables needed in simulations, such as the evapotranspiration and the overland and channel flow velocities. These correction factors are calibrated using a set of recorded events. The Nash-Sutcliffe efficiency coefficient (or R 2 ) is used to assess goodness-of-fit. In general, an acceptable value of R 2 is 0.6 and a value exceeding 0.8 is indicative of an excellent fit (Pappenberger and Beven, 2004). A unique set of correction factors was sought for 13 flood events (discussed in the next section). When this is achieved, assuming high overall R 2 values are maintained, the fitted model should be of practical value in flood studies. Furthermore, as the model is distributed and relies on three maps of parameters that have been estimated along the whole basin, discharges can be estimated at ungauged locations (Beven, 1985). Calibrating the model to heavier events helps ensure that a satisfactory fit is obtained to annual maximum discharges, which are important in the intended application. Although this may result in some reduced goodnessof-fit for smaller events, continuous daily simulations of 9 historical years gave R 2 values of 0.75, 0.69, and 0.80 for B1Z1, B1Z2, and B2Z1, respectively (similar to those obtained by Vélez et al., 2009), indicating an overall satisfactory fit (in particular, with respect to soil moisture conditions which are important in the simulation of peak discharges).

Fitted TETIS model
From 2001, there were 13 flood events in the Urola catchment, with data available from three gauges (Fig. 10). These events form the basis of the calibration used here; their main features, which are given in Table 6, cover a range of different antecedent conditions and rainfall intensities. Maximum and minimum temperatures are also available from those gauges, enabling the estimation of potential evapotranspiration by means of the Penman-Monteith equation and the simplifications suggested by the FAO (Allen et al., 1998). Starting with a multiple-event automatic calibration process based on the SCE-UA algorithm (Duan et al., 1994) and mean R 2 as an objective function, the nine correction factors in TETIS were adjusted to match the flows at each of the three sites, taking into account both R 2 and the absolute peak error (|sim − obs|/obs, where "sim" and "obs" are abbreviations for simulated and observed flows respectively). The Aiztu (B1Z1) and Ibai-Eder (B1Z2) sites were calibrated first and simulated series from these sites used to calibrate the inter-watershed between them and the Aizarnazabal (B2Z1) site. Based on the three optimum models, a regionalisation procedure was used to obtain a unique set of correction factors that could be applied to the whole catchment and set of observed events. An iterative process was used to reach optimum values of the correction factors for the whole basin to ensure an overall goodness-of-fit, without a significant reduction in R 2 . Hence, the correction factors used in this study differed from those in Vélez at al. (2009); for example, in this study a single regionalised correction factor for overland flow (0.06) was used for the whole basin, which differed from the approach in Vélez et al. (2009) in which model calibration was for the most downstream station, thus leading to poorer model performance at upstream sites. These differing approaches are due, in part, to different objectives: the objective here is to produce models for flood studies, whilst the objective in Vélez et al. (2009) is to produce models for water resources.
The absolute peak errors and R 2 values obtained in the model calibration are given at the bottom of  poor, probably due to a deficient rainfall representation, the overall fit is good with a median R 2 of 0.77, 0.78 and 0.64 at sites B1Z1, B2Z1 and B1Z2, respectively, which are above the acceptation threshold. The absolute peak errors are also satisfactory. In general, taking into account that 13 events have been used in calibration, the model can be regarded as suitable for representing the hydrological response of the catchment during major flood events (Brath et al., 2004).
As an example, Fig. 11 gives the simulated and observed hydrographs for the six largest flood events recorded at the Aizarnazabal (B2Z1) gauge.

Flow simulations and validation of spatiotemporal models
The Urola catchment overlaps the central and north-eastern homogeneous regions (Fig. 1). However, the majority of the catchment, and the upper part of the catchment that contributes the most to the flows, is contained in the central region. Hence, the parameter estimates for the central region were selected ( Table 2). The estimates of θ (Table 2; Sect. 3.3, Stage 3) were multiplied by 1.08, 1.04 and 1.08 for sites B1Z1, B2Z1 and B1Z2, respectively, to achieve exact fits to the annual mean rainfalls at each site (this is equivalent to just scaling simulated rainfall series at these sites by these factors). Five hundred years of hourly multisite rainfall and temperature data were simulated at the three sites. The simulated rainfall and temperature series were then used as input to the calibrated distributed catchment model to create a series of simulated discharges for each of the three sites. The annual maximum flows were extracted from the simulated series and compared with the annual maximum flows recorded at the gauges for each site. Summary statistics for the distributions of the maximum flows are given in Table 7 and quantile plots in Fig. 12.
In general, the distributions of the simulated annual maximum flows compare favourably to distributions of the observed maximum flows (Fig. 12). Some over-estimation in the mean (and median) maximum flow for site B1Z2 is evident in Fig. 12 and in Table 7, and a t-test of the difference in the mean values gave a possibly significant result with a p value of 0.08. A Kolmogorov-Smirnov test statistic for the differences in distributions of maximum flows was found for each pair of distributions. This provides a more formal test of the differences seen in Table 7 and Fig. 12. The test statistic and p values are given in Table 8. Table 8 shows that there is no statistical evidence of a difference in the distributions of the annual maximum flows for data at the same sites (shown in bold type). The table also shows that the test is capable of discerning differences in the distributions at different sites; e.g. compare the results for the observed flows at B1Z1 and B2Z1. Some of the differences are possibly significant; e.g. the test statistic for the observed flows at sites B1Z1 and B1Z2 had a p value of 0.088 (Table 8). This is probably because the sites are the closest and have the same scale factor and, hence, the same ensemble temporal rainfall properties. Whilst it is interesting to note variations in the distributions at different sites, the most important validation statistic is the measurement of the difference in simulated and observed distributions at the same site, for which there is no statistical evidence of a difference (Table 8).

Conclusions
In summary, spatiotemporal rainfall models were fitted to data from three homogeneous regions in the Basque Country. In general, good fits were obtained to sample properties of the observed rainfall series. The first principal component explained 92 % of the variability in daily temperatures and was hence used as predictor in the spatiotemporal temperature model. A distributed model was calibrated for the Urola catchment. Using the spatiotemporal models, hourly rainfall and temperature data were simulated for three sites in the catchment and flows generated using the simulated data as input for the catchment model. The distribution of observed annual maximum flows, taken from the site gauges, compared favourably to simulated maximum flows. The models are therefore validated for the catchment and can be used with confidence in further studies; these will include obtaining high return period discharges for river networks in the Basque Country that will be used to predict flood risks within the European Flood Directive framework.

Statistical properties of the rainfall model
In model fitting it is usually necessary to use equations for aggregated properties, because rainfall data are usually sampled over discrete time intervals (or aggregated to discrete time intervals in the case of tipping bucket data). Let (h) ij (x)} be the aggregated time series of rainfall due to type i storms at point x = (x, y) ∈ R 2 in the j -th time interval of width h, and let Y (h) j (x) be the total rainfall in the j -th interval due to the superposition of the n storm types. Then,

{Y
where Y i (x, t) is the rainfall intensity at point x and time t due to type i storms (i = 1, ..., n). Since the superposed processes are independent, statistical properties of the aggregated time series follow just by summing the various properties that were derived by Cowpertwait (1995Cowpertwait ( , 1998, Leonard et al. (2008), and Rodriguez-Iturbe et al. (1987), and are given below (Cowpertwait, 1995(Cowpertwait, , 1998. where Q 1 and Q 2 are high-order polynomials in η i and β i and are given in Cowpertwait (1998). For exponential cell intensities, E(X 2 ij k ) and E(X 3 ij k ) are replaced by 2 θ 2 i and 6 θ 3 i respectively. The probability that an arbitrary time interval [(j − 1)h, j h] is dry at a point is obtained by multiplying the probabilities of the independent processes and is given by the following: where (Cowpertwait, 1995;Eqs. 2.17 and 2.19).