Christophe Laplanche

Pedro M Leunda

Laurie Boithias

José Ardaíz

Francis Juanes

Advantages and insights from a hierarchical Bayesian growth and dynamics model based on salmonid electrofishing removal data

Keywords: Growth, Population dynamics, Salmo trutta, Depletion sampling, Iberian peninsula, Mesoscale

teaching and research institutions in France or abroad, or from public or private research centers.

Advantages and insights from a hierarchical Bayesian growth and dynamics model based on salmonid electrofishing removal data

Introduction

Growth is a fundamental ecological process of most organisms. This is especially true for fishes for three reasons. First, fish continue to grow though their lifetime, i.e. they have indeterminate growth, and body size can increase by several orders of magnitude (from an average size of 1 mm at the egg stage to several meters in the largest species) (Summerfeldt & Hall, 1987;[START_REF] Jobling | Life histories of fish[END_REF]. Second, growth rate is dynamic through the life history, typically high in early life and slower later in life, is the most variable component of fish energy budgets [START_REF] Jobling | Life histories of fish[END_REF], and can respond quickly in a compensatory fashion to changed conditions (Ali et al., 2003). Third, fish growth is driven by a variety of factors including genetics and both abiotic and biotic factors, as demonstrated both in the lab [START_REF] Brett | Growth rate and body composition of fingerling sockeye salmon, Oncorhynchus nerka, in relation to temperature and ration size[END_REF]Elliott, 1975a,b;Coleman & Fausch, 2007b) and in situ (Coleman & Fausch, 2007a;[START_REF] Robinson | Influence of variable interannual summer water temperatures on brook trout growth, consumption, reproduction, and mortality in an unstratified Adirondack lake[END_REF]Xu et al., 2010;[START_REF] Letcher | Robust estimates of environmental effects on population vital rates: an integrated capture-recapture model of seasonal brook trout growth, survival and movement in a stream network[END_REF].

Because most fish are ectothermic, their growth is especially sensitive to environmental variation, particularly temperature. But growth also responds in a context-specific way to interactions among multiple abiotic and biotic factors [START_REF] Klemetsen | Atlantic salmon Salmo salar L., brown trout Salmo trutta L. and Arctic charr Salvelinus alpinus (L.): a review of aspects of their life histories[END_REF] and as such is tightly related to population dynamics. In stream-dwelling salmonids, growth is sensitive to a wide variety of factors, including temperature, discharge, elevation, and conspecific density (Table 1, for a range of species) and is also strongly correlated with critical life history events [START_REF] Hutchings | Life histories of fish[END_REF][START_REF] Pepin | Reconsidering the impossible -linking environmental drivers to growth, mortality, and recruitment of fish[END_REF]. For example, growth can determine smolting age, size and age at migration, overwinter mortality, return timing, sexual maturity, success on the spawning grounds and emergence time of embryos, among others [START_REF] Quinn | The behavior and ecology of Pacific salmon and trout[END_REF][START_REF] Levings | Ecology of salmonids in estuaries around the world: Adaptations, habitats, and conservation[END_REF]. However, disentangling the relative importance of biotic and abiotic factors remains challenging as each is dynamic and either tightly related to climate and hydrology or to population dynamics, all of which have been shown to be sensitive to global climate change. The ability to accurately model fish growth and population dynamics thus becomes critical when making predictions about the future, e.g. effects of changes in both land use and climate [START_REF] Parra | Latitudinal and altitudinal growth patterns of brown trout Salmo trutta at different spatial scales[END_REF][START_REF] Parra | Unravelling the effects of water temperature and density dependence on the spatial variation of brown trout (Salmo trutta) body size[END_REF][START_REF] Boithias | Assessment of the water supply: DEmand ratios in a Mediterranean basin under different global change scenarios and mitigation alternatives[END_REF] on salmonid population dynamics [START_REF] Milner | The natural control of salmon and trout populations in streams[END_REF]Jonsson & Jonsson, 2009;[START_REF] Baumann | Reduced early life growth and survival in a fish in direct response to increased carbon dioxide[END_REF][START_REF] Martins | Climate effects on growth, phenology, and survival of sockeye salmon (Oncorhynchus nerka): a synthesis of the current state of knowledge and future research directions[END_REF][START_REF] Kovach | Impacts of climatic variation on trout: a global synthesis and path forward[END_REF][START_REF] Clavero | Historical citizen science to understand and predict climate-driven trout decline[END_REF].

Two main approaches are made available for researchers and managers to monitor growth and dynamics of salmonid populations in the field: Individual Tagging Methods (ITM, e.g., using Passive Integrated Transponder 'PIT' tags) and Removal Sampling by ElectroFishing (EFRS). ITM provides information about individuals but are expensive to operate at a large spatial scale (although this is possible; [START_REF] Marvin | The Success of the Columbia Basin Passive Integrated Transponder (PIT) Tag Information System[END_REF]). EFRS is less precise about some aspects, by providing information on open groups of individuals, but requires less sampling time (see below). In view of their relative advantages, ITM and EFRS have both been used to monitor growth and dynamics of freshwater salmonid populations in the field, although EFRS is more common when studying growth (Table 1).

Two main reasons explain the popularity of EFRS: the relatively short sampling time it requires to collect data and the ease and wide variety of methods that can be used with it to compute maximum likelihood estimates of population size (reviewed by [START_REF] Cowx | Review of the methods for estimating fish population size from survey removal data[END_REF]). Another option is to use a dedicated software for a wider choice of models (e.g., MARK, although its main use is for ITM data, [START_REF] White | Program MARK: survival estimation from populations of marked animals[END_REF]). The ease of monitoring fish populations with EFRS has led to uninterrupted series of long-term data over large spatial scales, usually collected for management perspectives and later used for research (see for instance [START_REF] Parra | Latitudinal and altitudinal growth patterns of brown trout Salmo trutta at different spatial scales[END_REF]; [START_REF] Filipe | Forecasting fish distribution along stream networks: brown trout (Salmo trutta) in Europe[END_REF]; [START_REF] Bergerot | Relating life-history traits, environmental constraints and local extinctions in river fish[END_REF]).

Monitoring fish populations with EFRS (e.g., to estimate recruitment or mortality rates) includes measuring fish age. Calcified structures -otoliths or scales-can be sampled (lethally or non-lethally, respectively) on collected individuals and used to estimate fish age [START_REF] Dortel | Accounting for Age Uncertainty in Growth Modeling, the Case Study of Yellowfin Tuna (Thunnus albacares) of the Indian Ocean[END_REF]. An alternative to otoliths and scales for fishes in temperate climates is to measure the length of collected individuals and infer population structure from the statistical distribution of length data. This is possible for stream-dwelling trout because the length distribution is multimodal, with one component per year of emergence ('cohort'). The main reason for the multimodality is that reproduction occurs during a short period in autumn/winter [START_REF] Isely | Age and growth[END_REF]. Many statistical methods are available to managers to easily separate overlapping length distributions across ages [START_REF] Pitcher | A bumpy old road: Sized-based methods in fisheries assessment[END_REF]. Individual fish length is consequently collected during EFRS surveys (later referred to as 'EFRS length data'), thus providing long-term data over large spatial scales, ideal for investigating effects of global change on the growth and dynamics of salmonid populations [START_REF] Naslund | Brown trout (Salmo trutta) habitat use and life history in Swedish streams: possible effects of biotic interactions[END_REF][START_REF] Parra | Latitudinal and altitudinal growth patterns of brown trout Salmo trutta at different spatial scales[END_REF][START_REF] Filipe | Forecasting fish distribution along stream networks: brown trout (Salmo trutta) in Europe[END_REF][START_REF] Bergerot | Relating life-history traits, environmental constraints and local extinctions in river fish[END_REF][START_REF] Kanno | Seasonal weather patterns drive population vital rates and persistence in a stream fish[END_REF].

More recently, Hierarchical Bayesian Modelling (HBM) has increased the interest in using EFRS data to study stream-dwelling salmonid ecology (also applies to ITM data, see [START_REF] Kéry | Bayesian population analysis using WinBUGS: a hierarchical perspective[END_REF]). One main reason for the renewed interest is that the HBM framework offers the ability to build observation models that are connected to ecological models, both possibly advanced (e.g., more than what dedicated tools such as MARK can offer), as integrated models [START_REF] Letcher | Robust estimates of environmental effects on population vital rates: an integrated capture-recapture model of seasonal brook trout growth, survival and movement in a stream network[END_REF]. More specifically, HBMs have proven to effectively model EFRS observations [START_REF] Rivot | Hierarchical Bayesian modelling with habitat and time covariates for estimating riverine fish population size by successive removal method[END_REF], multimodal length distributions [START_REF] Ruiz | A hierarchical model to estimate the abundance and biomass of salmonids by using removal sampling and biometric data from multiple locations[END_REF], growth [START_REF] He | Modeling annual growth variation using a hierarchical Bayesian approach and the von Bertalanffy growth function, with application to lake trout in southern Lake Huron[END_REF][START_REF] Bal | Effect of water temperature and density of juvenile salmonids on growth of young-of-the-year Atlantic salmon Salmo salar[END_REF][START_REF] Lecomte | A length-based hierarchical model of brown trout (Salmo trutta fario) growth and production[END_REF][START_REF] Sigourney | Combining a Bayesian nonparametric method with a hierarchical framework to estimate individual and temporal variation in growth[END_REF][START_REF] Dortel | Accounting for Age Uncertainty in Growth Modeling, the Case Study of Yellowfin Tuna (Thunnus albacares) of the Indian Ocean[END_REF], and population dynamics [START_REF] Kanno | Seasonal weather patterns drive population vital rates and persistence in a stream fish[END_REF][START_REF] Bret | Understanding inter-reach variation in brown trout (Salmo trutta) mortality rates using a hierarchical Bayesian state-space model[END_REF]. Other reasons for the growing popularity of HBMs include their ability to propagate uncertainty from observations to parameter estimates and to compare competing models to test hypotheses [START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF]. The HBM framework also allows the use of prior distribution with model parameters (e.g., based on earlier studies) and definition of a hierarchical structure that facilitates spatial inter-/extrapolation and forecasting [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF][START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF].

While raw EFRS length data have been used to infer somatic growth [START_REF] Lecomte | A length-based hierarchical model of brown trout (Salmo trutta fario) growth and production[END_REF], and pre-processed length data (into mean-length-at-age and density-atage estimates) have been used to model either growth or population dynamics separately (e.g., [START_REF] He | Modeling variation in mass-length relations and condition indices of lake trout and Chinook salmon in Lake Huron: a hierarchical Bayesian approach[END_REF]; [START_REF] Laplanche | Modeling the fish community population dynamics and forecasting the eradication success of an exotic fish from an alpine stream[END_REF]), to our knowledge, raw EFRS length data have never been used to infer growth and population dynamics at the same time as an integrated model. We thus present an integrated HBM that models observations, somatic growth, and basic population dynamics. We illustrate the capabilities of the modelling framework by applying it to long-term data collected over a large spatial scale (Salmo trutta; 20 years; 48 sites). There is an apparent wide diversity of factors that affect stream-dwelling salmonid growth in situ (Table 1), which results of our modelling approach help explain. We further highlight advantages of our integrated approach and suggestions for potentially rewarding model extensions.

Materials and methods

Growth, population dynamics, and observation models

As mentioned, reproduction of stream-dwelling salmonids follows a yearly pattern.

In contrast, growth is continuous. We thus needed to define two time structures: an index over years-of-emergence i.e. cohorts (y ∈ {1, . . . , Y }; y = 1 for the first modelled cohort; Y consecutive cohorts) and an additional continuous time variable (t, in days; t = 0 on January 1st of year y = 1). To simplify presentation, equations are presented below as if there were only one sampled/modelled location. The spatial dimension and the hierarchical structure of the model are presented later (section 2.2).

Modelling the time/size at emergence

Times of emergence are strongly year-dependent, because spawning is mainly triggered by a decline in photoperiod and temperature (Jonsson & Jonsson, 2009), and because development of trout eggs from spawning to emergence is mainly driven by water temperature [START_REF] Elliott | An individual-based model for predicting the emergence period of sea trout fry in a Lake District stream[END_REF][START_REF] Ojanguren | Thermal dependence of embryonic growth and development in brown trout[END_REF]Jonsson & Jonsson, 2009).

