Estimation of daily streamflow time series is of paramount importance for the design and implementation of river engineering and management projects (e.g., restoration, sediment-transport modelling, hydropower). Traditionally, indirect approaches combining stochastic simulation of rainfall with hydrological rainfall–runoff models are used. However, these are limited by uncertainties in model calibration and computational expense. Thus, this paper demonstrates an alternative, direct approach, for stochastic modelling of daily streamflow data, specifically seeking to address well-known deficiencies in model capability to capture extreme flow events in the simulated time series. Combinations of a hidden Markov model (HMM) with the generalised extreme value (HMM-GEV) and generalised Pareto (HMM-GP) distributions were tested for four hydrologically contrasting catchments in the UK (Rivers Dee, Falloch, Caldew and Lud), with results compared to recorded flow data and estimations obtained from a simpler autoregressive-moving-average (ARMA) model. Results show that the HMM-GP method is superior in performance over alternative approaches (relative mean absolute differences (RMAD) of <2% across all catchments), appropriately captures extreme events and is generically applicable across a range of hydrological regimes. In contrast, the ARMA model was unable to capture the flow regime successfully (average RMAD of 14% across all catchments).

## INTRODUCTION

In recent years, ensuring the sustainability of the water resources and aquatic habitats of our river systems has become integral to legislative, management and engineering practices (European Commission 2000). Implicitly, this requires analysis of the long-term natural fluctuations (seasonality, floods, droughts and climate change) in streamflow records via backward-facing (historical data) and forward-facing (modelling) evidence-based approaches. This permits, for example: decisions on abstraction licensing or design of freshet release strategies to ensure appropriate flows for ecosystem productivity and pollution assimilation capacity; dynamic simulation of sediment/morphological response to variable flows for scheduling dredging intervention associated with flood risk management; feasibility assessments for design and operation of hydropower schemes. As many UK rivers now offer flow gauge data over 50 years, informed sustainable river management is largely dependent on the development of robust tools which allow for the generation of synthetic streamflow sequences (past and projected) and capture the natural variability in flow. For many such projects, common practice is to provide stochastic representations of climate information (e.g., rainfall) and use a hydrological rainfall–runoff model to estimate the corresponding flows in the river channels. Although many of these models have been shown to produce reliable results with ever increasing computationally efficiency (e.g., Mellor *et al.* 2000; Dechmi *et al.* 2012; Zhang *et al.* 2013), depending on the catchment of interest, pre-existing models may be unavailable. Should this be the case, the computational time and data required for accurate calibration and subsequent estimation of (potentially multiple) flow sequences can be extensive and, therefore, costly. An alternative methodology combines stochastic models and gauged flow data to provide a more efficient procedure for the generation of multiple daily streamflow sequences. This requires less data, is much quicker to implement at multiple sites and is, therefore, receiving increased attention from the hydrological research community.

