MSCliGAN—A Structure-informed Generative Adversarial Model for Multi-site Statistical Downscaling of Extreme Precipitation Using Multi-model Ensemble

1. Introduction

As communities adapt to changes in local weather extremes brought on by climate change, accessing the best and latest information available at locally relevant scales is essential. However, climate models typically operate at spatial scales out of step with hydrological, ecological, economic decision-making needs and must be downscaled to higher spatial resolutions. The use of statistical downscaling is widespread for climate model outputs, however there remains a shortage of methods that can downscale precipitation from Atmosphere-Ocean General Circulation Model (AOGCM) simulated precipitation to regional-scale gridded precipitation (Tryhorn et. al. 2011, Chaudhuri et. al., 2017). Precipitation information, especially extremes, is critical for flood modelling and mitigation (Hirabayashi et al. 2013). There is considerable intrinsic complexity in precipitation compared to other climate variables such as temperature or pressure. Precipitation is more dispersed in space and associations at various atmospheric scales (local, meso, synoptic) are more apparent in observed precipitation trends. It is extremely challenging to model such patterns using continuous functions used in conventional statistical downscaling methods, incorporating nonlinearities into downscaling functions is seen as critically important if accurate representations are to be realized (Weichert and Burger 1998). The latest developments in machine learning (ML) approaches such as convolutional neural networks have started to tackle these long-standing problems (Shi et al. 2015). However, commonly used loss functions such as Mean Absolute Error (MAE) and Nash-Sutcliffe Loss (NS) optimize for overall simulation efficiency but neglect precipitation spatial structure; a crucial property if replicating observed local trends is a goal (Plouffe et al. 2015). Narrowly specified loss functions hinder the ability of machine learning models in downscaling, contributing to weak output of regional scale models (Zhang et al. 2006).

In the context of climate change science, severe precipitation can lead to flooding, changes in hydrological morphology, and contribute to water pollution (Carpenter et al., 2018; Eekhout et al., 2018; Raghavendra et al., 2010). In comparison to normal precipitation, extreme precipitation is unpredictable, ambiguously defined, and more fragmented in space, which is difficult to forecast with any certainty (Lee et al., 2017). Downscaling experiments have also found extreme precipitation to be challenging to reproduce at meaningful local scales (Castellano and DeGaetano, 2016). Anticipating changes in extreme precipitation remains an urgent climate information need in many regions of the world.

Global Climate Models (GCMs) have been widely used to obtain climate change data under various scenarios to determine the potential impacts of climate change on hydrological processes in nature (Cho et al., 2016; Li et al., 2017). However, the native resolution of GCMs are too coarse which makes them not suitable for hydrological studies if used directly (Xu, 1999). To mitigate this problem downscaling techniques have evolved where large-scale atmospheric data given by AOGCMs are used along with regional scale observations to generate high-resolution realizations of the variable of interest at over a specified area (Chen et al., 2012). There are two primary types of downscaling methods, statistical and dynamical downscaling (Yhang et al., 2017). Statistical downscaling is more common in climate change impact studies due to its clear design methodology, computational ease and the ability to generate synthetic datasets of any desired length (Hidalgo et al., 2008; Maraun, D., et al. 2010; Trzaska, Sylwia & Schnarr, Emilie. 2014; Widmann et al., 2003).

Popular statistical downscaling techniques include the transfer function method, weather pattern method, and stochastic weather generators (Kioutsioukis et al. 2008; Murphy et al. 2004; Wilby and Harris 2006). Linear regression was used frequently in the past to connect circulation factors with weather variables (Huth, 2002; Zorita and von Storch, 1999). However, these relationships are always complex and non-linear, making it hard to achieve acceptable downscaling results with linear relationships (Ghosh and Mujumdar, 2008). As a result, in recent years machine learning methods have been gaining popularity given their ability to model non-linear relationships between predictor and predictants (Chaudhary et al., 2019; Yi et al., 2018).

