Is there a trend in extremely high river temperature for the next decades? A case study for France

. After 2003’s summer heat wave, Electricit´e de France created a global plan called “heat wave-dryness”. In this context, the present study tries to estimate high river temperatures for the next decades, taking into account climatic and anthropogenic evolutions. To do it, a speciﬁc methodology based on Extreme Value Theory (EVT) is applied. In particular, a trend analysis of water temperature data is done and included in EVT used. The studied river temperatures consist of mean daily temperatures for 27 years measured near the French power plants (between 1977 and 2003), with four series for the Rhˆone river, four for the Loire river and a few for other rivers. There are also three series of mean daily temperatures computed by a numerical model. For each series, we have applied statistical extreme value modelling. Because of thermal inertia, the Generalized Extreme Value (GEV) distribution is corrected by the medium cluster length, which represents thermal inertia of water during extremely hot events. The µ and σ parameters of the GEV distributions are taken as polynomial or continuous piecewise linear functions of time. The best functions for µ and σ parameters are chosen using Akaike criterion based on likelihood and

HAL Id: hal-00330927 https://hal.archives-ouvertes.fr/hal-00330927Submitted on 15 Oct 2008 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not.The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

Introduction
This study is performed for applications concerning French nuclear thermal power plants cooling systems in a context of industrial safety checking.After the 2003 heat wave, it appeared necessary to re-evaluate the extreme hot temperatures, which the rivers water could reach in the next decades, by taking the climatic and anthropogenic changes into account.
Indeed, the river temperature fluctuations are forced by natural and anthropogenic factors, e.g.industrial thermal discharges or hydraulic installations that modify natural flow.
Among the natural causes, we can identify hydraulic forcing of tributaries or ground waters, which have their own water temperature and discharge.In a similar way, the atmospheric forcing of air temperature, moisture, precipitations, cloud cover and wind, take part in the river thermal balance.
However, it is now recognized from IPCC work that global mean annual temperature increased during the last century.Trends in summer mean air temperature in France have been identified too (Mestre, 2000 and IMFREX project).Similarly to the mean, hot extreme air temperatures appear to present positive trends with increasing heat waves events (Laurent et al., 2007).This atmospheric forcing influencing the hot river temperatures has a high impact as it generally occurs during dry periods of low river flows.
For this reason, it seems interesting to determine whether the river energy balance during heat wave episodes is a stationary process or not.This study examines this question, as efforts to assign the detected trends to natural or anthropogenic causes is a delicate issue.This point will be discussed more thoroughly later in the paper.Only few studies were carried out on this subject with applications to river ecosystems (e.g.Sinokrot et al., 1995or Gooseff et al., 2005) as well as for industrial facilities dimensioning.This trend analysis which is a goal of the paper, is then applied to evaluate the hot extreme temperatures in the future.In our study, as design values are usually defined as return levels, we propose a statistical model based on the extreme value theory.This approach has the advantage of providing a rigorous mathematical framework for extreme values with temporal variability (in the case of a non stationarity assumption) and the random fluctuations of these values.This theory introduced by Gumbel (1958) within a stationary framework, then largely enriched later on (Coles, 2001) was applied in hydrology (case of high water level; Katz et al., 2002) or in oceanography (high sea level close to ports, Coles, 2001).In a non stationary context, extreme value theory reveals itself as a powerful tool to detect changes and provides frequency analysis (Zhang et al., 2004).Recently, the extreme value theory was used to evaluate hot air temperatures in France with temporal variations (Parey et al., 2007).Moreover, this theory provides quantitative results in the form of return levels, which could be used for example for industrial work dimensioning.Extremes can be modelled through two approaches.The first one is based on the Generalized Extreme Value distribution (GEV), i.e. the limit distribution of the maxima of a sequence of independent and identically distributed random variables.The second one is the Peaks Over Threshold method (POT) which looks at values above a high threshold, the exceedances being fitted to a Generalized Pareto Distribution (GPD).
Several series of water temperatures near 12 power plants are available for our study.In Sect.2, we present the data used.Then in Sect. 3 we describe the applied methodology.In particular, we explain how the non stationarity assumption is introduced and how the GEV distribution is corrected to take the thermal inertia of the river into account.The results on the detected trends and the associated return levels are presented in Sect. 4. Finally, we conclude on the obtained results.Still, those must be considered with caution especially when they are applied in a context of public interest or safety.

