Performance assessment of spatio-temporal regression kriging with GAMLSS models as trends

: The main objective of this study is to propose different probabilistic models for adjusting the trend component, since it significantly influences the quality of the spatio-temporal interpolation of rainfalls. We used the monthly total precipitation data of the São Francisco River Basin (SFRB) for the period of 31 years, 1989–2019. The SFRB occupies 8% of the whole Brazilian territory, mostly located in the Northeast Brazilian region. For the trend component, we propose the fitted GAMLSS models by comparing different probability distribution families, which in most cases include the characteristics of these data. The results indicate the existence of a spatio-temporal pattern of the residues obtained from the adjustment of the trend with zero adjusted Gamma distribution for the accumulated monthly precipitation. The adjustment revealed a spatial dependence of up to 873 km between the pluviometric stations and temporal autocorrelation of approximately 1.6 months. The methodology used in this study enabled us to create rainfall maps, interpolating unobserved locations in differences years. The projection of these maps to the SFRB is considered extremely important for planning and implementing activities related to water resources across the river basin.


INTRODUCTION
In recent years, the increasing frequency and intensity of droughts have been the focus of several studies associated with climate risk, most of them provoked by the unprecedented water crisis.Some experts believe that the water crisis in the 21 st century is more a matter of management than the real crisis of scarcity and stress (Rogers et al. 2005).However, for other experts, this is the result of a series of problems that are aggravated by other issues related to economic and social development (Gleick & Cain 2004).In one of the most recent studies carried out by the Intergovernmental Panel on Climate Change (IPCC), published on August 9, 2021, confirms that the water cycle has been intensifying and will continue to intensify the planet warms, causing more intense rains and prolonged droughts all around the world (IPCC 2021).In the last 2,000 years, climate change, mainly caused by man, has intensified at an exponential rate, causing serious problems for humans and the planet, such as storms, droughts and heat waves (IPCC 2021).Therefore, there is a need to investigate how climate change has impacted the climate, especially in the spatial and temporal distribution of precipitation in the region.
Several studies have examined the spatial and temporal variability in rainfall precipitation (Foehn et al. 2018, Achite et al. 2017).Gathering precipitation information with measuring instruments placed at specific locations can be an easy task, but what we need is to obtain the precipitation estimates for locations not sampled.Several methodologies can be used to obtain these estimates, such as regression analysis (Stauffer et al. 2016) and geostatistical models that consider spatio-temporal interaction (Kilibarda et al. 2015, Martínez et al. 2017).
Once we find that the phenomenon under study reveals spatial and temporal variability, we can obtain better results with models to consider the space-time relationship as well as its spatial, temporal, and spatio-temporal component structure (Hu et al. 2017, Monteiro et al. 2017).
In order to apply the methodology presented in this paper, we use the precipitation data set for the São Francisco River Basin (SFRB) evaluated in space and time.The SFRB covers a drainage area of 639,19 km², corresponding to 8% of Brazil and an average flow of 2,850 m³/s, representing 2% of the entire Brazilian territory.Since 2012, the SFRB faces several adverse hydrological conditions, with below-average flows and precipitations reducing the volume of storage in their reservoirs.The 2014-2016 period recorded the lowest annual average natural outflows in the Sobradinho reservoir since 1931, the year the Operator of the National Electricity System (ONS) began recording observations (ANA 2017).The São Francisco River is of great importance to the semiarid Northeast Brazilian (NEB) region for the supply of water for human consumption and irrigation, generation of energy by the Paulo Afonso and Sobradinho hydroelectric power plant, tourism and fishing (Araujo & Celeste 2019, Silveira et al. 2016).
The frequency and intensity that climatic extremes increase also generates an increase in global warming and will continue to increase in a scenario of medium and high greenhouse gases emissions causing droughts and climate changes in rainy regions in the coming centuries, changing basins and slopes.One of the latest IPCC reports (IPCC 2014) predicts a drastic reduction of up to 22% in precipitation levels in NEB, impacting the distribution of hydraulic power, mainly due to an estimated decrease of 24.6% in São Francisco River flow.
The main objective of this study is to propose different probabilistic models for adjusting the trend component, since it significantly influences the quality of the spatio-temporal interpolation of precipitation.Precipitation maps in the SFRB is considered of extreme importance for the planning and execution of activities related to water resources, identifying priority regions for executing governmental actions in order to minimize the disasters caused by the scarcity of rain in the region, especially, the desertification process that has caused situations of socioenvironmental vulnerability in the region.