Machine learning and data science methods have significantly expanded in various fields in recent years leading to application of tools in new areas and development of new domain-specific variants of existing computational methods. Convolutional neural network (CNN) modelling has become common due to lower computational requirements than the complex dense networks (LeCun et. al. 1998). There are a wide variety of CNN implementations including satellite image change detection (Wang et al. 2019), image segmentation (Ronneberger et. al. 2015), image recognition (Krizhevsky et. al. 2017). The approach used by CNNs for increasing image resolution (Dong et. al. 2014) is an algorithm that is very close in some respects to downscaling of climate model outputs. The central question in these studies is simply, how do we detect and move representative information from one scale down to another scale?

CNNs, which have seen widespread adoption for image processing applications in a variety of fields, cannot recognize the variability present in the results when trained through the general loss functions. This leads to the degradation in performance of CNNs when the training loss function cannot provide significant gradients to the parameters of the network. We believe that the development of appropriate loss functions is integral to application of CNNs (Zhao et al. 2016). Goodfellow et al. (2014) recommended the Generative Adversarial Training (GAN), in which a zero-sum game is played between two networks and in recent years was commonly used in training deep neural networks. This approach provides superior network training and can yield results that appear superficially like the truth and intuitively tackle the issue of gradients.

In this paper we develop a new GAN downscaling approach which can downscale annual maximum rainfall from several AOGCMs to regional-scale gridded annual maximum rainfall at various sites. This has the potential therefore to develop regionally specific downscaling models incorporating local landcover and topography into the downscaling model. The following specific research objectives are tackled by this study;

• Develop a methodology to downscale large-scale precipitation, given by several AOGCMs, to regional-scale precipitation at multiple regions of interest by statistical downscaling using Convolution Neural Network and Generative Adversarial Training.
• Analyze the effect of input DEM and Land-use/ Land-cover on the simulation output.
• Evaluating the efficacy of the novel loss function that incorporates content loss, structural loss and adversarial loss that enhances the estimation of the downscaled precipitation’s global and regional quality.

2. Methods

2.1. Study Area and Datasets

The model is developed for seven cities located in different regions of Canada. The locations of these cities are given in Figure 1. Individual location plots are shown in Figure S1. The salient features of these regions are summarized in Table 1.

We used annual maximum daily precipitation as our target variable to establish a downscaling methodology. Annual maximum daily precipitation data from the NRCANMET daily gridded precipitation dataset (Hutchinson et al. 2009) were extracted for the period 1950-2010. Although grid data derived from Canada’s complex observation network density is not ideal validation data, these are the best available data for these regions and should reasonably capture regional trends. Wilby et al. (1999) and Wetterhall et al. (2005) addressed the properties of any predictor in statistical downscaling, (1) it should be simulated by GCMs with high reliability, (2) it can be easily accessed from the archives of GCM outputs, and (3) it has a strong correlation with the surface variables of interest (rainfall in the present case). GCMs simulated rainfall includes atmospheric dynamic information as well as knowledge about the impact of climate change on rainfall through various physical parameterizations. It encouraged us to choose the precipitation simulated by AOGCM as a predictor of the observed (i.e., gridded) precipitation.

Table 1. Salient geographical and climatic features of the regions of interest.

Name of City	Area (sq. Km)	Elevation (m)	Population	Temperature Range (°C)	Average Monthly Precipitation Range (mm)
Schefferville	39	513	213	−24.5 to 12.2	29.7 to 114.6
Goose Bay	305.69	12	8109	−17.6 to 15.5	56.8 to 121.3
Yellowknife	136.22	206	19,569	−25.6 to 17.0	11.3 to 40.8
Edmonton	767.85	645	932,546	−12.1 to 16.2	12.0 to 93.8
Calgary	825.56	1045	1,239,220	−7.1 to 16.2	9.4 to 94.0
Saskatoon	228.13	481.5	246,376	−15.5 to 18.5	8.8 to 65.8
Regina	179.97	577	215,106	−14.7 to 18.9	9.4 to 70.9

