# Bayesian network-informed conditional random forests for probabilistic multisite downscaling of precipitation occurrence

M.N. Legasa, R.E. Chandler, R. Manzanas

► To cite this version:

M.N. Legasa, R.E. Chandler, R. Manzanas. Bayesian network-informed conditional random forests for probabilistic multisite downscaling of precipitation occurrence. Environmental Modelling and Software, 2026, 199, pp.106907. 10.1016/j.envsoft.2026.106907. hal-05523996

HAL Id: hal-05523996

https://hal.science/hal-05523996v1

Submitted on 23 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

# Bayesian network-informed conditional random forests for probabilistic multisite downscaling of precipitation occurrence

M.N. Legasa<sup>a,b ID,*</sup>, R.E. Chandler<sup>c ID</sup>, R. Manzanas<sup>a ID</sup>

## ARTICLE INFO

**Keywords:**

Statistical downscaling

Bayesian network

Random forest

Multisite downscaling

Probabilistic downscaling

Spatio-temporal downscaling

## ABSTRACT

This work introduces Bayesian network-informed conditional random forests (BNICRF): a novel multiresponse classification method for downscaling the joint probability distribution of precipitation occurrence at multiple geographical locations from large-scale reanalysis predictors. BNICRFs combine a Bayesian network to model spatial dependencies and a set of random forests to predict local precipitation from large-scale inputs. Extending prior studies on Bayesian networks and random forests in climate downscaling, this method is validated under the experimental framework of the COST action VALUE (the largest, most exhaustive intercomparison study of statistical downscaling methods to date). Results demonstrate that BNICRFs effectively capture spatial relationships while maintaining single-site predictive performance comparable to single-site random forests and a well-performing generalized linear model in VALUE. Additionally, BNICRFs outperform a robust multisite approach proposed in Chandler (2020) in predictive capability while matching spatial performance. Incorporating temporal structures further enables BNICRFs to generate temporarily and spatially realistic precipitation occurrence fields.

## 1. Introduction

Global climate models (GCMs) are the main tool used nowadays to simulate the evolution of the climate system. However, their spatial resolution – typically around one hundred kilometers – is still insufficient for most practical applications (see e.g. Doblas-Reyes et al. 2013 and references therein). There are two main approaches aimed to tackle this issue. On the one hand, dynamical downscaling relies on regional climate models (RCMs, Giorgi 2019), which are nested into GCMs to solve the physical equations of the climate system over limited spatial domains, producing simulations at a much higher resolution. On the other hand, statistical downscaling (SD, von Storch et al. 1993) relies on the use of statistical and machine learning models to learn the empirical relationships that link a set of key large-scale predictors well reproduced by GCMs (like geopotential or synoptic winds) to the local surface predictand of interest (like precipitation or temperature). These two approaches come with their own advantages and disadvantages. RCMs are very costly in terms of computational and energy resources, typically requiring months to run over complex computing infrastructures. Differently, SD methods can be deployed on personal computers

and generally require at most a few days to run. Moreover, RCMs also tend to inherit part of the biases from their driving GCMs and introduce their own, thus requiring additional correction steps based on observed records. RCMs do, however, have the potential advantage that they are physically based and hence do not rely on an assumption, required in SD, that empirical relationships derived from historical data will remain valid into the future. In this article we focus on SD.

Under the perfect prognosis (PP) framework (Charles et al., 1999; Timbal et al., 2003; Bürger and Chen, 2005; Haylock et al., 2006; Fowler et al., 2007; Hertig and Jacobheit, 2008; Sauter and Venema, 2011; Gutiérrez et al., 2013), SD models are trained/calibrated based on observed data over a historical reference period, relying on some reanalysis (reanalyses combine various climate observations with numerical weather models to create a comprehensive and consistent gridded dataset of past climate conditions, see e.g. Dee et al. 2011) for the predictors. Typically, these trained/calibrated models are subsequently applied to downscale future GCM large-scale predictors in order to obtain the corresponding downscaled values. In this work we focus

<footnote><sup>a</sup> Departamento de Matemática Aplicada y Ciencias de la Computación (MACC), Universidad de Cantabria, Santander, Spain</footnote>

<footnote><sup>b</sup> Laboratoire des Sciences du Climat et de l'Environnement (LSCE-IPSL), CEA/CNRS/UVSQ, Université Paris Saclay, Centre d'Etudes de Saclay, Gif-sur-Yvette, France</footnote>

<footnote><sup>c</sup> Department of Statistical Science, University College London, London, United Kingdom</footnote>

<footnote>* Corresponding author at: Laboratoire des Sciences du Climat et de l'Environnement (LSCE-IPSL), CEA/CNRS/UVSQ, Université Paris Saclay, Centre d'Etudes de Saclay, Gif-sur-Yvette, France.</footnote>

<footnote>E-mail address: mikel.legasa@lsce.ipsl.fr (M.N. Legasa).</footnote>
exclusively on the training/calibration phase, i.e. we use reanalysis
predictors in a PP framework.

Daily precipitation, the focus of this study, is a key variable for
many applications in a wide range of sectors; however, it is notably
difficult to model, due to its positive mass at zero accounting for dry
days and a continuous non-normal intensity for wet days. Thus, most
of the currently available models deal separately with precipitation
occurrence (wet/dry) and precipitation amount (Coe and Stern, 1982;
Stern and Coe, 1984; Ambrosino et al., 2014; Zhou et al., 2019;
Chandler, 2020). This makes sense, since different physical mech-
anisms control the binary occurrence and the continuous amounts.
Historically, most SD methods have been developed for independently
downscaling at multiple locations, yielding predictions which fail to
take into account the relevant spatio-temporal relationships among
the downscaled locations. These relationships are key for producing
reliable multisite precipitation series, and are particularly important
in hydrological applications (Singh, 1997; Grayson and Blöschl, 2001).
Indeed, many applications require spatially consistent data, reflecting
not only the marginal statistics of the different sites, but also the
multisite (or spatial) statistics such as pairwise correlation (Fiener
and Auerswald, 2009; Haile et al., 2009). Streamflow also generally
depends strongly on the spatial precipitation structure (Segond et al.,
2007). As opposed to precipitation amount, which is often modeled via
copulas by using a transformation to marginal normality at each site
(which keeps the covariance matrix), modeling the spatial structure
of precipitation occurrence is notably problematic (Yang et al., 2005;
Chandler, 2020). Therefore, we focus in this work on precipitation
occurrence exclusively.

Several methodologies have been proposed in order to keep the
spatial dependence structure of downscaled precipitation occurrence.
Previous studies have tried to preserve this spatial structure by linking
the local climate covariance to the global circulation covariance, some-
thing often called expanded downscaling (Bürger, 1996). Within this
methodology, Cannon (2008) proposed an expanded artificial neural
network producing outputs that preserve the correlation matrix by
adding a penalization term in the cost function; and in Chandler (2020)
this was done by using latent Gaussian variables. In Yang et al. (2005),
the spatial structure was approximated by means of directly modeling
the number of wet sites using the beta-binomial distribution. Charles
et al. (1999) conditioned the downscaled series on a hidden weather
state variable. All these studies focused on particular statistics or aspects
of the joint probability distribution (JPD), instead of modeling the
JPD itself. Indeed, some approaches were designed for use in specific
situations, e.g. the beta-binomial model in Yang et al. (2005) was
designed for use in small areas where dependence is uniformly very
high, whereas the hidden Markov models of Charles et al. (1999) were
designed for use at large spatial scales at which days can be classified
into weather types with their own distinctive rainfall patterns.

To solve the limitations of the aforementioned studies, in this work
we propose a model that combines the explanatory power of Bayesian
networks (BN) (Pearl, 1988) and the predictive capability of random
forests (RF) (Breiman, 2001) to produce downscaled precipitation oc-
currence fields that keep the spatial JPD for all the target locations.
In essence, BNs provide an approximation of the multisite dependence
structure of the locations of interest (see Legasa and Gutiérrez 2020
for an in-depth description of BNs as weather generators), and RFs
extract the relevant information from the large scale predictors to
produce accurate multisite probabilistic predictions (see Legasa et al.
2022 for the application of random forest for probabilistic precipitation
downscaling). RFs have been extensively applied for modeling meteo-
rological variables, due to their skill handling non-linear, multi-class,
unbalanced, noisy and very high dimensional datasets (Hengl et al.,
2015; He et al., 2016; Fox et al., 2020; Georganos et al., 2021; Saha
et al., 2021). They also require very little to no pre-processing of the
data, and are robust to outliers (Belgiu and Drăgut, 2016; Rodriguez
et al., 2006; Wang et al., 2019). All this makes RFs, to date, the

most prevalent machine learning methodology for climate change risk
assessment studies (Zennaro et al., 2021).

In addition to the predictive capability of the models, in this work
we assess the reliability of the downscaled joint probability distribution
of precipitation occurrence. We do so by measuring the predictive
performance of the models and how accurately they represent sev-
eral multisite and single-site statistics when drawing samples from
the downscaled probability distributions. an approach often referred
to as stochastic downscaling or weather generation conditioned on
large-scale covariates.

The remainder of this article is structured as follows: Section 2 describes the experimental framework and data used. Section 3 explains the particularities of the different techniques considered, including a detailed description of the proposed methodology. The results obtained are presented through Section 4, with Section 5 gathering the main conclusions.

## 2. Experimental framework and data

The proposed methodology is evaluated following the approach
developed for Experiment 1c of the COST action VALUE (http://www.
value-cost.eu, Gutiérrez et al. (2019)), which was designed with the
goal of providing a common, transparent, reproducible and compre-
hensive framework for intercomparing SD methods for climate. The
experiment 1c (Widmann et al., 2019) considers a dense set of 53
stations within Germany and is specially designed to focus on the
spatial aspects of the predictions. For illustrative purposes, we test
our methodology on a smaller subset of 11 locations in southern
Germany, which are listed in Table 1, and are shown alongside with
the orography of the region in Fig. 1. As can be seen, the stations
are located in a region presenting a complex orography, with several
mountain ranges and valleys which yield a complex spatial distribution
of precipitation occurrence. This subset of 11 stations corresponds
to the same set of stations previously used in Legasa and Gutiérrez
(2020), thus allowing for a direct comparison of the results and for
better visualization of the in-depth analyses performed in this article,
which intends to illustrate the methodology. There is, in principle, no
constraint on the number of locations and, as discussed in Section 5,
we have also tested the methodology using the complete set of stations
(53), obtaining comparable results.