Presentation of water temperature time series
Water temperature is observed upstream and downstream from several measuring sites found on the river bank near nuclear power plants in France.For the purpose of our study, we will be interested only in water temperatures upstream of the power plants.

The observed temperature series
Series of daily averages were computed from hourly temperature measurements.In the presence of missing or erroneous measurements, temporal interpolation (linear or Fourier series) or regression starting from the nearest station were performed on each series to obtain complete chronology.Thus, we have complete series of daily average temperatures measured in water, for the period 1977-2003, except for Chooz whose data is far too incomplete and unreliable.For this last site, the water temperature sensor has been installed in an inappropriate manner to measure the upstream temperature accurately.For this reason for this power plant the series of measured temperatures is replaced by a series of computed temperatures, (see Sect. 2.2).Table 1 details the percentage of missing data that have been replaced by interpolation or by regression and the longest period of missing data, for each site.

The computed temperature series
The series of water temperatures for the Chooz site is computed and averaged over 24 h from 3-hourly data produced by a numerical model.The CALNAT model calculates the temperature in a point of the river using a deterministic way, integrating the equation of the temperature evolution (see Gras, 1969): Where, U river speed K thermal diffusivity along the river ρ mass of water per unit volume C specific heat capacity of water H depth of river thermal inertia.The five thermal fluxes SR, AR, WR, C and E are caused by solar radiation, atmospheric radiation, water radiation, wind convection, and evaporation, respectively.The flux of solar radiation depends on the albedo, the cloud cover, the location (latitude and longitude) and the date (day and hour of the year).The flux caused by the atmospheric radiation is computed using the air temperature and the cloud cover.
To calculate the water radiation, the knowledge of water temperature is needed.In our context, the convection describes the wind cooling on water surface.To estimate it, wind speed, water and air temperature must be known.That explains why Eq. ( 1) is solved implicitly.
Finally wind speed, air temperature, air humidity, atmospheric pressure and latent heat of vaporization are used to calculate the thermal evaporation flux.Atmospheric data are taken from the nearest Meteo-France observation station.
To sum up, the temperature evolution is forced by weather data (air temperature, air humidity cloud cover and wind) and needs a specific initial condition and boundary layer conditions.More details about the CALNAT model are given in Dupeyrat et al. (2006).
Similarly, series of water temperature are formed over the period 1949-2003, for the power stations of Belleville and Dampierre.
The validation of the model is described in Dupeyrat et al. (2006).It consists in comparing computed and measured temperatures when the latter are available.This validation is done over the periods 1978-1991and 1998-2004for the Chooz site, and 1977-2003 for the Belleville and Dampierre sites (see Fig. 1 for example).

Locations of the temperature series
Four series are located along the Rhône river (Bugey, Saint-Alban, Cruas and Tricastin), four along the Loire river (Belleville, Dampierre, Saint-Laurent and Chinon) and some others on other French rivers (Chooz, Cattenom, Fessenheim and Golfech).Figure 2 shows their locations.

Methodology
3.1 Extreme value method 3.1.1Choice of statistical extreme value method Some preliminary tests have been performed with the POT method.We need a time series of around hundred water temperatures above a given threshold in order to be able to fit a GPD model reliably.To obtain them, we have to take between 15% and 60% of summer data (depending on the river time series).Threshold specified in this way are thus far too low to select only extreme temperatures and GPD conditions will not be fulfilled.Therefore these pre-processing analyses showed that the GEV distribution is more suitable for the study of water temperatures.

The generalized extreme value distribution
For independent observations arising from a common continuous distribution, the probability for not exceeding a level z for a temperature T in any block of days can asymptotically be approximated by the cumulative distribution function of a GEV distribution. Where: The parameters µ and σ are proportional to the mean and to the variance of extremes respectively.The proportionality  coefficient depends on the shape parameter ξ .Of course in a non stationary case, these parameters depend on the considered block.This function is known as the generalized extreme value distribution, because it includes three types of distribution function, according to the value of ξ .The GEV distribution reduces to the Gumbel family for ξ =0.It belongs to the Fréchet family for ξ >0, and the Weibull family for ξ <0.
Parameter ξ characterizes the shape of the probability density function g.Indeed, according to its value ξ =0, ξ >0 and ξ <0, this function is not bounded, bounded from the left by z lim =µ− σ ξ and bounded from the right by z lim respectively (see Fig. 3).

