Monthly Rainfall Prediction at Catchment Level with the Facebook Prophet Model Using Observed and CMIP5 Decadal Data

Abstract

Early prediction of rainfall is important for the planning of agriculture, water infrastructure, and other socio-economic developments. The near-term prediction (e.g., 10 years) of hydrologic data is a recent development in GCM (General Circulation Model) simulations, e.g., the CMIP5 (Coupled Modelled Intercomparison Project Phase 5) decadal experiments. The prediction of monthly rainfall on a decadal time scale is an important step for catchment management. Previous studies have considered stochastic models using observed time series data only for rainfall prediction, but no studies have used GCM decadal data together with observed data at the catchment level. This study used the Facebook Prophet (FBP) model and six machine learning (ML) regression algorithms for the prediction of monthly rainfall on a decadal time scale for the Brisbane River catchment in Queensland, Australia. Monthly hindcast decadal precipitation data of eight GCMs (EC-EARTH MIROC4h, MRI-CGCM3, MPI-ESM-LR, MPI-ESM-MR, MIROC5, CanCM4, and CMCC-CM) were downloaded from the CMIP5 data portal, and the observed data were collected from the Australian Bureau of Meteorology. At first, the FBP model was used for predictions based on: (i) the observed data only; and (ii) a combination of observed and CMIP5 decadal data. In the next step, predictions were performed through ML regressions where CMIP5 decadal data were used as features and corresponding observed data were used as target variables. The prediction skills were assessed through several skill tests, including Pearson Correlation Coefficient (PCC), Anomaly Correlation Coefficient (ACC), Index of Agreement (IA), and Mean Absolute Error (MAE). Upon comparing the skills, this study found that predictions based on a combination of observed and CMIP5 decadal data through the FBP model provided better skills than the predictions based on the observed data only. The optimal performance of the FBP model, especially for the dry periods, was mainly due to its multiplicative seasonality function.

1. Introduction

Rainfall is a very important climate variable and precious natural resource which affects livelihood and agriculture in many dimensions. An early and accurate prediction of rainfall enables a more efficient management of floods, agriculture, water resources, power development, and the planning and development of infrastructure [1,2,3,4]. However, the prediction of this most important hydrological aspect has become a very challenging task in terms of accuracy due to its peculiar variation over time and space. Due to ongoing climate change, the temporal and spatial variations of rainfall have intensified in the past few decades. Over the past few years, rainfall prediction has become a greater concern to the climate research community [5,6,7,8,9,10,11]. Rainfall prediction approaches are broadly classified into two main categories: (i) knowledge-driven approaches; and (ii) data-driven approaches. Knowledge-driven approaches use scientific understanding, thermodynamic balance, and physical mechanisms of hydrological processes such as General Circulation Models (GCMs). GCMs predict climate variables of coarse spatial resolution on a global scale. However, the knowledge-driven approach needs extensive data and computational facilities that are sometimes unavailable [12]. The data-driven approach is the stochastic and/or empirical statistical modelling approach that is widely used in rainfall prediction at a local level based on the observational relationship of the predictand variable. Data-driven approaches have some limitations, and all approaches do not perform well in predicting longer time spans as they cannot capture the non-linearity and dynamic behaviour of rainfall over time [13,14]. Several statistical/stochastic methods have been used for rainfall prediction and most of them are based on regression analysis, such as simple regression analysis (SRA), exponential smoothing, decomposition, and auto-regressive integrated moving average (ARIMA). Every individual method has its strengths and weaknesses. For instance, ARIMA is a popular stochastic model for time series prediction with great flexibility. However, as a stochastic model, it needs the stationarity of data [15] and its presumed linear form of the associated data sometimes makes it inappropriate for complex nonlinear time series data, such as rainfall [14]. This is why a better output from ARIMA heavily depends on the expertise of the modeller [15]. Dastorani et al. [16] compared different forms of the ARIMA model and concluded that the model parameters need to be tuned to obtain a certain level of accuracy based on location and data type.