One of the simplest, and most common, methods for stochastic simulation of river flows is the use of an autoregressive-moving-average (ARMA) model. While these have typically been applied to generate sequences of monthly- and annually-averaged flow (O'Connell 1974; Mohammadi *et al.* 2006; Sadek 2006; Musa 2013; O'Connell & O'Donnell 2014), the potentially large day-to-day variation in flow values means that the use of ARMA models for simulating daily streamflow sequences is far less common. Although there have been recent studies (Martins 2011; Can *et al.* 2012) these were focussed on international rivers with large (>20,000 km^{2}) catchments and a much stronger seasonality (wet and dry seasons) than occurs in the UK. To the best of the authors' knowledge their use for smaller catchments or UK rivers has not been robustly tested.

A more complex approach, also employed in modelling daily streamflow, is the use of Markov models or variations (Yapo *et al.* 1993, Xu *et al.* 2001, 2003, Szilagyi *et al.* 2006, Augustin *et al.* 2008). In their case study of the River Tees, UK, Augustin *et al.* (2008) proposed the use of a Markov model within a wider, multinomial logit, modelling framework to generate daily streamflow. However, while the underlying transition of the flow values was described by a Markov type model, the estimation of state transition probabilities is complex (with a heavy dependence on a range of associated climate variables observed at time *t* and in the past) and aligns more with regression modelling rather than a traditional Markov model (MM). Also, as the precise flow rates within each state were estimated by random sampling of values from the associated continuous distribution, sampling bias towards values of high probability is possible; hence introducing errors and uncertainties into the approach.

To overcome these shortfalls, the current paper uses a hidden Markov model (HMM) based solely on statistical properties of daily averaged streamflow data. This approach to flow series simulation implicitly includes a more robust sampling methodology, via systemic selection of all values based on their probabilities. The flexible structure of the model permits comparison of a traditional HMM combined with the generalised extreme value GEV (HMM-GEV) or generalised Pareto (HMM-GP) distributions to better replicate the entire range of flows. As all previous studies have been restricted to single test catchments or comparison of sites with low variability in the flow regime, a further objective of the paper is to provide a robust assessment of the HMM's generic suitability to a number of catchments exhibiting a range of hydrological regimes. To achieve this, development of a meaningful, non-subjective and common method for partitioning the flow records also requires consideration.

## STUDY CATCHMENTS

River . | Record length (years) . | Base flow index . | Q_{max%} (m^{3}/s)
. | Q_{10%} (m^{3}/s)
. | Q_{50%} (m^{3}/s)
. | Q_{70%} (m^{3}/s)
. | Q_{95%} (m^{3}/s)
. |
---|---|---|---|---|---|---|---|

Dee | 83 | 0.53 | 648.50 | 73.91 | 26.52 | 18.49 | 8.29 |

Falloch | 42 | 0.19 | 123.60 | 16.22 | 2.17 | 1.00 | 0.27 |

Caldew | 32 | 0.49 | 96.39 | 10.10 | 2.69 | 1.66 | 0.75 |

Lud | 46 | 0.90 | 5.64 | 0.89 | 0.35 | 0.24 | 0.13 |

River . | Record length (years) . | Base flow index . | Q_{max%} (m^{3}/s)
. | Q_{10%} (m^{3}/s)
. | Q_{50%} (m^{3}/s)
. | Q_{70%} (m^{3}/s)
. | Q_{95%} (m^{3}/s)
. |
---|---|---|---|---|---|---|---|

Dee | 83 | 0.53 | 648.50 | 73.91 | 26.52 | 18.49 | 8.29 |

Falloch | 42 | 0.19 | 123.60 | 16.22 | 2.17 | 1.00 | 0.27 |

Caldew | 32 | 0.49 | 96.39 | 10.10 | 2.69 | 1.66 | 0.75 |

Lud | 46 | 0.90 | 5.64 | 0.89 | 0.35 | 0.24 | 0.13 |

### River Dee

The Dee originates in the Cairngorm mountain range at an elevation of 1,309 m above ordnance datum (AOD) and flows east through Braemar, Ballater, Aboyne, Banchory and Aberdeen before discharging into the North Sea. The size and location of the catchment means that the geography and land use varies considerably from source to mouth; with regions of mountains, moorland, forestry and arable farm land.

Like many rivers in this part of the UK, the Dee is well-known for salmon fishing. This attracts anglers from all over the world and provides a valuable source of income for the local economy. As such, it is important that the timing of spates, maintenance of flows and wider competition for water resources are managed in a manner to ensure long-term sustainability of the industry. The ability to generate long-term synthetic flow sequences, on rivers such as the Dee, can therefore be used in conjunction with habitat models to allow for the creation and implementation of future management plans.

The flow gauge is located at Woodend to the west of Banchory (OS grid reference NO634956) and has been active since 1929, with the data used in this study being recorded between 1929 and 2012. At this location the catchment area is 1,370 km^{2}. The Flood Estimation Handbook (FEH) catchment descriptors indicate that the mean annual rainfall is 1,108 mm and 0.05% of the catchment is considered urban.

### River Falloch

The Falloch is located near Crianlarich in Glen Falloch, Stirlingshire. It is the main river that drains the surrounding mountains (maximum elevation 1,120 m AOD), before discharging into Loch Lomond at its northern end. Owing to its location, the catchment is mainly mountainous characterised by moorland and small amounts of forestry at the lower levels.

With the requirement for sustainable energy sources becoming increasingly important, numerous hydropower schemes have been constructed within the West Highlands of the UK. Within Glen Falloch itself, there has been recent development of schemes on some of the smaller watercourses (Loch Lomond & Trossachs National Park 2012), with a run-of-the-river scheme proposed for the Falloch itself (The Scottish Government 2009); although this has yet to be commissioned. The generation of synthetic flow sequences, for rivers like the Falloch, can be used by engineers to provide a more detailed assessment of future hydropower schemes, ensuring their long-term sustainability.

The gauge is located in Glen Falloch (OS grid reference NN3121194) and has been active since 1970, with the flow data used in this study being recorded between 1970 and 2012. The FEH catchment descriptors show that the mean annual rainfall is 2,848 mm and 0.00% of the 80 km^{2} catchment is considered urban.

### River Caldew

The Caldew is located in Cumbria and has its source on Skiddaw Peak (926 m AOD) in the Lake District. The river then flows north into the town of Carlisle where it reaches its confluence with the River Eden. The catchment is characterised by heath and moorland in the headwaters, with both arable farming and urban centres in the lower reaches.

The Caldew is a relatively steep (channel bed gradient > 1:500) gravel bed river with a highly active sediment transport and morphological regime. Historically, it has been responsible for flooding parts of Carlisle and, as such, has been the subject of numerous flood risk studies (Neal *et al.* 2009, 2013; Mason *et al.* 2009; Horritt *et al.* 2010). After the floods in 2005, an alleviation scheme was constructed within Carlisle to protect the town from future flood events on the rivers Eden, Caldew and Petteril. Owing to the active morphological regime of the Caldew, a sediment transport modelling study was commissioned by the Environment Agency (EA) (Jacobs 2007), and some gravel bars are now regularly monitored to ensure levels do not undermine the flood defences in Carlisle. The ability to simulate long-term synthetic flow sequences, combined with up-to-date morphological modelling, would allow for more comprehensive sediment management plans to be implemented for rivers where sediment transport affects flood risk, such as the Caldew.

The gauge is located at Holm Hill approximately 10.5 km upstream of Carlisle (OS grid reference NY378468) and was active between 1968 and 2000. Since the closure of the gauge in 2000 a new one has now opened further downstream at Cummersdale (OS grid reference NY394527). As the historic gauge at Holm Hill has a longer flow record (31 years compared to 14 years) these measurements are used in the study. The FEH catchment descriptors show that the mean annual rainfall is 1,213 mm and 0.13% of the 147 km^{2} catchment is considered urban.

### River Lud

The Lud is located in Lincolnshire and flows east then northeast from its source in the Wolds hill range (maximum elevation of 153 m AOD). During the 1700s the river was canalised just downstream of the town of Louth, with the majority of the flow being diverted into this new canal system. Although the canal is no longer used for inland navigation, it is still an integral part of the drainage system for the region. As such, the EA maintains a gauging station to measure the discharge in the canal.

Apart from the town of Lough, the catchment is mainly rural, dominated by arable farming. It is predominantly (c. 70%) chalk, meaning that a high proportion of the river flow is generated from stored sources; resulting in a higher baseflow index (BFI) than the other rivers (0.90 compared to 0.53, 0.19 and 0.49). This means that, compared to the other catchments, a relatively high baseflow is maintained during dry periods. This, combined with the good chemical quality and low turbidity of the water, makes them notable for providing important habitats for a wide range of fish, invertebrate and plant species (Wilby *et al.* 1998; Mainstone 1999; Cranston & Darby 2004).

The gauge is located in the town of Louth (OS grid reference TF337879) and has been active since 1968, with the flow data used in this study being recorded between 1968 and 2013. The FEH catchment descriptors show that the mean annual rainfall is 698 mm and 2.56% of the 55 km^{2} catchment is considered urban.

### Site comparison

Comparing Figure 1(a)–1(d) it is evident that each catchment experiences the same seasonal variation of drier summer and wetter winter months; with the average summer (May–August) flow being 50%, 35% and 27% of that of the winter (November–February) for the Dee, Falloch and Caldew, respectively. Owing to the Falloch catchment being very mountainous, the hydrograph has a flashier nature than those of the other rivers; in evidence, Figure 1(b) exhibits sharper transitions (i.e., steeper rising limbs to peak flows) and appears most prone to summer convective storms (demonstrating a number of high flow events in mid-summer, June–August). This is also confirmed by the Falloch having the lowest BFI (0.19), which indicates little groundwater storage of the rainfall. In addition, the higher altitude and latitude of the Dee and Falloch catchments give rise to more notable snowmelt, in March and April specifically.

_{e}). This allows for easier visual comparison between the data sets. Comparing the probability densities for the recorded flows (Figure 2), the different hydrological nature of each catchment becomes evident. Owing to the predominant, rain inducing, low pressure weather systems originating in the Atlantic Ocean, there are considerable differences between the meteorology and climate in eastern and western regions (Kingston

*et al.*2006). This is visible in Figure 3, with the Dee (north east Scotland) and Lud (eastern England) having much narrower distributions than the Falloch and Caldew (south west Scotland and north west England, respectively). This indicates that the hydrographs for rivers in the west of the country are more variable in nature, which is confirmed by comparing the 1995 hydrographs in Figure 1.