Table 1. Salient geographical and climatic features of the regions of interest.

Name of City	Area (sq. Km)	Elevation (m)	Population	Temperature Range (°C)	Average Monthly Precipitation Range (mm)
Schefferville	39	513	213	−24.5 to 12.2	29.7 to 114.6
Goose Bay	305.69	12	8109	−17.6 to 15.5	56.8 to 121.3
Yellowknife	136.22	206	19,569	−25.6 to 17.0	11.3 to 40.8
Edmonton	767.85	645	932,546	−12.1 to 16.2	12.0 to 93.8
Calgary	825.56	1045	1,239,220	−7.1 to 16.2	9.4 to 94.0
Saskatoon	228.13	481.5	246,376	−15.5 to 18.5	8.8 to 65.8
Regina	179.97	577	215,106	−14.7 to 18.9	9.4 to 70.9

Figure 1. The locations of the cities are shown against Canada in red circles.

Figure 1. The locations of the cities are shown against Canada in red circles.

This research included nine AOGCMs of the CMIP6 database (Table 1). Different AOGCMs have different grids and many of them are in Gaussian grids, making it difficult to view them in a single mathematical framework. We have interpolated precipitation from different AOGCMs to 10 km resolution grids specified in output to address this obstacle. Further, we used GEBCO (GEBCO Compilation Group, 2019) land and ocean terrain model as topographic input and 2010 Canada Land cover (Canada Centre for Remote Sensing, 2019) as LULC input to our model. Topography was interpolated to the input grids using a bi-quadratic interpolation function. Figure 2 shows the topography over 7 regions of interest. The LULC was resampled using dominant categories over the input grid. Figure 3 shows the dominant LULC grid over the 7 cities where we are building the downscaling model.

Table 2. Features of different input AOGCMs.

AOGCM	Institution	Grid Type	Horizontal Dimension (Lon × Lat)	Vertical Levels
BCC ESM	Beijing Climate Center	T42	128 × 64	26
CAN ESM5	Canadian Centre for Climate Modelling and Analysis, Environment and Climate Change Canada	T63	128 × 64	49
CESM2	National Center for Atmospheric Research	0.9 × 1.25 finite volume grid	288 × 192	70
CNRM CM6.1	Centre National de Recherches Meteorologiques	T127	256 × 128	91
CNRM ESM2	Centre National de Recherches Meteorologiques	T127	256 × 128	91
GFDL CM4	Geophysical Fluid Dynamics Laboratory	C96	360 × 180	33
HAD GEM3	Met Office Hadley Centre	N96	192 × 144	85
MRI	Meteorological Research Institute, Tsukuba, Ibaraki 305-0052, Japan	TL159	320 × 160	80
UK ESM1	Met Office Hadley Centre	N96	192 × 144	85

Table 2. Features of different input AOGCMs.

AOGCM	Institution	Grid Type	Horizontal Dimension (Lon × Lat)	Vertical Levels
BCC ESM	Beijing Climate Center	T42	128 × 64	26
CAN ESM5	Canadian Centre for Climate Modelling and Analysis, Environment and Climate Change Canada	T63	128 × 64	49
CESM2	National Center for Atmospheric Research	0.9 × 1.25 finite volume grid	288 × 192	70
CNRM CM6.1	Centre National de Recherches Meteorologiques	T127	256 × 128	91
CNRM ESM2	Centre National de Recherches Meteorologiques	T127	256 × 128	91
GFDL CM4	Geophysical Fluid Dynamics Laboratory	C96	360 × 180	33
HAD GEM3	Met Office Hadley Centre	N96	192 × 144	85
MRI	Meteorological Research Institute, Tsukuba, Ibaraki 305-0052, Japan	TL159	320 × 160	80
UK ESM1	Met Office Hadley Centre	N96	192 × 144	85

Figure 2. Topography for (a) Schefferville, (b) Goose Bay, (c) Yellowknife, (d) Edmonton, (e) Calgary, (f) Saskatoon, (g) Regina.