We use observational records of precipitation from the European
Climate Assessment & Dataset project (ECA&D, Klein Tank et al. 2002).
Data spanning the period from 1979 to 2008 (30 years) have been used,
and with a threshold of 1 mm for the discretization of precipitation
into dry (lower than 1 mm) and wet (higher than or equal to 1 mm)
days. For the purpose of illustrating the methodology, only days with
complete observations (i.e. with no missing values at any of the 11
locations) have been used, leaving 10,636 training instances (days).
Algorithms exist for handling missing data for both Bayesian networks
and random forests without much computational overhead; however,
the implementation of these is out of scope of this work.

As predictors we consider the same large-scale variables used in
Gutiérrez et al. (2019), all of them coming from ERA-Interim reanal-
ysis (Dee et al., 2011) on a daily basis over a 2° resolution regular
grid (blue crosses in Fig. 1). Note, however, that we have also tested
the methodology using a different reanalysis dataset (ERA5, Hersbach
et al. (2020)), concluding that the choice of reanalysis does not make
any difference on the results. As in Gutiérrez et al. (2019), we rely
on local predictor information, using the 4 nearest grid points to
each target location. The list of predictors considered is shown in
Table 2. Throughout the manuscript we use X to refer to this set of
predictors and to $Y = (Y^1, ..., Y^{11})$ to refer to the downscaled variable,
i.e. absence/presence of precipitation.

To avoid overfitting the statistical models, we applied the same 5-
fold cross-validation employed in Gutiérrez et al. (2019), in which the
total period of study, 1979–2008, is divided into 5 non-overlapping
<caption>**Table 1**</caption>

List of the 11 stations used in this work, which have been previously used in Legasa and Gutiérrez (2020). Note that the code given in this work (1st column) is different, and corresponds to an ancestral order given by the directed acyclic graph in Fig. 1. We also show the ECA&D identifier (ID), altitude in meters (Alt), longitude (Lon), latitude (Lat), Wet-day Frequency (Wet Days) and the location name.

<table><caption>**Table 2**</caption><thead><tr><th>Code</th><th>ID</th><th>Alt</th><th>Lon</th><th>Lat</th><th>Wet days</th><th>Location name</th></tr></thead><tbody><tr><td>1</td><td>4572</td><td>415m</td><td>10.17</td><td>49.39</td><td>32.2%</td><td>ROTHENBURG OB DER TAUBER</td></tr><tr><td>2</td><td>4083</td><td>657m</td><td>11.84</td><td>49.98</td><td>40%</td><td>FICHTELBERG OBERFRANKEN</td></tr><tr><td>3</td><td>4472</td><td>435m</td><td>10.51</td><td>48.83</td><td>30.4%</td><td>REIMLINGEN</td></tr><tr><td>4</td><td>4004</td><td>365m</td><td>12.1</td><td>49.04</td><td>30.4%</td><td>REGENSBURG</td></tr><tr><td>5</td><td>4007</td><td>921m</td><td>9.94</td><td>50.5</td><td>42%</td><td>WASSERKUPPE</td></tr><tr><td>6</td><td>4617</td><td>937m</td><td>10.77</td><td>50.66</td><td>43.3%</td><td>SCHMUCKE</td></tr><tr><td>7</td><td>488</td><td>1213m</td><td>12.96</td><td>50.43</td><td>42.7%</td><td>FICHTELBERG</td></tr><tr><td>8</td><td>4954</td><td>418m</td><td>12.87</td><td>50.79</td><td>31.5%</td><td>CHEMNITZ</td></tr><tr><td>9</td><td>52</td><td>515m</td><td>11.54</td><td>48.16</td><td>34.7%</td><td>MUENCHEN</td></tr><tr><td>10</td><td>4079</td><td>472m</td><td>12.73</td><td>48.48</td><td>32.8%</td><td>KR.ROTTAL-INN FALKENBERG</td></tr><tr><td>11</td><td>483</td><td>227m</td><td>13.76</td><td>51.13</td><td>30%</td><td>DRESDEN-KLOTZSCHE</td></tr></tbody></table>

ERA-Interim (Dee et al., 2011) large-scale reanalysis variables used in this work, defined over the 2° regular grid shown in Fig. 1 as blue crosses. As in Gutiérrez et al. (2019), we use the 4 nearest gridboxes, and thus we have a total of 68 individual predictor variables for each location.

<table><thead><tr><th>Description</th><th>Height level</th></tr></thead><tbody><tr><td>Air pressure</td><td>Sea level</td></tr><tr><td>Geopotential</td><td>500,700,850 mbar</td></tr><tr><td>Temperature</td><td>2-meter, 500,700,850 mbar</td></tr><tr><td>Eastward wind</td><td>500,700,850 mbar</td></tr><tr><td>Northward wind</td><td>500,700,850 mbar</td></tr><tr><td>Specific humidity</td><td>500,700,850 mbar</td></tr></tbody></table>

folds {$F_1, ..., F_5$}, with each fold comprising 6 years: $F_1 = 1979 - 1984, F_2 = 1985 - 1990, F_3 = 1991 - 1996, F_4 = 1997 - 2002, F_5 = 2003 - 2008$. For each of the 5 folds, predictions are obtained by training the model using the remaining 4 folds. By doing so, we end up with an out-of-sample (i.e. based on previously *unseen* data) probabilistic predictions for the daily occurrence time series which covers the whole period, 1979–2008.

Then, the so-obtained out-of-sample downscaled series are validated in terms of the metrics shown in Table 3, extracted from the VALUE validation framework. Note that from now on we will refer to the down-scaled parameters of the JPD as the *predicted probability values* and to series drawn from these probabilities as *simulations*. The measures in Table 3 allow us to assess their quality in terms of their predictive power, using the area under the ROC curve (AUC, see Kharin and Zwiers (2003)); probabilistic calibration (BrierScore, see Weigel et al. (2007)); marginal distributional performance (R01); spatial consistency, as measured by the pairwise Pearson correlations (SPCor); autocorrelation (AutoCor); spell reproducibility (WetSpellMean, DrySpellMean); and extreme spell reproducibility (WetLongestMean, DryLongestMean).

Both AUC and Brier score are computed from the predicted probabilities. In order to validate the rest of the metrics, which we refer to as distributional metrics, we draw 250 simulated series from the predicted probabilities. Since the observed series provide a single realization of each quantity, we want to verify that these observed values can plausibly be considered as drawn from the corresponding simulated distributions. We show the simulated distribution of each metric using violin plots: if the simulations are adequate then the observed values should lie inside the range of the distribution in each plot. In the case of AUC (Brier score), a value closer to 1 (0) implies better performance.

### 3. Methods

Let $Y = (Y^1, ..., Y^D)$ be a vector of binary variables representing precipitation occurrence at a set of $D$ locations, with $Y^d$ taking the values 1 or 0 according to whether location $d$ is wet or dry. We aim to predict the joint probability distribution of $Y$, accounting for the spatial

dependence between its elements. Traditional single-site models that ignore this structure fail to reproduce the correlation exhibited among the downscaled locations.

In this work we introduce a novel model, which has not been, to our knowledge, proposed before: An ensemble of conditional random forests (CRFs) able to predict the conditional probability distributions (CPDs) relevant to capture the joint probability distribution (JPD) of the complete set of locations of interest. Which CPDs are relevant is informed by a Bayesian network (BN), which ensures they approximate the full JPD as accurately as possible.

The full specification of the proposed model, which we call Bayesian Network-Informed Conditional Random Forests, is explained in Section 3.3, with Sections 3.1 and 3.2 focusing on the most relevant aspects of RFs and BNs used in this work.

#### 3.1. Conditional random forests

For a dataset $\{(y_i, x_i) : i = 1, ..., N\}$, where each $y, x$ are considered as the realized values of random variable $Y$ (the predictand) and $X$ (the predictors) respectively, a random forest (Breiman, 2001) is a machine learning model built as an ensemble of $S$ individual classification and/or regression trees (Hastie et al., 2009, Chapter 15). For each predictor value $x$, the $s$th tree defines a set $t_s(x) \subseteq \{1, 2, ..., N\}$ based on a set of decision rules, and collects together the corresponding observations into $\{y_i : i \in t_s(x)\}$, a *leaf*, based on a set of decision rules. The decision rules are constructed by means of a *split function*, which is tasked with recursively splitting the space of $X$ in order to minimize a measure of dissimilarity in $Y$, typically entropy in the context of classification. Traditionally, a *summary measure* $\hat{y}_{t_s(x)}$ is computed from $\{y_i : i \in t_s(x)\}$ for each tree $t_s$, with the random forest yielding a prediction $\hat{y}(x)$ as the average output of the $S$ trees:

$$ \hat{y}(x) = \frac{1}{S} \sum_{s=1}^{S} \hat{y}_{t_s(x)} \quad (1) $$

In a random forest individual trees are not trained with the entire training dataset, but with a randomly selected bootstrap sample. Also, only a random subset of the available predictors $X$ is made available as candidate split each time a split is to be performed, potentially allowing trees to reach better solutions. These two mechanisms provide more robust and better-performing models (Hastie et al., 2009, Chapter 9).

Legasa et al. (2022) introduced a *posteriori* random forests for statistical downscaling of climate (AP-RFs hereafter), and showed how they are capable of capturing the probability distribution of rainfall while keeping full predictive power. As opposed to the traditional approach in Eq. (1) (referred to AVG-RFs hereafter), AP-RFs take advantage of the full range of values in the leaves (which we call *predictive sample*) to estimate the parameters of the probability distribution of $Y | X = x$ (see Figure 2 in Legasa et al. (2022)). Thus, Eq. (1) is replaced with

$$ \hat{y}(x) = \Theta \left( \bigcup_{s=1}^{S} \{y_i : i \in t_s(x)\} \right), \quad (2) $$
<caption>**Table 3**</caption>