## MODELLING TECHNIQUES

### ARMA modelling

*et al.*2008). The model is generally referred to as an ARMA(

*p*,

*q*) where

*p*and

*q*are the orders of the autoregressive and moving-average polynomials, respectively. The general form of an ARMA(

*p*,

*q*) model is described by Equation (2). where

*X*is a time series of estimated data;

_{t}*c*is a constant;

*ɛ*is a time series of white noise generated from a normal distribution (N(0,1)); and

_{t}*φ*and

*θ*are estimated parameters of the model.

Equation (2) allows for the generation of a synthetic time series from a fitted ARMA(*p*,*q*) model.

*Q*

_{deseason}is the de-seasonalised flow;

*Q*is the log-transformed recorded flow; and

*i*is the month of the year.

Upon de-seasonalisation of the recorded flows, the most suitable ARMA model is fitted to *Q*_{deseason} and the model parameters (*φ* and *θ*) estimated. The fitting and estimation procedure was conducted using the *R* statistical software package (www.r-project.org). The most suitable model was chosen following the method of Can *et al.* (2012), where a range of stationary ARMA(*p*,*q*) models are fitted to the data and the one that minimises the Akaike information criterion (AIC) identified as the most suitable. The AIC (Akaike 1974) is a widely used (Jones 1980; Findley 1985; Jones & Vecchia 1993; Ji & Chee 2011; Can *et al.* 2012) tool for model selection that identifies the most suitable parameter set by balancing the goodness-of-fit of the model with its complexity, to avoid over-parameterisation. The individual parameters of the model were then estimated using a maximum likelihood approach. Once the parameters of the model have been estimated, synthetic flow sequences can be generated using Equation (2). These flow sequences are then re-seasonalised to give physical values using the inverse of Equation (3). The order of the best-fit ARMA(*p,q*) models (Table 2) were determined to be an ARMA(4,1), ARMA(2,2), ARMA(3,2) and ARMA(2,2) for the Dee, Falloch, Caldew and Lud, respectively. Table 3 provides the parameter estimates for each model using a maximum likelihood estimation procedure.

