Improving precipitation interpolation using anisotropic variograms derived from convection-permitting regional climate model simulations

The consideration of the spatial variability of daily precipitation, assessed through spatial covariance, is crucial for hydrological modeling. Estimating this covariance is particularly challenging in regions with sparse rain gauge networks or limited radar coverage. To address this issue, this study explores the potential of Convection-Permitting Regional Climate Model (CP-RCM) simulations to estimate anisotropic variograms. We compare five approaches: (1) SPAZM, an interpolator based on local precipitation-altitude regressions, Trans-Gaussian Random Fields, differing by their covariance structure and data source with (2) isotropic covariance from rain gauges, (3) anisotropic covariance from rain gauges, (4) isotropic covariance from CP-RCM simulations, and (5) anisotropic covariance from CP-RCM simulations. The models are evaluated with cross-validation and spatial metrics using radar-derived analyses. Results demonstrate that Trans-Gaussian Random Fields outperform SPAZM. Anisotropic covariance models derived from CP-RCM simulations capture orography-induced directional precipitation structures more effectively than the other models, leading to improved interpolation accuracy and better representation of spatial variability. The generated ensemble of conditional simulations successfully reproduces intense precipitation events at the catchment scale, providing valuable uncertainty quantification. For a 17 km² catchment, mean catchment precipitation can range from 175–450 mm for a convective event, despite high rain gauge density. These findings highlight the benefits of using CP-RCM simulations to generate anisotropic variograms for probabilistic precipitation interpolation. This approach improves the spatial variability of precipitation, making it highly relevant for hydrological applications such as flood forecasting. Future work will explore the integration of these ensembles into probabilistic hydrological modeling.

1 Introduction

Gridded daily precipitation data is essential for various environmental applications, such as assessing flood and drought risks or modeling glacier mass balance. However, rain gauge stations provide sparse and irregular observations in space, necessitating spatial interpolation models to estimate precipitation fields. Common interpolation procedures include local regression (Daly et al., 1994; Gottardi, 2009; Verdin et al., 2016), data assimilation (Alpuim and Barbosa, 1999; Devers et al., 2021; Vernay et al., 2024), geostatistics (Goovaerts, 2000; Sideris et al., 2014; Guédé et al., 2024), and more recently machine learning (Hengl et al., 2018; Sekulić et al., 2020; Zandi et al., 2022) models. Geostatistical models frequently outperform other statistical models for precipitation interpolation (Haberlandt, 2007; Bostan et al., 2012; Masson and Frei, 2014). They are considered Best Linear Unbiased Estimator (BLUE) methods (Rao, 1973), minimizing error variance but often producing excessive spatial smoothing that underestimates high precipitation intensities (Hofstra et al., 2010; Hiebl and Frei, 2018). To mitigate this issue, researchers use ensembles of equiprobable fields, known as conditional simulations (Frei and Isotta, 2019; Yan et al., 2021), all consistent with the measurements and the observational variance while displaying distinct spatial patterns.

Daily precipitation exhibits spatial autocorrelation as neighboring rain gauges often record similar values (Tobler, 1970). Spatial autocorrelation is mathematically modeled through the variogram (Cressie, 1991), which links spatial distance to observational variability. Traditionally, variograms are assumed to be isotropic, ignoring directional variations due to the limited number of rain gauge stations or for modeling simplicity (Adhikary et al., 2017). However, in complex topographic regions, anisotropy can arise from interactions between atmospheric conditions and mountain ranges (Tobin et al., 2011), degrading interpolation quality with isotropic variograms. To improve variogram estimation, researchers have explored alternative sources of spatial information, such as gridded precipitation products. Radar data has been used for deriving parametric or non-parametric variograms (Velasco-Forero et al., 2009; Schiemann et al., 2011) but is often available for shorter timeframes than the multidecadal period required for the daily precipitation analyses. Numerical weather prediction ensembles have also been explored (Khedhaouiria et al., 2022) to infer background error covariances in data assimilation approaches. Alternatively, simulations from Convection-Permitting Regional Climate Models (CP-RCMs) (Rockel et al., 2008; Brousseau et al., 2016; Keuler et al., 2016; Gerber et al., 2018) are supplied on extended periods and deliver new possibilities to infer variogram estimation. Although daily CP-RCM precipitation fields are frequently biased (Caillaud et al., 2021), they might still provide relevant information about anisotropy structure. This study evaluates the accuracy of daily precipitation interpolation by comparing isotropic and anisotropic variograms derived from rain gauges. It also investigates the potential of CP-RCM simulations for deriving anisotropic variograms, a novel approach that has not been previously explored. Additionally, the study aims to quantify the uncertainty in spatial interpolation at the catchment scale, a critical factor for hydrological modeling, which is rarely found in precipitation analyses or reanalysis (Frei and Isotta, 2019; Devers et al., 2021, 2024).

To address these objectives, we evaluate the ensemble means (kriging) and ensemble spreads (conditional simulations) from geostatistical models in a cross-validation framework across 786 intense precipitation events. We compare four types of variogram: (1) isotropic variogram estimated with rain gauge stations, (2) anisotropic variogram estimated with rain gauge stations, (3) isotropic variogram estimated with CP-RCM simulations, and (4) anisotropic variogram estimated with CP-RCM simulations. The geostatistical models are also compared to SPAZM (Gottardi, 2009), an interpolator based on local precipitation-altitude regressions stratified by weather patterns. Beyond point-scaled validation, we assess spatial structure using radar-derived analyses as a reference and evaluate the ability of the ensemble spread to capture mean catchment precipitation, an essential factor for hydrological modeling. The analysis is conducted over a mountainous region near the French Mediterranean Sea. This area encounters heavy rainfall due to the interaction of the air masses with elevation, offering an ideal setting for the study. The paper is structured as follows: Section 2 describes the domain under study and the available precipitation datasets. Section 3 introduces the four covariance models used for spatial interpolation. Section 4 displays the cross-validation and spatial evaluation results. Section 5 compares the findings to the existing literature and discusses future perspectives. Finally, Sect. 6 outlines the conclusions.