Data description and hypothesis verification
We work with daily averages of water temperatures over the whole year (circle in Fig. 4) upstream from the power plants and over various periods.The temperature values are selected from 90-day hot seasons (square in Fig. 4), fixed between the 20 June and the 17 September.Thus, we suppress the seasonal cycle.From this hot season series noted (H n ), a sequence (T n ) of hot extremes is extracted by taking the temperature maximum per time block b of 30 consecutive days (triangle in Fig. 4).The extreme value theory (Coles, 2001) shows that the probability for not exceeding the value T =z is asymptotically well approximated by a Generalized Extreme Value (GEV) distribution.In order to apply the Generalized Extreme Value distribution for the sequence (T n ), two fundamental assumptions must be fulfilled: -The sequence (H n ) must have the same probability distribution along the fixed time block, but not necessarily for different blocks or years.The selection of the hot season has been made in such a way.
-The sequence (H n ) must consist of independent or weakly dependent values in order to apply GEV approximation.If they are not, some correction have to be done on the limit law using the clusterization index θ.In order to verify this point, the cluster lengths have been studied.

Episodes length
In order to study the independence of high values, the clusters are defined by a set of couples (date, temperature) whose temperature is higher than a given threshold T =u.Thus, two clusters are at least separated by one value (daily mean temperature) below the threshold.The cluster rate θ can then be estimated by the ratio of the number of clusters to the number of values above the threshold u.
Temperature is more clustered for weaker θ. θ can be regarded as an estimator of the thermal inertia of river water temperatures for the threshold u, estimated as the inverse of the average length of clusters.
The results in Table 2 are found by taking as thresholds the values of the 92nd percentile for all the hot seasons.This threshold has been found to represent a good compromise between a high enough level and a selection of a significant number of episodes for reliable statistics.
Maximum temperatures for each site show that 2003 is an unprecedented heat wave in intensity over the 1977-2003 period (see Table 2).It is also the case for the cluster lengths except at St. Laurent (19 days from 22/7/1994 to 9/8/1994), The results of Table 2 show an important clustering of the water temperature.However, the sites on the Loire river are less clustered than those positioned on the Rhône river.Cattenom and Chooz sites present an intermediate clusterization between that of the Loire and the Rhône.On average, Fessenheim and Golfech are the most clustered series.These differences are linked to different flow characteristics of the river.
Thus, observed river temperatures show a certain thermal inertia as the duration of a very hot temperature episode can reach almost a whole season.This is why it is necessary to correct the GEV distribution with the factor θ, (see Leadbetter et al., 1983):

Temporal evolution and trends of the parameters
To estimate the hot extreme temperature in the future, it is necessary to identify trends, which may be present in observational data of extreme water temperatures, and compute changes in the involved return levels.Therefore, the parameters of the GEV distribution are supposed time-dependent and the form of this evolution is modelled as polynomial or continuous piecewise linear functions.Only parameters µ and σ are considered to vary with time.We suppose that the shape parameter can reasonably be kept constant with time.
As a matter of fact, ξ is the most delicate parameter to estimate and as our series are quite short, trends could not easily be identified.In this non stationary context of time varying µ and σ parameters, the return levels need to be re-defined.

Definition of return level in a non stationary context
If we define the return level z as the level reached or exceeded on average once over the time period P ′ for block maxima, the knowledge of the quadruplet (µ (t) , σ (t) , ξ, θ ) is sufficient for the estimation of the return level.Indeed, z can be statistically defined as the level whose number of exceedances N has an expectancy E(N) equal to one, during the various hot seasons of P ′ .This is true for any block b belonging to t ′ 0 , t ′ 0 + P ′ with initial time t ′ 0 corresponding to the first future year after the end of the series.Thus, in our case we estimate the level reached or exceeded on average once between 2004 and 2004+ P ′ −1 .
But the number N of exceedances of level z in time's interval t ′ 0 , t ′ 0 + P ′ can be derived as: With m, the number of blocks corresponding to each hot season of P ′ , And Y b , the variable related to a given block b, defined as: Therefore, the expectancy of N is given by:  4a) and (4b), the expectancy of N equal to one can be written as: With the indicator function δ such as δ=1 if 1+ξ (z−µ) σ >0 δ=0 otherwise.
Therefore, the return level z can be computed resolving Eqs. ( 7a) and (7b)