River Dee . | River Falloch . | River Caldew . | River Lud . | ||||
---|---|---|---|---|---|---|---|

Model . | AIC . | Model . | AIC . | Model . | AIC . | Model . | AIC . |

ARMA(3,1) | 43,617 | ARMA(1,1) | 30,201 | ARMA(2,3) | 15,659 | ARMA(1,1) | 6,339 |

ARMA(4,0) | 43,656 | ARMA(1,2) | 30,075 | ARMA(3,3) | 15,675 | ARMA(1,2) | 4,898 |

ARMA(4,1) | 43,566 | ARMA(2,2) | 30,070 | ARMA(3,2) | 15,658 | ARMA(2,2) | 4,783 |

ARMA(4,2) | 43,619 | ARMA(2,3) | 30,072 | ARMA(4,4) | 15,659 | ARMA(2,3) | 4,786 |

ARMA(4,3) | 43,574 | ARMA(3,2) | 30,072 | ARMA(3,5) | 15,660 | ARMA(3,2) | 4,786 |

River Dee . | River Falloch . | River Caldew . | River Lud . | ||||
---|---|---|---|---|---|---|---|

Model . | AIC . | Model . | AIC . | Model . | AIC . | Model . | AIC . |

ARMA(3,1) | 43,617 | ARMA(1,1) | 30,201 | ARMA(2,3) | 15,659 | ARMA(1,1) | 6,339 |

ARMA(4,0) | 43,656 | ARMA(1,2) | 30,075 | ARMA(3,3) | 15,675 | ARMA(1,2) | 4,898 |

ARMA(4,1) | 43,566 | ARMA(2,2) | 30,070 | ARMA(3,2) | 15,658 | ARMA(2,2) | 4,783 |