2 Domain and data under study

2.1 Study domain

The present study focuses on the Cevennes region in the southern Central Massif (see Fig. 1). This mountainous study area is characterized by rugged terrain, including high plateaus and forested hills, with elevations reaching 1700 m. Annual precipitation experiences significant spatial variability. The higher elevations can collect more than 2000 mm while the lower foothills present a Mediterranean climate with annual precipitation ranging from 600–1200 mm (Canellas et al., 2014). Winter precipitation can fall as snow at high elevations, even if snow cover is intermittent. Autumn is the wettest season. The moisture-laden air from the warm Mediterranean Sea faces colder air from the north, leading to frequent extreme rainfall events of 300–500 mm and generating heavy floods in the valleys (Delrieu et al., 2005).

Figure 1 Study domain (1.7–5.5° E, 42.8–45.6° N), colored by the elevation (height above sea level in meters), borders of the Cevennes region (black lines), borders and ID of four catchments (red lines, used in Sects. 3–5), and location of rain gauge stations (dots, around 800 stations available per day). The colors of the stations (dots) are explained in Sect. 3. Table 1 gathers information about the catchments.

Table 1 ID, name, surface area in km², mean and maximum altitude in m of the four catchments. Figure 1 shows the locations of the catchments.

2.2 Rain gauge observations

To capture the spatial variability of precipitation induced by the complex topography, a dense network of 973 rain gauges is available, including 197 stations in the study domain and 776 additional stations near the borders of the study domain. Rain gauges are up to 1500 m, encompassing the full elevation distribution. The number of available rain gauges on a daily timescale varies over the years, ranging from a minimum of 625 to a maximum of 836. The observed data were not quality-controlled in this study. Because we do not work with climatological statistics, temporal discontinuity is not a major issue. However, outliers may still occur and could affect the interpolated fields. Such outliers may influence local results, but we do not expect them to alter the conclusions of this study. Before applying the method in an operational or real-time context, a quality control procedure is needed, including a spatial anomaly analysis to identify outliers that exceed physically plausible differences from nearby stations under similar terrain characteristics, followed by a homogeneity test such as the Standard Normal Homogeneity Test (SNHT, Alexandersson, 1986).

2.3 Meteorological data

2.3.1 COMEPHORE

The product called COmbinaison en vue de la Meilleure Estimation de la Precipitation HOraiRE (COMEPHORE) (Champeaux et al., 2009) is a high-resolution precipitation analysis product for France, developed at Meteo-France, and based on radar-gauge merging to generate spatial precipitation estimates at an hourly timescale with a spatial resolution of 1 km. COMEPHORE is assumed as a reference to assess the spatial variability of precipitation, taking advantage of the radar coverage. COMEPHORE covers the 1997–2025 period but presents a methodology change in 2007, which separates interpolation for convective and stratiform precipitation. It uses the radar corrected by rain gauges to estimate convective precipitation and performs ordinary kriging of the rain gauge amounts for stratiform precipitation using an isotropic variogram. Using radar data, COMEPHORE can accurately represent the size and direction of precipitation cells, and the very local storm cells that the rain gauge network could miss. For this study, we aggregated the COMEPHORE hourly precipitation fields at the daily timescale over the 2008–2017 period. Some studies use COMEPHORE as a reference to evaluate climate models (Caillaud et al., 2021), and others use it for hydrological modeling (Evin et al., 2024).

2.3.2 AROME

AROME simulations (Caillaud et al., 2021) are produced with the convection-permitting RCM AROME in its NWP configuration cycle 41t1, which uses 60 vertical levels from 10 m to 1 hPa, including 21 levels below 2000 m to better resolve the lower-tropospheric dynamics over complex Alpine terrain. In this CP-RCM configuration, deep convection is explicitly resolved, while only shallow convection remains parameterized. The AROME domain over the Alpine region lies approximately 300–400 km from the lateral boundaries, which are forced by hourly outputs from the CNRM-ALADIN RCM (Nabat et al., 2020). ALADIN uses 91 vertical levels together with spectral nudging to ensure consistency with the large-scale circulation imposed by the ERA-Interim reanalysis (Dee et al., 2011). AROME simulations are available at the hourly timescale for the Alpine region, as described in the Flagship Pilot Study of the Coordinated Regional Climate Downscaling Experiment (CORDEX-FPS, Fantini et al., 2018), at 2.5 km spatial resolution, and cover the 1982–2018 years. Hourly outputs are aggregated to a daily timescale. Previous studies (Ban et al., 2021; Caillaud et al., 2021; Monteiro et al., 2022) show that AROME provides a more realistic representation of intense precipitation than its driving model ALADIN, despite persistent biases. In a Lagrangian evaluation over the Mediterranean region, which includes our study domain, Caillaud et al. (2021) reports that AROME simulations reproduce well the location, intensity, frequency, and interannual variability of heavy precipitation events. Remaining biases are mostly due to the model, rather than insufficient constraints from ERA-Interim. In the AROME model, very intense daily amounts ($>\mathrm{200}\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$) tend to be underrepresented, the spatial extent of intense convective cells is overestimated, and their propagation speed is slightly too high. These biases could be reduced by further refining horizontal and vertical resolution and by improving the parameterization of residual shallow and dry convection. Unlike a short-term NWP forecast, AROME-climate simulations do not assimilate observations such as radar reflectivity; therefore, they are not expected to reproduce individual events exactly, but rather their typical spatial structures. Consequently, even if absolute precipitation amounts may be biased, the spatial organization and anisotropy of intense precipitation systems, key for informing the anisotropic covariogram, may be sufficiently well captured by AROME to support our interpolation framework.

2.3.3 SPAZM