Estimates of the GEV parameters
To take the climatic evolution into account, the parameters µ and σ are regarded as functions of time t b .
Several time models are then considered for µ and σ : -polynomial models for µ=µ (t b ) and σ =σ (t b ), i.e. as: Thus, these models are defined by the couple (i, j ) of the order of the polynomial equation for µ and σ respectively.The optimal polynomial could be selected by the maximum likelihood test whose asymptotic distribution is a chi square with a number of degrees of freedom defined as the difference in the number of parameters (Davison, 2003).
Still, there is a technical problem, because this test provides results only for nested models.For instance, it does not allow to choose between polynomials defined by degrees (1,2) and (2,1).But this situation does not occur for our data.
Another and perhaps simpler way to select the optimal polynomial (noted polyOPT) is to use the Akaike type criteria.Let L(n, i, j ) be the log likelihood, n being the sample size.
The penalized likelihood is A (n, i, j ) =L (n, i, j ) −ϕ (n) (i+j ).Then, we choose i and j maximising A(n, i, j ).A popular choice of the function ϕ(n) is ln(n), giving a consistent estimate for i and j (Davison, 2003).
-Another choice of parameterisation of µ consists of continuous piecewise linear models (CLM) with zero, one, two or three parts (respectively noted 0P, 1P, 2P et 3P).
They are defined such as: where a k are the angular points for i>k>0 and a 0 =t 0 (initial time of the series) and a i =t 0 +P with the series duration P .
To respect the continuity between two linear pieces, we have: The models in zero and one part correspond to the stationary (for µ and σ ) and linear (for µ) models respectively of the polynomial modelling.
Because of high computational time, the CLM were not studied for σ .Thus σ is supposed to be stationary, which is a plausible assumption according to the results of polynomial modelling.
The Akaike criterion (see above) and localization conditions of the angular points determine the optimal CLM model (noted clmOPT).Indeed, the classical theory of the likelihood tests can be applied only if separation constraints are imposed between the angular points.The same constraint needs to be applied to CLM edges as well.
These time dependent models will be used for prediction.More precisely, they will be extrapolated to predict the values of the parameters in the future and thus also to estimate the return levels.The validity of such an extrapolation is a strong hypothesis.Of course it requires choosing for these parameters an analytical form easy to estimate (as a polynomial form of low degree).Moreover, the resulting extrapolation should have physical sense.This approach seems to be the best way to perform this prediction.The main reason is that there does not exist any model (as a filter) giving scenarios for the variables (CO 2 concentration) linked to the greenhouse effect and its consequences for the rivers temperature.This problem is in a sense more complex than the one concerning air temperatures.

Model choice
Tables 3, 4 and 5 give the results obtained for both types of modelling using likelihood maximization.They also provide the model selected in fine.The models are defined by the couple (i, j ) of the order of the polynomial equation for µ and σ respectively.Continuous piecewise linear models with i parts are noted iP .
The two types of modelling show trends for all sites.With polynomials, the effect of variance appears only on St. Laurent and Cattenom with a negative trend.Thus, there does not seem to be generalized trend in variance as in the case for mean.Furthermore, as the samples are quite short, trends in variance are difficult to be accurately estimated and it then has been decided not to take them into account.Return levels are then evaluated using extrapolation of identified trends.Of course, for the location parameter µ trends have to be monotonous on a large time interval at the end of the observation period, and degrees above 2 cannot be retained.For piecewise linear functions, we also have to check that the last break point is not too close to the end of the series, which could reveal a strong influence of the last hot year 2003, leading to unrealistic extrapolations.Moreover, a coherence of the trends between close sites located on the same river and which do not have a distinctive hydrological feature is maintained.This coherence could be spread along the entire river, in the case of the Rhône.The cases of the Loire and Rhône rivers will be described in more details.