Validation metrics considered to assess the performance of the different SD methods used in this work. Note that, except for the AUC and BrierScore, we use the predicted/observed ratio for direct comparison among metrics. Excluding the Brier score, the closer to 1, the better. As for the Brier score, the lower, the better. The second column indicates how metrics are computed. Metrics measured on the prediction are computed over the downscaled probability values, whereas measures marked on simulations are computed from a binary simulation drawn from the downscaled prediction. In the latter case, 250 simulations are drawn, for a stable assessment of inter-realization variability.

<table><thead><tr><th>Code</th><th>Measured on</th><th>Description</th></tr></thead><tbody><tr><td>AUC</td><td>Prediction</td><td>Area under the ROC Curve</td></tr><tr><td>R01</td><td>Simulations</td><td>Number of wet days</td></tr><tr><td>AutoCor</td><td>Simulations</td><td>Lag-1 autocorrelation</td></tr><tr><td>SPCor</td><td>Simulations</td><td>Spatial pairwise correlation</td></tr><tr><td>WetSpellMean</td><td>Simulations</td><td>Mean duration of wet spells</td></tr><tr><td>DrySpellMean</td><td>Simulations</td><td>Mean duration of dry spells</td></tr><tr><td>WetLongestMean</td><td>Simulations</td><td>Mean of longest annual wet spells</td></tr><tr><td>DryLongestMean</td><td>Simulations</td><td>Mean of longest annual dry spells</td></tr><tr><td>BrierScore</td><td>Prediction</td><td>Σ|p - o|<sup>2</sup>, for the predicted probabilities p and observed values o</td></tr></tbody></table>

where $\Theta$ is a function which estimates the relevant parameters of the probability distribution of $Y$ | $X$. In addition, AP-RFs were shown to yield very robust results when used to downscale GCM projections for the future climate projections (Legasa et al., 2023).

In order to maximize the spatial coherence of our predictions, we propose in this work to use AP-RFs to predict, for each $X = x$, the JPD of $Y$ instead of just the marginal of each $Y^d$. For example, for two variables $Y^1, Y^2$, the random forest predicts the parameters of $P(Y^1)$ and $P(Y^2 | Y^1)$, fully determining the JPD of $Y$ for $X = x$ (thus $P(Y^2 | Y^1, X = x)$). We call these models Conditional Random Forests (CRFs).

Note, however, that this poses a hard practical limitation since, as we will see in the next section, the dimensionality of such a problem increases very rapidly, becoming intractable for relatively small values of $D$. Bayesian networks, explained in the next section, provide an efficient tool to tackle this issue by identifying conditional independence relationships among the elements of $Y$ that allow the problem to be solved using lower-dimensional conditional distributions. As we will explain in Section 3.3, BNs will allow us to build an ensemble of CRFs which can predict the relevant conditional probability distributions in $Y$, whilst keeping its spatial structure.

## 3.2. Bayesian networks

The joint probability distribution of any discrete multivariate $Y$ can be specified using the chain rule as a product of conditional distributions:

$$P(Y^1, \dots, Y^D) = P(Y^1) \prod_{d=2}^{D} P(Y^d | Y^1, \dots, Y^{d-1}). \quad (3)$$

The specification of the JPD for Bernoulli marginals using this formula requires the specification of $2^{(D-1)} + 2^{(D-2)} + \dots + 1 = 2^D - 1$ conditional probabilities, corresponding to each of the different configurations of the elements of $Y$: this quickly becomes intractable as $D$ increases. For instance, 11 stations require the specification of 2047 parameters.

A Bayesian network (Koller and Friedman, 2009; Scutari, 2010, BN) is a probabilistic graphical model which simplifies Eq. (3) by capturing conditional independencies in the dataset, thus yielding an approximation of the JPD. This approximation is determined by a directed acyclic graph (DAG, a directed graph with no cycles, see e.g. Fig. 1), and it is formalized through the concept of graphical separation (a BN is, by definition, an independence map, i.e. the graph implies conditional probabilistic independence, but not the other way around. See e.g. Scutari (2010) or Koller and Friedman (2009)). The number and configuration of arcs (directed edges) in the DAG determines the degree of approximation made by the BN, by lowering the threshold for probabilistic independence and thus placing fewer arcs. It should be noted that the direction of the arcs does not imply causality (nor dependence). We refer the reader to Koller and Friedman (2009) and Scutari

(2010) for further details on Bayesian network learning and inference algorithms and summarize here two properties of Bayesian networks that are particularly relevant for this work.

Firstly, each node $Y^d$ is conditionally independent of its non- descendants given its parents. As we will later see, this property is essential to build the ensemble of RFs proposed in this work. This is called the local Markov property. Secondly, the DAG of a BN always has an ancestral ordering of the nodes (not necessarily unique). That is, the nodes in a DAG can be arranged as $Y^1 \le \dots \le Y^D$, with $Y^i \le Y^j$ if and only if $Parents(Y^j) \subseteq \{Y^1, \dots, Y^i\}$, for all $i \ne j$. There is thus always at least one starting node, i.e. a node with no parents, $Y^1$. Note how the nodes in Fig. 1 have already been labeled following an ancestral order. As we will later see (Section 3.3), we use this property to simulate from the proposed model.

An immediate corollary can be extracted from these properties: the JPD described in Eq. (3) can be approximated by a Bayesian network as $P_{BN}$, with

$$P(Y^1, \dots, Y^D) \approx P_{BN}(Y^1, \dots, Y^D) = \prod_{d=1}^{D} P(Y^d | Parents_{BN}(Y^d)). \quad (4)$$

This simplification, together with the properties explained above, effectively allows us to divide the problem of predicting the JPD of $Y$ | $X$ into a set of $D$ tractable subproblems, as we explain in the next section. Each of the $D$ terms on the right hand side of Eq. (4) is called a conditional probability distribution (CPD).

To automatically learn the DAG from the observed occurrence precipitation records, score-based algorithms rely on the maximization of the Bayesian information criterion (BIC) over the space of all possible DAGs. BIC equals the log-likelihood plus a penalization term $k = \log(N)/2$ (where $N$ is the number of training instances) multiplied by the number of parameters of all the CPDs. Note that, as we will see in Section 4.3, even though BIC is strictly defined with $k = \log(N)/2$, we can fine-tune the density (number of arcs) of the DAG by varying $k$. Increasing $k$ results in sparser DAGs, whereas reducing it results in denser DAGs.

Since the number of possible DAGs that can be built for $D$ nodes is $2^{O(D^2)}$, heuristic algorithms have to be employed. Legasa and Gutiérrez (2020) showed how the Tabu search (Glover, 1990; Koller and Friedman, 2009) algorithm reliably captures the spatio-temporal links for precipitation occurrence, yielding robust approximations of the JPD for several locations. Tabu search is a standard heuristic search procedure for maximizing the BIC in the discrete space of possible DAGs, by essentially running several hill-climbing searches with a memory that stores a backtrack of already visited DAGs and performing additional searches prohibiting previously visited local maxima, performing permutations of the variables (note that we labeled locations according to the ancestral order, but that they were not input in that order) together with random restarts. It yields generally better results than the greedy hill-climbing. We thus employ Tabu search to learn the BN structures used in this work.
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** A directed acyclic graph (DAG) corresponding to a Bayesian network capturing the relevant independence structure for precipitation occurrence at the 11 stations described in Table 1. The DAG was learnt using Tabu search and BIC score (see Legasa and Gutiérrez (2020) for the details), and the first fold of the cross-validation procedure (check Section 2). The maximum number of parents has been set to 3. Note that the codes for the stations correspond to an ancestral ordering, given by the algorithm. We use these codes throughout all the article. ERA-Interim reanalysis grid cell centers (Dee et al., 2011) are shown as blue crosses. The list of predictors is shown in Table 2. Orography data was obtained from Yamazaki et al. (2017).</figcaption></figure>

### 3.3. Bayesian network-informed conditional random forests (BNICRFs)

The model we propose in this work aims to predict **Y** | **X** for a set of locations **Y** = (Y¹, ..., Y^d) while preserving the dependence structure in **Y**. Since, as explained before, the specification of the full JPD is intractable, the core idea is to use an ensemble of conditional random forests to link all the relevant conditional probability distributions within **Y** to the large scale predictors **X**, with a Bayesian network informing which CPDs are relevant. That is, thanks to the local Markov property and the factorization in Eq. (4) produced by the BN, we can use a CRF to predict the CPD of Y^d | Parents<sub>BN</sub>(Y^d) from the large scale predictors **X** for each location **d**. Once trained, these CPDs fully characterize **Y** | **X** and allow us to simulate values from the JPD of the predicted **Y**. We call this ensemble of CRFs a Bayesian network-informed conditional random forest (BNICRF).

In the following, we will explain how to train a BNICRF, using the DAG in Fig. 1 for illustration. The first step is to build an informative Bayesian network which approximates all the independence relationships of **Y**. As explained in Section 3.2, for each location Y^d the BN prescribes a CPD<sub>d</sub> = P(Y^d | Parents<sub>BN</sub>(Y^d)). For the DAG in Fig. 1 we thus have 11 CPDs. After the informative BN is built, an ensemble of D CRFs is trained to predict each CPD from each large-scale predictors' state **X** = **x**. Each CRF uses as predictors the 4 nearest gridboxes to each location including the parent locations, **X**, and {Y^d, Parents<sub>BN</sub>(Y^d)} as predictands, and is thus trained to reduce the entropy of the JPD of both the location of interest and its parents. In this regard, note that

BNICRFs assume that the relationships learnt by the BN continue to hold when we incorporate the large-scale predictors **X**, using the CRFs to predict P(Y^d | Parents<sub>BN</sub>(Y^d)) from each **X** = **x**, which we denote P(Y^d | Parents<sub>BN</sub>(Y^d), **X** = **x**). As we will later see in Section 4, for the case study considered here this is a reasonable assumption.

The first location, Y¹, has no parents, and thus the first CRF is actually a single-response RF model for Y¹, which is tasked with just predicting the probability distribution P(Y¹ | **X** = **x**). We then proceed to build the CRF for Y², which has as single parent Y¹, and thus this CRF learns the CPD₂(**X** = **x**) = P(Y² | Y¹, **X** = **x**). This process is iterated for each node/variable in **Y**, with each CRF using information from the child Y^d and its parents to predict the CPD<sub>d</sub>(**X** = **x**) = P(Y^d | Parents<sub>BN</sub>(Y^d), **X** = **x**) for each large-scale predictor configuration **X** = **x**. For example, the location Y¹¹, the last variable in the DAG in Fig. 1, requires the specification of CPD₁₁(**X** = **x**) = P(Y¹¹ | Y², Y⁷, Y⁸, **X** = **x**) for each **X** = **x**. The chart in Fig. 2 visually summarizes the training process for 5 locations and a simple illustrative DAG.