SPAZM (Gottardi, 2009) creates daily precipitation analyses at 1 km spatial resolution over a large part of France since 1948 and is continually updated by Electricité de France (EDF) for hydrological modeling. The daily precipitation fields constitute an adjustment of a climatological background field using the daily rain gauge amounts. A climatological background field represents the average daily structure of the precipitation field, where the influence of topography is easier to account for. The background fields incorporate local orographic effects conditioned by eight weather patterns. SPAZM leads to balanced water budgets at the annual scale, even in mountainous areas (Ruelland, 2020). However, the spatial variability of the daily fields is questionable, with an unrealistic correlation to the altitude. Moreover, SPAZM underestimates the spatial variance of intense precipitation, resulting in overly smoothed precipitation fields (Penot, 2014).

2.4 Study period

The study period ranges from 1982–2018, which corresponds to the availability of AROME simulations. 786 precipitation events (nearly 20 events per year), defined as the days with at least 50 mm recorded at a minimum of five rain gauges, are selected. Most of the 786 precipitation events arise from southerly atmospheric flow and central depression patterns, where strong orographic intensification is expected. Some winter events also include snowfall at high elevations.

3 Methods

3.1 Spatial Interpolation

This study develops a probabilistic geostatistical framework for daily precipitation interpolation, using four variograms derived from rain gauges and CP-RCM simulations. Unlike traditional approaches that rely solely on rain gauge data, our method exploits CP-RCM fields to estimate daily isotropic and anisotropic covariance structures. In the following, we detail the stochastic modeling and evaluation procedures.

3.1.1 Data transformation

Geostatistics involving kriging methods are built in the Gaussian paradigm (Diggle et al., 2003). However, daily precipitation data often exhibits a positively skewed distribution due to left truncation at zero precipitation and the presence of high values. To address this issue, we assume the positive precipitation follows a Gamma distribution and apply quantile-quantile transformation (Gyasi-Agyei, 2018) to normalize the data each day. We perform Gamma fitting using the maximum-likelihood approach from the fitdistr package (Delignette-Muller and Dutang, 2015). Spatial interpolation is then conducted in the Gaussian domain. Skewness and heteroscedasticity of the precipitation data are thus explicitly considered (Erdin et al., 2012), theoretically resulting in greater interpolation uncertainty in areas with high precipitation.

The modelization steps include:

1. Gamma to Normal transformation of positive precipitation:

   $$Y^{\ast }={\mathrm{\Phi }}^{-\mathrm{1}}[u_{\mathrm{0}}+(\mathrm{1}-u_{\mathrm{0}}\left)F_{\mathit{\gamma }}\right(Y\left)\right],$$

   where Y represents daily positive observed precipitations, F_γ the distribution function of a Gamma distribution, Φ⁻¹ is the quantile function of a Gaussian distribution with zero mean and unit variance, u₀ is the empirical probability of zero precipitation estimated by the number of rain gauges with less than 0.5 mm over the total number of rain gauges, and Y^∗ are the normalized precipitations.

2. Gamma to Normal transformation of zero precipitation:

   The zeros are transformed to ${\mathrm{\Phi }}^{-\mathrm{1}}\left[\mathcal{U}\right(\mathrm{0.01},u_{\mathrm{0}}\left)\right]$, 𝒰 is the uniform distribution.

3. Conditional simulations from geostatistical models (detailed in Sect. 3.1.2).