The Loire river
Figure 5 shows the trends computed for µ, in Belleville and Dampierre: the best polynomial obtained is (3,0) and the CLM model is 3P They correspond to the best statistical fit to describe the data, but regarding the shortness of the series, they cannot reasonably be used for extrapolation.As mentioned above, these models are very sensitive to the 2003 heat wave placed at the end of the sample.That remains true for all the sites.Thus, only extrapolation starting from linear and 2P CLM models can be considered as reasonable.The 2P model has the advantage of separating the first years of relatively cold summers from the last part of the trend, which can then more reasonably be extrapolated.Natural climatic variability is here sampled in an uncomfortable way for robust trend identification, since relatively cold years are frequent at the beginning of the sample and a particularly hot year occurs at the end of the observed period.In addition, as the water temperature adjustment with the atmospheric conditions is fast on this river, various trends on µ can be retained for sufficiently distant sites.The Rhône river Figure 6 presents the trends for µ, on various sections of Rhône river: Bugey, which is located rather upstream, and Cruas located rather downstream.On the two sites again, a polynomial extrapolation starting from models in (3,0) or in 3P CLM on µ, are not reasonable.These high order trends give a good fit with the observations but cannot be considered as representative of the global trend.A good example of a correct fit but with little physical sense is the characteristic model shape of the 3P.Even the 2P model shows a great influence of 2003 on the trend in Bugey, but not in Cruas.On the other hand, the linear trend seems to have a more reasonable physical meaning on the two sites.
For the Rhône, the river thermal behavior is different from the Loire.In particular, the river thermal balance does not depend only on atmospheric forcing but also on the glacial contributions (lake Léman and glacial tributaries) and on anthropogenic factors.In this complex and slow thermal balance, it appears necessary for the river to preserve the same modelling of µ for all the sites.

Discussion on trend models
According to the sample, the situation is very specific.CLM models are used in order to control and perhaps to improve polynomial models.But, the data present a very hot temperature at the end of the sample.Then, the last break point has a strong probability to be near to this end and therefore the last piece can have a very high slope.For a fit and descriptive aim, CLM is well adapted.But for a prediction by extrapolation, CLM could be lead to huge values which cannot easily interpretable.
To sum up, a positive trend for the evolution of µ is found for all the power plants.However, it is advisable to interpret this trend with great care, as the water temperature is the resultant of several causes, natural or not.
That is particularly true on the Rhône river, which is a socalled "transfer river".The heat contributions of industrial origin is not evacuated (or little) in the atmosphere and modify the water temperature all along the river downstream of the disposal (see F. Hendrickx et al., 2004).Thus, the detected trend includes both natural and industrial evolutions during the last 27 years.For the Loire, the trend is less complex because the adjustment of water temperature to the atmospheric conditions is faster.Thus, the river temperature is not influenced by the upstream thermal discharges.
It is also the case for Golfech, Chooz, Cattenom and Fessenheim, which do not have major industrial facilities carrying out upstream important thermal discharges.
On the other hand, the upstream hydraulic installations can influence the water temperature by modifying the flows.
For all these reasons, it is difficult to attribute trends only to climatic evolutions, even if they should have an influence.

Results on the GEV parameters
Figure 7 and Table 6 give the results of the GEV parameters for the chosen models.The value of the cluster rate θ is added too.
For each parameter the standard error is similar for all stations.On average, standard errors for ξ , the last slope of µ and σ are equal to 0.07, 0.27.10 −4 and 0.13 respectively.Thus, the 95% confidence intervals for these parameters are about equal to 0.14, 0.53.10 −4 and 0.25 around the estimated values.
With values of ξ ranging between −0.570 and −0.126, the distribution governing the extreme water temperatures belongs to a Weibull type for all the stations with a 70% confidence interval.Thus, the probability density of the extreme temperature levels is bounded, which seems physically reasonable in the case of hot extreme water temperature.
In the overwhelming majority of stations, no significant trends are detected in the σ parameter of the GEV, and because of the shortness of the samples, the identified ones must be carefully considered.For applications to return levels computing, it has been decided to not take them into account.For the parameter µ, a positive trend is detected for all stations, whatever the temporal evolution model chosen for this parameter.The slope coefficient of the last part of the model for µ lies between 2.1×10 −4• C/time and 4.75×10 −4• C/time.The 10-, 20-and 30-year return levels seem to be physically reasonable.Indeed, the 10-or 20-year return levels were already reached, sometimes exceeded in particular during the 2003 heat wave, on all sites.For Chooz, even the 30-year return level was approached.On the other hand, the 50-and 100-year return levels do not appear very reasonable in particular in comparison with the sample size (27 years, sometimes less when using piecewise linear functions) on which the extrapolation is based.
To evaluate the trend effects, return levels were computed with a stationary assumption.In this case, µ, σ and ξ parameters of GEV distribution are supposed to be constant.The results are given in Fig. 9 in the same form as in the non stationary case (see Fig. 8).
We also compute confidence intervals for the different return levels.The non stationarity hypothesis makes it impossible to use close asymptotics.Thus, we simulate 1000 trajectories of the estimated optimal model.Then, the return levels are computed for every trajectory.From these data, we obtain an empirical estimate of the return levels distribution and we choose as confidence interval an interquantile interval, for instance (q(0.05),q(0.95))for the risk level 0.1.We also check that the magnitude of these intervals does not depend on small variations of the estimated model parameters.Table 7 shows different confidence intervals computed using trajectories, for Saint Laurent.Results are almost the same for all stations, even if confidence intervals are narrower for the Rhône than for the Loire.Thus, we do not give the detailed results.The orders of magnitude are as follows: -The width of the confidence interval is 2 degrees at risk level 0.05 for a 10-year return level.
-It is divided by 2 for a risk level 0.3.
-It is multiplied by 1.8 for a 30-year return level.