Since the BN is an independence map (i.e. graphical separation implies independence, but not the other way around), it allows the CRF flexibility on the dependence relationships of **Y**. For example, a configuration in the form Y¹ → Y² implies that Y¹ and Y² are not assumed to be independent. However, the CRF can still decide the degree of dependence for each particular **X** = **x**, by setting P(Y¹ | Y¹, **X** = **x**) = P(Y² | **X** = **x**) to make Y² and Y¹ independent provided **X** = **x**. In this sense, the BN only forces relationships between different variables to be
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Diagram illustrating the BNICRF methodology. Once a Bayesian network (left) has been trained to capture the (in)dependence structure for the locations of interest *Y* (for simplicity, only 5 locations have been considered in this figure), a conditional random forest (CRF) is trained for each location. The aim of this CRF is to link the large scale predictors *X* to the target location and its parents informed by the Bayesian network. Using the a posteriori approach described in Legasa et al. (2022), each CRF predicts, for each *X* = *x*, the parameters of *Y*<sup>*d*</sup> | Parents(*Y*<sup>*d*</sup>), ending up with a daily JPD for each *Y* | *X* = *x*.</figcaption></figure>

approximated by others, but does not force a relationship to exist. The
precision and complexity of this approximation, in turn, can be decided
by controlling the complexity of the BN (as we will show in Fig. 6 in
Section 4.3).

Additionally, in the BNICRF formulation, situations can arise in
which the predictive sample (see Section 3.1) does not provide infor-
mation about some conditional distributions. For example, for node Y¹¹
and a particular predictor configuration X = x, the CRF may produce a

predictive sample with no data about a combination of its parents' state
(e.g. P(Y¹¹ | Y² = 1, Y⁷ = 0, Y⁸ = 1, X = x)). In these cases the locally
marginal probability of the child variable P(Y¹¹ | X = x) is used as an
approximation, and thus in this case P(Y¹¹ | Y² = 1, Y⁷ = 0, Y⁸ = 1, X =
x) = P(Y¹¹ | X = x). Since the CRF only uses the large-scale predictors,
we can interpret this situation as the CRF prescribing a conditional
independence, on top of the BN, given the large-scale predictors X.
This approximation has to be done in at least one configuration for
<figure><img src="image_5.png" /><figcaption>Fig. 3. Cross-validated results for the measures AUC, R01 and SPCor (see Table 3) obtained for different values of the minimum leaf size and number of trees parameters (axis x and y, respectively). R01 and SPCor correspond to the averaged ratio between 250 simulated and observed measures, whereas AUC is computed from the predicted probabilities. We show the average of the 11 stations for the AUC and R01 and 55 correlation pairs for SPCor. In all cases, the RFs were built using joint entropy as split function and the same informative BN configuration (with 3 parents maximum).</figcaption></figure>

approximately 10% of the training days for the reference model, which is explained and analyzed in Section 4. Note that, if the ensemble of CRFs is well calibrated (i.e., there is no overfitting), the probability of drawing that particular configuration of the parents will be virtually 0, as those configurations are very unlikely to occur on those days.

Once the ensemble of CRFs has been constructed, it can be used to make either deterministic predictions, using Eq. (4); or stochastic simulations for $Y | X$, simulating by iteratively drawing samples from the learnt CPDs following in the same ancestral order determined by the DAG, which guarantees we have already simulated parent nodes. We start, for each large-scale predictors state $X = x$, by simulating from the marginal distribution of $Y^1$, obtaining $Y^1 = y^1$. We then simulate from the probability distribution $P(Y^2 | Y^1 = y^1, X = x)$. We proceed iteratively by drawing samples from the CPDs until we obtain a complete simulated vector $Y$ of values at all locations.

The last issue to consider here is the split function, which, as mentioned in Section 3.1, is tasked with finding the optimal set of split points in the space of $X$ in order to reduce the entropy for each $P(Y^d | Parents_{BN}(Y^d))$. In this work we use, for each CRF, the joint information gain, i.e. the reduction of entropy of $(Y^d, Parents(Y^d))$. The formula for the entropy used by each CRF is, thus,

$$H(Y^d, Parents(Y^d)) = \sum_{y \in (Y^d, Parents(Y^d))} -P(y) \log(P(y)) \quad (5)$$

with the split function performing the splits that obtain the highest reduction of entropy. For example, the last CRF, which predicts $P(Y^{11} | Y^2, Y^7, Y^8)$, is trained to reduce the entropy of $(Y^2, Y^7, Y^8, Y^{11})$. Experiments carried out considering the entropy of the child (e.g. $Y^{11}$ for the last CRF) yielded worse results, particularly in terms of spatial consistency. The decision on the best parameters for the CRF ensemble, i.e. the number of trees and maximum depth of each tree, is tackled in Section 4.1.

## 4. Results

All the results shown in this section correspond to out-of-sample predictions/simulations, i.e. to the complete time-series downscaled by the models covering the period 1979–2008, obtained according to the 5-fold cross-validation scheme explained in Section 2. Recall that, in

order to assess the performance of the different models considered, we use the measures in Table 3, computed for the deterministic series (using the predicted probabilities) in the case of AUC and Brier score and on 250 stochastic simulations for the rest of the measures.

Note that in Sections 4.1–4.3 we aim to find the optimal model configuration, and thus, for brevity, we focus only on the predictive performance (AUC), the number of wet days (R01) and the spatial pairwise correlation (SPCor), which provide an overview of model performance. By contrast, in Sections 4.4 and 4.5 we use all the measures to provide a comprehensive assessment of the performance of the methodology.

### 4.1. Sensitivity to the random forest configuration

In this section we aim to find the optimal combination of parameters for the conditional random forest models. As explained in Section 3.1, two main parameters control the overall complexity of a RF: the number of trees and the minimum leaf size. Fig. 3 shows the cross-validated results obtained for the measures in Table 3. Ideally, we aim for a combination of parameters that attain the maximum AUC (i.e. the maximum predictive performance) while keeping a good distributional performance, as measured by the ratio in the R01 and the SPCor. For this analysis, the values of R01 and SPCor are summarized by averaging the ratios of simulated to observed measures: the variation between simulations (not shown) is relatively small so that this simple summary is considered adequate for the purpose of comparing parameter settings.

We conclude from Fig. 3 that there is a clear tradeoff between the predictive and distributional performance. Models with lower minimum leaf size values tend to obtain greater predictive performance (i.e. higher AUCs, but also slightly bias the distributional measures). From a minimum leaf size of 50 and upward the ensembles show no bias on distributional metrics but start gradually losing predictive performance. Therefore, we chose a minimum leaf size of 15 as a good compromise. In terms of the number of trees, models with just 10 trees tend to perform badly at both the distributional and predictive performance, with the difference being slim for more than 25 trees. We therefore choose CRFs with 100 trees and 15 minimum leaf size as the optimal configuration, and we call this model the reference BNICRF hereafter.
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Comparison of 3 multisite BNICRF models against 3 single-site models, for the measures AUC (panel a), R01 (panel b) and SPCor (panel c) (see Table 3). For the strip charts in panels a and b each dot represents each location. In panel c, each boxplot encompasses the correlation pairs. In all cases, measures have been averaged from 250 simulations. For all multisite models we use the same informative Bayesian network configuration. The BNICRF models correspond to the proposed approach, using a minimum leaf size 15 (red, AP-RF 15) and 50 trees (orange, AP-RF 50). We also include a BNICRF using the AVG-RF approach for the CRFs (brown, AVG-RF, check Legasa et al. (2022)). The single-site models correspond to a random forest using the AP-RF approach (cyan), a random forest using the AVG-RF approach (aquamarine), and the GLM described in Legasa et al. (2022) (blue).</figcaption></figure>

## 4.2. BNICRFs vs single-site models

Since capturing the spatial structure of Y requires additional model complexity with respect to predicting marginally at each site $Y^d$, a relevant question is whether our model introduces additional bias in the predictions for each individual station. This section addresses the trade-off between single-site predictive performance and multisite/spatial performance.

To answer this question we compare our proposed model against three single-site models: a generalized linear model (GLM) described in Legasa et al. (2022), a well-known robust downscaling technique. Specifically, this is a simple application of the approach, that was found to perform competitively compared with other methods (Gutiérrez et al., 2019); and 2 single-site random forests, one using the averaging approach (AVG-RF, standard classification random forest) and the other one using the AP-RF approach (Legasa et al., 2022), which employs the same configuration (the best performing models are employed). Results are shown in Fig. 4.

For the multisite BNICRFs, in addition to the reference model from Section 4.1 (red, AP-RF 15), we include for completeness two additional models. The first one corresponds to a random forest using the proposed methodology but with minimum leaf size of 50 (orange), which, as explained in the previous section, better captures the parameters of the

distribution (both multisite and singlesite), but slightly loses predictive performance. The second one uses the same parameters as the reference model but instead employs averaging (brown) instead of the AP-RF approach. All models use the default informative BN configuration (see Section 4.3).

For the spatial models, we conclude that using the AP-RF approach (red, orange) is essential to capture the relevant spatial relationships. Instead, traditional AVG-RF approach (brown) fails to reproduce pairwise spatial correlation (SPCor) and performs worse reproducing the number of wet days for individual stations (R01). This supports the results in Legasa et al. (2022), where it was concluded that the AP-RF is essential to capture the full probability distribution. In addition, when compared against single-site models (cyan, green, blue) which, as expected, fail to reproduce the relevant spatial relationships, both multiresponse BNICRFs (red, orange) slightly lose either some predictive performance (orange, AUC) or some distributional similarity (red, R01). Overall, the minor reduction in performance is due to BNICRFs having to estimate more parameters for the conditional probability distributions as the number of parents grows (see Section 4.3). Finally we also find a tradeoff for the two BNICRFs: as concluded in the previous section, those which better capture the distributional aspects tend to lose some predictive performance. This difference is small, with the median AUC falling from ~ 0.86 for the reference BNICRF to ~ 0.845
<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Simulated pairwise correlation values for the proposed Reference BNICRF (red) against a single-site GLM (blue, each $Y^p$ predicted independently), for the period and stations shown in Table 1.</figcaption></figure>