Applications of machine learning algorithms, including artificial neural networks (ANNs) of different forms of architecture, are popular for many time series predictions, including the time series of rainfall and enhancing prediction accuracy [6,8,9,17,18,19,20]. According to the level of complexity of the dataset, ANN can be combined with different types of algorithms due to its highly flexible characteristics. However, based on the needs and opportunities, different researchers have come up with different research interests and time scales with the application of ANN. For instance, Wu et al. [21] predicted monsoon rainfall in China more than 10 years ahead, whereas Chakraverty and Gupta [22] predicted Indian monsoon rainfall 6 years in advance. To predict the summer monsoon of India 1 year in advance, Chattopadhyay and Chattopadhyay [23] used 129 years of historical data. Though the ANN is good at capturing nonlinear relationships of data, the presence of outliers in time series data can critically affect the reliability of ANN, as it is a grey box model. Thus, ANN requires proper data pre-processing before its application, especially for climatic data [24,25]. Some other hybrid models have also come into existence and shown very good skills in rainfall prediction [14,26,27]. However, all of the above-mentioned data-driven approaches have used historically observed data and performed predictions for several years ahead based on historical relationships, assuming the climatic conditions will remain the same in the historical and prediction periods.

Compared to the other climatic variables, rainfall has been affected most due to ongoing climate change. Over the past few decades, temporal changes and shifts in rainfall patterns, extreme rainfall during wet periods, longer dry spells during dry periods, and an overall reduction in total precipitation amount have been very common phenomena around the globe. In recent decades, these changes have been intensified due to ongoing climate change [28]. Climate change will continue to change, but the rate of change may be higher in the future. The future higher rate of climate change may adversely affect future precipitation and its level of impact may be significant [28]. Therefore, researchers should not only rely on data-driven approaches (based on historical data only) for future rainfall predictions. For this reason, this study aimed to predict future rainfall for decadal time scales by combining the knowledge and data-driven approaches and employing both GCM-derived precipitation and historically observed data. To complete this, this study used the Facebook Prophet (description provided later) model, where historically observed rainfall was used as an input variable and GCM-derived precipitation data from the decadal experiment of Coupled Model Intercomparison Project Phase 5 (CMIP5) (description provided in Section 2.2) were used as additional regressors to guide the Prophet model in the prediction process. However, the application of the Prophet model in time series prediction is not new [29,30,31], but it is rarely found in the literature for predicting rainfall.

2. Study Area, Data and Methods

2.1. Study Area

The Brisbane River catchment in Queensland (Figure 1) was selected as the study area. It lies in the eastern states of Australia between the latitudes 26.50 S~28.150 S and the longitudes 151.70 E~153.150 E. It has an area of 13,549 square kilometres and a sub-tropical climate where maximum rainfall occurs during summer (December–January–February) and minimum rainfall in winter (June–July–August) [32]. Monthly observed rainfall (1911–2015) over the Brisbane River catchment varied from nil to 1360 mm with an annual average rainfall of 628 mm [33].

2.2. Data Collection