ARMA(4,2) | 43,619 | ARMA(2,3) | 30,072 | ARMA(4,4) | 15,659 | ARMA(2,3) | 4,786 |

ARMA(4,3) | 43,574 | ARMA(3,2) | 30,072 | ARMA(3,5) | 15,660 | ARMA(3,2) | 4,786 |

The bold values indicate the best-fit models for each data set.

River . | φ_{1}
. | φ_{2}
. | φ_{3}
. | φ_{4}
. | θ_{1}
. | θ_{2}
. |
---|---|---|---|---|---|---|

Dee | 1.787 ± 0.025 | −0.979 ± 0.026 | 0.256 ± 0.014 | 1.787 ± 0.025 | −0.078 ± 0.08 | NA |

Falloch | 1.067 ± 0.078 | −0.180 ± 0.064 | NA | NA | −0.320 ± 0.080 | −0.107 ± 0.012 |

Caldew | 1.146 ± 0.080 | −0.098 ± 0.135 | −0.090 ± 0.070 | NA | −0.318 ± 0.081 | −0.229 ± 0.081 |

Lud | 1.260 ± 0.023 | −0.263 ± 0.023 | NA | NA | −0.535 ± 0.022 | −0.212 ± 0.013 |

River . | φ_{1}
. | φ_{2}
. | φ_{3}
. | φ_{4}
. | θ_{1}
. | θ_{2}
. |
---|---|---|---|---|---|---|

Dee | 1.787 ± 0.025 | −0.979 ± 0.026 | 0.256 ± 0.014 | 1.787 ± 0.025 | −0.078 ± 0.08 | NA |

Falloch | 1.067 ± 0.078 | −0.180 ± 0.064 | NA | NA | −0.320 ± 0.080 | −0.107 ± 0.012 |

Caldew | 1.146 ± 0.080 | −0.098 ± 0.135 | −0.090 ± 0.070 | NA | −0.318 ± 0.081 | −0.229 ± 0.081 |

Lud | 1.260 ± 0.023 | −0.263 ± 0.023 | NA | NA | −0.535 ± 0.022 | −0.212 ± 0.013 |

### Hidden Markov modelling

HMM was first introduced by Baum & Petrie (1966) and, since then, has been applied successfully to a range of complex problems where the standard MM technique appears limited. For effective modelling of such processes, in addition to the time-evolution of observed states (as in a simple MM), a HMM allows the evolution of the process through underlying un-observed (hidden) states. Thus, HMM are theoretically more advanced.

The model is essentially a dynamic Bayesian network that allows for the modelling of time series using probability distributions (Ghahramani 2001). Statistically, the HMM is distinct from ARMA as it is based on accounting for evolution probabilities of the observed state at time *t* from one previous state (first order Markov model) at *t-1,* through a set of un-observed states; in comparison, ARMA models use only auto-correlation of observations at time *t*, with a few subsequent observations in the past. Thus, the general make up of HMMs is described by Jenkins *et al.* (2014) using the following five elements: (a) a set of distinct observed states; (b) a set of unobserved states within the observed states; (c) a matrix to describe the probability of transition between the states – *state transitional probability matrix*; (d) a matrix of the probabilities of unobserved states – *emission probability matrix*; and (e) a set of initial probabilities for the observed states. To conduct modelling of a time series using a HMM, these five elements have to be obtained from the recorded data.

The implementation of a HMM to simulate daily streamflow is summarised below:

Identify discrete states (

*S*) in the flow record.Define the set of

*N*unobserved states to account for all possible values between the discrete state limits. For example, if state A were between 5 and 15 m^{3}/s the unobserved states correspond to 5, 6, 7, 8, … 15 m^{3}/s and so on. Each discrete state has*N*unobserved states.Define the state transition probability matrix. For

*S*states this is an*S*×*S*matrix with the value in row*i*and column*j*corresponding to the probability of flow transition from state*S*to state_{i}*S*._{j}Define the emission probability matrix that contains the emission probabilities of all unobserved states that correspond to each observed state, i.e., an

*S*×*N*matrix.Define a set of

*S*initial probabilities for each discrete state following the method of Rabiner (1989).