Figure 2. Topography for (a) Schefferville, (b) Goose Bay, (c) Yellowknife, (d) Edmonton, (e) Calgary, (f) Saskatoon, (g) Regina.

Figure 3. Land-Use/Land-Cover for (a) Schefferville, (b) Goose Bay, (c) Yellowknife, (d) Edmonton, (e) Calgary, (f) Saskatoon, (g) Regina.

Figure 3. Land-Use/Land-Cover for (a) Schefferville, (b) Goose Bay, (c) Yellowknife, (d) Edmonton, (e) Calgary, (f) Saskatoon, (g) Regina.

2.2. Downscaling Method

We developed a generating network $G_{{\theta }_G}$ that can turn extreme precipitation AOGCM coarse resolution (P_GCM) into extreme precipitate (P_obs) fine resolution at the regional level. Increasing resolution, bias correction and regional precipitation feature corrections are the technical challenges of this problem. The number of years for training y = 1, …, n then the θ_G is obtained by minimizing the downscaling total loss function $l^D$;

$$\widehat{{\theta }_G}=\frac{argmin}{{\theta }_G} \frac{1}{n}{{\displaystyle \sum }}_{y=1}^n{l}^D\left(G_{{\theta }_G}\left(P_{GCM}^y\right),P_{obs}^y\right)$$

(1)

We developed a weighted combination of many components for downscaling total loss. These elements are discussed later in details.

2.2.1. Adversarial Training

As a discriminatory network (Goodfellow et al. 2014), Wasserstein GAN (WGAN) (Arjovsky et. al. 2017) is used in our context. The Earth-Mover (also called Wasserstein-1) distance W(q, p), which is informally defined as a minimum cost to transport mass to convert the distribution q into p (where the cost is mass times transportation distance), was approximated. Formally, the game between the generator $G_{{\theta }_G}$ and the discriminator $D_{{\theta }_D}$ is a minimax objective. The WGAN loss function is developed using the Kantorovich-Rubinstein [Villani 2008] duality to obtain;

$$\frac{min}{G_{{\theta }_G}}\frac{max}{D_{{\theta }_D}\in D}E \left[D_{{\theta }_D}\left(P_{obs}\right)\right]-E\left[D_{{\theta }_D}\left(G_{{\theta }_G}\left(P_{GCM}\right)\right)\right]$$

(2)

where, D is the set of 1-Lipschitz functions. To satisfy the Lipschitz constraint on the WGAN, Arjovsky et. al. (2017) propose clipping of the weights of the WGAN to ensure that the function lies within a compact space [−c, c]. The set of functions (discriminator under training) satisfying this constraint are subsets of the k-Lipschitz functions for some k which depends on c and the WGAN architecture.

2.2.3. Downscaling Total Loss

In our framework the Total loss is designed as a linear combination of 3-loss components: content loss, structural loss, adversarial loss.

$$Totalloss=contentloss+structuralloss+adversarialloss$$

Content Loss

The Nash–Sutcliffe model efficiency coefficient (NSE) is widely used loss estimate which assesses the predictive ability of any hydrology model. We have taken (1-NSE) as content loss.

$$Content Loss=\frac{{{\displaystyle \sum }}_{y=1}^n{\left(G_{{\theta }_G}\left(P_{GCM}^y\right)-P_{obs}^y\right)}^2}{{{\displaystyle \sum }}_{y=1}^n{\left(P_{obs}^y-\frac{1}{n}{{\displaystyle \sum }}_{y=1}^n{P}_{obs}^y\right)}^2}$$

(3)

This loss can take range 0 to $\infty$. A value of 0 loss means perfect match with the observation. The loss value 1 means that the model is as accurate as the average of the observed values. The loss greater than 1 signifies that the observed mean is a better predictor than the model.

Structural Loss