Monthly hindcast precipitation data at a decadal time scale from eight (EC-EARTH MIROC4h, MRI-CGCM3, MPI-ESM-LR, MPI-ESM-MR, MIROC5, CanCM4, and CMCC-CM) GCMs were downloaded from the CMIP5 data portal (https://esgf-node.llnl.gov/projects/cmip5/, accessed on 20 June 2018) for the period 1960–2005, and initialized at 1960, 1965, 1970, and so on, i.e., every five years up to 2005. The names of the models, spatial resolutions, and the available historical runs are given in

Table 1. List of models (GCMs) used as additional regressors in this study.

Modelling Centre (or Group)	Resolutions (Lon × Lat)	Initialization Year (1960–2005)
		60	65	70	75	80	85	90	95	00	05
		Number of Ensembles
EC-EARTH	(1.125 × 1.1215)	14	14	14	14	14	14	14	14	10	18
MRI-CGCM3	(1.125 × 1.1215)	06	08	09	09	06	09	09	09	09	06
MPI-ESM-LR	(1.875 × 1.865)	10	10	10	10	10	10	10	10	10	10
MPI-ESM-MR	(1.875 × 1.865)	03	03	03	03	03	03	03	03	03	03
MIROC4h	(0.5625 × 0.5616)	03	03	03	06	06	06	06	06	06	06
MIROC5	(1.4062 × 1.4007)	06	06	06	06	04	06	06	06	06	06
CanCM4	(2.8125 × 2.7905)	20	20	20	20	20	20	20	20	20	20
CMCC-CM	(0.75 × 0.748)	03	03	03	03	03	03	03	03	03	03

.

An observed monthly gridded rainfall of 0.05° × 0.05° (5 km × 5 km) spatial resolution for the entirety of Australia was collected from the Australian Bureau of Meteorology (BoM). The BoM produces the gridded data using the Australian Water Resources Assessment Landscape model (AWRA-L V5) [34].

2.3. Data Processing

In the first step, all the available ensembles of individual initializations were averaged to produce a single dataset, and then were subset for the Australian region. Secondly, the averaged ensembles were spatially interpolated, using the second-order conservative (SOC) method, into 0.05° × 0.05° spatial resolution, thus matching the grid used in the observed data. This study used the SOC method as it conserves the precipitation flux while sub-gridding the GCM data [35]. It is also marked as the most suitable spatial interpolation method, especially for the GCM-derived gridded dataset [36]. Then, both the models and observed datasets were subset for the Brisbane River catchment. Every initialization spanned a dataset of 10 years that overlapped 5 years with the dataset of the next initialization. In the third step, the last five years of each initialization, except 2005, were discarded and the first five years were combined to produce a single time series from 1961 to 2015. For the initialization year 2005 (2006–2015), the whole dataset was taken instead of the first five years to make the dataset longer.

2.4. Model Description

This study used the FBP model to predict monthly rainfall for a decade (2006 January–December 2015) and then the performance of the FBP’s predicted values was compared with the predictions from six different machine learning regression models: Multi-Layer Perceptron (MLP), Epsilon-Support Vector Regression (SVR), Light Gradient Boosting (LGB), Extreme Gradient Boosting (XGB), Random Forest (RDF), and a combination of these five models. The descriptions of all models are given below.

2.4.1. Facebook Prophet

FBP is a fully automatic open-sourced time series forecasting library developed by Facebook’s core data science team. Though Prophet was built for business purposes, it works for observed hourly, daily, weekly, and monthly time series data that has strong seasonality. It predicts time series as a generalized additive model combining the trend function, seasonality function, holiday effects, and an error term as given in Equation (1) below.

$$Y\left(t\right)=g\left(t\right)+s\left(t\right)+h\left(t\right)+{\in }_t$$

(1)

where $g\left(t\right)$ and $s\left(t\right)$ represent trend and seasonality, respectively, whilst $h\left(t\right)$ presents the holiday effect and ${\in }_t$ is the error term. As this study uses the monthly rainfall data as an input variable, therefore, holiday effect will be invalid here. FBP provides a decomposition regression model that is extendable and easy to use for time series forecasting with a wide range of tuneable parameters. It has functionality for cross-validation to measure forecasting errors and provision to include additional regressors and customize seasonality. The additional regressor feature enhances forecasting accuracy, makes the prediction process more transparent, and helps to tune the prediction process. The additional regressor must be a separately forecasted variable that should be available for both the training and prediction periods.

Prophet can handle outliers without any requirement for imputation or missing data, but the best way to handle outliers is to remove them. Taylor and Letham [37] described further information about Prophet for simulating historical forecasting. Compared to other data-driven approaches, Prophet has two main advantages: (i) Prophet automatically detects changes in trends by selecting change points from the historical data, and it is much more straightforward to create a reasonable, accurate forecast; (ii) its predictions are customizable in ways that are intuitive to non-expert users, and it does not need rigorous data pre-processing. It is easy to use and the components are easily explainable. Its predictions are decent; however, in some cases, certain parameters need to be tweaked compared to the default setting, but that can be easily completed.

2.4.2. Multi-Layer Perceptron (MLP) Regressor

MLP is a class of feedforward artificial neural network (ANN) that utilizes a supervised learning algorithm. It learns by training the dataset using backpropagation with no activation function in the output layer.

2.4.3. Epsilon-Support Vector Regression (SVR)

SVR is also a supervised learning algorithm that acknowledges the presence of non-linearity in the data and provides a proficient prediction using the similar principle of Support Vector Machines (SVMs). The basic idea of SVR is to find the best-fit line that has the maximum number of points. To fit the best line within a threshold, SVR tries to minimize the errors between the real and predicted values.

2.4.4. Gradient Boosting

Boosting is a strategy that combines several simple models into a composite single model. Gradient boosting is a type of boosting and a very popular supervised machine learning technique for regression problems. Light Gradient Boosting (LGB) uses histogram-based learning algorithms following a leaf-wise splitting approach whilst XGB (Extreme Gradient Boosting) uses a level-wise tree growth approach. XGB is a more regularized form of gradient boosting that delivers a more accurate prediction by using the strengths of the second-order derivative of the loss function.

2.4.5. Random Forest Regressor (RDF)

RDF is a supervised learning algorithm that uses the ensemble learning method for regression problems. RDF builds multiple decision trees during the training period and merges them to obtain a more stable and accurate prediction. To control the overfitting problem, a bootstrap technique is used in RDF.

In addition to the five regression models, another combined regression model was developed by stacking the above-mentioned five regression models (referred to as STC) and was used to predict the rainfall for the same period of 2006–2015.

For the predictions, at first, all models were trained from 1911/1961 to 2005 and then predicted for 2006–2015. This study considered two different cases: (i) only the observed data, to demonstrate the data-driven approach (Case-I); (ii) observed data along with GCM-derived precipitation data as the additional regressor (Case-II) to demonstrate the combination of the data and knowledge-driven approaches. For ease of comparison, each of these two cases was further divided into two different subcases; Case-I (a, b) and Case-II (a, b). For Case-I, the FBP model was trained from 1911 (referred to as Case-I (a)) and 1961 (referred as Case-I (b)). For Case-II, two different modes of additional regressors were added to FBP. The first additional regressor was the arithmetic mean of the best five GCMs: MIRCO4h, EC-EARTH, MRI-CGCM3, MPI-ESM-LR, and MPI-ESM-MR (henceforth referred to as MMEM) among the eight considered GCMs in Table 1. The arithmetic mean of these five models showed comparatively better performance compared to the arithmetic mean of the other combinations of the selected models [38]. FBP with MMEM as an additional regressor is referred to as Case-II (a), whilst FBP with all eight GCMs (in Table 1) as eight individual regressors (together) is referred to as Case-II (b). As the CMIP5 decadal data are available since 1961, Case-II and all regression models were trained from 1961 to 2005. For training the regression models, GCM-derived precipitation (MMEM) was used as the independent variable (feature) and corresponding observed data were used as the dependent variable (target variable). After training the regression model from 1961 to 2005, GCM-derived data from 2006 to 2015 was given to the trained regression models to predict the dependent variable (the observed data).

However, to train the model, the most important task was to optimize the model parameters. To optimize the models’ parameters, a wide range of parameter values were given and the best parameters combinations were chosen based on the minimum Mean Absolute Errors (FBP) using the Scikit-learn parameter grid function. For FBP, the optimization of the parameters’ values was performed at a single grid point (latitude 27.50 S and longitude 153.050 E, henceforth referred to as Point-I). The same parameter values were then applied to the two other points. The multiplicative seasonality function with cap and floor values of 600 mm and 0 mm, respectively, was used for PFB prediction. For the regression models, parameters were optimized for all three points and MMEM was used for training and prediction purposes. In the training process, the regression models eventually developed transfer functions (GCMs to observed values) using the best combination of the used parameters (see Table S1) based on the minimum mean squared errors after going through cross-validation 10 times. These transfer functions, obtained from the training period, were then used to transfer the GCM data (2006–2015) to the target variable (referred to as predated data).

As Prophet performs better without outliers, the monthly rainfall values above 250 mm in the observed dataset were set to 250 and zero values were replaced by 1.0 mm. Note that Prophet starts prediction from its training period and continues to predict the future values for the specified period. In addition to the predicted values, it also provides upper and lower limits of the predicted values along with other statistical parameters. The preliminary results revealed that FBP could not reproduce upper (in summer) and lower extreme rainfall (in winter). It also had a correction factor of 0.85 for July and August, 1.15 for December, and an average of the upper limit and the predicted values were employed for January and February. These correction factors were selected based on a trial-and-error basis as well as the performance of the FBP model with the historical prediction (training period) at different randomly selected grids within the Brisbane River catchment. The final predicted values were then examined using four different skills: Pearson Correlation Coefficient (PCC), Anomaly Correlation Coefficient (ACC), Index of Agreement (IA), and Mean Absolute Error (MAE). A brief description of the skill tests is given below, and the detailed descriptions can be found in [39].

2.5. Skill Tests

2.5.1. Pearson Correlation Coefficient (PCC)

PCC is a very commonly used performance metric that measures the linear correlation between two datasets. Here, it is used to measure the linear correlation between the predicted and observed values. Its value varies between −1 and 1 (perfect correlation).

$$\mathrm{PCC}=\frac{{{\displaystyle \sum }}_{t=1}^N\left(P_t-\overline{P}\right)\left(O_t-\overline{O}\right)}{\sqrt{{{\displaystyle \sum }}_{t=1}^N{\left(P_t-\overline{P}\right)}^2}\sqrt{{{\displaystyle \sum }}_{t=1}^N{\left(O_t-\overline{O}\right)}^2}}$$

(2)

where $P$ and $O$ present the predicted and observed values, respectively, and this notation will be the same for the following skill tests also. Bars over the predicted $\left(\overline{P}\right)$ and observed ($\overline{O})$ values represent the mean of the predicted and observed values, respectively. N is the maximum lead time (e.g., the maximum number of months—120).

2.5.2. Anomaly Correlation Coefficient (ACC)

ACC was suggested by Wilks [40] for measuring the correlation between the anomalies of two datasets. Here, ACC is used to measure the temporal anomaly correlation between the anomalies of the predicted and observed values. Anomalies are calculated by subtracting the mean ($C$, mean of the observed values over the entire prediction period) from both the predicted and corresponding observed values.

$$\mathrm{ACC}=\frac{\sum \left\{\left(P-C\right)-\left(\overline{P-C}\right)\right\}*\left\{\left(O-C\right)-\left(\overline{O-C}\right)\right\}}{\sqrt{\sum {\left(P-C\right)}^2}\sqrt{\sum {\left(O-C\right)}^2}}$$

(3)

ACC values range from zero to 1.0 and higher values of ACC do not represent the higher accuracy of the prediction values but the anomalies.

2.5.3. Index of Agreement (IA)

Wilmot [41] suggested the IA for measuring the accuracy of predicted data based on the corresponding observed values. IA values are bounded between 0 and 1 where the value closer to 1 presents the more efficient prediction.

$$\mathrm{IA}=1-\frac{{{\displaystyle \sum }}_{t=1}^N{\left(P_t-O_t\right)}^2}{{{\displaystyle \sum }}_{t=1}^N{(\left|P_t-O^{\prime }\right|+\left|O_t-O^{\prime }\right|)}^2}$$

(4)

Here, $O\prime$ presents the mean of every individual year of the predicted period.

2.5.4. Mean Absolute Error (MAE)

MAE measures the average magnitude of errors and the differences between the predicted and observed values. MAE values range from 0 to ∞ where the lower value indicates higher accuracy and vice versa.

$$\mathrm{MAE}=\frac{1}{N}{\displaystyle \sum }_{t=1}^N\left|P_t-O_t\right|$$

(5)

3. Results and Discussion

There are 496 grids of 5.0 × 5.0 km spatial resolution available in the Brisbane River catchment. This study predicted monthly rainfall for a decade at three different locations (see

Table 2. Comparison of skills and total rainfall prediction among the different cases of FBP models.

Location (Lon/Lat)	Cases	Skills				Under and Overestimation of Total Rainfall (%)
		MAE	PCC	ACC	IA	1Y	3Y	5Y	8Y
Point-I (153.05 E/27.50 S)	I-(a)	53.6	0.549	0.536	0.615	35.9	14.6	−7.2	−11.2
	I-(b)	55.9	0.526	0.418	0.491	11.9	−5.94	−24.6	−28.7
	II-(a)	54.9	0.533	0.517	0.622	33.5	15.1	−8.34	−12.8
	II-(b)	55.1	0.528	0.488	0.577	24.8	5.25	−16.0	−19.3
	MMEM	58.11	0.434	0.433	0.510	48.6	35.6	5.64	−3.1
Point-II (152.0 E/27.0 S)	I-(a)	40.8	0.497	0.496	0.603	50.2	26.5	2.12	−10.4
	I-(b)	41.0	0.484	0.484	0.581	50.7	27.1	2.46	−6.4
	II-(a)	40.8	0.489	0.486	0.593	47.3	26.3	0.53	−8.4
	II-(b)	39.8	0.519	0.517	0.611	38.5	22.7	−1.42	−8.2
	MMEM	41.4	0.494	0.493	0.612	58.2	39.3	13.8	−5.6
Point-III (152.05 E/27.30 S)	I-(a)	46.1	0..491	0.490	0.588	54.4	24.1	3.4	−6.7
	I-(b)	48.1	0.471	0.470	0.583	65.5	32.1	10.6	−0.15
	II-(a)	46.9	0.464	0.460	0.567	51.7	23.1	1.1	−9.9
	II-(b)	45.2	0.490	0.485	0.580	44.6	18.7	−1.6	−10.2
	MMEM	44.7	0.489	0.474	0.571	48.7	23.9	−0.8	−14.2

) closest to the automated weather stations operated by the Australian Bureau of Meteorology. Rainfall predictions for a few months, seasons, and to some extent, for a few years are commonly seen in the literature. This study predicted rainfall for a decade because of the additional regressor for a decadal time scale that was derived through the GCMs that contributed to the decadal experiment of CMIP5. In Section 1, FBP’s predictions of different cases are compared and assessed through different skill tests. In the second part, the performances of FBP for the monthly rainfall predictions are compared with the performance of six different regression models.

3.1. Prediction Using FBP

Figure 2 presents the comparison between observed and FBP-predicted monthly rainfall values of different cases at Point-I, where Case-I presents the data-driven approach, MMEM presents the knowledge-driven approach, and Case-II presents the combination of the knowledge and data-driven approaches. From the comparison, it is evident that FBP can reproduce seasonal variability with better performance when producing dry events, but none of the cases of FBP could reproduce the extreme peak values of the observed rainfall.

However, Case-II (a) showed comparatively better performance to catch the upper peaks, followed by Case-I (a). For the dry events, Case-II showed a considerably better resemblance with the observed values compared to Case-I. This means prediction skills improve, especially in dry events, when the combination of the knowledge and data-driven approaches is employed. These improvements were also observed in the skill tests of all considered cases. Table 2 presents the different skill test results and the percentage of over and under prediction of total rainfall, the cumulative sum over the different periods, and different cases of FBP along with MMEM at three selected points.

Comparing the predicted values and their skill tests at all three points, it is evident that Case-I (a) showed comparatively better skills and a lower percentage of under/over prediction of total rainfall than Case-I (b). This is due to the higher training period (1911–2015) of Case-I (a), which was about double the training period of all other cases. This means, with a higher training period, FBP can reach a better prediction performance. Among the cases of similar training periods, Case-II (a) showed comparatively better skill scores at Point-I where the FBP model parameter values were optimized. However, at points II and III, Case-II (b) showed higher skills than Case-II (a). The better performance of Case-II (b) may be due to the involvement of a range of climate responses provided by the different GCMs as individual regressors, or their higher skills than the MMEM. On the contrary, the better skills of Case-II (a) at Point-I may be either due to the better skills of the MMEM that guided FBP as an additional regressor, or the tuning of the FBP parameters (at the other two points, the FBP models’ parameter values were not tuned).

It is difficult to attach here any valid reason behind why Case-II (b) performed better in the other two points but not at Point-I, as no comparison was conducted during the performance of MMEM and all selected GCMs between the points. Anyway, Case-II (either a or b), as the combination of the knowledge and data-driven approaches, provided better prediction skills than only the data-driven approaches. From the predicted values and their skill comparisons at all three selected points, this study reveals that the skill improvement of FBP predictions with the combination of the knowledge and data-driven approaches was due to better capturing the dry events, and doing so even better than the MMEM (Figure 3). For reproducing the peak values (upper extremes), the knowledge-driven approach (MMEM) was found to be better than the FBP predictions. From Figure 3, it can be observed that, for the dry events (lower values), FBP predictions showed a comparatively better resemblance with the observed values and worse predictions for the upper peaks (wet events). For capturing the dry events, the combination of the knowledge and data-driven approaches was found to be better than any of the individual approaches; Case-II (b) was only slightly better than II (a) in terms of all tested skills, as well as total rainfall prediction over different time spans.

3.2. Prediction Using Regression Models

A comparison of monthly rainfall prediction at Point-I by five supervised machine learning regression algorithms and their stacking models (STC) using MMEM as the feature is presented in Figure 4. From the comparison, it can be observed that the regression models were also able to reproduce seasonal variations with very little improvement in reproducing the lower extremes compared to the MMEM. This study also used all selected models as independent variables, but the results were not as good as the MMEM. For this reason, only the skills of the regression models using MMEM as a feature are presented here. From the comparisons, it can be observed that all regression models, except RDF, showed similar skills in monthly rainfall prediction, among which MLP and SVR were comparatively better than the others (see

Table 3. Skills comparison of different regression models.

Models	Point-I				Point-II				Point-III
	MAE	PCC	ACC	IA	MAE	PCC	ACC	IA	MAE	PCC	ACC	IA
MLP	57.1	0.430	0.371	0.445	39.3	0.480	0.450	0.515	43.5	0.476	0.427	0.494
SVR	57.6	0.430	0.361	0.418	39.3	0.481	0.447	0.516	43.5	0.478	0.430	0.487
LGB	56.6	0.432	0.374	0.442	39.5	0.469	0.432	0.510	43.7	0.466	0.425	0.493
XGB	57.2	0.427	0.370	0.439	39.7	0.451	0.417	0.503	44.1	0.444	0.410	0.484
RDF	57.2	0.426	0.369	0.441	39.9	0.427	0.372	0.433	44.0	0.421	0.359	0.412
STC	57.1	0.434	0.365	0.435	39.1	0.483	0.425	0.475	43.4	0.464	0.405	0.464
FBP(II-a)	54.9	0.533	0.517	0.622	40.9	0.489	0.486	0.593	46.9	0.464	0.460	0.567
MMEM	58.1	0.434	0.433	0.510	41.4	0.494	0.493	0.612	44.7	0.489	0.474	0.571

). From the comparison of the predicted monthly rainfall values and the skills tests, this study finds that the regression models showed little improvement in reproducing the lower extremes compared to the MMEM, and the reverse in reproducing the peak values (upper extremes).

In reproducing the lower extremes, FBP showed better performance than the regression models, and this was the same even for the upper extremes. Note that the correction factors are employed in FBP predictions, but no such correction factors are considered for the regression models. However, the reason behind the better performance of FBP to catch the dry events could be due to employing the correction factor or due to choosing the multiplicative seasonality function, or a combination of both, enabling FBP to reproduce wider seasonal variation compared to the regression models [42]. The skill test results showed little weak prediction where IA values varied between 0.5 and 0.62, and MAE values varied between 40 and 50. The main reason behind these weak skills was that there were very frequent extreme peaks in the observed values (see Figure 2 or Figure 4) which were used to measure the skills. Another reason may be the comparatively shorter training period (1960–2005) where FBP was not familiar with the rainfall values above 250 mm for the target variable during the training process, as these were set to 250 mm to remove the outliers from the observed datasets.

Upon comparing Figure 3 and Figure 5, it can be seen that FBP and the regression models showed a comparatively better resemblance of the lower extreme values (dry events) than the MMEM. Note that GCMs are imperfect replicas of real-world phenomena and contain systematic biases [43]. They intend to overestimate the wet events and underestimate the dry events [44,45]. GCM outputs need rigorous correction before application [46,47,48,49]. This study reveals that rainfall prediction at the local level using FBP, using a combination of both GCMs and observed data, would enhance the overall prediction accuracy. A longer training period using a longer period of GCM-derived hindcast data may enable FBP to provide better skills in reproducing dry events, which was seen in the first case (Case-I). Note that this study used two different types of datasets (GCMs and observed) of the same time span and followed a supervised training approach where known observed values were used for both the training and the prediction period (for comparison with the predicted values). As no real future data were predicted here, one can consider the training period as verification and the assessed prediction skills as validation.

Every individual future prediction model, either an ANN or any ML algorithm, is different, and shows different performance depending on function, tuning parameters, and the variables considered for predictions. In addition, every individual precipitation time series is different at different geographical locations. For the prediction purpose, researchers used different models with different input variables and data pre-processing techniques. Then, performances were assessed based on the corresponding observed values. That is why it is difficult to compare and contrast different types of models for different regions and periods with the results obtained in this study.

The early prediction of upper and lower extremes can help in managing floods and droughts, respectively. However, in this study, it was observed that neither the FBP model nor the ML regression model could reproduce extreme wet events. Rather, they could reproduce the dry events considerably better than the wet events. Reproducing the dry events is also important to agriculture-dependent countries such as Australia, where the most variable climate exists. In Australia, a typical major drought in a season may reduce agricultural production by about 10% and gross national product by 1% [50]. This study will be beneficial for: water resource managers assessing future water availability; managing agriculture and agricultural-dependent businesses; and other water-related stakeholders planning and developing infrastructures.

4. Conclusions

Rainfall prediction is highly important from both social and economic perspectives. Predictions through GCMs and other time series models have been seen individually, but the application of GCM-derived precipitation together with observed values in a time series model on a decadal scale has not yet been found. On the contrary, as a relatively new time series forecasting library, the application of Facebook Prophet in rainfall prediction can hardly be found in the literature. For time series data, which have strong seasonality, Facebook Prophet works well for future prediction. As climate variables show seasonality over the cycle of a year, this study aimed to predict monthly rainfall using a combination of both GCM-derived and observed data through the Facebook Prophet model. In doing so, this study used historically observed monthly rainfall as the input and GCM-derived monthly precipitation from the CMIP5 decadal experiment as an additional regressor. Multiple additional regressors were implemented and compared with the performance of a single additional regressor. A few correction factors were introduced for the predicted values of different months that enabled FBP to provide better prediction accuracy. From the comparison of skills, this study finds that the combination of both GCM-derived and observed values gives better prediction accuracy compared to predictions based on the observed data only. Using GCM-derived data as an additional regressor guided the FBP model for future prediction. GCM-derived data includes not only scientific understanding, but also the historical records that guided FBP to come up with higher prediction accuracy. Based on the outlined skills assessments, the following conclusions are drawn.

(i)

FBP can reproduce dry events considerably better than wet events. This may be due to a better understanding of FBP of dry periods through the training and its multiplicative seasonality function;

(ii)

Following the combination of GCM-derived data (as an additional regressor) and the corresponding observed values, FBP should be able to reproduce future rainfall with higher prediction accuracy than the predictions based on the observed values only;

(iii)

A higher number of regressors will provide comparatively better prediction accuracy than a single additional regressor. In this case, a longer period of GCM hindcast data would elicit a higher prediction accuracy.

However, this study compared the performance of FBP with six regression models for the same places and same datasets, and found that FBP outperformed them. This study highly encourages the cross-validation of a similar approach by using different forms and architecture of deep ANNs that may increase prediction accuracy by utilizing different tuneable features.