Study area and data set
In this study, we used the monthly total precipitation data of the SFRB provided by the Agência Nacional das Águas (ANA) for the period of 31 years, 1989-2019 (Figure 1).We considered only the pluviometric stations that presented 10 years or more of information, resulting in 341 stations.Consequently, several pluviometric stations kept of this study presented incomplete data, however, the methodology proposed in this research, spatio-temporal geostatistics, incorporates this type of question, taking into account the neighboring stations for precipitation interpolation in non-sampled locations.
The main river in the SFRB is the São Francisco River.This river originates in Serra da Canastra, the State of Minas Gerais, and covers a drainage area of 639,219 km² in an extension of 2,863 km, corresponding to 8% of Brazil.On its natural course to the ocean, the river passes through the states of Minas Gerais, Bahia, Pernambuco, Sergipe, and Alagoas, before flowing into the Atlantic.For planning, the SFRB is divided into four physiographic regions: Upper, Middle, Lower-Middle, and Lower São Francisco.Figure 1 shows the SFRB representation identifying these four regions.
According to the classification of Köppen-Geiger (Alvares et al. 2013), the SFRB includes five different climatic regions: tropical climate (Awwith a summer rainy season; As -with a winter rainy season), semi-arid climate (BSh -low latitude and altitude), and humid subtropical (CWa -with dry winter and hot summer; CWb -with dry winter and temperate summer).Consequently, the different ecosystems throughout the region comprise the Atlantic Rainforest (SFRB Headland), Cerrado (Upper and Middle São Francisco), Caatinga (East and Middle-Lower São Francisco), Caatinga, Atlantic forest, mangroves, and coastal vegetation (Lower São Francisco) (Bezerra et al. 2019).
Because of its extension, several atmospheric systems act in the distribution of rainfall on the SFRB: the Intertropical Convergence Zone (ITCZ) and High-Level Cyclonic Vortices (HLCV) act in Lower São Francisco mainly in the first months of the year; the South Atlantic Convergence Zone (SACZ) and frontal systems (cold fronts) act in the Upper, Middle, and Lower-Middle São Francisco during the austral summer; the East Undulatory Disorders (DOL), Instability Lines, and breezes act throughout the year (Cavalcanti 2016).
The El Niño phenomenon, which designates the significant changes in temperature on the surface of waters in the Pacific Ocean, has negatively influenced the rainfall regime in regions of Northeastern Brazil.Costa (2012) found that the strongest El Niño events imply longer droughts in northeastern Brazil, resulting in cycles that alter weather and climate conditions and negatively impact rainfall distribution across the region.Gurjão et al. (2012) found that El Niño and the action of the HLCV inhibited the development of clouds and contributed to the reduction of rainfall in several locations in the SFRB.Bezerra et al. (2019) stated that the Lower-Middle region of São Francisco is the most vulnerable to severe droughts during periods of strong El Niño.