Additional inter-individual differences in spawning times (spawning lasts for several weeks around a peak of activity; e.g., [START_REF] Riedl | Timing of brown trout spawning in Alpine rivers with special consideration of egg burial depth[END_REF]; [START_REF] Isely | Age and growth[END_REF]) are magnified by inter-individual differences in the development of eggs and parr, causing inter-individual differences in emergence times for a given year. We define d emerg y (day) as the (year-dependent) median time of emergence, i.e. the day of year y when half the fry have emerged.

Many studies identified additional inter-individual variation in size at emergence. As it was not possible to disentangle variation in both time and size at emergence with EFRS length data alone, we model inter-individual variation in both the time of and size at emergence as a single source of variation, in the form of the distribution of theoretical fish length at time d emerg y . While we consider d emerg y an unknown parameter in the model, we assume that mean trout length at emergence is known and constant. We denote L emerg (in mm) this quantity, i.e. mean length at time d emerg y of the cohort which emerged in year y.

Modelling cohort growth

The growth model which follows is the consequence of a similar model working at the individual fish level with random variation of growth parameters among individuals. We present the individual growth model in Appendix A and keep to the cohort level in the following, which is of greater interest given the available data (e.g., EFRS length data).

One central assumption of the cohort growth model is that fish of a given cohort grow under similar environmental conditions. This is the case for stream-dwelling salmonids, due to limited movement, which includes long-distance return migration for reproduction, meso-habitat movement as habitat needs change through their life-time, and daily microhabitat movements [START_REF] Schlosser | Stream fish ecology: a landscape perspective[END_REF][START_REF] Gido | Community ecology of stream fishes: Synthesis and future direction[END_REF][START_REF] White | A 'behaviorscape' perspective on stream fish ecology and conservation: linking fish behavior to riverscapes[END_REF][START_REF] Matthews | Stream fish community dynamics: A critical synthesis[END_REF][START_REF] Laplanche | Modeling the fish community population dynamics and forecasting the eradication success of an exotic fish from an alpine stream[END_REF].

The mean length at time t of the cohort which emerged on year y is denoted µ y (t), where t highlights the fact that growth is time-dependent, and subscript y specifies that mean length is also year-of-emergence dependent, since several cohorts exist at the same time. Cohort growth is modelled as

dµ y (t) dt = H y (t) for t ≥ d emerg y , (1)
starting from µ y (d emerg y ) = L emerg , where H y (t) is daily length increase, which is also year-of-emergence dependent and time dependent.

Growth rate of stream-dwelling salmonids decreases with fish age, which can be appropriately modelled using empirical, concave growth functions (von Bertalanffy, Gompertz, etc.). We follow [START_REF] Elliott | A new, improved growth model for brown trout, Salmo trutta[END_REF], who modelled growth as linear for certains powers of weight, and retain some of their notations to further facilitate comparison of results.

Hence, growth rate is expressed as

d(W y (t) b ) dt = b G y (t) 100 for t ≥ d emerg y , (2)
where W y (t) is the mean weight at time t of the cohort which emerged in year y, b is the power when weight raised to this power grows linearly, and G y (t) is a year-of-emergence dependent and time-dependent parameter. In the case of a one-to-one length-weight relationship (W y (t) = a w µ y (t) bw ; see data section), the model is equivalent to having daily length increase proportional to the power of length

dµ y (t) dt = 1 b w a 1/bw w G y (t) 100 µ y (t) 1-bbw . (3)
Growth curves covered by this 'power growth model' are illustrated in Appendix S1.

High correlation between growth parameters b and G y (t), due to multiplying them in eq.

(2), and as illustrated in Appendix S1, compels us to set one of the two parameters as a constant, in our case b, and be more flexible on the other, G y (t).

The effect of water temperature on growth rate is introduced into the model by defining

G y (t) = X(T w (t))G y (t), (4)
where X(T w (t)) models the effect of temperature on growth and G y (t) is a random effect (defined later). The function X ∈ [0, 1] defines the suitability of water temperature for growth, equaling to 0 below a minimum (T min ) and above a maximum water temperature (T max ), and reaching 1 at an optimal temperature (T opt ). We chose a rational function (e.g., [START_REF] Mallet | Growth modelling in accordance with daily water temperature in European grayling (Thymallus thymallus L.)[END_REF]).

In sum, daily growth rate is the product of 3 terms: µ y (t) 1-bbw /100b w a 1/bw w , which models a decrease in growth rate with increasing age; X(T w (t)), which models the suitability of water temperature for growth; and G y (t), which accounts for other sources of variation.

Modelling growth dispersion

Differences in growth trajectories among individual fish led us to model distribution of trout length at any time of a given cohort with a normal distribution (Appendix A). The individual growth model also led us to express the standard deviation of length within a cohort (denoted σ y (t)) as proportional to its mean, thus modelling the spread of the length distribution of cohorts over time, as follows

σ y (t) = νµ y (t), (5)
where the coefficient of variation (CV) ν is a direct measure of the variation in growth rates among individual fish (Appendix A).

Modelling the distribution of fish length

Because of the normal distribution of trout length within a cohort, trout length from all the included cohorts is modelled as a mixture of Gaussian distributions, one component per age (Figure 1). The theoretical probability density function of trout length at time t is consequently

f (t, x) = K k=1 λ y(t)-k+1 (t) λ(t) 1 √ 2πσ y(t)-k+1 (t) exp - (x -µ y(t)-k+1 (t)) 2 2σ y(t)-k+1 (t) 2 , ( 6
)
where

(t) (m -2 )
is the density at time t of the cohort that emerged in year y; λ(t) = K k=1 λ y(t)-k+1 (t) is the overall trout density at time t.