4. Back-Transformation of precipitation predictions into Gamma distribution:

   $$\widehat{Y}=\left\{\begin{array}{ll}{F}_{\mathit{\gamma }}^{-\mathrm{1}}\left\{\frac{\mathrm{\Phi }\left(Y^{\ast }\right)-u_{\mathrm{0}}}{\mathrm{1}-u_{\mathrm{0}}}\right\},& \mathrm{\Phi }\left(Y^{\ast }\right)>u_{\mathrm{0}},\\ \mathrm{0},& \text{otherwise}.\end{array}\right.$$

Once the precipitation data is normalized, we model its spatial structure using geostatistical models, focusing on estimating covariance structures through variogram fitting.

3.1.2 Geostatistic models

Precipitation should be centered in the first instance to estimate its covariance structure accurately (Schabenberger and Gotway, 2005). To handle the strong non-Gaussianity of daily precipitation (positivity and skewness), we assume that the rainfall field can be represented as a Trans-Gaussian Random Field. That is, there exists a transformation such that the transformed field follows a second-order stationary Gaussian random field. The expectation of the Trans-Gaussian Random Field is modeled in the Gaussian domain as a function of geographical predictors (longitude, latitude, and altitude) and a seasonal climatological background field. The literature commonly employs geographical predictors such as altitude in kriging with external drift model (Bárdossy and Pegram, 2013). The seasonal climatological background fields have been built through Random Forest modelization using CP-RCM simulations in (Dura et al., 2024). Their construction is based on the idea that the fine-resolution simulations of climate models should summarize the topographical influence on precipitation. Climatological background fields have been used broadly (see e.g. Hunter and Meentemeyer, 2005; Gottardi, 2009; Masson and Frei, 2014) to summarize long-term precipitation patterns influenced by orography, improving the robustness of daily precipitation interpolation. We successfully check the absence of multicollinearity by computing the correlation matrix among the predictors.

Spatial interpolation requires the estimation of a spatial covariance. We fit each day four different exponential variograms, all including a nugget effect:

• an isotropic variogram estimated with rain gauge observations (rgISO),

• an anisotropic variogram estimated with rain gauge observations (rgANISO),

• an isotropic variogram derived from daily AROME precipitation field (arISO),

• an anisotropic variogram derived from daily AROME precipitation field (arANISO).

The nugget parameter corresponds to micro-scale precipitation variability and measurement errors. The exponential choice is standard in numerous geostatistical studies (Masson and Frei, 2014; Frei and Isotta, 2019) because it effectively captures the gradual decline in autocorrelation as the separation distance increases. The weighted least-square estimation of the variograms is done in the Gaussian domain for both rain gauge and AROME precipitations. The weights are proportional to the number of pairs of observations and inversely proportional to the squared average distance between them. We selected only 25 % of AROME grid cells to compute variograms to balance spatial representativeness and computational cost. These selected cells were treated as virtual gauges when computing empirical variograms. To confirm robustness, we performed a sensitivity test using 100 % of cells for a subset of events, which produced similar anisotropy parameter estimates (not shown). We note η and θ as the estimated anisotropy ratio and angle. Schiemann et al. (2011) describes the methodology in more detail, with an ordinary least square estimation in contrast to a weighted least square estimation. The variograms are estimated with the gstat packages (Pebesma, 2004).

3.1.3 Conditional simulations

Conditional simulation of Trans-Gaussian Random Field is conducted using sequential Gaussian simulations implemented in the gstat package (Pebesma, 2004) to create each day an ensemble of 100 simulations. Gyasi-Agyei (2018) states that the conditional simulations obtained with gstat are too granular. They decided to apply a 3×3 window smoothing technique. We do not apply the same post-processing step because we will, in future developments, employ mean catchment precipitation, which is not sensitive to the granularity. While this ensemble size is sufficient for representing precipitation variability (Frei and Isotta, 2019), it does not account for all sources of uncertainty, such as covariance parameters inference or precipitation undercatch (see Sect. 5.3).

3.2 Evaluation

In our analysis, the ensemble mean of the conditional simulations is equivalent to the result obtained using deterministic kriging with external drift (Schabenberger and Gotway, 2005) (KED). Later in the article, we will use the terms “kriging” and “ensemble mean” interchangeably. While kriging captures the most likely precipitation value, the conditional simulations capture the entire ensemble of possible realizations, providing spatial uncertainty. To assess the performance of the spatial interpolation models, we evaluate both deterministic and probabilistic interpolated precipitation fields. The evaluation focuses on three key aspects:

• point-scale accuracy using cross-validation on rain gauge stations;

• spatial structure by comparing interpolated fields with radar-derived precipitation analyses;

• uncertainty quantification by analyzing the ensemble spread of conditional simulations at the catchment scale.

The models are evaluated on 786 precipitation events (nearly 20 events per year), defined as the days with at least 50 mm recorded at a minimum of five rain gauges.

3.2.1 Cross-validation of ensemble means

The deterministic performance of the spatial interpolation models is evaluated through a leave-cluster-out cross-validation scheme. Unlike traditional leave-one-out cross-validation, this procedure removes entire clusters of neighboring stations to better mimic ungauged conditions and accentuate performance differences between the models. An accurate covariance estimation is crucial for sparse rain gauge station networks. Sequentially, a cluster of neighboring stations is removed from the training dataset and then predicted. Ten groups of stations are created by the K-means algorithm (MacQueen, 1967). The clustering variables are the longitude and the latitude of the stations. Figure 1 indicates the rain gauge station clusters.

The evaluated models include (1) SPAZM as the deterministic interpolation model (2) rgISO: Kriging with an external drift with an isotropic variogram estimated from rain gauge data, (3) rgANISO: Kriging with an external drift with an anisotropic variogram estimated from rain gauge data, (4) arISO: Kriging with an external drift with an isotropic variogram derived from CP-RCM precipitation fields, (5) arANISO: Kriging with an external drift with an anisotropic variogram derived from CP-RCM precipitation fields. Mean error (ME) measures bias, indicating systematic over- or underestimation. Mean absolute error (MAE) quantifies overall prediction accuracy.

3.2.2 Spatial Structure Evaluation with Radar Analyses

Beyond point-scale accuracy, we assess whether interpolated precipitation fields reproduce the true observed spatial variability. The Teweles–Wobus Score (TWS, Teweles and Wobus, 1954) is used to compare the spatial gradients of interpolated precipitation fields with those from the radar-based COMEPHORE fields. Section 2.3.1 describes the COMEPHORE data. We consider COMEPHORE fields as a reference in the study domain because of the good radar coverage, allowing an accurate representation of the size and direction of the heavy precipitation cells. It should be noted that for stratiform events, COMEPHORE could favor geostatistic models with isotropic covariance by construction.

TWS is defined as:

$$\text{TWS}={\displaystyle \frac{{\sum }_{(s,s^{\prime })\in \text{Adj}}\left|(p_s-p_{s^{\prime }})-({\widehat{p}}_s-{\widehat{p}}_{s^{\prime }})\right|}{\mathrm{2}{\sum }_{(s,s^{\prime })\in \text{Adj}}max\left(|p_s-p_{s^{\prime }}|,|{\widehat{p}}_s-{\widehat{p}}_{s^{\prime }}|\right)}},$$

where Adj is the set of adjacent grid points $(s,s^{\prime })$ in the northern-southern and eastern-western direction within the study domain, p_s (${\widehat{p}}_s$) is the COMEPHORE (predicted) daily precipitation amount at the grid point s divided by the maximum grid amount. In this part, the interpolated fields are obtained using all rain gauges in the study domain. This metric is particularly relevant for evaluating the performance of variograms in capturing directional precipitation patterns. Lower values of TWS indicate better agreement with radar-based spatial structures.

Additionally, spectral analysis allows comparison of spatial variability across scales and is well-suited to assess whether predictions preserve the multiscale structure of precipitation. The two-dimensional Fourier power spectrum was computed for the reference (COMEPHORE) and predicted fields, and averaged radially in wave-number. For scaling processes, the radius-averaged spectrum follows a power-law relationship: $E\left(k\right)\propto k^{-\mathit{\beta }}$, where k is the wave number and β is the spectral slope. The slope β, estimated from a linear fit in log–log space, was used as the comparison metric. Similar values of β indicate similar spatial variability between reference and predicted fields. We compare the slopes using ME, RMSE, and slope best-fit line metrics. Whereas TWS evaluates short-range gradient image similarities, this comparison evaluates spatial multi-scale dependencies of the daily precipitation fields.

3.2.3 Probabilistic Performances and Uncertainty Quantification

The ensemble spread of conditional simulations provides uncertainty estimates for interpolated precipitation fields. We evaluate whether these uncertainty estimates are both statistically reliable and hydrologically relevant.

We evaluate probabilistic performances from the geostatistical models (rgISO, rgANISO, arISO, arANISO) using Continuous Ranked Probability Score (Matheson and Winkler, 1976) (CRPS) metrics on the 786 precipitation events using the same leave-one-cluster out procedure. A lower CRPS indicates better probabilistic distribution. To assess whether interpolation uncertainty translates into hydrologically realistic precipitation estimates, we compare, for the 20 most intense events, the ensemble spread of mean catchment precipitation (from conditional simulations) to the radar-derived mean catchment precipitation (from COMEPHORE). A reliable spatial interpolation model should ensure that the observed catchment precipitation falls within the simulated uncertainty.

Table 2 summarizes the considered metrics.

Table 2 Metrics computed for evaluation.

3.2.4 Case studies

We illustrate applications of the spatial interpolation models to provide concrete examples of how variograms affect precipitation interpolation. In Sect. 4.4, we show ensemble means, a sample of conditional simulations, and fitted anisotropic variograms for the 26 May 2008 and the 18 September 2014 days. The first event is an example of a South–North anisotropic event. The second is a Cevenol episode, with stationary storm cells, resulting in 300–400 mm precipitation amount in the first foothills of the Cevennes region.

4 Results

We conduct cross-validation at rain gauge stations and evaluate spatial structure using a radar-based metric to assess the performance of deterministic and stochastic models. The reliability of ensemble spreads is also examined. The ensemble of simulations provides uncertainty quantification, which is crucial for hydrological applications. Additionally, we illustrate the methodology with two case studies of daily precipitation fields, showing ensemble means and conditional simulations.

4.1 Evaluation of ensemble means

Leave-one-cluster-out cross-validation evaluates the accuracy of ensemble means for all rain gauge stations in the Cévennes region. We compute the mean error and mean absolute error on intense precipitation (≥50 mm) for the 786 events.

Figure 2 Cross-validation mean error (ME, (a)) and mean absolute error (MAE, (b)) for rain gauge stations in the Cevennes region, considering days with more than 50 mm. Results are shown for deterministic models (SPAZM) and stochastic models (rgISO, rgANISO, arISO, arANISO). For stochastic models, the cross-validation is conducted on the mean of 100 ensemble members from conditional simulations. Black contours indicate the Cevennes region.

SPAZM has lower spatial interpolation skills than rgISO, rgANISO, arISO, and arANISO (Fig. 2). Notably, SPAZM systematically underestimates intense precipitation in the northernmost high-altitude areas, with mean errors ranging from −30 % to −50 %, whereas other models (rgISO, rgANISO, arISO, arANISO) show less severe underestimation (−20 % to −30 %). SPAZM simultaneously underestimates and overestimates precipitation at neighboring stations, raising concerns about its robustness. rgISO, rgANISO, arISO, and arANISO all underestimate intense precipitation on average. The anisotropy contribution is only visible in the mean absolute error for the most northeasterly rain gauge stations.

4.2 Spatial evaluation using radar-analysis

To complement cross-validation, the TWS investigates the image gradient similarities, taking COMEPHORE fields as references. Lower TWS values indicate a better resemblance to COMEPHORE. Figure 3 illustrates the distribution of TWS values for three event classes: (1) Low anisotropy events ($\mathrm{0.75}\le \mathit{\eta }\le \mathrm{1}$), (2) Medium anisotropy events ($\mathrm{0.5}\le \mathit{\eta }\le \mathrm{0.75}$), (3) Strong anisotropy events ($\mathrm{0}\le \mathit{\eta }\le \mathrm{0.5}$). η results from the anisotropy ratio estimation with the daily AROME fields. TWS scores show that SPAZM exhibits lower image gradient similarity to COMEPHORE than all geostatistical approaches, regardless of the event type, indicating a poorer representation of spatial gradients. For both isotropic and anisotropic configurations, AROME-based methods (arISO and arANISO) consistently outperform their counterparts that use only rain gage data (rgISO and rgANISO). This systematic improvement suggests that the use of AROME fields allows for a more accurate estimation of the precipitation covariance structure and, consequently, a better reproduction of spatial gradients. Using only rain gauges, rgANISO outperforms rgISO for strongly anisotropic events, whereas the opposite behavior is observed for weakly anisotropic events. This indicates that accounting for anisotropy can improve gradient similarity when anisotropy is pronounced, but the limited robustness of anisotropic variogram estimation from rain gauges alone prevents drawing firm conclusions for weakly anisotropic cases. In contrast, arANISO consistently outperforms arISO, with a marked improvement for strongly anisotropic events and a slight but systematic improvement for weakly anisotropic ones. As no robustness issues arise in the estimation of anisotropic variograms when using AROME fields, these results provide stronger evidence of the added value of explicitly accounting for anisotropy in the representation of the spatial variability of precipitation. COMEPHORE relies on kriging for stratiform precipitation and radar observations for convective precipitation. While COMEPHORE uses isotropic covariance in the kriging step, anisotropy is introduced through the radar component. This suggests that anisotropic variograms, particularly when informed by AROME fields, help better reproduce the spatial structure of radar-derived precipitation patterns.

Figure 3 Distribution of TWS scores for deterministic models (SPAZM) and stochastic models (rgISO, rgANISO, arISO, arANISO). Low values of TWS scores indicate better image gradients, using COMEPHORE fields as references. The distributions of TWS scores are plotted for low estimated anisotropy ratio ($\mathrm{0.75}\le \mathit{\eta }\le \mathrm{1}$, left), medium estimated anisotropy ratio ($\mathrm{0.5}\le \mathit{\eta }\le \mathrm{0.75}$, center), and strong estimated anisotropy ratio ($\mathrm{0}\le \mathit{\eta }\le \mathrm{0.55}$, right). The number of events is displayed per event class. For the stochastic models, the evaluation is conducted over the mean of 100 ensemble members of conditional simulations.

Figure 4 presents the spectral slopes of the reference and predicted precipitation fields. SPAZM systematically underestimates the highest slopes, indicating an overestimation of fine-scale spatial variability. Compared to rgISO, rgANISO degrades the representation of spatial variability for weakly anisotropic events and provides only limited improvement for strongly anisotropic events. In contrast, arISO and arANISO improve the representation of spatial variability, particularly for strongly anisotropic events. Among them, arANISO better reproduces spatial variability, as the highest spectral slopes are no longer underestimated.

Figure 4 Scatter plots of the slopes of the radius averaged spectral density, using COMEPHORE data (reference) and the other gridded precipitation products models (SPAZM, rgISO, rgANISO, arISO, arANISO). The scatter plots of slopes are plotted for low estimated anisotropy ratio ($\mathrm{0.75}\le \mathit{\eta }\le \mathrm{1}$, top), medium estimated anisotropy ratio ($\mathrm{0.5}\le \mathit{\eta }\le \mathrm{0.75}$, center), and strong estimated anisotropy ratio ($\mathrm{0}\le \mathit{\eta }\le \mathrm{0.55}$, bottom). The dashed line represents the diagonal. Best model are highlighted per anisotropy class and metrics (ME, RMSE, best-fit line). For the stochastic models, the slopes are computed over the mean of 100 ensemble members of conditional simulations.

4.3 Probabilistic performances and Uncertainty Quantification

To assess the reliability of ensemble simulations, we apply the same cross-validation approach used for evaluating ensemble means (see Sect. 4.1). The Continuous Ranked Probability Score (CRPS) is computed to quantify probabilistic prediction accuracy, where lower values indicate better performance. Distribution of CRPS scores reveals similar performances for rgISO, rgANISO, and arISO but a notable improvement with arANISO (Fig. 5). The error reduction is substantial for strong anisotropy events. This suggests that integrating anisotropy into the covariance structure using AROME enhances the ability of the ensemble to capture observed daily precipitation. Figure S1 in the Supplement demonstrates that this result holds across all weather regimes. Figure S2 in the Supplement further shows that the result also applies to lower precipitation intensity ranges (1–20 and 20–50 mm). For this reason, we later quantify precipitation uncertainty using conditional simulations from anisotropic covariance, derived from AROME simulations.

Figure 5 Distribution of cross-validation CRPS scores for stochastic models (rgISO, rgANISO, arISO, arANISO). 1 mean value of CRPS is obtained for each station within the Cevennes region for the days with more than 50 mm. Lower CRPS values indicate better spatial interpolation. The distributions of CRPS scores are plotted for low estimated anisotropy ratio ($\mathrm{0.75}\le \mathit{\eta }\le \mathrm{1}$, left), medium estimated anisotropy ratio ($\mathrm{0.5}\le \mathit{\eta }\le \mathrm{0.75}$, center), and strong estimated anisotropy ratio ($\mathrm{0}\le \mathit{\eta }\le \mathrm{0.5}$, right).

For the 20 most intense events, the mean catchment precipitations derived from the ensembles of 100-member arANISO from conditional simulations closely agree with COMEPHORE estimates (Fig. 6) with a high rain gauge density. The interquartile range of the simulations frequently contains the mean COMEPHORE precipitation for the Chassezac at Sainte-Marguerite and the Thines at Pont de Gournier catchments. For the Chassezac at Pont du Mas catchment, the interquartile range does not necessarily encompass COMEPHORE precipitation, but the full range of the ensemble does encompass it. This suggests the reliability of the uncertainty estimates provided by the ensemble. arANISO is robust, as conditional simulations nearly agree with COMEPHORE estimates in an ungauged scenario. The artificial decrease of rain gauge density causes a slight negative bias in the conditional simulations. The latter frequently contain COMEPHORE precipitations through wider confidence intervals.

Figure 6 Scatter plot between mean catchment COMEPHORE (x axis) and simulated (y axis) precipitations in two rain gauge density scenarios: using all rain gauges available (Gauged, left panels), removing the rain gauge of the main catchment (Ungauged, right panels). The 20 d with the highest COMEPHORE precipitation are considered for each catchment represented in the vertical panels. The simulations are derived from conditional simulations with anisotropic covariance (arANISO). The vertical lines in red (grey) illustrate the 25–75 (0–100) % range of simulations. Black points indicate the ensemble means. The dashed line represents the diagonal.

As expected, uncertainty increases with precipitation intensity, decreases with catchment size and rain gauge density, and provides asymmetrical distributions on mean precipitations. For instance, on 18 September 2014, the mean rainfall on the small Thines at Pont de Gournier catchment ranges from 175 mm to more than 450 mm with the complete set of rain gauges, representing a factor higher than 2.5. In contrast, the uncertainty is lower (factor of 1.33) for larger catchments such as the Chassezac at Sainte-Marguerite.

4.4 Case studies

To further illustrate model performance, we analyze precipitation fields for two selected events, using COMEPHORE as the reference. Section 2.3.1 describes the COMEPHORE data. Figure 7 presents daily precipitation maps for four catchments. While all models broadly reproduce the spatial structure of COMEPHORE fields, important differences emerge. SPAZM exhibits both local overestimation and underestimation, producing overly broad precipitation cells. rgISO also yields overly smooth fields, leading to high-intensity underestimation. arISO slightly mitigates this issue, as seen in the 18 September 2014 case, where it captures more intense precipitation than rgISO. rgANISO captures the radar anisotropy for the 26 May 2008 events but produces an inaccurate spatial structure for the 18 September 2014 event. arANISO provides the most realistic interpolation, reducing both excessive smoothing and underestimation of high intensities. It is the only stochastic model able to capture the highest precipitation values (388–444 mm) in the 18 September 2014 event, demonstrating the advantage of incorporating directional effects.

Figure 7 Daily precipitation field illustrations, with a focus on the four catchments for the 26 May 2008 (a) and 18 September 2014 (b) days. Daily precipitation fields are obtained with AROME, deterministic models (COMEPHORE, SPAZM) and stochastic models (rgISO, rgANISO, arISO, arANISO). For the stochastic models, the fields correspond to the mean of 100 ensemble members derived from conditional simulations. TWS (COMEPHORE as reference) are displayed for the gridded precipitation products. Precipitation is expressed in mm. Figure 1 shows the locations of the catchments.

A single precipitation field does not reflect the interpolation uncertainty, leading to unwanted smoothing spatial patterns. Conditional simulations provide an attractive alternative to derive a set of equiprobably plausible fields. Figure 8 illustrates five conditional simulations with anisotropic covariance (arANISO) for the 26 May 2008 and 18 September 2014 days, revealing different spatial patterns and precipitation intensities. The Thines at Pont de Gournier catchment, located near the center of the strong cell intensities, has a substantial mean precipitation uncertainty on 18 September 2014. Simulation 55 gives rise to nearly 200 mm mean catchment precipitation, compared to more than 500 mm for simulation 95.

Figure 8 Anisotropic (arANISO) illustration of conditional simulations of daily precipitation fields, with a focus on the four catchments for the 26 May 2008 (a) and 18 September 2014 (c) days. Five members are randomly sampled from a 100 ensemble, and the distribution of their precipitation values are shown ((b), (d)) for each catchment . Precipitation is expressed in mm. Figure 1 shows the locations of the catchments.

Figure 9 displays the AROME fields, as well as the isotropic and anisotropic fitted variograms using rain gauges or the corresponding AROME fields, for the two case studies. Although the AROME fields do not show the same spatial variability and precipitation intensity as COMEPHORE (Fig. 7), they include nearly the same anisotropy. The 26 May 2008 event shows a higher spatial correlation in the north–south direction. The 18 September 2014 one is characteristic of a Cevenol episode with a higher spatial correlation in the southwest-to-northeast direction. Figure 10 summarizes the anisotropy parameters η (anisotropy ratio) and θ (anisotropy angle) estimated from AROME and COMEPHORE. Both datasets show similar preferred anisotropy directions, predominantly oriented south–north (S–N) and southwest–northeast (SW–NE), with AROME exhibiting a slightly stronger S–N component. In both cases, the anisotropy is generally more pronounced along the SW–NE direction.

Figure 9 Examples of two AROME daily precipitation fields used to derive anisotropic variograms for the 26 May 2008 (a) and 18 September 2014 (c) days. Precipitation is expressed in mm. Black contours represent the Cevennes region. In the 2D variograms (b, d), low (high) gamma values indicate strong (weak) spatial autocorrelation between observations, spaced from dx(dy) in horizontal (vertical) km scale.

Figure 10 Anisotropy parameters estimated with daily AROME and COMEPHORE precipitation fields, namely the angle (θ) of dominant spatial autocorrelation and the anisotropy ratio (η) with the corresponding orthogonal direction.

5 Discussion

In this study, we derive anisotropic variograms from CP-RCM simulations for spatial interpolation of precipitation, an approach that, to our knowledge, has not been previously explored. Estimating anisotropy parameters based solely on rain gauge stations is challenging due to the sparse and irregular station networks. However, CP-RCM simulations provide a way of characterizing anisotropy. We demonstrate that anisotropic variograms improve probabilistic precipitation analysis compared to isotropic ones through both point-based and spatial evaluations. In the following, we discuss the advantages and limitations of our approach and outline the further developments that are needed.

5.1 Ensemble from conditional simulations versus local regression approach

This study highlights that Trans-Gaussian Random Fields (kriging with external drift) outperform local regression models in robustness and spatial variability. SPAZM, a local regression interpolator, excessively spreads precipitation cells and overcorrelates precipitation with altitude. This overcorrelation results in poor performance for convective events, leading to regional inconsistencies where SPAZM significantly underestimates or overestimates intense precipitation in localized areas. In contrast, the ensemble mean from Trans-Gaussian Random Fields shows better regional coherence, despite the global underestimation of intense precipitation. This bias is well-known in kriging studies (Hiebl and Frei, 2018). Conditional simulations from Trans-Gaussian Random Fields provide an alternative to generate more realistic precipitation. This approach accurately captures large-scale convective events. Further evaluation is needed in ungauged areas with localized convection cells and complex topography, where orographic effects drive uncertainty.

5.2 Covariance estimation

Estimating anisotropic parameters typically requires a high density of rain gauge stations, a condition rarely met in practice. Our results show that CP-RCM simulations, despite their imperfections, effectively capture anisotropic precipitation structures, making them valuable for informing variogram parameters in regions with limited observational data.

In this study, we estimate four covariance models using weighted least squares: (1) isotropic covariance from rain gauge stations (rgISO), (2) anisotropic covariance from rain gauge stations (rgANISO), (3) isotropic covariance from CP-RCM simulations (arISO), and (4) anisotropic covariance from CP-RCM simulations (arANISO). Considering anisotropy in the covariance structure using rain gauges (rgANISO) does not necessarily improve interpolation performances compared to isotropic covariance (rgISO). More specifically, for both point cross-validation and spatial variability evaluation, rgANISO outperforms rgISO during strongly anisotropic events, whereas the opposite behavior is observed for weakly anisotropic events. This indicates that anisotropy can be significant but remains difficult to robustly estimate using a rain-gauge network alone, even when the network is dense. This result agrees with the few studies on covariance comparisons (Haberlandt, 2007; Shi et al., 2007; Haylock et al., 2008). Using a limited number of rain gauges generates a loss of robustness in estimating the covariance parameters despite the actual anisotropy presented in the daily precipitation fields. arISO shows similar performances to rgISO. This suggests that the rain gauge density is high enough to capture spatial autocorrelation in the study region. This also implies that CP-RCM simulations could serve as a viable alternative for estimating spatial autocorrelation in areas with sparse rain gauge station coverage. Estimating anisotropic parameters usually requires a large number of rain gauge stations, which is rarely met in practice. arANISO outperforms the other stochastic models. This approach aligns with previous studies that have used radar data to infer anisotropy through nonparametric estimation (Velasco-Forero et al., 2009; Schiemann et al., 2011; Gyasi-Agyei, 2016). High-quality CP-RCM simulations, increasingly available worldwide, might correctly carry anisotropic information.

The covariance modeling in this study has two main limitations. Firstly, CP-RCM simulations may miss or inaccurately generate precipitation cells, leading to biases in the estimation of covariance parameters such as the range, the anisotropy ratio, and the anisotropy angle. Secondly, the assumption of second-order stationarity may not hold on large spatial domains. Topographical effects might induce spatial non-stationarity in the covariance structure. The windward mountainsides might encounter more spatial variability than the leeward ones. The study region has daily oceanic or Mediterranean influences, but not both simultaneously, limiting spatial non-stationarity. Several authors implement spatial non-stationarity of the covariance (Paciorek and Schervish, 2006; Risser and Calder, 2015; Risser et al., 2019). For computational reasons, we decided not to include non-stationarity. Extending the method to larger and topographically complex domains would require non-stationary covariance. A practical option is to partition the region into climatologically homogeneous sub-regions, ideally preserving major watershed boundaries to maintain hydrological consistency. Alternatively, a solution is to incorporate non-stationary covariance structures, for example, through geographical coordinate deformation (Youngman, 2023) or locally stationary covariance models (Paciorek and Schervish, 2006; Risser and Calder, 2015), which would allow spatial dependence to evolve smoothly across the domain. These approaches would make the method suitable for operational applications over larger domains.

Despite those limitations, we are confident in the study's transposability to other regions of similar size with high-quality CP-RCM simulations. Furthermore, the scientific community should explore whether the same methodology can be applied using RCM simulations in regions where CP-RCM simulations are unavailable. The methodology could also be extended to real-time, sub-daily interpolation. CP-RCM simulations are not continuously updated, so their replacement by hourly numerical weather forecasts (NWP) should be investigated. NWP assimilate past radar reflectivity and should therefore display a higher correlation with observations than CP-RCM simulations. As a result, NWP could allow us to extract precipitation intensity, spatial patterns, and spatial variability, while quantifying interpolation uncertainty through conditional simulations and the use of ensemble NWP forecasts. A natural follow-up would be to use NWP forecasts as both drift and covariance structures within a kriging-with-external-drift framework (Velasco-Forero et al., 2009; Schiemann et al., 2011). At the daily timescale, timesteps are typically considered independent, but this assumption no longer holds at the hourly scale. To address this, temporal dependence should be incorporated into the model, as done in Sideris et al. (2014); Frey and Frei (2025). A major limitation is the limited availability of sub-daily rain gauges. One potential solution to bypass this shortcoming is to disaggregate daily interpolated fields using radar data or numerical weather forecasts.

5.3 Quantification of interpolation uncertainty

We assess precipitation interpolation uncertainty at the catchment scale, which is relevant for hydrological applications. For very small catchments (17 km²), the simulated mean catchment precipitation can vary widely by a factor of 2.5, highlighting substantial uncertainty despite high rain gauge density. Using a gamma distribution to normalize daily observations creates an asymmetrical distribution of mean catchment precipitation, which appears to be physically reasonable. Other authors (Erdin et al., 2012; Frei and Isotta, 2019) frequently used the Box–Cox transformation to normalize data, causing uncertainties in the back-transformation step. Box–Cox transformations can lead to over- or under-normalization, artificially inflating or deflating the uncertainty in spatial interpolation. Standard square root transformation, a special case of Box–Cox (parameter equal to 0.5), often leads to flattening high precipitation values (Erdin et al., 2012), causing underestimation of spatial interpolation uncertainty.

Uncertainty in precipitation analysis mainly arises from precipitation undercatch and spatial interpolation methods.

In this study, rainfall prevails over snowfall, and precipitation undercatch might be considered negligible. However, in snow-dominated regions, precipitation undercatch is significant (Sevruk et al., 2009; Pollock et al., 2018) under windy conditions. We propose fitting an asymmetrical gamma distribution to the catch ratio to handle this issue. The parameters of this distribution would depend on several factors, including the precipitation phase (rain or snow), precipitation intensity, rain gauge type, its exposure, and wind speed.

While we modeled the primary sources of uncertainty related to spatial interpolation, additional uncertainty remains unaddressed:

• We assume fixed parameters for the gamma distribution used to normalize the data. A natural extension would be to generate multiple sets of gamma parameters and run conditional simulations for each.

• The seasonal climatological background fields used in this study are treated as deterministic despite inherent uncertainties. Further developments should explore probabilistic background fields to propagate uncertainty better.

• In addition, we did not consider the uncertainty related to the estimation of the covariance parameters. Bayesian inference to estimate covariance parameters (Frei and Isotta, 2019), using informative priors derived from CP-RCM simulations, is a logical extension.

Bayesian hierarchical models present an ideal framework to deal with this combination of interdependent uncertainties that can be propagated into hydrological models.

6 Conclusion

This study compares different covariance estimations for conditional simulations of daily precipitation using rain gauge stations. The evaluated models include (1) a local regression interpolator called SPAZM; Trans-Gaussian Random Fields with (2) isotropic covariance estimated with rain gauges, (3) anisotropic covariance estimated with rain gauges, (4) isotropic, and (5) anisotropic covariances estimated from daily CP-RCM simulations. The expectation of the Random Fields is a function of topographical variables and a seasonal climatological background field. We assess interpolation performance through cross-validation and a spatial metric, using radar-derived analysis as a reference.

Results indicate that Trans-Gaussian Random Fields modeling, whatever the covariance used, consistently outperforms the local regression interpolator SPAZM in cross-validation and spatial metrics. The geostatistical models are more robust to low rain gauge density and resemble radar fields more. Among the covariance models tested, anisotropic covariance derived from CP-RCM simulations better captures directional precipitation structures observed in radar data and shows superior cross-validation scores for both ensemble mean and spread. By bypassing the lack of robustness of anisotropy estimation using sparse rain gauge networks, this approach reveals the clear value of incorporating anisotropy into the spatial interpolation of daily precipitation. The obtained ensemble of 100 mean catchment precipitation simulations successfully encompasses radar-based estimates for intense precipitation events, highlighting its potential for uncertainty quantification, a key consideration for hydrological modeling.

These findings suggest that deriving anisotropic variograms from high-quality CP-RCM simulations is a promising approach for probabilistic precipitation interpolation. A natural extension of this work is to integrate the generated precipitation ensembles into probabilistic hydrological modeling, further improving flood risk assessment and water resource management.