The overall low error cannot maintain the regional information common to precipitation at high resolution. We implemented a structural loss function to tackle this problem and preserve regional details. The luminance, contrast, and structure of two 2D fields are compared by the Structural Similarity Index (SSIM) (Wang et al. 2003). This matric was compared to preserve the regional precipitation structure on a fine scale. Our loss was described by (1-SSIM)/2 based on structural dissimilarity.

$$Structural Loss=\frac{1}{2}\left(1-SSIM\right)=\frac{1}{2}\left(1-l_i^{{\alpha }_i}{c_i}^{{\beta }_i}{s}_i^{{\gamma }_i}\right)$$

(4)

${\alpha }_i$, ${\beta }_i$, and ${\gamma }_i$ are weights of luminance, contrast, and structure components respectively at scale i. We used α = β = γ = 1. For any given scale the luminance comparison is given by;

$$l=\frac{2{\mu }_{obs}{\mu }_G+c_1}{{\mu }_{obs}^2+{\mu }_G^2+c_1}$$

(5)

The contrast comparison is given by;

$$c=\frac{2{\sigma }_{obs}{\sigma }_G+c_2}{{\sigma }_{obs}^2+{\sigma }_G^2+c_2}$$

(6)

and the structure comparison is given by;

$$s=\frac{{\sigma }_{obs,G}+c_3}{{\sigma }_{obs}{\sigma }_G+c_3}$$

(7)

where ${\mu }_{obs}$ and ${\mu }_G$ are the mean of observation and generated values respectively, ${\sigma }_{obs}$ and ${\sigma }_G$ are standard deviation of observation and generated values respectively, and ${\sigma }_{obs,G}$ is the co-variance between observation and generated values. $c_1$, $c_2$, and $c_3$ are small constants to stabilize the divisions with weak denominator. We used $c_1={\left(k_{1}L\right)}^2$, $c_2={\left(k_{2}L\right)}^2$, $c_3=\frac{c_2}{2}$ where $k_1=0.01$, $k_2=0.03$ and $L$ is the dynamic range of the input which we have taken as 100. We used a scale of 3x3 windows to define SSIM.

Adversarial Loss

Adversarial loss is given by the approximate earth moving distance between the observation and generator output. It is given by;

$$Adversarial Loss=D_{{\theta }_D}\left(G_{{\theta }_G}\left(P_{GCM}\right)\right)$$

(8)

2.2.4. Networks

Generative network is developed into a U-net (Ronneberger et al. 2015) convolution-style encoder decoder with skipped link. LeakyReLu with momentum 0.2 was used for activating the hidden layers and for downsampling. In the last generator layer, ReLu has been used to provide positive output which is important for precipitation simulation.

A convolution encoder with a strided convolution downsampler is used as the discriminative network. We used LeakyReLu with momentum 0.2 as the activation function in the hidden layers and output layer is activated by a linear function. The approximate earth moving distance (Wasserstein distance) between the observed and simulated precipitation is calculated by the discriminator. We have built a compact support for layer weights by cutting them between [-0.1, 0.1] following recommended parameterization for WGAN models.

The schematic for the networks and training procedure are shown in Figure 4. The black arrows indicate the flow of data and red arrows indicate the feedback of objective functions to the models in terms of optimized parameters.

2.2.5. Training Details

The training was done on the 4000 GPU NVIDIA Quadro. Input / output data rescaling were avoided in any way and the model was able to learn and evolve the climate signal present in the various model datasets and observations. Optimization of the generator was done using the L_inf norm (Adamax) (Ruder et. al. 2016) Adam gradient-base optimization, with initial learning rate 0.02 and decay 0.5. We wanted the discriminator to calculate the gradient in the current iteration based on an adaptive window and did not want to use the previous gradient data generator results. We used the Adam (Adadelta) moving window update (Zeiler et. al. 2012) to train the Wasserstein distance computer, discriminator. The model was trained for 15,000 iterations which are more than sufficient to simulate high-resolution realistic precipitation patterns. With just 10,000 iterations, even the error stabilized.