Influence of θ on return level
The difference of return levels according to whether the series are considered stationary or not, increases with return   level duration.On average on all the stations, this difference is between 1.8 • C for 10 years and 3.5 • C for 30 years.For 30-year return level, the temperature variation lies between 2.5 • C and 5 • C according to the stations.
Figure 10 shows the differences between the 30-year level with or without θ in the stationary and non stationary contexts.Globally, these differences decrease with θ as ex-pected.Indeed, the probability (1-G) of exceeding a fixed return level is a decreasing function of θ .
For all return levels, the differences vary between 0.2 • C and 2 • C. Finally, if we want to take temporal evolution entirely into account, the θ parameter should be a function of time.But, in that case the problem becomes difficult to solve, as our series are much too short to allow a reliable study of evolution of clustering.This issue could guide a future development of our methodology.

Results for the long time series
In order to analyse the effect of the sample size, two long temperature time series computed with a thermal model of rivers behaviour, are used (see paragraph 2.2).For important computing time reasons, only the polynomial modelling of the parameters µ and σ has been tested with these 55-year series.The main results are given in Table 6.
The optimal polynomial models are (1,1).The variance effect becomes important for these series and tends to increase the return levels compared to the model (1,0).
For the same modelling of parameter µ, the return levels estimated from 1977-2003 and 1949-2003 series show significant differences, which increase with the return period.This sensitivity test shows quite well that there is a strong sampling effect on the return level calculation.
However, Table 9 shows that the shape parameter ξ is of the same order of magnitude for the 27-year and the 55-year series, which supports the idea of the time independence of ξ .
The choice of the computation period, although crucial for the extreme temperature evaluation, is delicate.Indeed, it is difficult, or even impossible, to determine with exactitude the beginning of the climatic change effects on water temperatures.
However as we said before for the Loire, the thermal equilibrium with the atmospheric forcing is fast because of low flow and small depth.Thus, the thermal influence of hot water rejected by upstream power plant is weak.The analysis of monthly mean temperature of the Loire river using a series since 1976 and a series since 1881 shows a temperature increase, with a significant acceleration since the late 1980s (Moatar, 2006) due to the rise in air temperatures and to lower discharge rates.
Besides, the trends observed for µ for the average and extreme air temperature, show an acceleration from the years 1970 (Parey et al., 2007).
For these reasons, to calculate short term return levels, it could seem wiser to carry out calculations by extrapolating the most recent trends and thus starting from the series based on the last 30 years.