Fish length range is divided into L intervals of width ∆x (mm), from 0 to maximum length x max = L∆x (mm) (class centers are denoted x l = (l -1/2)∆x in mm; l ∈ {1, . . . , L} is an index over length classes). The expected density of fish of size class l at time t is therefore λ l (t) = λ(t) l∆x (l-1)∆x f (t, x)dx m -2 .

Modelling observations

The number of fish actually present is modelled as a Poisson variate (Wyatt, 2002)

N l (t) ∼ Poisson(Aλ l (t)), (7)
where A (m 2 ) is the area which is sampled by EFRS. The Poisson distribution models stochasticity of fish presence and assumes that the distributions of individuals for a given size class are independent of one another and are not spatially structured, e.g., via physical habitat characteristics [START_REF] Peterson | On the estimation of detection probabilities for sampling streamdwelling fishes[END_REF].

The number of fish of size class l caught at time t by electrofishing the area during removal r, observations of which were referred to as EFRS length data, is modelled as a binomial variate (Wyatt, 2002;[START_REF] Kanno | Seasonal weather patterns drive population vital rates and persistence in a stream fish[END_REF] C l,r (t) ∼ Binomial(R l,r (t), p l,r (t)),

where R l,1 (t) = N l (t) and R l,r (t) = R l,r-1 (t) -C l,r-1 (t) (r ≥ 2) is the stock left before removal r. The binomial distribution models stochasticity of fish capture by assuming that capture of fish of a given size class in the sampled area is independent with the same probability. Capture probability p l,r (t) increases with increasing fish size, which is modelled as logit(p l,r (t)) = αx l /1000 + β. Equations and values of the variables of the HBM are shown in Table 2. Relationships between HBM variables of the growth and population dynamics models are illustrated with a Directed Acyclic Graph (DAG, Figure 2). HBM variables, either measured or unknown, may be scalar, vectors, or multi-dimensional, as indicated by their subscript(s).

As an illustration, (known) times of EFRS surveys are grouped together in the variable d s,y (days), which has 2 dimensions: site and year.

Details due to the study of Salmo trutta

Based on other studies, we chose L emerg = 30 mm for the size of emergence [START_REF] Nika | Change in allometric length-weight relationship of Salmo trutta at emergence from the redd[END_REF]. We used parameter values published by [START_REF] Elliott | A new, improved growth model for brown trout, Salmo trutta[END_REF] for T min , T max , and T opt to calculate the temperature-dependent growth rate (Appendix S1). Brown trout is an autumn spawner, which makes the length distribution of trout of age 0+ observable in summer, when our sampling took place.

2.2.2. Details due to the temporal structure of the sampling scheme

The data sampling scheme (uninterrupted series of one EFRS survey in summer each year, see section 2.4) influenced details of the temporal structure of the HBM. The mean, standard deviation, and density of each component that defines the multimodal distribution of trout length at survey times (eq. ( 6)) are denoted µ s,y,k (mm), σ s,y,k (mm), and λ s,y,k (m -2 ), respectively. These parameters play a special role in the HBM by being directly connected to the growth model (in the case of µ s,y,k and σ s,y,k ), to the population dynamics model (λ s,y,k ), and to the observation model (µ s,y,k , σ s,y,k , and λ s,y,k ), as highlighted in the DAG (Figure 2). The resulting expected number of fish in each size class present in sampled area A s,y at survey times is denoted E(N s,y,l ) = A s,y λ s,y f s,y,l , where f s,y,l is found by integrating eq. ( 6) over size class l, and λ s,y = k λ s,y,k denotes overall trout density. The probability of capturing fish during EFRS surveys is denoted p s,y,l,r . EFRS length data, for each site, year, size class, and removal, are gathered into a 4-dimensional contingency table, denoted C s,y,l,r .

Our sampling scheme also implies that 'only' one observation of the multimodal distribution of trout length is available each year. As a result, we defined the random effect in eq. ( 4) in our HBM as site-and year-dependent (G s,y ). The population dynamics model reduces to a Markov process, with site-, year-, and age-dependent apparent survival rates (denoted Sur s,y,k ) between subsequent survey times. The abundance of trout of age 0+ at survey times is, using the notation defined earlier, λ s,y,1 .

Continuous variable t used in the growth, population dynamics, and observation mod- 3) was thus integrated at a daily time step using Euler's forward method. In this case, we approximated predicted mean lengths at survey times as

µ s,y,k = (L emerg ) bbw + b a b w ΣG s,y,k 100 1/bbw , (10)
where ΣG s,y,k is the cumulative sum of G s,y,k,d from emergence to observation.

Trout of age 1 and older on the year of the first EFRS survey emerged K -1 years before this year. As a result, we modelled growth and population dynamics K -1 years before the year of the first EFRS survey. The lack of need of backcasting/forecasting in our case study led us to model cohorts from this point to the year of the last EFRS survey.

Index y, defined earlier, thus still represents cohorts in the HBM, while the first EFRS survey corresponds to y = K and the last one to y = Y .

Spatial structure of the HBM

Some quantities defined during model presentation become spatially dependent, which we highlighted with variable subscript s (Table 2 andFigure 2). Prospective spatial similarities of growth rates (G s,y ), recruitment (λ s,y,1 ), and apparent survival rates (Sur s,y,k ) are modelled as follows

         G s,y ∼ Lognormal(log(G s ), σ 2 G s ) λ s,y,1 ∼ Lognormal(log(λ 1 ), σ 2 λ 1 ) Sur s,y,k ∼ Lognormal(log(Sur k ), σ 2 Sur,k ) . ( 11
)
We chose log-normal distributions to model multiplicative errors for G s,y and Sur s,y,k and to model variation in animal density [START_REF] Limpert | Log-normal Distributions across the Sciences: Keys and Clues[END_REF]. Hyperparameters log(λ 1 ), log(Sur k ), log(G s ) and σ 2 λ 1 , σ 2 Sur,k , σ 2 G s are regional means and variances of log(λ s,y,1 ), log(Sur s,y,k ), and log(G s,y ), respectively.

Priors

All model parameters were provided with vague uniform priors (Table 2), between 0 and 1 for λ 1 , Sur k , b, and ν s , and between 0 and 10 for σ λ 1 , σ Sur,k , G s , and σ G s . We provided time of emergence (denoted d emerg s,y

) with a uniform prior of a 6-month amplitude (±120 days) around a site-dependent, known value.

2.3.

Computations and measure of model fit

Simulating samples of the posterior distribution of the HBM

Equations provided in the previous section can be combined to express the posterior distribution of model parameters. The HBM is, however, too complex for such a distribution to be analytically tractable. Samples from the posterior distribution can be relatively easily simulated via Markov chain Monte Carlo (MCMC), and we used OpenBUGS for this purpose [START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF]. The code of our HBM and a tutorial are both available as Appendix S2. Data pre-processing and output post-processing were implemented in R (R Core Team, 2014). MCMC convergence was assessed by computing inter-chain variances of simulated latent variable samples across 3 chains; initializations were computed using 5, 50, and 95% prior marginal quantiles. After convergence, 5,000 samples were simulated. Only effective sample sizes (ESS) are reported. Point estimates are posterior means.

Absolute measure of model fit

We assessed model fit by comparing the EFRS data collected in the field to their respective values simulated by the model. Observed catch, all removals pooled together, is denoted C s,y,l = r C s,y,l,r (Figure 1; sum of the stacked bins for each size class). The distributions of the observed (C s,y,l ) and the expected modelled (C pred s,y,l = E(N s,y,l ) r p s,y,l,r (1p s,y,l,r ) r-1 ) distributions of the catches were compared using standard quantile-quantile (Q-Q) plots. A finer comparison of the distributions with a measure of the level of significance of the potentially under-and over-estimated values (for each site, year, and length class) were obtained in a Bayesian framework by using posterior predictive p-values [START_REF] Gelman | Bayesian data analysis[END_REF][START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF][START_REF] Chambert | Use of posterior predictive checks as an inferential tool for investigating individual heterogeneity in animal population vital rates[END_REF]. For this purpose, replicated data (C rep s,y,l ) were simulated by the fitted model, which is relatively easy to perform with BUGS [START_REF] Ntzoufras | Bayesian modeling using WinBUGS[END_REF]; [START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF]; Appendix S2). The scope of our model checking is to evaluate the fitness of the survey layer (i.e. Poisson-layer and removal sampling) given the estimated length distribution, and consequently given values for model hyperparameters such as growth and survival rates. For that reason, replicated data were simulated using the Poisson and Binomial models (eqs ( 7)-( 8)). The desired p-value, as the probability that the replicated data (of size class l, site s, year y) could be more extreme than the observed data, is to recreational angling except from some reaches that have preserved sections. Stocking in upper and middle reaches of the study area stopped in 1992 but continues in lower reaches, where brown trout is not the dominant species. Human population density is low in the study area (<10 inhabitants/km 2 ), and rivers are not degraded by anthropogenic land uses or pollution so their ecological status is good or very good (see internal reports commissioned by the Department of the Environment of the regional Government of Navarra (DEGN) considering physical-chemical water parameters and biological wa-ter quality). Agricultural land use, hydroelectric power stations, and dams are the main human pressures in the study area [START_REF] Parra | Latitudinal and altitudinal growth patterns of brown trout Salmo trutta at different spatial scales[END_REF].

