Accepted Manuscript Exploiting the spatial pattern of daily precipitation in the analog method for regional temporal disaggregation Julie Carreau, Nada Ben Mhenni, Frédéric Huard, Luc Neppel PII: DOI: Reference:
S0022-1694(18)30878-3 https://doi.org/10.1016/j.jhydrol.2018.11.023 HYDROL 23264
To appear in:
Journal of Hydrology
Please cite this article as: Carreau, J., Mhenni, N.B., Huard, F., Neppel, L., Exploiting the spatial pattern of daily precipitation in the analog method for regional temporal disaggregation, Journal of Hydrology (2018), doi: https:// doi.org/10.1016/j.jhydrol.2018.11.023
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
1
2
Exploiting the spatial pattern of daily precipitation in the analog method for regional temporal disaggregation Julie Carreaua,∗, Nada Ben Mhennib , Frédéric Huardc , Luc Neppela
3
4
a HydroSciences
5 6 7
8
Montpellier (CNRS/IRD/UM), Université de Montpellier - Case 17, 163 rue Auguste Broussonet 34090 Montpellier, France b Nagoya University, Nagoya, Japan c AgroCLIM, INRA, Avignon, France
Abstract
9
The analog method relies on the subdaily distribution of analog days - also called fragments - to perform
10
temporal disaggregation. The features that serve to define analog days are generally at the daily scale
11
although they are assumed to convey information on the subdaily scale. In most cases, the only feature used
12
is the daily total at the location of interest. Features pertaining to the synoptic circulation are also often
13
used although with the purpose of increasing the spatial rather than the temporal resolution. Precipitation is
14
conditioned not only on synoptic circulation but also on regional features such as orography. This combined
15
effect is reflected on the spatial pattern of daily precipitation. We propose to rely on descriptions of the spatial
16
pattern to select analog days in order to perform temporal disaggregation. Three analog method variants
17
based on different sets of spatial pattern features, a variant that rely on synoptic circulation alone and a
18
variant that combines synoptic circulation with spatial pattern features are evaluated and compared. All the
19
features considered characterize a region therefore the same analog day is employed at all the locations within
20
the region. Inter-site coherence is thus enforced. A leave-one-out procedure is carried out on hourly radar
21
precipitation reanalyses covering a mountainous watershed in the French Mediterranean on a 1 km2 grid
22
over the period 1997-2006. Uncertainty is taken into account by generating five realizations of each variant.
23
Performance is measured in terms of return levels (1 to 20 years) at various temporal scales (1 h to 12 h) and
24
at two spatial scales - grid box and watershed. The proposed regional approach to temporal disaggregation
25
globally preserves inter-site coherence. This can be seen from the relatively good reproduction of the regional
26
pattern of return levels of the study area by all variants. Synoptic circulation conveys less information on the
27
subdaily scale than spatial pattern features. The variant that depends on synoptic circulation alone largely
28
overestimates return levels especially at lower temporal scales (1 h and 3 h). Combining synoptic circulation
29
with spatial pattern features yielded no significant improvement.
30
Keywords: analog method, temporal disaggregation, spatial pattern of daily precipitation, weather types,
31
spatial and temporal coherence, French Mediterranean
∗ Corresponding
author Email addresses:
[email protected] (Julie Carreau),
[email protected] (Nada Ben Mhenni), Preprint submitted to Journal of Hydrology November 2, 2018
[email protected] (Frédéric Huard),
[email protected] (Luc Neppel)
1
1. Introduction
2
In many study areas, daily precipitation data are available from relatively dense networks with compara-
3
tively long observation periods. In contrast, subdaily data is often scarce owing to the recent installation of
4
rain gauges with subdaily resolution (there is about 10 years of hourly data in the Cevennes-Vivarais area in
5
the French Mediterranean, see http://ohmcv.osug.fr/). However, subdaily precipitation is the main input
6
to several impact models such as rainfall-runoff models. These models are designed to reproduce flash floods,
7
a typical phenomenon of the Mediterranean region that occurs within a few hours or less [7, 4]. To augment
8
subdaily data bases, temporal disaggregation methods can be applied to daily data bases in order to increase
9
the temporal resolution.
10
Temporal disaggregation methods can be divided into three main categories. The first category encom-
11
passes random cascade models and approaches based on scale invariance properties [26, 13]. Stochastic
12
weather generators constitute the second category [17, 1, 12]. Analog methods (sometimes called methods
13
of fragments) form the third category [18, 32, 16]. In this work, we focus on the last category. Despite
14
being conceptually simple, analog approaches have been proven to be competitive [23, 2]. These methods are
15
non-parametric and rely on the core concept of analog days.
16
To perform temporal disaggregation, the subdaily distribution of analog days, also known as fragments,
17
is employed to distribute daily precipitation totals into subdaily intervals. Analog days are defined as days
18
that are similar in terms of a number of features. They are selected for each day of the daily data base
19
among the days of the subdaily data base. Often the only feature used to determine analog days is the daily
20
precipitation total at a given location [16, 2]. In some cases, this is supplemented with the information on the
21
wetness states (whether it rained or not) of the days preceding and following the day of interest [32, 23]. The
22
features employed to select analog days are generally at the daily scale. In spite of this, they are implicitly
23
assumed to convey information on the subdaily scale.
24
In operational forecasts or in climate change impact studies, analog approaches rely on the premise that
25
precipitation depends on synoptic circulation to increase the spatial resolution and correct biases, a process
26
commonly referred to as statistical downscaling [21, 3, 24, 33]. Similarity between days is based on the
27
synoptic circulation predicted or simulated by numerical weather prediction models (NWPMs) or general
28
circulation models (GCMs). NWPMs and GCMs typically operate at a low spatial resolution with respect
29
to the needs of impact studies. The increase in spatial resolution is achieved by selecting analog days from
30
a precipitation data base at high resolution. An additional benefit of the approach is to bypass the bias
31
in the precipitation fields predicted or simulated by NWPMs and GCMs. As stated in Obled et al. [21],
32
although precipitation is conditioned in part on synoptic circulation, it also depends on regional features such
2
1
as orography.
2
The spatial pattern of precipitation in the area of interest (where and how much precipitation occurred)
3
can be thought of as translating the combined effect of regional features and synoptic circulation. To our
4
knowledge, there are very few attempts at including features of the spatial pattern of precipitation in statistical
5
downscaling or in temporal disaggregation approaches. One such attempt can be found in Mezghani &
6
Hingray [18]. A single feature, the daily spatial average, is computed on both temperature and precipitation.
7
First, the daily spatial averages are estimated from the synoptic circulation - a statistical downscaling step.
8
They are then employed to define analog days to obtain subdaily distributions - a temporal disaggregation
9
step. We propose to describe further the spatial pattern of daily precipitation and to consider several of its
10
features.
11
In this study, our aim is twofold. First, we seek to assess the ability of features of the spatial pattern
12
of daily precipitation at conveying information on the subdaily distribution of precipitation. This ability
13
is evaluated through the performance of temporal disaggregation with the analog method. Second, we
14
seek to evaluate whether synoptic circulation, although conventionally used in statistical downscaling, can
15
bring useful information in temporal disaggregation. Since spatial pattern features and synoptic circulation
16
characterize a region, temporal disaggregation is performed simultaneously for all the sites in the region of
17
interest. More precisely, the same analog day is used at all the sites thereby enforcing inter-site coherence.
18
Our analyses are based on the evaluation and comparison of five analog method variants. We set a variant
19
based on synoptic circulation as summarized by weather types [10]. We define three variants that rely on
20
different sets of spatial pattern features. The last variant combines synoptic circulation and spatial pattern
21
features. The analyses are carried out on radar precipitation reanalyses [28]. The selected 1 km2 grid covers
22
a mountainous watershed in the French Mediterranean prone to flash floods. This is meant as a perfect case
23
study in which the definition of features of the spatial pattern of daily precipitation is not affected by the
24
density of the rain gauge network. The analog method variants are evaluated with a leave-one-out procedure.
25
To assess the temporal and spatial coherence of the variants, their performance is measured in terms of their
26
ability to reproduce return levels (1 to 20 years) at various temporal scales (1 h to 12 h) and at two spatial
27
scales : (1) the grid box scale and (2) at the watershed scale.
28
2. Study area and data
29
The study area is located in the Cevennes range (highest peak about 1700 m), a mountainous area of the
30
French Mediterranean that sits at the southern end of the Massif Central (Fig. 1a). The study area has a
31
rectangular shape of about 60 km by 30 km (1800 km2 ) and encompasses the Gardon at Anduze watershed
32
whose area is approximately 545 km2 (Fig. 1b). The watershed features a steep orography with narrow 3
1
valleys separated by ridges. Two rivers, the Gardon at Saint-Jean and Gardon at Mialet flows from north-
2
west to south-east. The study area is prone to rare but catastrophic flash floods like those which occurred
3
on September 2002 [7].
(b)
(a)
(c) Figure 1: (a) Topographic map of France with the location of the study area (dashed white rectangle). (b) Zoom on the study area, with extended Lambert II coordinates (km), encompassing the Gardon at Anduze watershed (dotted black line). (c) Two sub-regions defined by clustering the grid boxes with K-means according to their annual totals - in white (resp. in grey), the mean annual total is 1427 mm (resp. 1057 mm).
4
Météo-France, the French weather service, provided us with a 10 year (1997-2006) radar reanalysis of
5
hourly precipitation data over a 1 km2 grid [28]. The reanalysis is obtained by merging 5 min radar reflectivity
6
images with hourly and daily raingauge data. The grid covering the study area in Fig. 1b contains 2016 boxes,
7
63 in longitude and 32 in latitude. We computed daily precipitation totals by aggregating hourly precipitation
8
in the 24 time intervals preceding 7:00 am.
9
Two subregions, that correspond approximately to the mountainous area and the valley, are obtained by
10
clustering the grid boxes with K-means according to their annual totals for each of the 10 year (1997-2006)
11
period, see Fig. 1c. As expected, the mean annual total in the mountainous area (in white in Fig. 1c) is
12
higher than in the valley (in grey in Fig. 1c), 1427 mm in the former as oppposed to 1057 mm in the latter.
13
This partition into two subregions is used in § 3.1 to define an analog method variant.
14
We also make use of the daily classification into eight weather types of Garavaglia et al. [10]. The
15
classification is obtained by first identifying typical spatial patterns of daily precipitation over France. These
16
patterns are then associated to synoptic circulation as described by geopotential height fields measured at
17
pressure levels 700 and 1000 hPa and at times 0 h and 24 h, see Garavaglia et al. [10] for all the details. 4
Figure 2: Hietograph of the daily maximum precipitation over the study area (see Fig. 1) for the 10 year (19972006) period; annual maxima are depicted with diamond shapes, the highest one being the maximum of the year 2002 (665.4 mm/day) of the catastrophic flash flood event analysed in Delrieu et al. [7]. The colors correspond to the weather type classification proposed in Garavaglia et al. [10].
1
The hietograph of daily maximum precipitation over the study area for the 10 year period is shown in
2
Fig. 2. We can observe that intense precipitation events, that can be defined as daily totals greater than
3
100 mm as in Tramblay et al. [29], tend to occur mostly in autumn (from September to December). In
4
Fig. 2, each color corresponds to a weather type. While the 4th weather type called south circulation has
5
been connected with the occurrence of heavy rainfall [30], we see from Fig. 2 that the 7th weather type called
6
central depression is also associated with heavy rainfall, in particular, with the event in September 2002
7
which yielded catastrophic flash floods [7].
8
For the subsequent analyses, we selected 1113 days considered as rainy. The following criteria were
9
applied : at least one grid box has at least 4 mm/day and at least 100 grid boxes have at least 1 mm/day.
10
Besides, as the reanalyses contain a lot of small values, we use a threshold of 1 mm to define positive
11
precipitation totals.
12
3. Analog method for regional temporal disaggregation
13
Let us call target day a day d∗ for which only the daily precipitation total Rd∗ is assumed to be known.
14
Moreover, let us call analog day a day d which is similar to the target day d∗ with respect to some features.
15
Similarly, we denote Rd and Rd,h with 1 ≤ h ≤ 24, the daily and hourly totals respectively of the analog day.
16
The hourly profile Pd (also known as fragments) of day d describes its subdaily distribution. It represents
17
the proportion of daily precipitation that fell into each hour and is computed as follows
5
Pd =
Rd,24 Rd,1 ,..., Rd Rd
.
(1)
1
Temporal disaggregation of the target day d∗ is achieved by applying the hourly profile of day d to the daily
2
total of day d∗ yielding
ˆ d∗ ,1 , . . . , R ˆ d∗ ,24 ) = Pd × Rd∗ . (R
(2)
3
The observed daily total of the target day is automatically preserved and is distributed into time intervals of
4
1 h following the hourly profile ot the analog day.
5
When performing temporal disaggregation at several locations simultaneously, the same analog day serves
6
for all the locations thereby ensuring spatial coherence (see Domínguez-Mora et al. [9] for a similar way to
7
extract a spatial shape in a resampling method). If it occurs that Rd = 0 for some locations, then the
8
corresponding hourly profiles are set to Pd = (1/24, . . . , 1/24).
9
3.1. Selection of analog days
10
We define five analog method variants that differ in the way analog days are selected. One variant selects
11
randomly an analog day that shares the same weather type and occurred during the same month as the target
12
day. When no days fulfilled these two requirements, the selection is made among days with the same weather
13
type only. Fig. 3 illustrates the mean daily precipitation intensities per weather type for the reanalysis data
14
over the study area. The spatial pattern changes, although moderately so, with the weather type.
15
Three variants are based on different sets of features of the spatial pattern of daily precipitation. In these
16
variants, an analog day is selected so as to maximize the similarity with the target day in terms of the spatial
17
pattern features. To achieve this, standardization is performed to bring all the features to a comparable
18
range of values by computing the standard Normal quantiles associated to the empirical frequencies of the
19
features. More precisely, let x be the value taken by a given feature, let Fˆn (·) be the associated empirical
20
cumulative distribution function and Φ(·) be the standard Normal cumulative distribution function. Then
21
the standardized feature value is given by Φ−1 (Fˆn (x)). Analog days are selected as days whose standardized
22
features minimize the Euclidean distance with the standardized features of the target day.
23
The first variant relies on the following three features : (1) the spatial average, (2) the spatial coefficient
24
of variation and (3) the spatial proportion of daily precipitation intensities [15]. Let Rd (s) be the daily total
25
at location s on day d. These three features are computed from Sd = {1 ≤ s ≤ 2016 : Rd (s) > 1}, the set of
26
locations with positive daily precipitation totals on day d :
6
(a) atlantic wave (16.1%)
(b) steady oceanic (19.5%)
(c) south-west circulation (7.55%)
(d) south circulation (18.4%)
(e) north-east flux (9.97%)
(f ) east return (8.45%)
(g) central depression (7.37%) (h) anticyclonic circulation (12.7%) Figure 3: Average daily precipitation intensities of the reanalysis data for each weather type over the 1113 days analyzed in this work together with the percentage of occurrence of the weather type indicated in parentheses.
7
(1)
Spatial average
x ¯d =
(2)
Spatial coefficient of variation
(3)
Spatial proportion
q
1 |Sd | P
1/(|S |−1) d
P
s∈Sd
Rd (s),
xd ) s∈Sd (Rd (s)−¯
2
/x¯d ,
|Sd |/2016.
1
Fig. 4 illustrates the aforementioned three spatial features for each weather type. There are some, albeit
2
mild, changes in the pattern made by these features with the weather type, e.g. between central depression
3
in Fig. 4g and anticyclonic circulation in Fig. 4h.
4
The second and third variants seek to include information on the location of the rain shower. The
5
second variant contains, in addition to the three features of the first variant, the x- and y-coordinates of the
6
barycenter of the rain shower. These are computed as the average of the x- and y-coordinates weighted by
7
the precipitation intensities considering only values above the spatial third quartile of the daily precipitation
8
totals. More precisely, let S˜d be the set of grid boxes whose daily totals are among the 25% largest on day d.
9
Moreover, for s ∈ S˜d , let x(s) and y(s) be the coordinates of the grid box and let wd (s) = Rd (s)/
10
P
˜ s∈S d
Rd (s).
The coordinates of the barycenter on day d are then defined as
X
˜d s∈S
X
wd (s)x(s),
˜d s∈S
wd (s)y(s) .
11
An illustration of a barycenter with the corresponding rain shower shown as the map of its daily totals
12
is given in Fig. 5 for September 9th, 2002 - the catastrophic flash flood event analysed in Delrieu et al. [7].
13
In Fig. 6, the barycenters are presented for each weather type. There is a tendency for the barycenters to be
14
located over the mountain range, particularly for the steady oceanic and south circulation weather types.
15
In the third variant, information on the location of the rain shower is included by computing the three
16
features of the first variant - the spatial average, the spatial coefficient of variation and the spatial proportion
17
- on the two subregions of the study area shown in Fig. 1c. The representation of the three features restricted
18
on each subregion is very similar to the representation in Fig. 4 for the whole study area and is therefore not
19
shown. Variants 1, 2 and 3 have respectively 3, 5 and 6 features.
20
As weather types seem to bring a complementary information to the three features of the first variant,
21
see Fig. 4, we consider a fourth variant which combines the first variant and the weather type classification.
22
An analog day is selected among the days sharing the same weather type as the target day by maximizing
23
the similarity in terms of the three features of the first variant.
8
(a) atlantic wave
(b) steady oceanic
(c) south-west circulation
(d) south circulation
(e) north-east flux
(f ) east return
(g) central depression (h) anticyclonic circulation Figure 4: For each weather type defined in Garavaglia et al. [10], the spatial average (x-axis in logarithmic scale), the spatial coefficient of variation (y-axis) and the spatial proportion (size of plotting symbol) of daily precipitation intensities for the selected 1113 days are represented.
9
Figure 5: Barycenter (white dot) of the rain shower of September 9th, 2002 computed as the average x- and y-coordinates weighted by the precipitation totals for the grid boxes whose daily totals are among the 25% largest for that day. Table 1: Five variants to select analog days : variants 1 to 3 are described in terms of features of the spatial pattern of daily precipitation which are used to define the similarity with target days, see § 3.1.
WT Variant Variant Variant Variant
1
1 2 3 4
random day from the same month with the same weather type spatial average, coefficient of variation and proportion of daily precipitation intensities variant 1 + barycenter of rain shower variant 1 on two sub-regions variant 1 among days with the same weather type
4. Return level estimation
2
To evaluate and compare the five analog method variants described in Table 1, we employ performance
3
criteria based on return level estimates. For a given temporal scale, say n-hourly, a return level of period
4
T years is the n-hourly total level which is expected to be exceeded on average once every T years. It
5
corresponds to a quantile of the distribution of n-hourly totals.
6
To estimate return levels, we make use of the peaks-over-threshold approach and rely on the generalized
7
Pareto (GP) distribution. Indeed, for almost all distributions, the GP can appproximate the distribution of
8
the excesses over a large threshold arbitrarily well [22]. Let Y be the random variable that represents the
9
excesses of n-hourly totals over u, a sufficiently high threshold. Then Y is approximately GP distributed and
10
return levels of period T years for n-hourly totals can be estimated as follows u + σ ((T Nexc )ξ − 1) ξ u + σ log(T Nexc )
if ξ 6= 0
(3)
if ξ = 0
11
where σ > 0 and ξ ∈
are the scale and shape parameters respectively of the GP and Nexc corresponds to
12
the average number of excesses per year. This can be computed as the product of the probability of exceeding
13
the threshold u and 365.25 × 24/n, the number of observations per year.
14
4.1. Regional approach
15
As the radar reanalyses only extend over a 10 year period, we employ a regional approach to increase
16
the sample size at each grid box. This is particularly useful to reduce the uncertainty linked to the shape 10
(a) atlantic wave
(b) steady oceanic
(c) south-west circulation
(d) south circulation
(e) north-east flux
(f ) east return
(g) central depression (h) anticyclonic circulation Figure 6: For each weather type defined in Garavaglia et al. [10], barycenters of rain showers (average of the x- and y-coordinates weighted by the precipitation totals for grid boxes whose daily totals are among the 25% largest) for the selected 1113 days.
11
1
parameter estimate [14]. We follow the regional approach described in Carreau et al. [5]. The main assump-
2
tion is that the region of interest can be partitioned into subregions with a constant shape parameter. The
3
approach is summarized below.
4
In our case, a subregion is a set of M grid boxes. Let us assume that the shape parameter of the GP is
5
constant over these M grid boxes, i.e. ξi = ξ R ∀i, with ξ R a regional shape parameter. Let µi be the mean
6
of the excesses Yi over a threshold ui . Then the distribution of the scaled excesses
7
1 ≤ i ≤ M and only depends on the regional shape parameter [5].
Yi/µi
is identical for all
8
The inference strategy of the regional approach relies, in addition to µi , on the second probability weighted
9
moment (PWM) of the scaled excesses [8]. First, µi is approximated with kernel regression [20, 31] and
10
the scaled excesses are computed. Second, a partition into subregions is performed automatically with K-
11
means [25] based on the second PWM of the scaled excesses. Finally, ξ R is estimated for each subregion from
12
the second PWM of the scaled excesses and σi is determined by combining the regional shape estimates with
13
µi estimates.
14
For this case study, the threshold is set to the 95% quantile of precipitation intensities, i.e. daily totals
15
greater than 1 mm. See Ceresetti et al. [6] for instance for a justification of this choice of threshold in a
16
study area that encompasses the present one. Moreover, we fix the partition to contain three subregions in
17
order to capture the upper tail behaviour in the study area.
18
5. Evaluation and comparison of analog method variants
19
To evaluate and compare the five analog method variants (see Table 1), we performed a leave-one-out
20
analysis (also known as jackknife) over the 1113 wet days of the 10 year observation period. We applied an
21
ensemble approach and selected five analog days for each target day : the five days that are most similar in
22
terms of spatial pattern features for variants 1 to 4 and five random days for the weather type variant. For
23
each variant, we thus obtained an ensemble of five realizations of temporally disaggregated precipitation.
24
5.1. Return levels at the grid box scale
25
Temporal coherence is assessed by comparing return levels at several temporal scales. From the hourly
26
precipitation of the reanalysis and the disaggregated data from each analog variant, we estimated return
27
levels for several periods (1, 2, 5, 10 and 20 years) at four temporal scales (1 h, 3 h, 6 h and 12 h) at each of
28
the grid boxes. For each temporal scale, we employed the regional peaks-over-threshold approach explained
29
in § 4.
30
The maps of the return levels of period 1 year and 20 years for two temporal scales, 1 h and 12 h, are
31
provided in Figs. 7-10. For these maps, we picked a single realization of the ensemble. For variants 1 to 4, it 12
1
corresponds to the realization provided by the most similar analog. For longer return periods that correspond
2
to more intense events and/or short temporal scales (see Figs. 7a, 9a and 10a), precipitation intensities are
3
higher at the foothills than on the mountain range. This is typical of this area as the orography tends to have
4
a blocking effect [7] (see also Gardes & Girard [11] who obtained a map of hourly return levels of period 10
5
years estimated from rain gauge data with a similar pattern). This pattern is relatively well captured by the
6
analog method variants. However, the return levels provided by the WT variant are largely overestimated at
7
the temporal scale 1 h.
(a) Reanalysis
(b) Variant 1
(c) Variant 2
(d) Variant 3
(e) Variant 4 (f ) Variant WT Figure 7: Return level maps of period 1 year for a temporal scale of 1 h estimated with a regional peaks-over-threshold approach with a partition into three subregions [5]. A single realization is presented for each analog method variant, it corresponds to the most similar analog for variants 1 to 4. See Table 1 for a description of the five variants.
8
The root-mean-square error (RMSE) can serve to evaluate quantitatively the ability of each analog method
9
variant at reproducing the return levels estimated from the reanalysis data at each of the grid boxes. We use
10
a normalized RMSE so that the error for each temporal scale is comparable. It is defined as the ratio of the
13
(a) Reanalysis
(b) Variant 1
(c) Variant 2
(d) Variant 3
(e) Variant 4 (f ) Variant WT Figure 8: Return level maps of period 1 year for a temporal scale of 12 h estimated with a regional peaks-over-threshold approach with a partition into three subregions [5]. A single realization is presented for each analog method variant, it corresponds to the most similar analog for variants 1 to 4. See Table 1 for a description of the five variants.
14
(a) Reanalysis
(b) Variant 1
(c) Variant 2
(d) Variant 3
(e) Variant 4 (f ) Variant WT Figure 9: Return level maps of period 20 years for a temporal scale of 1 h estimated with a regional peaks-over-threshold approach with a partition into three subregions [5]. A single realization is presented for each analog method variant, it corresponds to the most similar analog for variants 1 to 4. See Table 1 for a description of the five variants.
15
(a) Reanalysis
(b) Variant 1
(c) Variant 2
(d) Variant 3
(e) Variant 4 (f ) Variant WT Figure 10: Return level maps of period 20 years for a temporal scale of 12 h estimated with a regional peaks-over-threshold approach with a partition into three subregions [5]. A single realization is presented for each analog method variant, it corresponds to the most similar analog for variants 1 to 4. See Table 1 for a description of the five variants.
16
1
RMSE and of the standard deviation of the return levels estimated from the reanalysis data :
q
1 2016
P2016 s=1
(zsdisag −zsrea )2
/σrea .
(4)
2
In Eq. (4) above, zsrea and zsdisag are respectively the return levels estimated from the reanalysis data and
3
from an analog method at each grid box 1 ≤ s ≤ 2016 and σ rea is the spatial standard deviation of the return
4
levels estimated from the reanalysis data for a given temporal scale.
5
In Fig. 11, each analog variant is shown as a band whose lower (resp. upper) value is the minimum (resp.
6
maximum) of the normalized RMSE of the ensemble of five realizations. The WT variant has overall the
7
worst performance, especially for lower return periods - 5 years or less - and for the shorter time scales - 6
8
hours or less. Variants 2 and 3 that include information on the location of the rain shower do not perform
9
systematically better than variant 1 that uses only three features. Similarly, variant 4 which combines
10
variant 1 with weather types did not yield a performance significantly different from variant 1. For all
11
variants, the error decreases in magnitude and in spread with increasing temporal scale. This is expected
12
since, by construction, the daily totals are preserved.
13
5.2. Return levels at the watershed scale
14
Although the return level maps provide information on the spatial coherence of the analog method variants,
15
further assessment is made by considering areal rainfall (see Müller & Haberlandt [19] for other criteria
16
pertaining to spatial coherence). We estimated return levels of several return periods - 1, 2, 5, 10, and
17
20 years - at several temporal scales - 1 h, 3 h, 6 h and 12 h - for the areal rainfall over the Gardon at
18
Anduze watershed. To this end, the mean of hourly precipitation totals is computed for the the grid boxes
19
belonging to the watershed to obtain hourly areal totals. For each temporal scale, hourly areal precipitation
20
is aggregated with a sliding window. To remove temporal dependence, a declustering procedure is applied
21
(see the function declust in the R package from Southworth et al. [27] and the references therein). We then
22
estimated return levels for each temporal scale thanks to the GP approximation, see § 4. For the analog
23
methods, the estimation is carried out for each member of the ensemble of five realizations.
24
In Fig. 12, return levels of rainfall intensities are shown per time scale - 1, 3, 6 and 12 hours. The
25
reanalyses’ levels are represented by a black curve while the analog variants’ levels are shown as coloured
26
bands - the minimum and the maximum of the ensemble of five realizations. The bands from the WT variant
27
overestimates the reanalyses’ curve at all temporal scales for all return periods albeit they get closer at time
28
scale 12 h. In contrast, the bands from variants 1 to 4 include most of the times the reanalyses’ curve.
29
Variants 2 and 3, which contain information on the location of the rain shower and have more features, tend
30
to have larger bands at time scales 1 h and 3 h. Variant 4 does not yield a significant improvement over 17
(a) 1 h
(b) 3 h
(c) 6 h (d) 12 h Figure 11: Bands of normalized RMSE (min and max of ensemble of five disaggregated realizations) for return levels of periods ranging from 1 to 20 years (in logarithmic scale) at four temporal scales for the five analog variants (see Table 1 for the definition of the variants). Lower values mean better performance.
18
1
variant 1 in this respect. Overall, the analog variants’ bands are closer to the reanalyses’ curve and are
2
narrower as the temporal scale increases. As was mentioned for the normalized RMSE in Fig. 11, this is
3
expected since daily totals are preserved.
(a) 1 h
(b) 3 h
(c) 6 h (d) 12 h Figure 12: Areal return level intensity curves for four time scales (in logarithmic scale) with respect to the return period in years (in logarithmic scale). The reanalyses’ levels are shown as a black curve while the analog variants’ are shown as coloured bands - min and max of ensemble of five disaggregated realizations (see Table 1 for the definition of the variants).
4
6. Summary and outlook
5
We evaluated and compared five analog method variants to temporally disaggregate daily precipitation
6
totals simultaneously at several sites. The variants differ in the sets of features used to select analog days. We
7
considered various sets of features of the spatial pattern of daily precipitation together with information on
8
the synoptic circulation in the form of weather types (see Table 1). Our case study relies on radar reanalyses
9
of precipitation which is meant as a perfect case study. We selected a 63 by 32 grid with a spatial resolution
10
of 1 km2 . The grid covers the Gardon at Anduze watershed located in the Cevennes mountain range in the 19
1
French Mediterranean. We carried out a leave-one-out procedure taking in turn as target day each of 1113
2
days defined as rainy. Performance assessment is based on return levels at various temporal scales and at two
3
spatial scales - grid box and watershed - thereby assessing the temporal and spatial coherence of the analog
4
method variants.
5
The conclusions drawn from our case study are twofold. First, the ability of the spatial pattern features at
6
conveying information on the subdaily distribution of precipitation is assessed. This can be seen through the
7
good performance of the analog method variants 1, 2, and 3. Second, although synoptic circulation is often
8
used in statistical downscaling, we found that it brings only limited information in temporal disaggregation.
9
More specifically, our analyses based on Figs. 7 - 12 lead to the following findings.
10
• Information on the location of rain showers, as defined in this work, does not bring significant additional
11
useful information on the subdaily distribution of precipitation. This is deduced from the comparison
12
of variants 2 and 3 with variant 1.
13
• Synoptic circulation alone as summarized by weather types is not sufficient to recover accurate subdaily
14
precipitation distribution. This is particularly striking at lower temporal scales (1 h and 3 h). This
15
conclusion is drawn from the poor performance of variant WT.
16
• The combination of spatial pattern features with synoptic circulation - variant 4 - did not yield signif-
17
icant improvement over spatial pattern features alone. This can be seen by comparing variants 1 and
18
4.
19
In this work, we exploited the spatial pattern of daily precipitation to select analog days. As it is
20
well known, synoptic circulation influences precipitation and is typically employed in statistical downscaling
21
approaches [21, 3, 24, 33]. Precipitation is also conditioned by regional features such as orography [21].
22
The spatial pattern of daily precipitation can be thought of as a way to capture the interaction of synoptic
23
circulation and regional features.
24
In many studies, the selection of analog days is based on the daily total precipitation value at a given
25
location on the target day [16, 2]. To our knowledge, only Mezghani & Hingray [18] made use of the spatial
26
pattern of daily precipitation in the framework of temporal disaggregation with the analog method. In their
27
work, a single feature, the spatial average, was considered. In this work, we built on the three features of the
28
spatial pattern defined in Leblois [15] in order to identify precipitation regimes. We proposed two variants
29
that include information on the location of the rain shower.
30
As opposed to daily totals which are specific to each location, spatial features contain information on the
31
whole area. Therefore, temporal disaggregation can be performed simultaneously and straightforwardly at
20
1
several sites by selecting a unique analog day for all the sites. On the other hand, the computation of the
2
spatial pattern features relies on the availability of daily totals at a significant number of sites in the area of
3
interest.
4
Potential next steps to evaluate the performance of spatial pattern features in the analog method for
5
temporal disaggregation are as follows. It could be applied on real raingauged networks or by sub-sampling
6
the radar reanalysis grid to assess the sensitivity of the spatial pattern features to the density of gauged sites.
7
It would also be interesting to apply the proposed analog approach in other areas in which precipitation is
8
driven by different types of regional features. Besides, inter-day coherence was not assessed and it is probably
9
under-estimated by the proposed analog variants. Following Westra et al. [32], Pui et al. [23], the wetness
10
states of the days preceding and following the target day could be considered when selecting analog days.
11
Acknowledgements
12
This work was partially supported by the Agence Nationale de la Recherche (French Research Agency)
13
with the projects Almira and Starmip. Part of this has also been supported by the LEFE-INSU-Cerise
14
project. We are grateful to Météo-France for providing the radar precipitation reanalyses. We thank our two
15
reviewers, H. Müller-Thomy and J. Bliefernicht, for their valuable comments, which greatly helped improve
16
the quality of the paper.
21
1
[1] Allard, D., & Bourotte, M. 2015. Disaggregating daily precipitations into hourly values with a transformed
2
censored latent Gaussian process. Stochastic Environmental Research and Risk Assessment, 29(2), 453–
3
462.
4
5
6
7
[2] Bárdossy, A., & Pegram, G. 2016. Space-time conditional disaggregation of precipitation at high resolution via simulation. Water Resources Research, 52(2), 920–937. [3] Bliefernicht, J., & Bárdossy, A. 2007. Probabilistic forecast of daily areal precipitation focusing on extreme events. Natural Hazards and Earth System Science, 7(2), 263–269.
8
[4] Braud, I., Ayral, P.-A., Bouvier, C., Branger, F., Delrieu, G., Le Coz, J., Nord, G., Vandervaere, J.-P.,
9
Anquetin, S., Adamovic, M., Andrieu, J., Batiot, C., Boudevillain, B., Brunet, P., Carreau, J., Confoland,
10
A., Didon-Lescot, J.-F., Domergue, J.-M., Douvinet, J., Dramais, G., Freydier, R., Gérard, S., Huza,
11
J., Leblois, E., Le Bourgeois, O., Le Boursicaud, R., Marchand, P., Martin, P., Nottale, L., Patris, N.,
12
Renard, B., Seidel, J.-L., Taupin, J.-D., Vannier, O., Vincendon, B., & Wijbrans, A. 2014. Multi-scale
13
hydrometeorological observation and modelling for flash flood understanding. Hydrology and Earth System
14
Sciences, 18(9), 3733–3761.
15
16
[5] Carreau, J., Naveau, P., & Neppel, L. 2017. Partitioning into hazard subregions for regional peaks-overthreshold modeling of heavy precipitation. Water Resources Research, 53(5), 4407–4426.
17
[6] Ceresetti, D., Ursu, E., Carreau, J., Anquetin, S., Creutin, J.-D., Gardes, L., Girard, S., & Molinie, G.
18
2012. Evaluation of classical spatial-analysis schemes of extreme rainfall. Natural hazards and earth system
19
sciences, 12, 3229–3240.
20
[7] Delrieu, G., Nicol, J., Yates, E., Kirstetter, P.-E., Creutin, J.-D., Anquetin, S., Obled, C., Saulnier,
21
G.-M., Ducrocq, V., Gaume, E., Payrastre, O., Andrieu, H., Ayral, P.-A., Bouvier, C., Neppel, L., Livet,
22
M., Lang, M., du Châtelet, J. P., Walpersdorf, A., & Wobrock, W. 2005. The Catastrophic Flash-Flood
23
Event of 8-9 September 2002 in the Gard Region, France: A First Case Study for the Cévennes-Vivarais
24
Mediterranean Hydrometeorological Observatory. Journal of Hydrometeorology, 6(1), 34–52.
25
[8] Diebolt, J., Guillou, A., & Rached, I. 2007. Approximation of the distribution of excesses through a
26
generalized probability-weighted moments method. Journal of Statistical Planning and Inference, 137(3),
27
841–857.
28
[9] Domínguez-Mora, R., Arganis-Juárez, M.L., Mendoza-Reséndiz, A., Carrizosa-Elizondo, E., Guzmán-
29
García, H., Echavarría-Soto, B., & Baños-Martínez, J.J. 2014. Time and spatial synthetic hourly rainfall
30
generation in the Basin of Mexico. International Journal of River Basin Management, 12(4), 367–375. 22
1
[10] Garavaglia, F., Gailhard, J., Paquet, E., Lang, M., Garçon, R., Bernardara, P., et al. . 2010. Introducing
2
a rainfall compound distribution model based on weather patterns sub-sampling. Hydrology and Earth
3
System Sciences Discussions, 14.
4
5
6
7
8
9
10
11
12
13
14
15
16
17
[11] Gardes, L., & Girard, S. 2010. Conditional extremes from heavy-tailed distributions: An application to the estimation of extreme rainfall return levels. Extremes, 13(2), 177–204. [12] Koutsoyiannis, D., & Onof, C. 2001. Rainfall disaggregation using adjusting procedures on a Poisson cluster model. Journal of Hydrology, 246(1), 109–122. [13] Kundu, P. K., & Siddani, R. K. 2011. Scale dependence of spatiotemporal intermittence of rain. Water Resources Research, 47(8). [14] Kyselý, J., Gaál, L., & J., Picek. 2011. Comparison of regional and at-site approaches to modelling probabilities of heavy precipitation. International Journal of Climatology, 31, 1457–1472. [15] Leblois, E. 2012. Le bassin versant, système spatialement structuré et soumis au climat. Mémoire présenté en vue de l’obtention du diplôme d’Habilitation à Diriger les Recherches, Université de Grenoble-INPG. [16] Lee, T., & Jeong, C. 2014. Nonparametric statistical temporal downscaling of daily precipitation to hourly precipitation and implications for climate change scenarios. Journal of Hydrology, 510, 182–196. [17] Marty, R., Zin, I., & Obled, C. 2008. On adapting PQPFs to fit hydrological needs: the case of flash flood forecasting. Atmospheric Science Letters, 9(2), 73–79.
18
[18] Mezghani, A., & Hingray, B. 2009. A combined downscaling-disaggregation weather generator for
19
stochastic generation of multisite hourly weather variables over complex terrain: Development and multi-
20
scale validation for the Upper Rhone River basin. Journal of Hydrology, 377(3), 245–260.
21
22
[19] Müller, H., & Haberlandt, U. 2015. Temporal rainfall disaggregation with a cascade model: from singlestation disaggregation to spatial rainfall. Journal of Hydrologic Engineering, 20(11), 04015026.
23
[20] Nadaraya, E. A. 1964. On estimating regression. Theory of Probability & Its Applications, 9(1), 141–142.
24
[21] Obled, C., Bontron, G., & Garçon, R. 2002. Quantitative precipitation forecasts: a statistical adaptation
25
of model outputs through an analogues sorting approach. Atmospheric research, 63(3), 303–324.
26
[22] Pickands, J. 1975. Statistical inference using extreme order statistics. Annals of Statistics, 3, 119–131.
27
[23] Pui, A., Sharma, A., Mehrotra, R., Sivakumar, B., & Jeremiah, E. 2012. A comparison of alternatives
28
for daily to sub-daily rainfall disaggregation. Journal of Hydrology, 470-471, 138–157. 23
1
[24] Radanovics, S., Vidal, J.-P., Sauquet, E., Ben Daoud, A., & Bontron, G. 2013. Optimising predictor
2
domains for spatially coherent precipitation downscaling. Hydrology and Earth System Sciences, 17, p–
3
4189.
4
[25] Ripley, B. D. 1996. Pattern recognition and neural networks. Cambridge university press.
5
[26] Serinaldi, F. 2010. Multifractality, imperfect scaling and hydrological properties of rainfall time series
6
simulated by continuous universal multifractal and discrete random cascade models. Nonlinear Processes
7
in Geophysics, 17(6), 697–714.
8
9
[27] Southworth, H., Heffernan, J. E., & Metcalfe, P. D. 2017. texmex: Statistical modelling of extreme values. R package version 2.4.
10
[28] Tabary, P., Dupuy, P., L’Henaff, G., Guenguen, C., Moulin, L., Laurantin, O., Merlier, C., & Souber-
11
oux, J.-M. 2012. A 10-year (1997-2006) reanalysis of Quantitative Precipitation Estimation over France:
12
methodology and first results. IAHS-AISH publication, 255–260.
13
14
15
16
17
18
19
20
21
22
[29] Tramblay, Y., Neppel, L., & Carreau, J. 2011. Climatic covariates for the frequency analysis of heavy rainfall in the Mediterranean region. Natural Hazards and Earth System Sciences, 11(9), 2463. [30] Tramblay, Y., Neppel, L., Carreau, J., & Najib, K. 2013. Non-stationary frequency analysis of heavy rainfall events in southern France. Hydrological Sciences Journal, 58(2), 280–294. [31] Watson, G. S. 1964. Smooth regression analysis. Sankhy¯ a: The Indian Journal of Statistics, Series A, 359–372. [32] Westra, S., Mehrotra, R., Sharma, A., & Srikanthan, R. 2012. Continuous rainfall simulation: 1. A regionalized subdaily disaggregation approach. Water Resources Research, 48(1). [33] Yiou, P. 2014. Anawege: a weather generator based on analogues of atmospheric circulation. Geoscientific Model Development, 7(2), 531–543.
24