As one of the objectives of this paper is to assess the general applicability of a HMM approach to UK catchments, it is important that the partitioning of the recorded flow time series into distinct states is consistent across all catchments; this process is well-known to suffer from subjectivity (Augustin *et al.* 2008) and hence standard hydrological practices of using percentiles of flow is employed (Table 4). Each flow record is divided into 11 states, using increments of 10%, which ensures that each is sampled approximately the same number of times throughout the simulation and provides a transferrable methodology between test catchments. This is important, so as to ensure that the underlying flow values are well represented, allowing for an accurate emission probability matrix to be formulated. Owing to the difference in magnitude of the flows in the rivers, all recorded flows were rounded to zero decimal places for the Dee, one decimal place for the Falloch and Caldew and two decimal places for the Lud.

State . | Percentile range (%) . | Dee Q range (m^{3}/s)
. | Falloch Q range (m^{3}/s)
. | Caldew Q range (m^{3}/s)
. | Lud Q range (m^{3}/s)
. |
---|---|---|---|---|---|

A | 0–10 | <11.0 | <0.4 | <1.0 | <0.15 |

B | 10–20 | 11.1–15.0 | 0.5–0.7 | 1.1–1.3 | 0.16–0.19 |

C | 20–30 | 15.1–19.0 | 0.8–1.0 | 1.4–1.7 | 0.20–0.24 |

D | 30–40 | 19.1–22.0 | 1.1–1.5 | 1.8–2.1 | 0.25–0.29 |

E | 40–50 | 22.1–27.0 | 1.6–2.2 | 2.2–2.7 | 0.30–0.35 |

F | 50–60 | 27.1–32.0 | 2.3–3.3 | 2.8–3.5 | 0.36–0.44 |

G | 60–70 | 32.1–40.0 | 3.4–5.3 | 3.6–4.6 | 0.45–0.54 |

H | 70–80 | 40.1–51.0 | 5.4–9.0 | 4.7–6.5 | 0.55–0.68 |

I | 80–90 | 51.1–74.0 | 9.1–16.2 | 6.6–10.1 | 0.69–0.89 |

J | 90–99 | 74.1–182.0 | 16.3–45.0 | 10.2–27.4 | 0.90–1.52 |

K | 99–100 | >182.0 | >45.0 | >27.4 | >1.52 |

State . | Percentile range (%) . | Dee Q range (m^{3}/s)
. | Falloch Q range (m^{3}/s)
. | Caldew Q range (m^{3}/s)
. | Lud Q range (m^{3}/s)
. |
---|---|---|---|---|---|

A | 0–10 | <11.0 | <0.4 | <1.0 | <0.15 |

B | 10–20 | 11.1–15.0 | 0.5–0.7 | 1.1–1.3 | 0.16–0.19 |

C | 20–30 | 15.1–19.0 | 0.8–1.0 | 1.4–1.7 | 0.20–0.24 |

D | 30–40 | 19.1–22.0 | 1.1–1.5 | 1.8–2.1 | 0.25–0.29 |

E | 40–50 | 22.1–27.0 | 1.6–2.2 | 2.2–2.7 | 0.30–0.35 |

F | 50–60 | 27.1–32.0 | 2.3–3.3 | 2.8–3.5 | 0.36–0.44 |

G | 60–70 | 32.1–40.0 | 3.4–5.3 | 3.6–4.6 | 0.45–0.54 |

H | 70–80 | 40.1–51.0 | 5.4–9.0 | 4.7–6.5 | 0.55–0.68 |

I | 80–90 | 51.1–74.0 | 9.1–16.2 | 6.6–10.1 | 0.69–0.89 |

J | 90–99 | 74.1–182.0 | 16.3–45.0 | 10.2–27.4 | 0.90–1.52 |

K | 99–100 | >182.0 | >45.0 | >27.4 | >1.52 |

To give the best representation of the flow regime, it is preferable that the values within each state have a continuous distribution. The location of the first 10 states (A–J) ensured this for all catchments but this was not possible for the upper state (K) due to the sparseness of data at the extreme limit. To ensure that the extreme limits of the synthetic sequences are not constrained by the measured data, a continuous distribution is required to describe the flows within this state. Although any number of probability distributions could be sampled from, due to the extreme nature of the data within this state, it makes sense that an extreme value distribution is used. To overcome this limitation, Augustin *et al.* (2008) used the GEV distribution to estimate flows in their upper state. As the GEV model combines the traditional extreme value families (Weibull, Fréchet and Gumbel) it is considered the most suitable for inclusion (Coles 2001). Thus, a similar approach has been adopted here (HMM-GEV).

*μ*is the location parameter;