for the BNICRF with 50 as minimum leaf size, but this is a symptom of some potential underfitting.

In addition to the SPCor measure, we show in Fig. 5 all the observed and simulated pairwise correlation values for the single-site GLM and the BNICRF. It can be clearly seen how the proposed model accurately captures the spatial dependence structure. However, and although the predictors (recall that the 4 nearest gridpoints to the locations of interest are used) do carry some spatial information, the single-site GLMs are unable to capture the full dependence, greatly underestimating the correlation between locations. This effect is particularly strong for the higher correlations: the maximum observed correlation is 0.75, well simulated as 0.761 by the BNICRF but the correlation falls to 0.46 by the GLM. Instead, the BNICRF tends to slightly underestimate the lower – and thus less relevant – correlation values.

Following from these results, from now on we stick to the reference BNICRF, which provides a tradeoff between predictive performance and distributional similarity.

### 4.3. Sensitivity to the Bayesian network configuration

As explained in Section 3.2, structure learning algorithms maximize the log-likelihood in combination with a complexity parameter $k$ multiplied by the number of parameters of the model. For the BIC score, $k$ is set to $k = \log(N)/2$ (where $N$ is the number of training instances), but this parameter can be decreased/increased for denser/sparser DAGs. Although we found in Legasa and Gutiérrez (2020) the parameter $k = 1$ to be the best choice for modeling the dependence among the stations used in this study, the CRFs introduce an additional layer of complexity,

and thus we reassess in this section the choice of the Bayesian network complexity parameter $k$.

Fig. 6a shows the pairwise correlation for different values of $k$. We stick to the CRF parameters selected in the previous section (100 trees and minimum leaf size 15). Note that an additional way to control the complexity is by limiting the number of parents that each node in the DAG can have, and thus we show, in Fig. 6b, the same analysis for different limitations on the number of parents. Since no significant difference is found in single-site performance by varying the Bayesian network complexity, we focus this section on the spatial (pairwise correlation) aspect.

From Fig. 6a, we decide to stick to the default penalization value corresponding to the BIC score: $k = 4.64 \approx \log(N)/2$, where $N$ is the number of training instances in the complete dataset. Decreasing the penalization results in too dense DAGs that could potentially bias the final results and require the CRFs having insufficient information to estimate that many CPDs. It should be noted that using sparser DAGs, e.g. with $k = 25$ (with 23 arcs for this example) may be an option in certain applications, such as when having computational constraints, as the difference in spatial performance is relatively small. We also conclude that the improvement when allowing 4 parents or more (no more than 4 parents are used for the default $k$) is marginal, and thus it does not compensate for the additional computational cost and overfit potential.

Therefore we fully determine the specification of the BNICRF model for the final complete experiment: an ensemble of CRFs with 100 trees and minimum leaf size 15 each, informed by a BN learnt using Tabu search with BIC score and 3 parents maximum.
<figure><img src="image_8.png" /><figcaption>Fig. 6. Simulated pairwise correlation values for different values of the parameter k (panel a) and maximum allowed number of parents using the default k = log(N)/2 (panel b). In all cases the CRFs have 100 trees and minimum leaf size 15. Topleft legend shows the number of arcs and the log likelihood (divided by the number of days i.e. training instances) of the DAGs.</figcaption></figure>

## 4.4. Validation vs a multisite GLM

So far, different variants of the BNICRF model have been compared with each other and with simple single-site downscaling procedures to illustrate the effects of the various implementation choices. In this section we present a tougher test, comparing the BNICRF performance with that of an existing multisite downscaling approach. This approach is based on multisite GLMs as implemented in the Rglimclim software package (Chandler, 2020). Like the single-site GLMs considered above, Rglimclim handles precipitation occurrence using logistic regression models, but includes a variety of options for inducing inter-site dependence in the simulated time series of daily precipitation occurrence. This methodology has been chosen for comparison against a BNICRF because it has been shown to perform competitively by comparison with other multisite downscaling techniques (e.g., Frost et al. 2011, Vaittinada Ayar et al. 2016).

Although the 5-fold cross-validation described in Section 2 was used to develop and train and test the Rglimclim models, a key difference between the two approaches is that a BNICRF is trained automatically whereas Rglimclim models must be developed manually, using diagnostics and expert judgement to identify an appropriate model structure. Due to the element of judgement involved, it is hard to develop completely independent Rglimclim models for each fold. Acknowledging this feature of statistical downscaling models in general, Maraun et al. (2015) suggested that a pragmatic approach to cross-validation for such models is to identify the model structure using the complete training data set and then, conditional on this structure, to estimate the model parameters separately for each fold. Although this suggestion potentially saves time (it typically takes an experienced Rglimclim user two or three hours to develop a precipitation occurrence model for a multisite data set of the size considered here), we have carried out both the model identification and estimation separately for each fold to ensure that the Rglimclim and BNICRF modeling results are as comparable as possible. To minimize the possibility that the chosen Rglimclim model structure from one fold will influence the decisions made for subsequent folds, a structured model development protocol was followed throughout: this is described below.

The manual nature of the Rglimclim GLM development also has implications for the number of predictors that can be considered. In particular, it is not feasible to consider all of the 68 atmospheric predictors listed in Table 2 (17 separate predictors at each of 4 grid

nodes) for each location. The predictor data were simplified for the Rglimclim models therefore, by averaging each predictor over all of the 12 ERA-Interim grid nodes shown in Fig. 1. This approach is standard in other applications of GLMs for downscaling.

In the logistic regression models used in Rglimclim, the probabilities of rain at each location are linked explicitly to the values of multiple covariates. Covariates include the original variables described in Section 2, and additional ones representing the effects of seasonality, topography, temporal dependence, large-scale atmospheric structure and interactions. This contrasts with the BNICRF approach where the effects of, for example, seasonality and topography are not represented explicitly but, instead, are reflected in the (automatically) chosen predictor configurations. For covariates $X_1^d, ..., X_r^d$, the logistic regression model specifies the probability of a wet day at location d as

$$ \log \left( \frac{P(Y^d = 1)}{1 - P(Y^d = 1)} \right) = \beta_0 + \beta_1 X_1^d + \dots + \beta_r X_r^d, \quad (6) $$

where the coefficient vector $\beta = (\beta_0 \ \beta_1 \ \dots \ \beta_r)^T$ is estimated by maximizing a log-likelihood defined as though all locations are conditionally independent given the covariates. Standard errors and confidence intervals for the coefficients, together with hypothesis tests allowing the formal comparison of nested models, are adjusted for potential dependence between locations (Chandler et al., 2006; Chandler and Bate, 2007; Chandler, 2020). When simulating from the fitted models, several mechanisms are available for inducing realistic levels of inter-site dependence: that used here is based on a correlated latent Gaussian field as introduced by Ambrosino et al. (2014), which is appropriate in situations where (as confirmed during the model development for each fold) there is a clear decay of inter-site dependence with distance.

The protocol for model development was, first, to define a candidate set of covariates for each primary source of variation (seasonality, topography and so forth); then, for each fold, to consider each source of variation in turn and examine alternative representations of its effects from among the candidates. For example, the candidates to represent seasonality were $\cos(2\pi kt/365.25)$ and $\sin(2\pi kt/365.25)$ for $k = 1, 2, 3$ where t denotes time in days; those representing temporal dependence were various functions of rainfall occurrence at all locations on each of the four days preceding the time point of interest; those representing topographic variation included a variety of measures such as altitude and aspect as in Ambrosino et al. (2014); and those representing large-scale atmospheric structure were a subset of the spatially averaged ERA-Interim predictors (several predictors with highly significant
BNICRFs for Probabilistic Multisite Downscaling of Precipitation Occurrence

<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Comparison of the proposed model (Reference BNICRF, using minimum leaf size 15 and 100 trees), and the Rglimclim GLM, for the measures listed in Table 3. Panel a shows the AUC, BrierScore, as well as the average (for all the simulations) distributional performance metrics for each location. For the distributional performance metrics, we show the relative bias (%), with stripcharts showing each location and boxplots encompassing station pairs. Panel b shows the distributional measures (columns) for each location (rows), including the variability of the models: each violin plot encompasses 250 simulations, with dashed lines displaying the observed value.</figcaption></figure>

time trends were excluded from consideration, because there were no
significant time trends in the precipitation occurrences). Throughout,
modeling decisions were based on a combination of hypothesis tests,
diagnostic plots and monitoring of standard errors as protection against
overfitting (Chandler, 2020, Section 4). The full process is documented