Spatio-temporal model
Observations of rainfall precipitation were modeled following the spatio-temporal stochastic process, considering the random function this is d=2, corresponding to the spatial domain dimension where the geographical coordinates (latitude and longitude) of each sampled point was used.The Z(s,t) is decomposed into the sum of the trend component m(s,t) and stochastic residual ε(s,t): ( In Equation 1, m(s,t) = E[Z(s,t)] (s,t) gives the expected value of the process Z(s,t), representing the large-scale variation in the process.

Trend component
In most of the observed phenomena, component cannot be modeled by a deterministic function because of the randomness in the process.Thus, we need to model this component through stochastic functions designed to model data patterns observed in space and time (Monteiro et al. 2017).The modeling of this component included the covariates collected along with the response variable.In this study, we used each sample point's geographical coordinates and annual effects inserted through a time index.To model this regression, we considered a Generalized Additive Models for Location, Scale, and Shape (GAMLSS) since it allows for greater flexibility in the regression fitted model.For comparison purposes, for the fitted model, we used the normal Linear Model (LM) and the GAMLSS with Normal (NO), Zero-Adjusted Gamma distribution (ZAGA) and Generalized Gama (GG) models.Let , , , =  T n z z z z be a vector for response variable with n observations.Also, assume that observations ( ) ( ) In Equation ( 2) θ k and η are vectors of size n; β k is a vector of size parameters J' k ; X k and Y jk are design matrices of dimension n×J' k and n×q ij , respectively; and ϕ jk is a q ij -dimensional random.
If Y jk = I n , where I n is an identity matrix of dimensions n×n e ϕ jk = h jk (x jk ) for all combinations of j and k, Equation (2) takes the form: where x jk are vectors of size n, h jk (x jk ) is a vector function of covariates x jk or evaluating the function h jk by x jk .Equation ( 3) is called a semiparametric GAMLSS model, because it contains parametric and semi-parametric terms.In this study, when GAMLSS models were fitted, we used the most current version of Smoothing Cubic Spline functions with singular value decomposition (Stasinopoulos et al. 2017).
In the regression fitted with covariates latitude (x 11 ), longitude (x 12 ) and temporal index (x 13 ).The temporal index was constructed as follows: number 1 was assigned to Jan/1989, number 2 to Feb/1989, and so on.Thus, the models fitted for comparison with the trend component are expressed as: In Equation 4, β 10 , β 11 , β 12 , β 13 , β 20 and β 30 represent the parameters to be estimated in the model; μ, σ and ν represent the probability distribution parameters related to location, scale, and shape in the GAMLSS models, respectively.The function g k ( .) is defined as a monotonous link function relating θ k to the explanatory variables and random effects through an additive model.Note that with a Normal distribution in the GAMLSS regression fitted model, we need not include the ν parameter because this distribution includes only location and scale parameters.A major advantage of using GAMLSS regression in the trend component is that it allows to insert different link functions into the linear predictor, and it is not necessary to apply Box-Cox transformations to the response variable.
We used in the temporal index the Smoothing Cubic Spline functions.These functions are extremely useful for modeling the relationships between linear predictors and the response variable, especially when these relationships are of high degree or are non-linear polynomials.Figure 2 shows the behavior of precipitation between the covariates longitude, latitude and the temporal index in the Basin.There is an effect of precipitation as a function of geographic coordinates and a seasonal behavior, justifying the use of spline functions in the modeling.
In this paper, the fit for Equation 1 is performed separately, taking a two-step approach.First, we estimate the trend component and then model the spatio-temporal residuals having spatiotemporal dependence.In the next section, we give the details of modeling these residuals.

Stochastic residual
For a regression fitted for the trend component, we use the maximum likelihood method.It was assumed that the errors are uncorrelated.However, we know that the residuals produced by this fitted model present a spatio-temporal autocorrelation.Thus, to model these residuals, we proposed the fitted model for a valid space-time variogram, considering the spatial, temporal, and spatio-temporal dependence structures of the residuals.The theoretical variogram model was fitted to the empirical variogram model as follows: (5) In Equation ( 5), |L(r s , r t )| is the cardinality for the set L(r s , r t ), and r s and r t are the spatial and temporal distances, respectively.
Several methods are presented in the literature for incorporating a spatio-temporal dependence structure of the stochastic residual (Cressie & Huang 1999, Iaco et al. 2001, Gneiting 2002, Snepvangers et al. 2003).The main difference between space and time covariance models is in forming their structures since some consider space and time to act separately, while others consider them non-separable (Iaco 2010).In the separable models class, the space-time covariance functions can be written as the sum or product of purely spatial and purely temporal components.Although such models present facilities in the computational aspect, it does not allow for space-time interaction, which is the main disadvantage of these models (Montero et al. 2015).Thus, this study proposes a nonseparable type of model approach using the generalized sum-product type spatio-temporal covariance model (Gräler et al. 2016, Medeiros et al. 2019).The generalized sum-product model (γ st (h s ,h t )) it expressed as follows: (6) In Equation ( 6), γ st (h s ,0) and γ st (0,h t ) represent the marginal variogram models of the spatial and temporal effects, respectively.The parameter k represents the magnitude of the effect of the spatio-temporal interaction.

Regression kriging
The choice of the variogram model for each model component in spatio-temporal kriging is of fundamental importance to obtain the best unbiased predictor.This predictor is given by: ( ) where ( ) ε s t is the value obtained using the kriging method for the space-time point ( ) 0 0 , s t ; 0 T c is the vector for size n of the covariance between the residuals of the estimated and observed points; C n is the covariance matrix for the residuals of dimension n × n; and ε is the vector for the residuals of fitting Model (Equation 7).The final predictor (Equation 8) of the variable Z for location ( ) 0 0 , s t is given by: ( ) ( ) ( ) where ( ) ˆ, m t is the trend components estimated value for location ( ) 0 0 , s t .Note that regression kriging consists of the separately fitted trend component and residuals.
By fitting a regression model, we obtain the predicted values for the response variable and the residuals for all the points sampled in spacetime.We also obtain the empirical variogram for these residuals and for the theoretical fitted model.The residuals are then interpolated and added to the estimated trend.In recent years, regression kriging has made theoretical advances, with several applications in real phenomena (Hengl et al. 2012, Heuvelink et al. 2012, Kilibarda et al. 2014).

Validation test
The best structure was determined using the validation test.The method consisted of randomly removing 10 pluviometric stations and then making the prediction of the entire time series of these 10 stations that had 2315 observations.Finally, the data set will have the observed and predicted values.For select the best spatio-temporal variogram model we calculate the Root Mean Square Error (RMSE), Normalized Root Mean Square Error (NRMSE), Nash-Sutcliffe Efficiency (NSE) and Determination Coefficient (R²).After selecting the best spatio-temporal variogram model, we built scatter plots between observed and predicted precipitation values, considering the trend component, spatio-temporal regression kriging and data estimated by the NASA POWER project (Sparks 2018).

Statistical software
We performed all analyses in this research using the software R (R Core Team 2020).For the fitted trend, we used the library GAMLSS (Stasinopoulos et al. 2017).For spatio-temporal regression kriging we used the library gstat (Gräler et al. 2016) and for the construction of the maps the sp library was used (Bivand et al. 2013).

Exploratory analysis
Figure 3 shows the histogram for the variable monthly total precipitation observed in all the 341 points sampled over 31 years.This variable showed a positive skewness, with 21% of the observations equal to zero.The median value of precipitation was 41 mm/month with a standard deviation of 110 mm/month.Additionally, 25% of the observations showed less than 2 mm/month of precipitation.These climatic characteristics of SFRB are associated with the performance and behavior of many atmospheric teleconnection systems and standards influencing the annual average precipitation ranging from 350 mm (Lower-Middle SFRB) to 1500 mm (Upper SFRB), characterizing high spatial climatic variability of the region (Bezerra et al. 2019), as may be seen in Figure 4.
Figure 4 highlights the maps of the median values obtained from the 1989-2019 historical precipitation series for all the 12 months.Generally, the highest rainfall indices are found in the Upper São Francisco region, followed by the Middle region, with lower rainfall indices in the Lower and Lower-Middle São Francisco regions.
The hydrographic regions that divide the basin show different behavior of the spatial variability of precipitation.One justification for this is the climatic differences, especially in  the region of the mouth of the São Francisco River, due to its proximity to the Atlantic Ocean.Upper São Francisco, mainly in the region of the state of Minas Gerais, has the highest altitudes and a climate classified as humid.Middle São Francisco, having the largest area among the hydrographic regions, is located in the semiarid region, presenting lower precipitation indices in relation to the Alto São Francisco and predominating the semiarid climate.Lower-middle is known for having the lowest precipitation indices and with an average annual temperature of 27 ºC, the climate in this region is classified as semi-arid and arid.With its proximity to the Atlantic Ocean, the Lower São Francisco has a milder climate, ranging from semi-arid to humid tropical (Hermuche 2002).
In addition to spatial variability, the SFRB indicates high temporal variability, justifying the application of a model considering spatial, temporal, and space-time interactions.Furthermore, the distribution of rainfall in certain SFRB areas shows monthly variability.

Spatio-temporal geostatistics
Cross-validation was used to choose the best structure, combining the different regression models for the trend with the different theoretical variogram models for the stochastic residual (Table I).In all, regression kriging fitted nine spatio-temporal structures for the residual component and using four probability density functions in the trend component, resulting in 36 adjusted spatio-temporal models.
According to the values of the RMSE, NRMSE(%), NSE and R² statistics, we have the choice of the ZAGA probability distribution for the trend component.For the stochastic residue, the selected adjustment to the generalized product-sum variogram occurred with the Spherical models for the spatial and temporal components.
Table II shows the estimates of the parameters of the generalized product-sum models adjusted to the residues of the GAMLSS model with ZAGA distribution.
The results indicate the existence of a spatio-temporal pattern of the residues obtained from the adjustment of the trend with ZAGA distribution for the accumulated monthly precipitation (Table II).It is noted through the sill parameters that the temporal component is responsible for greater variability compared to the spatial component.In addition, the adjustment revealed a spatial dependence of up to 873 km between the pluviometric stations and temporal autocorrelation of approximately 1.6 months.
The spatio-temporal variogram equation based on Table II is .
Figure 5 shows the empirical variogram ("sample") and the adjustment of the generalized product-sum model ("productSum") for the monthly accumulated precipitation after removing the trend.It is noticed that the residues present a correlation in space and also in time, justifying once again the use of the methodology proposed in this research.
In Figure 6 we can visualize the marginal behavior in relation to space (Figure 6a) and time (Figure 6b) for the empirical spatiotemporal variogram (Figure 5a).Note that the highest values of variances are present in the temporal structure.As expected, when regression kriging is used, the predicted values are close to those observed, providing evidence that the proposed methodology performed well, even surpassing the data estimated by NASA POWER.

Spatio-temporal prediction
As an illustration of regression kriging, Figure 8 shows the predictions of SFRB precipitation values for all months of 2019.This type of approach is essential to detect the possible changes in rainfall distribution in this region, enabling preventive and mitigation activities (Araujo & Celeste 2019).Figure 8. Maps of estimated precipitation in SFRB for all months of the year 2019.
In this study, it was not intended to go deeper into the results of spatio-temporal interpolation, but rather to present an alternative to trend component modeling through GAMLSS models, which presented good results for Regression Kriging.Thus, for future work, it is necessary to present an analysis of constructions of precipitation scenarios in the SFRB basin, since there are studies showing a reduction in precipitation in the basin in the coming years.Marengo et al. (2011) show decreasing rainfall projections, with a 15% of annual precipitation by 2040 in the SFRB, which could go up to 20% in the rainy season (Bezerra et al. 2019, Silveira et al. 2016, Santos et al. 2018).This means a further decrease in volume of recharge of the São Francisco River, leading to climate risks very Table II.Estimates of parameters (Nugget, Sill, Range and k) in the generalized product-sum model and the spatial and temporal dependency index (DI) for the residuals of accumulated monthly rainfall of the São Francisco River Basin, Brazil.The parameter k represents the effect of the space-time interaction and contemplates the global sill.

Model performance
Figure 7 shows the difference between the observed and estimated precipitation with the trend component (Figure 7a), spatio-temporal regression kriging (Figure 7b) and data estimated by the NASA POWER project (Figure 7c).Note the absence of a good fit when only the trend component is used to predict precipitation.shortly, with even more problems in the country's economy due to lack of water resources, mainly for human consumption, agriculture, and energy production in the NEB.

DISCUSSION
Knowledge about the distribution of precipitation in space and time is of utmost importance for scientific knowledge, which is used to understand regional and global changes to processes better involving water and its associated materials.Efficient monitoring of a region's water resources, which takes into account rainfall distribution, can deal with natural disasters such as floods, prolonged periods of drought.In addition, rainfall estimation is of fundamental importance during hydrological modeling, models that have a high level of reliability of these estimates are required, as these can be used to prevent flooding and to assess runoff of water quality and ecosystem health in urbanized areas (Knight et al. 2005).So that we can be aware of the distribution of rainfall throughout an area, it is necessary to use tools capable of estimating values with high quality and resolution, mainly for urban interpolation and in catchment-scale (Hu et al. 2019).Additionally, the spatial variability of rainfall can significantly influences the flow generation mechanisms of urban basins (Segond et al. 2007).
The geostatistical methodology presented in this study can, in principle, be used for urban and catchment interpolation.However, when using analysis for the spatial distribution of precipitation data on a daily scale, geostatistical methods present several difficulties, mainly in the presence of large temporal variability (Javari 2017) and spatial variability, as it increases the complexity of estimates on the uptake scale (Knight et al. 2005).However, once considered that the variable presents an irregular behavior, these difficulties can be considered in the spatio-temporal geostatistics, because this methodology, it is possible to enter autocorrelation structures space-time capable, such as those used in this manuscript, to incorporate modeling the occurrence of high rainfall variability.
In this study, we present a novel method to model the trend component in the space-time stochastic process.GAMLSS regression models (Stasinopoulos et al. 2017) have become a powerful tool in modeling this trend since they allow for the presence of non-parametric terms and use of different probability distribution families in regression adjustment.The modeling of this component is of paramount importance because it is responsible for the large variation in the phenomenon studied.The results demonstrated that GAMLSS regression models are suitable and preferable in regression-fitted models.Additionally, the approach taken in this manuscript for precipitation value interpolation allows for flexible modeling for the trend component, since it is possible to determine different link functions, in particular to the location parameter (mean), no longer need to apply the boxcox type changes in the response variable.In the literature, some studies have used methods involving GLM for trend adjustment, such as (Pebesma et al. 2005) performed a purely spatial interpolation, modeling the trend with Poisson regression as a function of water depth and distance to the coast and residues by ordinary kriging.
Recent studies have suggested the adjustment of these models in phenomena having characteristics that do not follow a normal distribution, such as discrete quantitative data, censored data, and data with excessive zeros in continuous quantitative variables (Stauffer et al. 2016, Sá et al. 2018).However, when modeling the trend component with GAMLSS, independence was assumed between the observations of the response variable for the different spatio-temporal points, this assumption was necessary to obtain the estimates by the maximum penalized likelihood method.However, this drawback was circumvented by using regression kriging, which by adjusting the co-variance function took into account for the autocorrelation of the residuals.It is noteworthy that the regression kriging interpolated only for the location parameter and that in the modeling of the residuals, it was considered a Gaussian random field.A solution to these weaknesses pointed out in our modeling would be to adjust in one-step a GAMLSS regression and consider in the covariance matrix spatiotemporal autocorrelation structures.This type of solution has already been presented for area data, considering regression model with spatial effects (Bastiani et al. 2018).
To apply the methodology proposed in this paper, we used the set of monthly total precipitation data evaluated in the SFRB region space-time.A descriptive analysis of these data found that the precipitation presented a positive skewness distribution, with excessive zeros and irregular distribution of the monthly median value in space and time, indicating the need for spatio-temporal modeling for these data.For these reasons, it was hypothesized that the normal regression model might not be the most adequate one for the trend component in the multiple regression fitted model, thus suggesting the GAMLSS regression modeling with ZAGA, GG, and Normal distributions.The normal regression model has been widely used in the trend component fitted models in indexed space-time phenomena (Gräler et al. 2012, Kilibarda et al. 2014, Hu et al. 2015, Iaco et al. 2015).However, recent studies have suggested adjustment of generalized linear models to evaluate the variables with skewed distributions (Menezes et al. 2016, Monteiro et al. 2017).Thus, this paper's proposal is to consider GAMLSS models with different probability distributions in the trend component fitted model present in the spatio-temporal stochastic process.

CONCLUSIONS
Spatial-temporal regression kriging was used to analyze the SFRB precipitation data, combining the kriging on the residuals with GAMLSS regression, generating better results in comparison with the regression with Normal distribution.A weakness of this research consisted in the non-construction of precipitation projections in the SFRB, which is considered a possibility for works that will follow this research.

Figure 1 .
Figure 1.Map of São Francisco River Basin (SFRB) and distribution of the 341 pluviometric stations.
where p is the number of parameters.Now, function k g (•) is defined as a monotonous link relating θ k to the explanatory variables and random effects through an additive model expressed by:

Figure 2 .
Figure 2. Relation of precipitation observed between longitude a), latitude b) and time c) covariates in the SFRB for the 1989-2019 period.

Figure 4 .
Figure 4. Map of the median value of precipitation in the SFRB for the 1989-2019 period.

Figure 3 .
Figure 3. Histogram of precipitation.The colors in the figure represent the percentages of values equal to or greater than zero in the data set.

Figure 6 .
Figure 6.Empirical marginal variograms for space a) and time b) of precipitation residues in SFRB.

Figure 7 .
Figure 7.Comparison between the observed and predicted precipitation with the trend component a), regression kriging b) and NASA POWER c), in the randomly removing 10 pluviometric stations.

Figure 8 .
Figure 8. Maps of estimated precipitation in SFRB for all months of the year 2019.

Table I .
Results of cross-validation considering the theoretical models of Spherical, Gaussian and Exponential variogram for the Spatial and Temporal components, and the different probability density (PDF) functions for the trend component, being them Generalized Gamma (GG), Zero Adjusted Gamma (ZAGA), Normal GAMLSS (NO) and Linear Model Normal (LM).Metrics to evaluate the error were Root Mean Square Error (RMSE), Normalized Root Mean Square Error (NRMSE), Nash-Sutcliffe Efficiency (NSE) and Determination Coefficient (R²).