*σ*is the scale parameter; and

*ξ*is the shape parameter.

River . | μ
. | σ
. | ξ
. |
---|---|---|---|

Dee | 214.30 | 33.23 | 0.43 |

Falloch | 51.33 | 6.79 | 0.62 |

Caldew | 31.49 | 4.49 | 0.59 |

Lud | 1.61 | 0.11 | 0.84 |

River . | μ
. | σ
. | ξ
. |
---|---|---|---|

Dee | 214.30 | 33.23 | 0.43 |

Falloch | 51.33 | 6.79 | 0.62 |

Caldew | 31.49 | 4.49 | 0.59 |

Lud | 1.61 | 0.11 | 0.84 |

*et al.*2007; Callaghan

*et al.*2008; Pender & Karunarathna 2013; Karunarathna

*et al.*2014). Although this GP model only fits the distribution to the extreme tail of the data, the HMM-GP framework was implemented using the same >99th percentile fitting approach as the HMM-GEV. The selection of the thresholds was made from exploratory analysis of the data following the procedure recommended by Coles (2001). The cumulative distribution function of the GP is provided in Equation (5) with details of the parameter estimates for each catchment given in Table 6. Figure 5 provides the quantile plots in order to assess the goodness-of-fit of the model, for the estimated parameters. where

*u*is the threshold level;

*σ*is the scale parameter; and

*ξ*is the shape parameter;

*ζ*= Pr{

_{u}*X*>

*u*}.

River . | u
. | σ
. | ξ
. |
---|---|---|---|

Dee | 250 | 62.52 | 0.04 |

Falloch | 50 | 15.95 | 0.06 |

Caldew | 35 | 10.88 | 0.07 |

Lud | 2 | 0.50 | 0.26 |

River . | u
. | σ
. | ξ
. |
---|---|---|---|

Dee | 250 | 62.52 | 0.04 |

Falloch | 50 | 15.95 | 0.06 |

Caldew | 35 | 10.88 | 0.07 |

Lud | 2 | 0.50 | 0.26 |

Comparison of the quantile plots for both fitted models (Figures 4 and 5) shows that the use of the GP provides a better representation of the extreme tail of the flow data, for all catchments. This indicates that, for these data records, there are too many data above the 99th percentile to ensure a good fit with the GEV. As such, it is expected that, although improving on the ARMA model, the HMM-GEV will provide less realistic estimations of the upper extremes, compared to the HMM-GP.

As with the ARMA modelling, both HMM approaches were implemented using the *R* statistical software package to allow for efficient generation of synthetic flow sequences for each site.

## RESULTS AND DISCUSSION

### Analysis methods

*Q*,

_{p}_{synthetic}and

*Q*,

_{p}_{recorded}are the synthetic and recorded flow values at percentile

*p*, respectively.

*Q*

_{p}_{,synthetic}and

*Q*

_{p,}_{recorded}are the synthetic and recorded flow values at percentile

*p*, respectively.

### ARMA model

The results from the ARMA modelling show that generating synthetic flow sequences is more successful on the Dee, Caldew and Lud (RMAD 5.6%, 6.0% and 5.1%) compared to the Falloch (RMAD 38.5%). This reflects differences in the geographic and catchment characteristics, indicating poor model performance for the highly variable flow regime of the runoff-dominated (Figure 1(b)) mountainous catchment of the Falloch. In explanation, the ARMA model is based on a random white noise component; this fails to account for the probability of change in the flow regime and therefore induces poorer model performance where flows exhibit a rapid change in state.

Thus, only for the more attenuated hydrological regimes of the Dee, Caldew and Lud can the ARMA model capture and replicate the flow regime within the ±10% range of acceptability, between approximately the 5th and 95th percentiles (Figure 6(m)–6(p)). At <5th and >95th percentiles the ARMA model estimates fail to adequately simulate flow, indicating weakness in the modelling approach at extreme flows. On average, sequences for the Dee and Lud underestimate the most extreme values, with the percentage difference at the peak flow (100th percentile) being −11.1% and −25.1%, respectively. Although the ARMA model has been fitted to the data to give the most appropriate parameter set-up, it is important to note that this fitting only provides the best representation of the average changes in the flow regime, not extremes. This explains its inability to effectively capture the extreme flows.

### HMM extreme value models