in the notebook accompanying this paper (Rglimclim Multi-Site GLM -
Model Development and Reproducibility of Results.ipynb, located in https:
//github.com/MNLR/BNICRF).

Under this protocol, most modeling decisions were clear-cut and led
to similar model structures across all folds. In the limited number of
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Comparison between the observed and simulated pairwise correlation values for the proposed Reference BNICRF (red) and a multisite GLM (Blue) constructed using Rglimclim, for the 11 stations in Table 1.</figcaption></figure>

exceptions, the final decision was typically adopted the more complex of the alternatives under consideration. The justification is that the increased complexity typically involved one or two additional parameters and, with more than 90 000 observations in total, the performance of a model is unlikely to be affected adversely overall by a handful of potentially redundant terms. The final models selected for each of the five folds have 65, 65, 71, 66 and 65 coefficients respectively.

The results of the comparison between the BNICRF and Rglimclim are shown in Figs. 7 and 8. Fig. 7 shows that the BNICRF performs better than the Rglimclim model on measures that describe predictive performance: thus it has a higher median AUC (0.87, compared with around 0.81 for the Rglimclim model) and a lower Brier score. Overall, both models accurately reproduce the proportion of wet days (thus the marginal Bernoulli probability distribution) and they also capture the inter-site dependence structure accurately as represented by pairwise correlations between locations (see Fig. 8). By contrast, Rglimclim slightly outperforms the BNICRF capturing autocorrelation (AutoCor) for some locations, which translates into slightly better reproduction of the spell measures (WetSpellMean, DrySpellMean, WetLongestMean and DryLongestMean).

Note that we have also carried out this same analysis separately for each season for the BNICRF, by dividing the simulations into the four boreal seasons and computing the measures in Table 3, excluding the

spells for simplicity. We show the results in Figure S1 (supplementary information file). We observe that, even though the model has not been trained separately for each season, BNICRFs adequately represent seasonality, with no significant differences found in performance across seasons: the maximum difference is found between summer (JJA, red) and winter (DJF, cyan), with AUC falling from around 0.89 in winter to around 0.84 in summer. This effect is also seen for the Brier score, and coincides with results found in Legasa et al. (2022), in which AP-RFs were shown to perform worse in summer, probably due to the locations of interest having significantly less wet days in summer than in winter. In addition, we also observe this for the Rglimclim GLM, hinting at a systematic, rather than model-dependent, higher bias in summer.

These results may be summarized as: both BNICRFs and multisite Rglimclim GLMs perform well overall in the distributional measures, but BNICRFs exhibit a better predictive performance. This is due to its greater flexibility to represent nonlinear predictor–predictand relationships. The representation of autocorrelation and spell lengths in the Rglimclim model follows from the explicit inclusion of previous days’ rainfall occurrences in the model structure. In the next section, we show how a simple extension of BNICRFs can vastly improve the reproduction of all temporal aspects.
<figure><img src="image_11.png" /><figcaption>**Fig. 9.** As Fig. 7, for the T-BNICRF (red), that takes into account information from each location in the previous day (purple), by training the CRFs to predict $Y_i^d | \text{Parents}_{BN,i}(Y^d), Y_{i-1}^d$, with $Y_{i-1}^d$ and $Y_i^d$ representing absence/presence of precipitation at location $Y^d$ on the previous and present day, respectively. $\text{Parents}_{BN,i}(Y^d)$ represents the parents of $Y^d$ in the DAG for the present day.</figcaption></figure>

### 4.5. Accounting for the temporal component in BNICRFs

Up to this point, this work has been devoted to building a model capable of capturing the multisite dependence structure at the down-scaled locations. With the BNICRF relying only on the large-scale predictors in Table 2 and no explicit term accounting for the temporal dependence structure, it is outperformed by the Rglimclim GLM in

terms of the autocorrelation of the series, underestimating both the wet and dry day spell length, as can be seen in Fig. 7.

A straightforward way to fix this is to use CRFs to learn the relationship of each location with its past in the same way we learn the inter-site relationships. This approach was proposed in Legasa et al. (2023) and was shown to yield very good results. Let us denote $Y_i^d$ and $Y_{i-1}^d$ absence/presence of precipitation at location $i$ on the present and
previous day, respectively, and Parents<sub>BN,t</sub>(Y<sup>d</sup>) the parents determined by the DAG in the present day. We propose the CRFs to predict Y<sup>d</sup><sub>t</sub> | Parents<sub>BN,t</sub>(Y<sup>d</sup>), Y<sup>d</sup><sub>t-1</sub>, which is roughly equivalent to adding 11 extra nodes (the 11 locations considered in the previous days) in the BN and forcing an arc in the direction Y<sup>d</sup><sub>t-1</sub> → Y<sup>d</sup><sub>t</sub>. We call this model T-BNICRF.

Fig. 9 shows the results for T-BNICRF using the reference configuration (100 trees and minimum leaf size 15) in purple, put in context with the Rglimclim GLM of the previous section. With no loss of predictive performance (as measured by the AUC and the Brier score) or spatial consistency (SPCor), a clear improvement can be seen in all measures dealing with temporal aspects: the T-BNICRF model is able to properly capture the lag-1 autocorrelations (Autocor) for all locations. The new models is also able to produce realistic spell lengths, both dry and wet (WetSpellMean, DrySpellMean, WetLongestMean and DryLongestMean). We conclude that the T-BNICRF model performs very well in all metrics and is able to keep predictive capability, leading to comparatively better performance than the Rglimclim benchmark.

## 5. Conclusions

We introduced in this work Bayesian network-Informed Conditional Random Forests (BNICRFs), a novel machine learning model capable of accurately downscaling precipitation occurrence at several locations from large-scale predictors whilst preserving the multisite dependence structure. With no loss of predictive performance with respect to well-performing single-site downscaling models (two single-site random forests and a GLM), our BNICRFs provide a full and explicit description of the downscaled joint probability distribution of precipitation occurrence for each individual day. We use this description to provide spatially consistent downscaled precipitation occurrence fields at 11 locations in Germany, and thoroughly assess them using several diagnostic indices used in the COST action VALUE validation framework. We also put BNICRFs in context with other state-of-the-art single-site and multisite downscaling methodologies.

The multisite model used as benchmark corresponds to a GLM implemented in the Rglimclim software package, which has been shown to perform competitively with respect to other multisite downscaling techniques. With similar performance capturing the multisite structure, we observe that BNICRF exhibits better predictive performance, probably due to its greater flexibility to represent nonlinear relationships. This highlights how the BNICRF framework enables us to explore the relationship between precipitation and daily atmospheric predictors in a lot of detail, albeit focusing on predictive performance rather than on interpretation since the BNICRF does not provide a compact set of interpretable coefficients analogous to the vector $\beta$ in Eq. (6). For example, it would not be feasible to build a Rglimclim GLM using the individual grid-cell values of all the predictors, so we just spatially-averaged them, sacrificing thus some potentially relevant information in the predictor field. This is partly due to the difficulty in incorporating the individual predictors instead of their spatial averages, a limitation of Rglimclim with respect to BNICRFs. Nevertheless, we consider that this work has highlighted the difference between the two approaches in a very helpful way, a comparison which has not been properly attempted before in the downscaling literature.

Since the predictor fields chosen in this work alone do not carry enough information to capture the temporal structure of daily precipitation, in Section 4.5 we include lagged precipitation values to reproduce the temporal dependence structure with the BNICRF, leading to the T-BNICRF model. T-BNICRFs are able to produce spatially and temporarily consistent simulations of precipitation occurrence fields, that account for autocorrelation and realistically simulate both wet and dry spells, even reproducing extreme spells, and performing better than the Rglimclim GLM used as benchmark.

T-BNICRFs also highlight how the methodology can incorporate expert knowledge into the model easily: in this work we forced the model

to include lagged information, but other relationships are certainly an option. In this respect, another alternative could be considering as nodes of the BN the set {$Y_t^1, ..., Y_t^D, Y_{t-1}^1, ..., Y_{t-1}^D$} and letting the structure learning algorithm to learn all the arcs. This could potentially help the proposed model to capture spatio-temporal aspects like e.g. day-1 lagged correlations (Legasa and Gutiérrez, 2020). A thorough study should be carried out to assess the added value and potential downsides of this approach.

A clear advantage of the BNICRF methodology introduced in this paper is that there is no need to select predictors. This avoids the task of manually selecting the most adequate large-scale variables, which currently relies on human expertise, is very time-consuming (see e.g., Manzanas et al. (2020) and the effort devoted to the Rglimclim GLMs in Section 4.4), and potentially subjective. Also, our BNICRF selects the local predictor information, also avoiding the need to manually search the optimal geographical domain of interest (Gutiérrez et al., 2013).

Another relevant advantage of BNICRFs is the fact that they do not impose the observed correlation structure over the downscaled series. In this sense, the BNICRF model only forces some relationships to be approximated by others, while still letting the model learn the strength of these relationships from the large-scale predictors. This implies that the multisite structure can also change according to the predictors' state. In addition, this new approach does not require adding an additional penalization term in the cost or split function, nor any correction done on the final downscaled probabilities to force the observed correlations into the downscaled values, like in Cannon (2008). All this implies that BNICRFs can potentially provide adaptability to the potential variations in the multisite structure that may occur under climate change conditions. The adaptability to climate change conditions of random forests was already proven promising in Legasa et al. (2023). The suitability of BNICRFs to downscale climate change projections is something we are exploring on a follow-up article.

In terms of computational complexity, BNICRFs are more time-consuming than both single and multisite GLMs. Yet, for the configuration used in this work (100 trees and minimum leaf size 15), a BNICRF for 11 locations can be trained and used to produce 250 simulations in around 30 min using an Intel® Core™ i5-8500 @ 3.00 GHz Processor. Note that half of the computing time is spent on building the informative Bayesian network using Tabu search, and that we can reduce the complexity of Tabu search if required by time-constraints by trading off some spatial performance. In this regard, and to show the scalability of the approach, we have also tested the methodology on the complete set of 53 ECA&D stations, and included a notebook showing the process (check the notebook in *Illustrative Example for 53 Locations for Bayesian Network-Informed Conditional Random Forests.ipynb* in https://github.com/MNLR/BNICRF). For 53 locations, the whole process (a 5-fold validation procedure including training the BNICRF, predicting the probabilities and simulating out-of-sample) can be undertaken in less than 5 h in a personal computer, obtaining comparable results. Additional details are given in the notebook. We conclude, then, that GLMs are quicker to run than BNICRFs, but at the cost of higher expert time, which has to be devoted to training and calibration.

Finally, a relevant open question in which we are currently working on is the extension of the approach to continuous and hybrid (semi-continuous) probability distributions. In a follow-up paper we plan to extend this to the full distribution of precipitation, i.e. to predict both wet/dry days and rainfall on wet days. There are several issues to consider here, since Bayesian networks have several limitations on the probability distribution that can be used on the nodes. The most straightforward way to extend the methodology is by performing a sensible discretization of the variables of interest for the informative BN. The extension of discrete Bayesian networks to continuous and semi-continuous distributions is still an active field of research. In particular, Gaussian Bayesian networks only allow for normally distributed nodes
and are limited to linear relationships among them. There are exten-
sions proposed in the literature, including non-parametric (Friedman
and Nachman, 2000; Schwaighofer et al., 2007) and semi-parametric
Bayesian networks (Atienza et al., 2021). An approach we are currently
exploring is Copula Bayesian networks (Elidan, 2010), which provide
a promising extension to discrete BNs to model continuous multivari-
ate distributions via copulas utilizing the independence assumptions
encoded by a DAG.

CRediT authorship contribution statement

**M.N. Legasa:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **R.E. Chandler:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Formal analysis, Data curation, Conceptualization. **R. Manzanas:** Writing – review & editing, Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization.

Software and data availability

The software required is listed next:

1. For the Bayesian networks:

• Name of software: bnlearn

• Developers: Marco Scutari

• Contact: marco.scutari@supsi.ch

• Software required: None

• Program language: R

• Source code at: https://www.bnlearn.com

* Documentation: Detailed documentation for application installation, testing, and deployment can be found at https: //www.bnlearn.com/documentation/

2. For the random forests:

• Name of software: RandomForestDist

• Developers: Mikel N. Legasa

• Contact: mikel.legasa@lsce.ipsl.fr

• Software required: None

• Program language: R

• Source code at: https://github.com/MNLR/RandomForestDist

* Documentation: Detailed documentation for application installation, testing, and deployment can be found at http s://github.com/MNLR/RandomForestDist/blob/master/W orkedExample.ipynb and Legasa et al. (2022).

3. For the multisite generalized linear models:

• Name of software: Rglimclim

• Developers: Richard E. Chandler

• Contact: r.chandler@ucl.ac.uk

• Software required: None

• Program language: R

• Source code at: https://github.com/Richard-Chandler/Rglimclim

* Documentation: Detailed documentation for application installation, testing, and deployment can be found at https://github.com/Richard-Chandler/Rglimclim and Chandler (2020).

In addition, all the data used in this work, both observed rainfall
from ECA&D and ERA-Interim predictors, are publicly available and
can be downloaded from http://www.value-cost.eu/data. For repro-
ducibility of the results, we have created a Github directory (https:
//github.com/MNLR/BNICRF), containing the functions and data re-
quired to build the BNICRF and Rglimclim models used in this work.
The interested reader can find two Jupyter notebooks, one for repro-
ducing the results of the BNICRF model (BNICRF: *Worked Example and
Reproducibility of Results.ipynb*), and the other one explaining in detail
the Rglimclim GLM model development process (Rglimclim Multi-Site
GLM - Model Development and Reproducibility of Results.ipynb).

Declaration of competing interest

The authors declare no conflict of interest.

Acknowledgments

This study is part of the R&D project "Eventos extremos compuestos para la evaluación de los impactos del cambio climático en la agricultura" (COMPOUND: TED2021-131334A-100) funded by MCIN/AEI/10.13039/501100011033 and by the European Union NextGenerationEU/PRTR. It has also received funding from Agence Nationale de la Recherche - France 2030 as part of the PEPR TRACCS programme under grant numbers ANR-22-EXTR-0005 (EXTENDING project) and ANR-22-EXTR-0011 (LOCALISING project).

**Appendix A. Supplementary data**

Supplementary material related to this article can be found online
at https://doi.org/10.1016/j.envsoft.2026.106907.

Data availability

All data and code used and produced for this study is publicly
available and the results can be reproduced easily. Please refer to the
Section “Software and Data Availability”.

References

Ambrosino, C., Chandler, R.E., Todd, M.C., 2014. Rainfall-derived growing season characteristics for agricultural impact assessments in South Africa. Theor. Appl. Climatol. 115, 411–426. http://dx.doi.org/10.1007/s00704-013-0896-y.

Atienza, D., Zielza, C., Larrañaga, P., 2021. Semiparametric Bayesian networks. http://dx.doi.org/10.48550/ARXIV.2109.03008.

Belgiu, M., Drăgut, L., 2016. Random forest in remote sensing: A review of applications and future directions. ISPRS J. Photogramm. Remote Sens. 114, 24–31. http://dx.doi.org/10.1016/j.isprsjprs.2016.01.011.

Breiman, L., 2001. Random Forests. Mach. Learn. 45 (1), 5–32. http://dx.doi.org/10.1023/A:1010933404324.

Bürger, G., 1996. Expanded downscaling for generating local weather scenarios. Clim. Res. 7 (2), 111–128. http://dx.doi.org/10.2307/24864541.

Bürger, G., Chen, Y., 2005. Regression-based downscaling of spatial variability for hydrologic applications. J. Hydrol. 311 (1), 299–317. http://dx.doi.org/10.1016/j.jhydrol.2005.01.025.

Cannon, A.J., 2008. Probabilistic multisite precipitation downscaling by an expanded Bernoulli-Gamma density network. J. Hydrometeorol. 9 (6), 1284–1300. http://dx.doi.org/10.1175/2008JHM960.1.

Chandler, R.E., 2020. Multisite, multivariate weather generation based on generalised linear models. Environ. Model. Softw. 134, 104867. http://dx.doi.org/10.1016/j.envsoft.2020.104867.

Chandler, R.E., Bate, S.M., 2007. Inference for clustered data using the independence log-likelihood. Biometrika 94, 167–183.

Chandler, R.E., Isham, V., Bellone, E., Yang, C., Northrop, P.J., 2006. Space-time modelling of rainfall for continuous simulation. In: B. Finkenstadt, L.H., Isham, V. (Eds.), Statistics of Spatial-Temporal Systems. CRC Press, pp. 169–207.

Charles, S.P., Bates, B.C., Hughes, J.P., 1999. A spatiotemporal model for downscaling precipitation occurrence and amounts. J. Geophys. Res.: Atmospheres 104 (D24), 31657–31669. http://dx.doi.org/10.1029/1999JD900119.

Coe, R., Stern, R., 1982. Fitting models to daily rainfall data (Jordan, niger, botswana, Sri Lanka). J. Appl. Meteorol. 21 (7), 1024–1031.
Dee, D.P., Uppala, S.M., Simmons, A.J., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M.A., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A.C.M., Berg, L.v.d., Bidlot, J., Bormann, N., Dolsel, C., Dragani, R., Fuentes, M., Geer, A.J., Haimberger, L., Healy, S.B., Hersbach, H., Hólm, E.V., Isaksen, L., K° allberg, P., Köhler, M., Matricardi, M., McNally, A.P., Monge-Sanz, B.M., Morcrette, J.-J., Park, B.-K., Peubey, C., Rosnay, P.d., Tavolato, C., Thépaut, J.-N., Vitart, F., 2011. The ERA-Interim reanalysis: configuration and performance of the data assimilation system. Q. J. R. Meteorol. Soc. 137 (656), 553-597. http://dx.doi.org/10.1002/qj.828.

Doblas-Reyes, F.J., García-Serrano, J., Lienert, F., Biescas, A.P., Rodrigues, L.R.L., 2013. Seasonal climate predictability and forecasting: status and prospects. WIREs Clim. Chang. 4 (4), 245-268. http://dx.doi.org/10.1002/wcc.217.

Elidan, G., 2010. Copula Bayesian networks. In: Lafferty, J., Williams, C., Shaw-Taylor, J., Zemel, R., Culotta, A. (Eds.), In: Advances in Neural Information Processing Systems, vol. 23, Curran Associates, Inc.

Fiener, P., Auerswald, K., 2009. Spatial variability of rainfall on a sub-kilometre scale. Earth Surf. Process. Landf. 34 (6), 848-859. http://dx.doi.org/10.1002/esp.1779.

Fowler, H.J., Blenkinsop, S., Tebaldi, C., 2007. Linking climate change modelling to impacts studies: Recent advances in downscaling techniques for hydrological modelling. Int. J. Climatol. 27 (12), 1547-1578. http://dx.doi.org/10.1002/joc.1556.

Fox, E.W., Ver Hoef, J.M., Olsen, A.R., 2020. Comparing spatial regression to random forests for large environmental data sets. PLoS One 15 (3), 1-22. http://dx.doi.org/10.1371/journal.pone.0229509.

Friedman, N., Nachman, I., 2000. Gaussian process networks. In: Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence. UAI '00, Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, pp. 211-219.

Frost, A.J., Charles, S.P., Timbal, B., Chiew, F.H.S., Mehrotra, R., Nguyen, K.C., Chandler, R.E., McGregor, J.L., Fu, G., Kirono, D.G.C., Fernandez, E., Kent, D.M., 2011. A comparison of multi-site daily rainfall downscaling techniques under Australian conditions. J. Hydrol. 408, 1-18, doi:10.1016/j.jhydrol.2011.06.021.

Georganos, S., Grippa, T., Gadiaga, A.N., Linard, C., Lennert, M., Vanhuysse, S., Mboga, N., Wolff, E., Kalogirou, S., 2021. Geographical random forests: a spatial extension of the random forest algorithm to address spatial heterogeneity in remote sensing and population modelling. Geocarto Int. 36 (2), 121-136. http://dx.doi.org/10.1080/10106049.2019.1595177.

Giorgi, F., 2019. Thirty years of regional climate modeling: Where are we and where are we going next? J. Geophys. Res.: Atmospheres 124 (11), 5696-5723. http://dx.doi.org/10.1029/2018JD030094.

Glover, F., 1990. Tabu Search: A Tutorial, vol. 20, http://dx.doi.org/10.1287/inte.20.4.74.

Grayson, R.B., Blöschl, G., 2001. Spatial Patterns in Catchment Hydrology. Cambridge University Press, http://dx.doi.org/10.1002/hyp.5215.

Gutiérrez, J.M., Maraña, D., Widmann, M., Huth, R., Hertig, E., Benestad, R., Roessler, O., Wibig, J., Wilcke, R., Kotlarski, S., Martín, D.S., Herrera, S., Bedia, J., Casanueva, A., Manzanas, R., Iturbide, M., Vrac, M., Dubrovsky, M., Ribalaygua, J., Pórtoles, J., Räty, O., Räisänen, D., Casado, M.J., Ramos, P., Zerenner, T., Turco, M., Bosshard, T., Štěpánek, P., Bartholy, J., Pongracz, R., Keller, D.E., Fischer, A.M., Cardoso, R.M., Soares, P.M.M., Czennecki, B., Pagé, C., 2019. An intercomparison of a large ensemble of statistical downscaling methods over Europe: Results from the VALUE perfect predictor cross-validation experiment. Int. J. Climatol. 39 (9), 3750-3785. http://dx.doi.org/10.1002/joc.5462.

Gutiérrez, J.M., San-Martín, D., Brands, S., Manzanas, R., Herrera, S., 2013. Reassessing statistical downscaling techniques for their robust application under climate change conditions. J. Clim. 26 (1), 171-188. http://dx.doi.org/10.1175/JCLI-D-11-00687.1.

Haile, A.T., Rientjes, T., Gieske, A., Gebremichael, M., 2009. Rainfall Variability over Mountainous and Adjacent Lake Areas: The Case of Lake Tana Basin at the Source of the Blue Nile River. J. Appl. Meteorol. Clim. 48 (8), 1696-1717. http://dx.doi.org/10.1175/2009JAMC2092.1.

Hastie, T., Tibshirani, R., Friedman, J., 2009. The Elements of Statistical Learning. In: Springer Series in Statistics, Springer New York, New York, NY, http://dx.doi.org/10.1007/978-0-387-84858-7.

Haylock, M.R., Cawley, G.C., Harpham, C., Wilby, R.L., Goodess, C.M., 2006. Downscaling heavy precipitation over the United Kingdom: A comparison of dynamical and statistical methods and their future scenarios. Int. J. Climatol. 26 (10), 1397-1415. http://dx.doi.org/10.1002/joc.1318.

He, X., Chaney, N.W., Schleiss, M., Sheffield, J., 2016. Spatial downscaling of precipitation using adaptable random forests. Water Resour. Res. 52 (10), 8217-8237. http://dx.doi.org/10.1002/2016WR019034.

Hengl, T., Heuvelink, G.B.M., Kempen, B., Leenaars, J.G.B., Walsh, M.G., Shepherd, K.D., Sila, A., MacMillan, R.A., Mendes de Jesus, J., Tamene, L., Tondoh, J.E., 2015. Mapping soil properties of afrika at 250 m resolution: Random forests significantly improve current predictions. PLoS One 10 (6), 1-26. http://dx.doi.org/10.1371/journal.pone.0125814.

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Flemming, J.,

Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R.J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., Thépaut, J.-N., 2020. The ERA5 global reanalysis. Q. J. R. Meteorol. Soc. 146 (730), 1999-2049. http://dx.doi.org/10.1002/qj.3803.

Hertig, E., Jacobeit, J., 2008. Assessments of mediterranean precipitation changes for the 21st century using statistical downscaling techniques. Int. J. Climatol. 28 (8), 1025-1045. http://dx.doi.org/10.1002/joc.1597.

Kharin, V.V., Zwiers, F.W., 2003. On the ROC score of probability forecasts. J. Clim. 16 (24), 4145-4150. http://dx.doi.org/10.1175/1520-0442(2003)016<4145:OTSROP>2.0.CO;2.

Klein Tank, A.M.G., Wijngaard, J.B., Können, G.P., Böhm, R., Demarée, G., Gocheva, A., Mileta, M., Pashiardis, S., Hejkrlik, L., Kern-Hansen, C., Heino, R., Bessemoulin, P., Müller-Westermeier, G., Tzanakou, M., Szalai, S., Pálsdóttir, T., Fitzgerald, D., Rubin, S., Capaldo, M., Maugeri, M., Leitass, A., Bukantis, A., Aberfeld, R., van Engelen, A.F.V., Forland, E., Mietus, M., Coelho, F., Mares, C., Razuvaev, V., Nieplova, E., Cegnar, T., Antonio López, J., Dahlström, B., Moberg, A., Kirchhofer, W., Ceylan, A., Pachaliuk, O., Alexander, L.V., Petrovic, P., 2002. Daily dataset of 20th-century surface air temperature and precipitation series for the European climate assessment. Int. J. Climatol. 22 (12), 1441-1453. http://dx.doi.org/10.1002/joc.773.

Koller, D., Friedman, N., 2009. Probabilistic Graphical Models, Principles and Techniques. The MIT Press.

Legasa, M.N., Gutiérrez, J.M., 2020. Multisite Weather Generators Using Bayesian Networks: An Illustrative Case Study for Precipitation Occurrence. Water Resour. Res. 56 (7), http://dx.doi.org/10.1029/2019WR026416, e2019WR026416.

Legasa, M.N., Manzanas, R., Calviño, A., Gutiérrez, J.M., 2022. A posteriori random forests for stochastic downscaling of precipitation predicting probability distributions. Water Resour. Res. 58 (4), http://dx.doi.org/10.1029/2021WR030272, e2021WR030272.

Legasa, M.N., Thao, S., Vrac, M., Manzanas, R., 2023. Assessing three perfect prognosis methods for statistical downscaling of climate change precipitation scenarios. Geophys. Res. Lett. 50 (9), http://dx.doi.org/10.1029/2022GL012525, e2022GL012525.

Manzanas, R., Fiwa, L., Vanya, C., Kanamaru, H., Gutiérrez, J.M., 2020. Statistical downscaling or bias adjustment? A case study involving implausible climate change projections of precipitation in Malawi. Clim. Change 162 (3), 1437-1453. http://dx.doi.org/10.1007/s10584-020-02867-3.

Maraun, D., Widmann, M., Gutiérrez, J.M., Kotlarski, S., Chandler, R.E., Hertig, E., Wibig, J., Huth, R., Wilcke, R.A.I., 2015. VALUE: A framework to validate downscaling approaches for climate change studies. Earth's Futur. 3 (1), 1-14. http://dx.doi.org/10.1002/2014EF000259.

Pearl, J., 1988. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Elsevier, http://dx.doi.org/10.1016/C2009-0-27609-4.

Rodriguez, J., Kuncheva, L., Alonso, C., 2006. Rotation forest: A new classifier ensemble method. IEEE Trans. Pattern Anal. Mach. Intell. 28 (10), 1619-1630. http://dx.doi.org/10.1109/TPAMI.2006.211.

Saha, A., Basu, S., Datta, A., 2021. Random forests for spatially dependent data. J. Amer. Statist. Assoc. 1-19. http://dx.doi.org/10.1080/01621459.2021.1950003.

Sauter, T., Venema, V., 2011. Natural three-dimensional predictor domains for statistical precipitation downscaling. J. Clim. 24 (23), 6132-6145. http://dx.doi.org/10.1175/2011JCLI4155.1.

Schwaighofer, A., Dejori, M., Tresp, V., Stetter, M., 2007. Structure learning with nonparametric decomposable models. In: de Sá, J.M., Alexandre, L.A., Duch, W., Mandic, D. (Eds.), Artificial Neural Networks – ICANN 2007. Springer Berlin Heidelberg, Berlin, Heidelberg, pp. 119-128. http://dx.doi.org/10.1007/978-3-540-74690-4_13.

Scutari, M., 2010. Learning Bayesian networks with the bnlearn r package. J. Stat. Softw. 35 (3), 1-22. http://dx.doi.org/10.18637/jss.v035.i03.

Segond, M.-L., Wheater, H.S., Onof, C., 2007. The significance of spatial rainfall representation for flood runoff estimation: A numerical evaluation based on the lee catchment. UK. J. Hydrol. 347 (1), 116-131. http://dx.doi.org/10.1016/j.jhydrol.2007.09.040.

Singh, V.P., 1997. Effect of spatial and temporal variability in rainfall and watershed characteristics on stream flow hydrograph. Hydrol. Process. 11 (12), 1649-1669. http://dx.doi.org/10.1002/(SICI)1099-1085(19971015)11:12<1649::AID-HYP495>3.0.CO;2-1.

Stern, R.D., Coe, R., 1984. A model fitting analysis of daily rainfall data. J. R. Stat. Soc.: Ser. A (General) 147 (1), 1-18. http://dx.doi.org/10.2307/2981736.

von Storch, H., Zorita, E., Cubasch, U., 1993. Downscaling of Global Climate Change Estimates to Regional Scales: An Application to Iberian Rainfall in Wintertime. Journal Climate 6 (6), 1161-1171. http://dx.doi.org/10.1175/1520-0442(1993)006<1161:DOGCC>E>2.0.CO;2.

Timbal, B., Dufour, A., McAvaney, B., 2003. An estimate of future climate change for western France using a statistical downscaling technique. Clim. Dyn. 20 (7-8), 807-823. http://dx.doi.org/10.1007/s00382-002-0298-9.

Vaittinada Ayar, P., Vrac, M., Bastin, S., Carreau, J., Deque, M., Gallardo, C., 2016. Intercomparison of statistical and dynamical downscaling models under the EURO- and MED-CORDEX initiative framework: present climate evaluations. Clim. Dyn. 46 (3-4), 1301-1329.
Wang, N., Cheng, W., Zhao, M., Liu, Q., Wang, J., 2019. Identification of the debris flow process types within catchments of Beijing mountainous area. Water 11 (4), http://dx.doi.org/10.3390/w11040638.

Weigel, A.P., Liniger, M.A., Appenzeller, C., 2007. The discrete brier and ranked probability skill scores. Monthly Weather. Review 135 (1), 118–124. http://dx.doi.org/10.1175/MWR3280.1.

Widmann, M., Bedia, J., Gutiérrez, J.M., Bosshard, T., Hertig, E., Maraun, D., Casado, M.J., Ramos, P., Cardoso, R.M., Soares, P.M.M., Ribalaygua, J., Pagé, C., Fischer, A.M., Herrera, S., Huth, R., 2019. Validation of spatial variability in downscaling results from the VALUE perfect predictor experiment. Int. J. Climatol. http://dx.doi.org/10.1002/joc.6024.

Yamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O'Loughlin, F., Neal, J.C., Sampson, C.C., Kanae, S., Bates, P.D., 2017. A high-accuracy map of global terrain elevations. Geophys. Res. Lett. 44 (11), 5844–5853. http://dx.doi.org/10.1002/2017GL072874.

Yang, C., Chandler, R.E., Isham, V.S., Wheater, H.S., 2005. Spatial-temporal rainfall simulation using generalized linear models. Water Resour. Res. 41 (11), http://dx.doi.org/10.1029/2004WR003739.

Zennaro, F., Furlan, E., Simeoni, C., Torresan, S., Aslan, S., Critto, A., Marcomini, A., 2021. Exploring machine learning potential for climate change risk assessment. Earth-Sci. Rev. 220, 103752. http://dx.doi.org/10.1016/j.earscirev.2021.103752.

Zhou, L., Meng, Y., Abbaspour, K.C., 2019. A new framework for multi-site stochastic rainfall generator based on empirical orthogonal function analysis and Hilbert-Huang transform. J. Hydrol. 575, 730–742. http://dx.doi.org/10.1016/j.jhydrol.2019.05.047.