p B s,y,l = Pr(C rep s,y,l ≥ C s,y,l ). ( 12

EFRS survey network

Electrofishing data were collected by the Fish and Game section of the DEGN (http://cazaypesca.navarra

The survey network is composed of 61 sampling sites (Figure 3) which are located in every river in upper, middle, and lower reaches and some scattered tributaries. Streams were surveyed once a year every summer (July-September) in 1992-2014 with the exception of 5 sites, where sampling started later (2 in 1997; 1999; 2000; 2005). Surveys suspected of being influenced by the presence of stocked individuals were excluded: (1) surveys before 1995 and (2) sites in lower reaches where stocking continues (Figure 3). As a result, the survey network we used is composed of 48 sites sampled for 19.5±1.7 consecutive years for a total of 48 × 20 -23 = 937 EFRS surveys. Catchment areas of upstream sampling sites ranged from 9.2 to 614.5 km 2 (mean: 87.9 km 2 ), and slopes at sampling sites ranged from 0.27 to 7.68% (mean: 1.47%). Sampled area differed among sites and years depending on stream width (8.2±3.6 m) and reach length (105.1±35.3 m). The sampling time required to survey 1000 m 2 ranged from 23 to 127 minutes (mean: 55); variability depending on habitat heterogeneity and fish density.

Fish assessment

One-to three-pass depletion electrofishing was performed, with the two-pass design being the most frequent (not sampled due to surveys started later than 1995: 2.4%; 1 removal: 7.0%; 2: 89.2%; 3: 1.5%). Each captured individual was measured for fork length (±1 mm) before being released, for a total of 189,533 fish-length data samples.

Modelling drove us to code fish length data by 10-mm length class. Trout are relatively small in the study area with a short life-span, and thus the maximum size class was set to 400 mm (99% quantile of trout size is 298 mm, length of 146 individuals > 400 mm were truncated to 400 mm) with a maximum age of 4 (K = 4). Raw length data collected during EFRS surveys were therefore turned into a 4-dimensional table with the number of trout caught at each site, year, length class, and removal (later refered to as EFRS length data), with missing values depending on when surveys started and how many fish were removed. A total of 14,296 brown trout collected in the study area were also measured for weight. Predicted weights (g) were modelled as W = a w L bw , where L is length (mm), and a w = 1.09 10 -5 and b w = 3.010 are the scaling coefficient and exponent respectively (R 2 > 0.99; linear regression on log-transformed variables). This relationship was used to convert between length and weight in the growth model (e.g., eq. ( 10)).

Water temperature

As presented below and in additional detail in Appendix S3, we used air temperature as a proxy for water temperature to calculate values of water temperature at EFRS sites for the entire period (as the model requires). Air temperature, water temperature, and EFRS length data were collected using three independent survey networks, at distinct sites and for different time periods (Figure 3). In short, daily air temperature was spatially interpolated by universal kriging using elevation as a linear predictor with day-dependent regression coefficients. Monthly water temperature was linearly related to monthly air temperature using site-dependent regression coefficients. The seasonal trend in water temperature at EFRS sites was obtained by using the air-water temperature relationship of the nearest water temperature gauge with the seasonal trend in daily air temperatures at EFRS sites as inputs. Simulated water temperatures (range 0.7-23.2 • C) covered the range of temperature for brown trout growth (Appendix S1).

Results

Fish length distribution

The multimodal fish length distribution was predicted for each of the 937 EFRS surveys (Figure 1). QQ-plots of observed vs. expected modelled fish length distributions indicate that the observed and expected modelled distributions of the catch are fairly similar (Figure 4). The results of the posterior predictive checks highlight size classes of underand over-estimated catch (Figure 1). The combined results of the posterior predictive checks indicate that the model fits correctly, to the exception of site 1620, due to some over-estimated values, and to the exception of sites 1010 and 1830, due to scarce data (Figure 5). From these results, we conclude that the model provides a picture of the distribution of fish length data that is well supported by the data, both in terms of expected values (mixture of multimodal distribution) and dispersion around these values (Poisson and binomial model), thus allowing us to produce and interpret estimates for model's hyperparameters, e.g. related to growth and population dynamics.

Growth and population dynamics

The model simulates seasonal and interannual growth variation (Figure 6) as well as differences in growth profiles between sites (see below). The estimated value for parameter b = 0.525 indicates that growth curves are moderately concave in the study area. Estimated values for the growth dispersal parameter (ν s ) ranging from 0.09 to 0.17 (Appendix S4) indicate that growth dispersion is strongly site-dependent.

The model also simulates decrease in the apparent survival rate with increasing age (Figure 6), as indicated by estimated values of apparent survival rates at the regional Marginal posterior distributions for model parameters (λ 1 , σ λ 1 , Sur k , and σ Sur,k ) are illustrated in Appendix S4.

Hierarchical structure of growth rate

Daily growth rate was modelled as the product of three terms: size-dependence, temperature-dependence, and other sources of variation. We investigated further the relative contribution of the two latter terms to variation in the daily growth rate. For that purpose, we calculated the variance of the log of the product X s,d G s,y (0.36), which sums up into the variances of log(X s,d ) (0.24) and of log(G s,y ) (0.14) plus twice their covariance (-0.05). These results indicate that (1) the deterministic, temperature-dependent term X s,d and the random term G s,y are weakly correlated with each other (Pearson's r = -0.13) and that (2) they contribute respectively and approximately 2/3 and 1/3 of the variance in daily growth rate, size-dependence excluded, on a log-scale.

We investigated further the relative contribution of the two latter terms to the modelled spatio-temporal variation in daily growth rate. We thus calculated ANOVA sums of squares (SSQ) of both log(X s,d ) (using site, year, and month as factors) and log(G s,y ) (using site and year). Percentages of SSQs for each term are shown in Table 3. Results show that one major source of variation in the deterministic, temperature-dependent term is seasonal (month: 30.6%). Two other main contributions highlight the importance of the site-specificity of temporal variation driven by temperature (site*month: 24.7%; site*day:

31.7%). Variation in the random term is mainly spatial (site : 73.0%). This analysis also indicates that there is a negligible global trend in the growth rate over the last 20 years (year: 1.5-2.5%), although there is a considerable site-dependent, yearly trend for G s,y

(site*year: 24.5%). We did not detect any systematic increase in water temperature over time from 1995-2014, either at each site or at the regional scale (linear regression using mean water temperature).

We investigated further reasons for the spatial trend in the random term. We found that G s (mean value of G s,y at each site) increases in streams in the downstream direction (Figure 7). We found that G s was weakly correlated with mean water temperature (T w s ; log-transformed; r = -0.13), indicating that the variation of G s in the downstream direction is not due to water warming while flowing downstream. We instead found that G s was positively correlated with catchment area (denoted wsa s ; log-transformed;

Pearson's r = 0.73) leading towards other possible explanations for spatial variation in growth rate, as discussed later. Parameter G s was also negatively correlated with stream slope (slope s ; log-transformed; r = -0.59). Further linear regression analysis showed that a combination of catchment area and stream slope predicts the random term well, leading to the relationship log 10 (G s ) = 0.175 log 10 (wsa s /slope s ), which explained up to 64% of its variability (R 2 = 0.64).

Discussion

Strength of the approach

Collecting fish length data by electrofishing has been used widely for several decades to monitor riverine fish populations. This practice has resulted in long-term monitoring over large spatial scales. Researchers have taken advantage of it to measure, e.g., effects of global change on fish populations [START_REF] Naslund | Brown trout (Salmo trutta) habitat use and life history in Swedish streams: possible effects of biotic interactions[END_REF][START_REF] Parra | Latitudinal and altitudinal growth patterns of brown trout Salmo trutta at different spatial scales[END_REF][START_REF] Filipe | Forecasting fish distribution along stream networks: brown trout (Salmo trutta) in Europe[END_REF][START_REF] Bergerot | Relating life-history traits, environmental constraints and local extinctions in river fish[END_REF]. The HBM framework has increased the value of EFRS data by allowing the investigation of more challenging scientific questions and the refinement of data use (e.g., large scale data analysis by [START_REF] Kanno | Seasonal weather patterns drive population vital rates and persistence in a stream fish[END_REF]). Our study illustrates the ability of the HBM framework to consider a relatively complex model for EFRS length observations (mixture of distributions, capture probability increasing with increasing fish size, and Poisson/binomial errors), a growth model with a relatively complex hierarchical structure, and a population dynamics model (relatively simple in our case, as a Markov process with a simple hierarchical structure). There are three major advantages in connecting these three models together, or more generally in connecting observation and ecological models together, as an integrated model. The first is to combine diverse datasets by connecting two observation models together [START_REF] Myers | Stock and recruitment: Generalizations about maximum reproductive rate, density dependence, and variability using meta-analytic approaches[END_REF][START_REF] Kéry | Bayesian population analysis using WinBUGS: a hierarchical perspective[END_REF]. The second is to share parameters both ways by connecting two ecological models together, in our case to model density-dependent somatic growth or size-dependent mortality, which both turned out to be important processes to consider [START_REF] Sogard | Size-Selective Mortality in the Juvenile Stage of Teleost Fishes: a Review[END_REF][START_REF] Imre | Density-dependent growth of young-of-the-year Atlantic salmon Salmo salar in Catamaran Brook, New Brunswick[END_REF][START_REF] Myrvold | Density dependence and its impact on individual growth rates in an age-structured stream salmonid population[END_REF]. The third is to infer ecological processes from the data, as a result of model calibration, by connecting an ecological model to an observation model (e.g., [START_REF] Laplanche | Modeling the fish community population dynamics and forecasting the eradication success of an exotic fish from an alpine stream[END_REF]). Another major advantage in connecting an ecological model to an observation model is to 'enlighten' data processing with the knowledge brought by the ecological model, in a theory-guided data science paradigm [START_REF] Karpatne | Theory-guided data science: A new paradigm for scientific discovery from data[END_REF]. In our case, consideration of the observation and ecological models as an integrated model allowed us to separate age classes from length frequency data as a function of the ecological processes (growth, mortality). Although the length distribution of trout of age 0+ usually clearly stand out from the rest [START_REF] Crozier | Interacting effects of density and temperature on body size in multiple populations of Chinook salmon[END_REF]Xu et al., 2010;[START_REF] Logez | Variation of brown trout Salmo trutta young-of-the-year growth along environmental gradients in Europe[END_REF], the distributions of older fish usually overlap due to growth dispersion and growth rate decreasing with fish age, making 'blind' separation more challenging, possibly leading to misclassification [START_REF] Pitcher | A bumpy old road: Sized-based methods in fisheries assessment[END_REF]FAO Fisheries and Aquaculture Department, 2013). The subsequent loss of fit caused by constraining the observation model (e.g., via the ecological model) is an opportunity to measure the discrepancy between observations and the assumed model and how much better/worse the constrained model is than the null model, thus serving as a rational guide for model improvement [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF][START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF].

We chose to conduct Bayesian posterior predictive checks for model evaluation and model checking, because they are effective at identifying poorly fitted models without requiring further data, although conservative when parameters are not estimated accurately [START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF]. Other options [START_REF] Conn | A guide to Bayesian model checking for ecologists[END_REF] for a review) include crossvalidation, still possible in the case of models and data with a hierarchical structure, e.g. using leave-one-out cross-validatory assessment or blocking, at a cost in terms of computational requirements [START_REF] Marshall | Approximate cross-validatory predictive checks in disease mapping models[END_REF][START_REF] Roberts | Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure[END_REF]. We directly compared catch (eq. ( 12)), although it would have been possible to compare the fit for each size class using some measure involving both the data and the model, e.g., Pearson χ 2 statistic, thus measuring the dissimilarity between the observed and modelled distributions in a manner close to a frequentist χ 2 test. The approach we used has the advantage of telling about the direction of the poor fits (under-and over-estimates) while the χ 2 statistic allows to measure the fit at different scales (by summing the χ 2 statistics, e.g., over size classes to get a measure at the survey level). A χ 2 type statistic can also be used to measure model fit for a wider range of models [START_REF] Gelman | Bayesian data analysis[END_REF][START_REF] Ntzoufras | Bayesian modeling using WinBUGS[END_REF][START_REF] Bal | A hierarchical Bayesian model to quantify uncertainty of stream water temperature forecasts[END_REF]. We did not use the deviance statistic as a metric for model fit [START_REF] Ntzoufras | Bayesian modeling using WinBUGS[END_REF] either, due to the externalized computations it requires in the case of a model with multiple error terms in the survey layer (in our case: Poisson and binomial). The scope of our model checking was to evaluate the fitness of the survey layer given the estimated length distribution. The model has five remaining random effects (see below), which were consequently not subject to examination in our case, although this would still be possible using replicated data and associated Bayesian p-values. The hierarchical framework offers other options than the measure of absolute fit that can serve as a rational guide for model improvement. We deem important mentionning the model comparison approach [START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF], either related to information theory such as the AIC [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF], which seeks to identify which model would be the most efficient data compression algorithm for the observed data, or fully Bayesian approaches (e.g., product space search;

see [START_REF] Tenan | Bayesian model selection: The steepest mountain to climb[END_REF]), which seeks to identify which of the alternative models has the highest relative credibility of being the true model, considering given data. In any case, measuring the absolute fit as well as model selection benefits from an integrated approach, by requiring formulation of the likelihood of model parameters given the data, which is a direct result of expressing the observation process ('external' errors) in the model.

The HBM framework also allows for 'internal' errors or random effects, which account for additional sources of variability. We did not use internal, additive errors for the cohort mean sizes predicted by the growth model, as Lecomte & Laplanche (2012) did.

Our model thus represents an appreciable advance compared to theirs. We still considered five random effects in the model, however: times of emergence (d emerg s,y

), density of trout of age 0+ (λ s,y,1 ), apparent survival rate (Sur s,y,k ), growth rate (G s ), and growth dispersion rate (ν s ). All of these terms have an ecological meaning, the estimated values of which are of great value by themselves. The drawback of having these random effects is that they prevent using the model in its current state for inter-or extra-polation, either temporal (e.g., forecasting) or spatial (e.g., to the stream continuum). All of these random effects, however, open the possibility of adding a connection to covariates, since the framework offers the ability to input spatio-temporal series of forcing variables into the model (e.g., water temperature).

We set parameters T min , T opt , and T max as known and constant, using values from laboratory experiments. The HBM framework makes it relatively easy to update the model to adjust and estimate values of additional unknown parameters, since data are informative, simply by defining these parameters as stochastic in the model [START_REF] Ntzoufras | Bayesian modeling using WinBUGS[END_REF][START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF]. The temperature range covered in our study area (0.7-23.2 • C) would make estimating these three parameters conceivable. We chose to use low informative priors for all of our parameters to facilitate a posteriori comparison of our estimates to those of other studies. Another option is to tighten parameter priors by using results of other studies (e.g., from [START_REF] Forseth | Thermal growth performance of juvenile brown trout \textit{Salmo trutta}: no support for thermal adaptation hypotheses[END_REF] for T min , T opt , and T max ). In any case, the flexibility offered by the HBM framework allows researchers to adjust and update their model as a function of the knowledge available. It also allows modelling parameters as functions of environmental covariates and random effects as residuals, using hyperparameters that are considered to be prefectly known, partially known and defined with narrow priors, or unknown and defined with vague priors.

Ecological results and discussion of main assumptions

We modelled the distribution of fish length as a mixture of Gaussian components, which is the usual option [START_REF] Pitcher | A bumpy old road: Sized-based methods in fisheries assessment[END_REF]. We have shown that the Gaussian mixture model directly resulted from our growth model, which assumed that there are inter-individual differences in growth rates and that fish individuals keep their advantages and disadvantages over their life-time (Appendix A). Other studies have highlighted the importance of inter-individual differences in growth traits and possible relative superiority or inferiority among individual fish [START_REF] Juanes | Ecology of stream fish: insights gained from an individual-based approach to juvenile Atlantic salmon[END_REF][START_REF] Peck | Inter-individual differences in rates of routine energy loss and growth in young-of-the-year juvenile Atlantic cod[END_REF][START_REF] Biro | Individual and sex-specific differences in intrinsic growth rate covary with consistent individual differences in behaviour[END_REF].

We chose to model inter-individual differences in growth traits as random variation in growth parameters, which is a standard approach (e.g., see [START_REF] Sainsbury | Effect of individual variability on the von Bertalanffy growth equation[END_REF] or [START_REF] Tang | A hierarchical Bayesian approach for estimating freshwater mussel growth based on tag-recapture data[END_REF] with respect to the von Bertalanffy model). Our model led to a theoretical relationship between the mean and standard deviation of the length distribution of each cohort (i.e. proportional) that fitted our data well. Other authors have considered a connection between the mean and standard deviation of length distributions, with CVs with a range similar to ours [START_REF] Lobón-Cerviá | Field assessment of the influence of temperature on growth rate in a brown trout population[END_REF][START_REF] Lobón-Cerviá | Density dependence constrains mean growth rate while enhancing individual size variation in stream salmonids[END_REF]. The latter authors, however, considered year-dependence, while we modelled site-dependence. Our assumption of site-dependence could thus be relaxed into some spatio-temporal variation (ν s,y ), offering the ability to further study growth dispersion.

We modelled catchability as a logistic function of fish size and site (as in [START_REF] Ruiz | A hierarchical model to estimate the abundance and biomass of salmonids by using removal sampling and biometric data from multiple locations[END_REF]). The reasons for this choice included consideration of the increase in catchability with increasing fish volume [START_REF] Dolan | Immobilization thresholds of electrofishing relative to fish size[END_REF] as well as a dependence on physical habitat. Some authors have considered other important covariates, such as time (e.g., due to discharge), stream width [START_REF] Letcher | Robust estimates of environmental effects on population vital rates: an integrated capture-recapture model of seasonal brook trout growth, survival and movement in a stream network[END_REF][START_REF] Kanno | Seasonal weather patterns drive population vital rates and persistence in a stream fish[END_REF], or removal [START_REF] Vøllestad | Growth-rate variation in brown trout in small neighbouring streams: evidence for density-dependence?[END_REF][START_REF] Laplanche | A hierarchical model to estimate fish abundance in alpine streams by using removal sampling data from multiple locations[END_REF]. Although of little impact regarding inference on growth, the relevance of the catchability model becomes crucial when interpreting estimates or when further modelling fish density (e.g., recruitment or mortality rates).

We considered a Poisson model for fish dispersion, which seems acceptable in view of the results of the posterior predictive checks, not issuing warnings with high numbers of under-and over-estimates. The model assumes that the distribution of the fish of a given size class in proximity of a sampled site is not spatially structured [START_REF] Peterson | On the estimation of detection probabilities for sampling streamdwelling fishes[END_REF]. While riverine salmonid do not show gregarious behavior, due to strong intraspecific competition, salmonid distribution can be patchy as a consequence of a spatial structure of the physical habitat. In such cases, it would be necessary to consider another statistical model for dispersion, e.g. using the negative binomial distribution, or expressing dispersion as a function of physical habitat covariates.

We chose to model growth rate with a power function instead of the more widely used von Bertalanffy growth function [START_REF] He | Modeling annual growth variation using a hierarchical Bayesian approach and the von Bertalanffy growth function, with application to lake trout in southern Lake Huron[END_REF][START_REF] Lecomte | A length-based hierarchical model of brown trout (Salmo trutta fario) growth and production[END_REF]. [START_REF] Elliott | Validation and implications of a growth model for brown trout, Salmo trutta, using long-term data from a small stream in north-west England[END_REF] have suggested that salmonid growth is not asymptotic and that non-asymptotic models should be used instead. Some studies have modelled Salmo trutta growth close to linear (power growth model with b = 0.31; [START_REF] Elliott | A new, improved growth model for brown trout, Salmo trutta[END_REF]; [START_REF] Jensen | Atlantic salmon (Salmo salar ) in the regulated River Alta: Effects 33 of altered water temperature on parr growth[END_REF]; [START_REF] Elliott | Validation and implications of a growth model for brown trout, Salmo trutta, using long-term data from a small stream in north-west England[END_REF]; [START_REF] Forseth | Thermal growth performance of juvenile brown trout \textit{Salmo trutta}: no support for thermal adaptation hypotheses[END_REF]), where we found a stronger curvature in this case study (b = 0.525). To our point of view, however, EFRS data may not be the most appropriate tool for comparing growth models and to investigate whether salmonid growth is asymptotic or not. The reason for this is that trout of age 2+ and older have relatively low density and their length distributions overlap [START_REF] Pitcher | A bumpy old road: Sized-based methods in fisheries assessment[END_REF], thus providing a low amount of information on trout growth. Information on individuals, either from laboratory experiments or in situ via capture-recapture [START_REF] Tang | A hierarchical Bayesian approach for estimating freshwater mussel growth based on tag-recapture data[END_REF], seems more appropriate. The choice of the empirical growth function becomes critical if estimating the time of recruitment from length data becomes a priority, however. Age-dependence of the growth rate due to gonad maturation and periodical changes in growth trajectories could be approached with a biphasic growth model, still applicable in a HBM framework [START_REF] Quince | Biphasic growth in fish I: Theoretical foundations[END_REF][START_REF] Dortel | Accounting for Age Uncertainty in Growth Modeling, the Case Study of Yellowfin Tuna (Thunnus albacares) of the Indian Ocean[END_REF]Armstrong & Brooks, 2013;[START_REF] Higgins | Modelling growth in fish with complex life histories[END_REF].

Finally, an alternative to empirical growth functions is to use a mechanistic, bioenergetic model, such as Net Rate of Energy Intake (NREI) models, which simulate growth of drift-feeding salmonids [START_REF] Piccolo | Food and space revisited: The role of drift-feeding theory in predicting the distribution, growth, and abundance of stream salmonids[END_REF][START_REF] Weber | Estimation of salmonid habitat growth potential through measurements of invertebrate food abundance and temperature[END_REF]. NREI models are, however, not applicable at large spatial scales due to their considerable data and computationnal requirements (e.g., [START_REF] Hayes | Process-based modelling of invertebrate drift transport, net energy intake and reach carrying capacity for drift-feeding salmonids[END_REF]; [START_REF] Urabe | Evaluation of habitat quality for stream salmonids based on a bioenergetics model[END_REF]). Computational requirements would also make such models impractical in a Bayesian framework, which requires simulation of the ecological model at each iteration of the MCMC sampler.

We found that the density of trout of age 0+ varied greatly among sites and among years, as usual with salmonids [START_REF] Milner | The natural control of salmon and trout populations in streams[END_REF][START_REF] Lobón-Cerviá | Spatial and temporal variation in the influence of density dependence on growth of stream-living brown trout (Salmo trutta)[END_REF][START_REF] Vøllestad | Non-additive effects of density-dependent and densityindependent factors on brown trout vital rates[END_REF]. The number of trout of age 0+ present at survey times results from the combination of three ecological processes: spawning success the year before, survival between spawning and emergence, and survival/movement (apparent survival) between emergence and survey times. Survival between spawning and emergence may be related to environ-mental conditions, due to high discharge damaging gravel redds before emergence [START_REF] Kanno | Seasonal weather patterns drive population vital rates and persistence in a stream fish[END_REF]. The same applies to apparent survival between emergence and survey times, due to high discharge flushing parr following emergence, thus causing mortality and downstream movement [START_REF] Jensen | The functional relationship between peak spring floods and survival and growth of juvenile Atlantic Salmon (Salmo salar ) and Brown Trout (Salmo trutta)[END_REF][START_REF] Lobón-Cerviá | Numerical changes in stream-resident brown trout (Salmo trutta): uncovering the roles of density-dependent and density-independent factors across space and time[END_REF][START_REF] Nislow | Towards a life-history-based management framework for the effects of flow on juvenile salmonids in streams and rivers[END_REF]. Spatio-temporal variations in stream water temperature (via snowmelt) and in precipitation (and resulting stream discharge) can thus cause large variations in the density of trout of age 0+ at survey times. Another reason is that adult spawning success and competition of trout of age 0+ following emergence may be density-dependent [START_REF] Milner | The natural control of salmon and trout populations in streams[END_REF][START_REF] Liermann | Using accessible watershed size to predict management parameters for Chinook salmon, Oncorhynchus tshawytscha, populations with little or no spawner-recruit data: a Bayesian hierarchical modelling approach[END_REF]. Further variation in apparent survival strongly depends on age and on season [START_REF] Lobón-Cerviá | Environmental determinants of recruitment and their influence on the population dynamics of stream-living brown trout Salmo trutta[END_REF]. Survival can decrease significantly with increasing age in some rivers as a result of angling pressure, thus decreasing the apparent survival rate with increasing age, as we observed. Our population dynamics model is basic in its current version and does not consider spawning success nor stage-dependent survival rates. However, the HBM framework offers the possibility to model population dynamics [START_REF] Bret | Understanding inter-reach variation in brown trout (Salmo trutta) mortality rates using a hierarchical Bayesian state-space model[END_REF].

We found that a non-negligible portion of the growth rate was unrelated to water temperature. Several studies have illustrated the ineffectiveness of water temperature alone to predict salmonid growth rate in the wild (Table 1). A more recent meta-analysis [START_REF] Kovach | Impacts of climatic variation on trout: a global synthesis and path forward[END_REF] showed that temperature 'was rarely related to growth'. On the other hand, laboratory studies have reported a strong relationship between observed temperature and growth (e.g., pseudo-R 2 > 0.99 in [START_REF] Elliott | A new, improved growth model for brown trout, Salmo trutta[END_REF]). The latter study, however, considered fish fed to satiation and consequently did not consider growth limitation due to lack of food. The effect of temperature on growth is strongly mediated by food consumption; the optimum temperature for growth varies considerably depending on ration size [START_REF] Brett | Growth rate and body composition of fingerling sockeye salmon, Oncorhynchus nerka, in relation to temperature and ration size[END_REF]Elliott, 1975a,b;[START_REF] Piccolo | Food and space revisited: The role of drift-feeding theory in predicting the distribution, growth, and abundance of stream salmonids[END_REF]. Since salmonids are highly territorial and juvenile and subadult fish feed on benthic and drifting macroinvertebrates [START_REF] Oscoz | Diet of 0+ brown trout (Salmo trutta L., 1758) from the river Erro (Navarra, north of Spain)[END_REF][START_REF] Johnson | Diel Resource Partitioning among Juvenile Atlantic Salmon, Brown Trout, and Rainbow Trout during Summer[END_REF][START_REF] Johnson | Comparative diets and foraging strategies of subyearling Atlantic salmon, brown trout, and rainbow trout during winter[END_REF], in situ food ration is likely related to intrinsic (conspecific density) and extrinsic factors (discharge, macroinvertebrate density, etc.), which are the factors that were found to be connected to salmonid growth (Table 1). Moreover, we found a strong connection between growth rate and catchment area, which is considered as an integrated metric of habitat capacity and incorporates habitat requirements at multiple scales [START_REF] Rosenfeld | Assessing the habitat requirements of stream fishes: An overview and evaluation of different approaches[END_REF][START_REF] Liermann | Using accessible watershed size to predict management parameters for Chinook salmon, Oncorhynchus tshawytscha, populations with little or no spawner-recruit data: a Bayesian hierarchical modelling approach[END_REF][START_REF] Ayllón | Modelling carrying capacity dynamics for the conservation and management of territorial salmonids[END_REF]. Our results thus corroborate those of other studies indicating that growth rate is related to both water temperature and food ration, and that this relationship remains predominant at a large spatial scale. We also conclude that further modelling of the growth rate should include a relationship to environmental variables related to daily food ration (e.g., discharge and macro-invertebrate density).

Incorporating other kind of data/information

Water temperature data are generally not available at EFRS sites, especially in the case of long-term historical EFRS surveys initiated for management purposes. One possible consequence of using a model to predict water temperature is integrating finer scale variation in water temperature and not detecting the resulting effect on growth, such as thermally heterogeneous stream waters that fish exploit [START_REF] Ruff | Temperature-associated population diversity in salmon confers benefits to mobile consumers[END_REF]Armstrong et al., 2013;[START_REF] Kanno | Paired stream-air temperature measurements reveal fine-scale thermal heterogeneity within headwater brook trout stream networks[END_REF]. We predicted water temperature from air temperature using a mixed-effects linear model, which is effective using monthly data [START_REF] Caissie | The thermal regime of rivers: a review[END_REF].

Making predictions at a finer temporal scale requires either a non-linear empirical model (e.g., 4-parameter logistic, see [START_REF] Mohseni | A nonlinear regression model for weekly stream temperatures[END_REF] or [START_REF] Baerum | Interacting effects of temperature and density on individual growth performance in a wild population of brown trout[END_REF]) or a mechanistic model, which requires more data [START_REF] Caissie | The thermal regime of rivers: a review[END_REF]. We therefore kept a linear model, in view of the large spatial (∼10-km between sites) and temporal (∼1 year) sampling frequency of our EFRS survey, which prevented us from evaluating the importance of the effect of temperature on growth at a smaller spatiotemporal scale and removed the need for water temperature values at a finer resolution. Moreover, our growth model numerically integrates water temperature values (eq. ( 10)), which simulates the ecological growth process in which fish integrate environmental conditions, thus decreasing the effect of daily temperature variation on growth. Other studies have shown that modelling water temperature as linear is a poor choice to capture a gradual shift in water temperature (e.g., [START_REF] Bal | A hierarchical Bayesian model to quantify uncertainty of stream water temperature forecasts[END_REF]). We detected a spatial structure of the growth rate in the downstream direction that we attributed to food availability rather than to temperature increase, which was consistent with results of other studies (as discussed). Our thermal model included the gradual shift in air temperature with elevation as well as a streamdependent thermal regime. For this reason, while we do not exclude the possibility that our model missed a portion of the gradual shift in water temperature in the downstream direction, it seems highly unlikely that the spatio-temporal variation of the growth rate that we observed was due exclusively to water temperature. In any case, studying the effect of fine-scale thermal heterogeneity requires considering temperature data at a similar scale [START_REF] Ruff | Temperature-associated population diversity in salmon confers benefits to mobile consumers[END_REF][START_REF] Kanno | Paired stream-air temperature measurements reveal fine-scale thermal heterogeneity within headwater brook trout stream networks[END_REF].

We focused on juvenile and subadult stream-dwelling salmonids in upper streams, with movements limited to switches between micro-and macro-habitat [START_REF] Schlosser | Stream fish ecology: a landscape perspective[END_REF][START_REF] Gido | Community ecology of stream fishes: Synthesis and future direction[END_REF][START_REF] White | A 'behaviorscape' perspective on stream fish ecology and conservation: linking fish behavior to riverscapes[END_REF][START_REF] Matthews | Stream fish community dynamics: A critical synthesis[END_REF][START_REF] Laplanche | Modeling the fish community population dynamics and forecasting the eradication success of an exotic fish from an alpine stream[END_REF].

This considerably reduces the impact of water temperature heterogeneity on growth.

While brown trout has high rates of site fidelity [START_REF] Budy | Exploring the relative influence of biotic interactions and environmental conditions on the abundance and distribution of exotic brown trout (Salmo trutta) in a high mountain stream[END_REF], other stream-dwelling salmonids can move over larger distances (e.g., cutthroat trout, see [START_REF] Hilderbrand | Movement patterns of stream-resident cutthroat trout in Beaver Creek, IdahoUtah[END_REF]). In such cases, a movement model should be considered, which would require specific data such a fish's successive locations using ITM [START_REF] Hilderbrand | Movement patterns of stream-resident cutthroat trout in Beaver Creek, IdahoUtah[END_REF][START_REF] Marvin | The Success of the Columbia Basin Passive Integrated Transponder (PIT) Tag Information System[END_REF].

As indicated, salmonid growth may be predicted more accurately by modelling the daily ration. Ration size depends mainly on fish size, macro-invertebrate drift density, conspecific density, and discharge [START_REF] Serchuk | Rainbow trout: a population simulation based on individual responses to varying environmental and demographic parameters[END_REF][START_REF] Hughes | Physiological ecology meets the ideal-free distribution: Predicting the distribution of size-structured fish populations across temperature gradients[END_REF][START_REF] Weber | Estimation of salmonid habitat growth potential through measurements of invertebrate food abundance and temperature[END_REF]. Fish size and conspecific density are intrinsic model variables, which are directly available in an integrated model. Discharge can be either directly measured in situ, or predicted using hydrological models (e.g., the catchment-scale SWAT; [START_REF] Arnold | Large area hydrologic modeling and assessment-Part 1: model development[END_REF]). In the latter case, land use, soil type, topography and climate data including precipitation are required as inputs, and measurement of discharge at the catchment outlet is needed for discharge calibration. On the other hand, macro-invertebrate drift density can be effectively sampled in rivers (Allan, 1987;[START_REF] Boyero | Spatial and temporal variation of macroinvertebrate drift in two neotropical streams[END_REF][START_REF] Hay | Macroinvertebrate drift density in relation to abiotic factors in the Missouri River[END_REF].

Salmonid spawning is triggered mainly by photoperiod and temperature (Jonsson & Jonsson, 2009) and is directly observable (timing and intensity) by monitoring spawning grounds [START_REF] Gallagher | Redd Counts[END_REF]. The time required for development of trout eggs from spawning to emergence is driven mainly by water temperature [START_REF] Ojanguren | Thermal dependence of embryonic growth and development in brown trout[END_REF]Jonsson & Jonsson, 2009), which can be used to make accurate predictions of the time of emergence (e.g., [START_REF] Elliott | An individual-based model for predicting the emergence period of sea trout fry in a Lake District stream[END_REF]). However, time of emergence is more difficult to observe in situ, due to the small size of the emerged fry. Predictions and observations of the times of spawning and emergence could be still incorporated into the model, in a HBM framework (e.g., [START_REF] Lecomte | A length-based hierarchical model of brown trout (Salmo trutta fario) growth and production[END_REF]).

We considered the parameters of the length-weight relationship as constant (R 2 > 0.99). Modelling seasonal growth variation, reproduction, or temporal variation in the food ration might additionally require considering time-dependent length-weight parameters [START_REF] Kimmerer | Variability in lengthweight relationships used to estimate biomass of estuarine fish from survey data[END_REF][START_REF] Froese | Cube law, condition factor and weight-length relationships: History, meta-analysis and recommendations[END_REF]. Moreover, studies in larger or more heterogeneous areas are expected to show larger spatial variation in length-weight parameters [START_REF] Froese | Cube law, condition factor and weight-length relationships: History, meta-analysis and recommendations[END_REF]. Spatial variation in length-weight parameters could be considered in this case, even in a Bayesian framework [START_REF] He | Modeling variation in mass-length relations and condition indices of lake trout and Chinook salmon in Lake Huron: a hierarchical Bayesian approach[END_REF], which would require additional measurement of fish weight.

As stated, growth is strongly correlated with critical life history events. Development of a growth model would serve the development of a population dynamics model by providing access to growth-related variables (e.g., size, spawning time, emergence time).

As an illustration, mortality from angling could be modelled with an exponential decay (see [START_REF] Serchuk | Rainbow trout: a population simulation based on individual responses to varying environmental and demographic parameters[END_REF]; Lobón-Cerviá et al. ( 2012) for natural mortality), which would be activated only during the angling season and to length classes which are above the minimum legal capturable size.

Conclusion

Presentation of our model and results in their current form, and even more our suggestions of model updates, both from growth and population dynamics perspectives, illustrate that using a HBM allows for (1) modelling of ecological processes, (2) quantification of measurement errors, and (3) links to covariates, resulting in (1) an increased range of ecological applications, (2) improved hypothesis testing, and (3) increased predictive power, which would allow researchers and managers to better understand a variety of salmonid ecology issues at large spatio-temporal scales. Coupling our modelling approach to a basin-scale hydrological model will expand the range of application of this HBM frame- Observed quantile (mm) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q Observed quantiles q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q site 110 site 1010 site 1620 site 1830 remaining sites 95% CI (n=375) sizes and densities are used to plot the modelled distribution of fish length at survey times (e.g., Figure 1 for 1997). Site 110 was chosen here as an example.

ν s λ 1 b Sur k σ Sur k G s,y,k X s,d T w s,d T min T opt T max d emerg s,y σ G s G s

  els becomes a daily time step in the HBM, indexed with d ∈ {1, . . . , D} spanning the Y years that are considered in the model. Daily mean water temperatures are denoted T w s,d in the HBM and growth parameter G y (t) becomes G s,y,k,d . The ordinary differential equation (

  ) P-values lower than 0.05 highlight underestimated values and values greater than 0.95 highlight overestimates. P-values of the posterior predictive checks are uniformly distributed if the model fits correctly[START_REF] Marshall | Approximate cross-validatory predictive checks in disease mapping models[END_REF]. We thus compared with another series of Q-Q plots the distributions of the p-values for each site to their theoretical uniform (0,1) distribution. Computation of the posterior predictive p-values is also relatively easy to perform with BUGS[START_REF] Ntzoufras | Bayesian modeling using WinBUGS[END_REF];[START_REF] Lunn | The BUGS Book: A practical introduction to Bayesian analysis[END_REF]; Appendix S2).2.4. Study area and data sets2.4.1. Study areaThe study area represents the natural range of brown trout (Salmo trutta) distribution in the region of Navarra (northern Spain, 0 • 43'-2 • 29' W, 41 • 54'-43 • 19' N). This area of nearly 6,420 km 2 is geologically heterogeneous (http://geologia.navarra.es) and drains northwards into the Bay of Biscay and southwards to the Mediterranean sea via the Ebro river. The Mediterranean drainage of the study area can be further split into two subbasins that show a clear east-west altitudinal gradient (Figure3). Elevation in the study area ranges from 0 to 2,444 m.a.s.l.Brown trout is the dominant fish species throughout the study area, and its populations consists exclusively of resident individuals (except for lower reaches of the Atlantic Bidasoa basin, where anadromous individuals exist at low densities). Rivers are open

Figure 2 :

 2 Figure 2: Directed acyclic graph (DAG) of the growth and population dynamics components of the HBM. Forcing variables and known parameters (in magenta) are connected to top-level parameters (yellow) via intermediate nodes (white). Variables are either deterministic (rectangles) or stochastic (ellipses) expressions. Variable equations/values are gathered together in Table2. Most variables are multi-dimensional, as indicated by their subscripts and overlapping colour frames, one frame per index: site s (red); year-of-emergence y (orange); age k (blue); day d (green); variables outside frames are scalar.Brown trout length is modelled as a mixture of Gaussian distributions, which parameters (means µ s,y,k , standard deviations σ s,y,k , and contributions λ s,y,k of each Gaussian component; rectangles with thick borders) are connected together with the growth and the population dynamics model (as represented) as well as to observations (not shown in this DAG, see text).

Figure 3 :

 3 Figure 3: Sampling design and study area. The study area (6,420 km 2 ) is the northern section of Navarra (thick black line: regional border) at the most western part of the Pyrenees mountain range (dotted area) in northern Spain (thick grey line: national border). The study area separates into 3 regions (grey areas): Atlantic to the Ega River; West Mediterranean to the Arga river; East Mediterranean to the Aragón river. A total of 61 sites (0110, ...) have been sampled by electrofishing once a year from 1992 to 2014. Only data from rivers unaffected by the presence of stocked individuals(48 sites; 1995-2014) were considered in the present analysis (black dots). Water (green dots) and air (orange) temperatures were also sampled in the area.

Figure 4 :

 4 Figure4: QQ-plot of the predicted and observed distributions of fish length. The quantiles of the predicted (x-axis) and observed (y-axis) multimodal distributions of fish length data (Figure1) are plotted against each other for a graphical comparison of the two probability distributions. Quantiles were computed for probability values between 0.05 and 0.95 every 0.05 (one dot per probability interval), for each year (colour; see right legend), and for each site. Dot size (see left legend) is proportional to the total number of fish that were caught during surveys, thus highlighting scarce data. Site 110 was here chosen as an example.

Figure 5 :Figure 6 :

 56 Figure5: QQ-plots of the Bayesian p-values of the posterior predictive checks. The quantiles of the theoretical (x-axis) and observed (y-axis) distribution of the p-values of the posterior predictive checks are plotted against each other for a graphical comparison of the two probability distributions. P-values of the posterior predictive checks are uniformly distributed if the model fits correctly; the theoretical quantiles are thus those of the uniform (0,1) distribution. Quantiles were computed for probability values between 0.05 and 0.95 every 0.05 (one dot per probability interval), for each site (black lines; site 110 chosen as an example in Figures1, 4, and 6 highlighted in orange). The number of p-values used to compute the QQ-plots ranges from 183 to 523 (375 ± 72.6), depending on the site. The model fits well, to the exception of site 1620 (green) with a larger number of significant overestimates (Bayesian p-values > 0.95), and sites 1010 (red) and 1830 (blue) due to extremely low fish density. The 95 % confidence interval on the estimates of the quantiles of a theoretically uniform (0,1) distribution using n = 375 samples is plotted (magenta) as an indication of the variation in the QQ-plots due to random sampling.

Table 2 :

 2 Variables of the HBM. Most variables are multi-dimensional, as indicated by their subscript (s, y, k, d, l, r; Figure 2). Deterministic variables (=) are either measured without errors, known constant, or deterministic expressions from upper nodes. Stochastic variables (∼) are either stochastic expressions from upper nodes or top-level random variables, in the latter case priors are indicated.

	work, including the assessment of potential global change impacts on fish population
	dynamics.

Acknowledgements

This study was undertaken under the research collaboration agreement between DEGN and Ecolab. Authors are grateful to all the DEGN staff involved in the collection and management of electrofishing data since 1992, with special thanks to Javier Álvarez. Access to the HPC resources of CALMIP was granted under the allocation P1113. The authors acknowledge both anonymous reviewers for their insight and constructive comments which considerably improved the quality of this paper.

Appendix A: Derivation of the growth model at the population level

The growth of a fish individual is modelled as a consecutive length increase since its emergence

where i is an index on fish individuals, L i (t) is the length of individual i at time t, H s,y (t)(1 + i ) is its growth rate, L emerg i is its length at emergence that took place at time

. We assume that all the fish follow the same growth scheme, that is to say H s,y (t) depends on the site s, on the year-of-emergence y, and on time t but not directly on i. On the other hand, we allow individuals to have superior/inferior growth rate with respect to each other and we assume that individuals keep their advantage ( i > 0) or disadvantage ( i < 0) over their life-time (similar to [START_REF] Sainsbury | Effect of individual variability on the von Bertalanffy growth equation[END_REF]

As a consequence, under such assumptions, fish length can be modelled as a mixture of Gaussian distributions (eq. ( 6)), the mean size of each cohort (µ s,y (t)) is dictated by as similar growth model (eq. ( 1)), and the standard deviation of the length of each cohort is proportional to its mean (eq. ( 5)). The coefficient of variation (CV) is in the latter case the standard deviation of the variate i .

Supporting information

Additional Supporting Information may be found in the online version of this article.

• Appendix S1. Growth curves and temperature-dependence.

• Appendix S2. Model code and tutorial.

• Appendix S3. Water temperature. ). The ANOVA sums of squares (SSQ; here in %) of log(X s,d ) (using site, year, and month as factors; site*day refers to residuals) and log(G s,y ) (site and year; site*year refers to residuals) highlights the major sources of variability of these terms. Spatial variation of G s,y is represented in Figure 7.

E N W S 0 10 20 km q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 3015 3155 q q q q G s 1 2 3 4 5 6 Figure 7: Spatial variation in the growth rate. Daily growth rate was modelled as the product of 3 terms: size-dependence, temperature-dependence, and other sources of variation (G s,y ). The mean (denoted G s ) at each site of G s,y is represented on the map of the study area; G s increases in streams in the downstream direction. Results showed that G s is strongly correlated with catchment area (r = 0.73).