The probability density and percentage difference plots for the HMM-GEV model are provided in Figure 6(e)–6(h) and 6(m)–6(p), respectively. The modelling resulted in a RMAD of 1.9%, 3.8%, 4.4% and 2.9% for the Dee, Falloch, Caldew and Lud, respectively. Equivalent plots for the HMM-GP model are provided in Figure 6(i)–6(l) and 6(m)–6(p), with the RMAD determined to be 1.4%, 1.5%, 1.9% and 1.1% for the Dee, Falloch, Caldew and Lud, respectively.

The HMM-GEV model is shown to produce more realistic results than that of the ARMA model. For the HMM-GEV, all data between the 0 and 95th percentiles lie within the ±10% acceptability limit, with the majority falling within ±5% (Figure 6(m)–6(p)) and the overall mean RMAD are <5% for all catchments. These results show significant improvement compared to the simple ARMA model, particularly in simulating the flashier hydrograph of the Falloch, where the RMAD is now only 3.8%. However, examining the fit to the upper extremes (Figure 6(m)–6(p)) indicates the inadequacy of the GEV method for >99th percentile. While these extreme events on the highly attenuated flow regime of the River Dee do fall within the acceptability range (±10%), the flashier regimes of the Caldew and Falloch are very poorly represented. Reviewing the GEV methodology, this is due to there being too many data above the 99th percentile to result in an accurate fit for the tail of the distribution, as discussed previously and demonstrated by Figures 4 and 5.

The use of the alternative HMM-GP approach is shown to produce a more realistic representation of the entire flow regime, including the upper extremes. From Figure 6(m)–6(p) it is clear that all simulated flows fall within the acceptability limits, with the GP method yielding the lowest overall RMAD of <2% for all catchments. All flow extremes are fully captured by this approach, with all high flows simulated well within the ±10% acceptability range, such that the HMM-GP approach appears generically applicable for simulating a wide range of hydrographic regimes across the geographic regions of the UK. The main reason behind the improvement, compared to the HMM-GEV model, is the improved representation of the extreme tail of the flow data, as evident from Figures 4 and 5.

## CONCLUSIONS

This paper has compared and contrasted a simple ARMA approach to more complex HMM methods for stochastic generation of daily streamflow sequences. In developing an appropriate HMM methodology, the use of GEV and GP extreme value distributions has been analysed and a generic methodology for the partitioning of the flow record has been developed in line with standard flow percentile descriptors. All models were tested across four UK catchments of geographic diversity and varying hydrologic regime.

The ARMA model has the advantages of: being much simpler to implement; being unconstrained by estimates to a number of decimal places; and having no requirement for an emission probability matrix. These results showed that the ARMA model is reasonably capable of producing acceptable (within ±10% of recorded values) realisations within the 5th to 95th percentile range, but unable to effectively capture extreme values. This limits the use of an ARMA model to catchments with flow regimes with highly attenuated hydrographs and, therefore, it has been determined that ARMA is unsuitable for accurate representation of the entire flow regime for the range of UK catchments modelled.

As the alternative HMM method makes use of the actual probabilities of transitioning between flow states and probabilities of occurrence within these states, a HMM has been shown to be more suitable for modelling a range of flow time series with different characteristics. As the traditional HMM is limited for streamflow data sets due to the long data tail and sparseness of recorded values in this part of the distribution, it has been proven that accurate simulation of extreme flows requires the HMM to be combined with a GP extreme value distribution, to sample values above the 99th percentile. The HMM-GP model output accurately (RMAD < 2%) simulated the full range of flows for all hydrological regimes of the UK case studies tested, showing robustness of methodology and confidence in generic applicability. It is therefore advocated that, provided that a sufficient historic data set exists, use of a HMM-GP approach, combined with a consistent percentile-based method for partitioning of the flow record, makes for an accessible and straightforward modelling tool that can be readily applied in practice. The use of the outputs of such an approach, in conjunction with more detailed hydraulic models, will significantly benefit a range of detailed sustainability assessments paramount to the future security of water resource, flood risk and ecosystem management in UK river systems.

## ACKNOWLEDGEMENTS

This research was carried out as part of the Engineering and Physical Sciences Research Council funded Flood MEMORY grant EP/K013513/1 held by HH. The authors would like to thank the National River Flow Archive for the provision of the gauged data that enabled the work to be undertaken. Finally, the authors would like to thank the reviewers, whose constructive comments helped improve the original manuscript.