Conclusions
This paper proposes a methodology to identify if trends in extremely high river temperatures exist for the next decades.
Several models has been studied to estimate the trends.As application, one of these trend models is used to evaluate the return levels of extreme hot water temperature.We studied the trends of extremely hot water temperatures at 12 French stations, based on measured or computed (for one station) series over the period 1977-2003, using statistical extreme value theory, and particularly the GEV distribution.The GEV used for the return levels calculation takes the thermal water inertia into account via the θ corrective factor.The distribution modelling the extreme series is found to be of Weibull type, with a negative shape parameter.Thus, there is an upper bound for the extreme temperatures, which seems physically reasonable.
All these stations show an increasing trend in the mean level of extremes.On the other hand, we do not identify a reliable trend in the dispersion of the extreme values.This identified trend takes at the same time climatic change and evolutions of human induced activities along the rivers into account, particularly for the Rhône river.Analyses of longer temperature time series  confirm a sample size effect in the trend calculations.This is why it is advisable to take these results with care especially if they must be applied with the aim of public interests and safety.1949-2003 and 1977-2003.Belleville (1,0) Dampierre (1,0) Series period 1949-2003 1977-2003 1949-2003 1977-2003 ξ −0.377 −0.377 −0.357 −0.244 This method provides 10-, 20-and 30-year return levels, which seem physically reasonable.Indeed, the 10-or 20-year return levels were reached and even exceeded in particular during the 2003 heat wave on all sites.For Chooz even the 30-year return level was approached.Extrapolation to longer return periods seems hazardous compared with the sample size of temporal series of water temperatures.
In addition, return levels are subject to an important sample size effect.Indeed, all the observed series present generally the same behaviour: they start from one period of several years with rather cold summers to finish with the 2003 heat wave.Moreover, sensitivity studies show that the initial cold period has more influence in trend calculations than the final period including 2003.This behaviour, due to natural variability, is not entirely attributable to anthropogenic climate change and leads to stronger trends.This is why, it is necessary to confront the results with those of a physical and deterministic modelling of the rivers thermal behaviour.The physical model (e.g.CALNAT) (see Gras, 1969 andDupeyrat et al., 2006) could be used with meteorological parameters (air temperature, air humidity, wind, . . . ) scenarios taking the climatic change into account.
Finally, these results should be moderated for two reasons: -They suppose that the thirty last years trend will continue in the thirty next years.This working hypothesis cannot be checked.
-In addition, they also suppose that current forcing by the atmospheric conditions will have the same effect on water temperature in thirty years.However a high air temperature modifies the steam pressure of water surface and then enhances its evaporation, a phenomenon that cools water surface (Mohseni et al., 1999).The relationship between air temperature and water temperature is complex, and could be modified in case of very high air temperatures.
Nevertheless, the linear extrapolation of water temperature trend remains the most reasonable assumption in the short term to estimate the return levels of hot extreme temperatures.

Fig. 2 .
Fig. 2. Power stations site and related water temperature series.

4. 2
Figures 8 and 9 give the 10-, 20-and 30-year return levels while taking (right figure) or not (left figure) series clusterization into account.The 10-, 20-and 30-year return levels seem to be physically reasonable.Indeed, the 10-or 20-year return levels were already reached, sometimes exceeded in particular during the 2003 heat wave, on all sites.For Chooz, even the 30-year return level was approached.On the other hand, the 50-and 100-year return levels do not appear very reasonable in particular in comparison with the sample size (27 years, sometimes less when using piecewise linear functions) on which the extrapolation is based.To evaluate the trend effects, return levels were computed with a stationary assumption.In this case, µ, σ and ξ parameters of GEV distribution are supposed to be constant.The results are given in Fig.9in the same form as in the non stationary case (see Fig.8).We also compute confidence intervals for the different return levels.The non stationarity hypothesis makes it impossible to use close asymptotics.Thus, we simulate 1000 trajectories of the estimated optimal model.Then, the return levels are computed for every trajectory.From these data, we obtain an empirical estimate of the return levels distribution and we choose as confidence interval an interquantile interval, for instance (q(0.05),q(0.95))for the risk level 0.1.We also check that the magnitude of these intervals does not depend on small variations of the estimated model parameters.Table7shows different confidence intervals computed using trajectories, for Saint Laurent.

Fig. 8 .
Fig. 8.Return levels (10, 20 and 30 years) by sites with non stationary hypothesis (without θ on left and with θ on right).

Fig. 9 .
Fig. 9.Return levels (10, 20 and 30 years) with stationary hypothesis for all the sites (without θ on left and with θ on right).

Fig. 10 .
Fig. 10.Difference of 30-year return level with or without θ in non stationary and stationary cases, for all sites.

Table 1 .
Percentage and longest period of missing data.

Table 2 .
Characteristic size of clusters during heat wave period (the number of days written in bold corresponds to the maximum of the clusters size over the 1977-2003 period).

Table 3 .
Optimal trend computed for GEV parameters (Loire sites).

Table 4 .
Optimal trend computed for GEV parameters (Rhône sites).Linear on µ Linear on µ Linear on µ

Table 5 .
Optimal trend computed for GEV parameters (others sites).

Table 6 .
Last slope value of µ parameter (in C degrees/time).

Table 7 .
Cconfidence intervals for return levels (in • C degrees) at Saint-Laurent.

Table 8 .
Return levels computed from 1949-2003 series.Return levels Without θ With θ Without θ With θ Without θ With θ Without θ With θ

Table 9 .
ξ value for the two periods