Downscaling Daily Satellite-Based Precipitation Estimates Using MODIS Cloud Optical and Microphysical Properties in Machine-Learning Models

: This study proposes a method for downscaling the spatial resolution of daily satellite-based precipitation estimates (SPEs) from 10 km to 1 km. The method deliberates a set of variables that have close relationships with daily precipitation events in a Random Forest (RF) regression model. The considered variables include cloud optical thickness (COT), cloud effective radius (CER) an cloud water path (CWP), derived from MODIS, along with maximum and minimum temperature (Tx, Tn), derived from CHIRTS. Additionally, topographic features derived from ALOS-DEM are also investigated to improve the downscaling procedure. The approach consists of two main steps: firstly, the RF model training at the native 10 km spatial resolution of the studied SPEs (i.e., IMERG) using rain gauge observations as targets; secondly, the application of the trained RF model at a 1 km spatial resolution to downscale IMERG from 10 km to 1 km over a one-year period. To assess the reliability of the method, the RF model outcomes were compared with the rain gauge records not considered in the RF model training. Before the downscaling process, the CC, MAE and RMSE metrics were 0.32, 1.16 mm and 6.60 mm, respectively, and improved to 0.48, 0.99 mm and 4.68 mm after the downscaling process. This corresponds to improvements of 50%, 15% and 29%, respectively. Therefore, the method not only improves the spatial resolution of IMERG, but also its accuracy.


Introduction
In the context of climate change, precipitation patterns around the world have suffered spatial and temporal changes due to the rising temperatures [1][2][3]. In order to assess the impacts of changes in precipitation on the dynamic processes of water resources, ecosystems and society, it is desirable to obtain accurate precipitation estimations [4][5][6]. Traditionally, these estimations have been obtained from rain gauge records. However, rain gauge measurements provide a limited characterization of precipitation, since each station monitors only a small area of the studied region. Therefore, in remote areas, where the rain gauge networks are sparse and unevenly distributed [7,8], their application in the analysis of the spatio-temporal heterogeneity of precipitation entails a certain degree of uncertainty [9,10].
Alternatively, satellite-based precipitation estimates (SPEs) have the capacity to estimate precipitation continuously around the globe on regular grids. In this context, SPEs have become an alternative source of information for various water-related applications, particularly for regions encompassing poorly monitored catchments [11][12][13]. For example, SPEs have been used for flash flood monitoring [14][15][16], streamflow simulations [13,17,18], drought monitoring [19,20] and precipitation monitoring [21][22][23]. The aforementioned investigations demonstrated that, due to their coarse spatial resolution (i.e., several kilometers), SPEs do not provide an effective source of data for local-scale applications, but do for large-scale applications. Indeed, SPEs' low spatial resolution prevents an optimal consideration of climate and/or topographic characteristics influencing precipitation at the local scale.
To overcome this limitation, different studies have reported spatial downscaling procedures. Two techniques have been developed to yield better spatial resolution in SPEs: (i) dynamical and (ii) statistical downscaling.
The dynamic downscaling approach is based on regional climate models, which consist of high-resolution models nested in low-resolution global climate models. This technique relies on deterministic weather models that provide high-resolution precipitation estimates and other climate variables by simulating physical processes of the land-atmosphere system. This approach is known to be limited by the substantial computational resources it requires, as well as its sensitivity to the boundary conditions obtained from global climate models [24][25][26].
The statistical downscaling approach considers a regression analysis between environmental variables and precipitation estimates. The environmental variables involve high-resolution auxiliary datasets as the independent variables of the regression analysis. The precipitation data (i.e., the dependent variable in the regression analysis) are characterized by their low resolution. The resolution of the precipitation dataset is downscaled to a finer resolution using the auxiliary data products through the regression procedure. This data-driven downscaling method requires a large amount of data in order to establish statistical relationships and, unlike the dynamical methods, has lower computational demands [12,[26][27][28].
In the statistical downscaling approach, nonlinear regression methods such as artificial neural networks (ANNs), random forests (RF) and support vector machine (SVM) have become feasible alternatives for downscaling low-resolution datasets due to their high potential for complex and nonlinear input-output mapping [12,24].
Easily accessible through satellite datasets, topography features, land surface temperature (LST), and normalized difference vegetation index (NDVI) are generally used as auxiliary data in statistical downscaling techniques [29,30]. The use of these auxiliary data in statistical downscaling techniques often outperforms simpler interpolation methods (e.g., bilinear interpolation) [31].
Nevertheless, at fine temporal scales (i.e., daily), the use of NDVI can become a source of uncertainty, since the interaction between precipitation and this environmental variable is weak. Similarly, LST estimates are generally unavailable for the days with precipitation events, as cloud cover prevents LST measurement being performed from satellites. Therefore, these approaches are mainly used to downscale precipitation at the monthly time step [30][31][32].
One alternative to overcome the limitations of using auxiliary data that have weak interactions with precipitation at finer temporal scales (i.e., daily) is the use of cloud optical and microphysical property datasets. Indeed, clouds are a prerequisite for the occurrence of precipitation [8], and precipitation is a key factor that controls the amount of clouds on the globe [33]. Campos Braga et al. [34], and Kobayashi and Masuda [33,35], studied the relationship between cloud optical and microphysical properties, such as cloud optical thickness (COT) and cloud effective radius (CER), and precipitation. They found that the CER was a crucial parameter for establishing the initiation and amount of precipitation at the top of convective clouds. Moreover, the changes in COT resulted from changes in the size distribution of cloud droplets which are associated with precipitation.
In light of the context described above, in this study, the potential of cloud optical and microphysical properties, as observed by MODIS, to enhance the spatial resolution of SPEs across the Ramis River basin, the most important tributary of Lake Titicaca, located in the South American Andean Plateau, is investigated. To overcome the computational demand limitations of dynamic downscaling methods, an RF regression model is selected to implement a statistical downscaling approach. MODIS cloud optical and microphysical properties have been previously considered for SPE downscaling for isolated extreme daily precipitation events in central Europe [8]. Therefore, this study will explore the possibility of downscaling continuous daily precipitation for an entire year (i.e., 2012) to provide a step forward toward the construction of continuous daily precipitation observations at high spatial resolution.

Study Area
The Ramis catchment is located in Peru to the northwest of Lake Titicaca, between latitude 14°03′ S and 15°27′ S and longitude 71°07′ W and 69°25′ W, and is the largest tributary of Lake Titicaca [36] (Figure 1). The catchment area is estimated at 15,671 km 2 with an elevation ranging between 3810 m and 5740 m, and a mean streamflow of 73 m 3 ·s −1 [37]. The catchment has a semi-arid-sub-humid climate, with a mean annual precipitation of about 700 mm [38]. The hydrological year is divided in two seasons: (i) the wet season, which is about 5 months long (November to March); and (ii) the dry season, which is about 7 months long (April to October) [6].
The Ramis catchment is the largest tributary to Lake Titicaca [36] which, as the largest freshwater lake in South America, represent a crucial water resource in the central Andes [39]. Along with climate warming, anthropogenic pressure on Ramis water resources has increased in recent decades due to population growth and increased evapotranspiration losses due to agricultural activity [40,41]. Furthermore, the recent extreme drought in 2022-23 in the Lake Titicaca basin also highlights the need for accurate precipitation information for sustainable water resources management. Moreover, the Ramis catchment has a rain gauge network that is denser than that of most of the Plateau, facilitating the calibration and validation of SPEs in order to extend the methodology and results to neighboring regions with similar climate contexts.

Rain Gauges
The Servicio Nacional de Meteorología e Hidrología (SENAMHI) in Peru is in charge of the meteorological network in the study region. Fifty-eight ground stations, collecting data on a daily basis from 1975 to 2017, were available in the study area and subjected to quality control [37]. However, many stations suffer from missing values in their temporal series from few days to complete years. In this context, we pre-selected stations with daily records for the period 2000-present (n = 31), which is the observation period of the considered SPE (IMERG). Then, an analysis of the precipitation records identified 2012 as the year with the highest number of operational stations (n = 28) with more than 80% of daily records available. Sixteen out of the twenty-eight stations were located outside the Ramis catchment, but were used to increase the reference observations, thus extending the study area ( Figure 1). Based on the 28 rain gauges, the regional annual precipitation in 2012 was determined to be 822 mm, with more than 60% of the precipitation amount falling in January, February and December.

Satellite-Based Precipitation Estimates: IMERG
The Global Precipitation Measurement (GPM) mission was launched in February, 2014, as the successor of the Tropical Rainfall Measuring Mission (TRMM) [31,32]. The GPM improves the spatial resolution (from 0.25° to 0.1°), coverage (from 50° N-50° S to 60° N-60° S) and the sensitivity when detecting light and solid precipitation with respect to the TRMM [42]. The Integrated Multi-satellitE Retrievals for GPM (IMERG) datasets are available in three different products: IMERG Early-Run, Late-Run and Final-Run. At the beginning, IMERG precipitation datasets were only derived from GPM satellite observations, and were therefore only available for the 2014-present period. However, in 2020, the algorithm used to retrieve IMERG precipitation estimates was updated to fuse observations from both the TRMM (2000-2015) and GPM (2014-present) satellites in order to obtain long-term precipitation time series from 2000 to present.
IMERG Final-Run includes a gauge-based adjustment that generally provides more reliable precipitation than its Early-and Late-Run counterparts [43]. Therefore, in the current investigation, we adopted the daily product of the IMERG Final-Run, which will hereafter referred to as IMERG. It is worth mentioning here that IMERG daily records are obtained as the sum of 48 half-hour records for the considered day (00:00 to 24:00, GMT). IMERG data were downloaded from NASA (2022) ( Table 1).

Satellite Cloud Cover Properties: MODIS
The Moderate-Resolution Imaging Spectroradiometer (MODIS) instruments are part of the Terra and Aqua platforms. Both platforms have a sun-synchronous near-polar orbit and have a global coverage of every 1 to 2 days [44]. Having twin MODIS instruments makes it possible to improve cloud imaging during the morning and afternoon [8]. Near the equator, Terra acquires cloud data at about 10:30 and 22:30 (local solar time), whereas Aqua acquires cloud data at about 15:30 and 01:30 (local solar time) during a five-minute period. Their Level 2 products provide cloud optical and microphysical properties with spatial resolutions of 1 and 5 km [45], which are observed by Terra (MOD06) and Aqua (MYD06) during both daytime (10:30 and 15:30 for Terra and Aqua, respectively) and nighttime (22:30 and 01:30 for Terra and Aqua, respectively). The cloud optical and microphysical property variables used in the downscaling methods are (i) cloud optical thickness (COT), (ii) cloud effective radius (CER) and (iii) cloud water path (CWP), at a spatial resolution of 1 km. The MOD06 and MYD06 products were pre-processed using the software HEGTOOL [46] in order to obtain COT, CER and CWP datasets with a spatial resolution of 1 km in a suitable computer file format that stores raster graphics (e.g., GeoTIFF) for the year 2012.
For each day, we calculated the mean COT, CER and CWP, as observed by MOD06 and MYD06, before using them in the downscaling procedure ( Table 1). As most of the MOD06 and MYD06 nighttime records are missing, we only consider daytime observations in the daily average process to ensure consistency in the temporal series. As both convective cloudiness and the associated rainy episodes in the Andean Plateau occur mostly in the afternoon and evening [47], MOD06 and MYD06 daytime acquisitions look very adequate for estimating precipitation in this region.

Satellite-Based Temperature Estimates: CHIRTS
The Climate Hazards Group InfraRed Temperature with Station (CHIRTS), developed by the Climate Hazards Center (CHC) at the University of California Santa Barbara and the US Geological Survey (USGS) [48], is a satellite-based product with daily air temperature data. Its coverage ranges from 60° S to 70° N with a spatial resolution of 0.05° × 0.05° (approximately 5 km × 5 km). CHIRTS daily maximum and minimum 2 m air temperature (Tx, Tn, respectively) for 2012 were used. CHIRTS relies on temperature estimates derived from (1) thermal infra-red observations (TIR from GridSat), (2) reanalysis datasets (MERRA2 and ERA5), and (3) meteorological stations. The presence of clouds inevitably affects TIR observations. However, the TIR observations used to produce CHIRTS were previously cloud-screened to limit cloud interference [49]. Moreover, the consideration of other temperature datasets that are not directly influenced by clouds (i.e., reanalysis and meteorological stations) tend to decrease CHIRTS sensitivity to cloud cover interference even more on TIR observations. Finally, it is worth mention here that in comparison to other gridded temperature datasets (CPC, EWEMBI, ERA5, PGF and CPC), CHIRTS provided the most reliable daily temperature estimates across the study region [50]. For more information about CHIRTS daily temperature data, please refer to Funk et al. and Verdin et al. [49,51].
The use of this satellite-based temperature information could be advantageous in downscaling SPEs since the variability of the daily air temperature is linked to precipitation events. For instance, the moisture held by air increases with rising temperatures. As a result, precipitation events may become more intense [52]. This trend is in accordance with the Clausius-Clapeyron relation, which states that if relative humidity remains constant, atmospheric humidity will increase at a rate that follows the dependency of vapor saturation on temperature. This is an important consideration for understanding changes in rainfall events [53].

Digital Elevation Model: ALOS-DEM
Topography is another variable that could have an influence on the atmospheric circulation and precipitation patterns. Precipitation gradients with elevation are common, but they are usually nonlinear (e.g., [54]). The impact of topography on precipitation is greater if a source of moisture (e.g., the Titicaca Lake) is located near the analyzed region [55]. Furthermore, the relationship between precipitation and topography is frequently employed in the statistical downscaling of SPEs [31]. We used the Advanced Land Observing Satellite (ALOS) Digital Elevation Model (DEM) with a resolution of 30 m × 30 m [56] to retrieve elevation and terrain slope (Table 1).

Machine Learning Model: Random Forest
As stated in the introduction, there are two approaches for downscaling SPEs: dynamic and statistical. Due to the substantial computational demand required by dynamic techniques [24,25], the application of a statistical approach (e.g., random forest models-RF) was deemed appropriate to downscale the daily SPEs (i.e., IMERG) in the current study.
RF is an ensemble-learning algorithm composed of several tree predictors. Each tree predictor or decision tree is constructed randomly from a collection of training data [57]. RF models have the advantage of having low rates of overfitting problems, because the prediction is carried out by averaging all decision trees [9]. Moreover, if a large amount of data is to be processed, RF has the advantage of being able to train the model fast and efficiently [26]. In addition, the RF's meta-parameters can be modified in order to optimize the model outcomes. In this study, we used the scikit-learn ensemble RF regressor and the scikit-learn Hyper-parameters tuning estimator [58] to train and optimize the RF models, respectively. Figure 2 illustrates the two-stage procedure followed to train and apply the RF model for downscaling IMERG at 1 km resolution.

Downscaling Setup
The first stage consists of training the RF model at the native spatial resolution of IMERG (i.e., 10 km). To do so, the model inputs (or regressors) IMERG, CHIRTS (Tx, Tn), ALOS-DEM (elevation, slope) and MODIS (COT, CER, CWP) were upscaled to 10 km using a bilinear resampling method. The RF model targets entail the rainfall observations of 20 of the 28 rain gauges randomly selected (Figure 1). During RF model training, the model inputs were shaped only in the 20 cells corresponding to the locations of the 20 rain gauges. The precipitation records of the remaining 8 stations were applied to verify the downscaling outcomes.
The second stage consists of the application of the trained RF model at 1 km spatial resolution. In order to obtain the corresponding outcomes, the model inputs IMERG, CHIRTS (Tx, Tn) and ALOS-DEM (elevation, slope) were resampled to 1 km using bilinear resampling methods. No resampling process was applied to MODIS (COT, CER, and CWP), as its native resolution is 1 km. The model outcomes of the second stage consist of 366 matrices of cells (i.e., the 366 days of the year 2012) that cover the study area and have a spatial resolution of 1 km.
In order to assess the benefits of the different regressors, three RF models were contemplated, which will be called scenarios from now on. Each scenario uses different model inputs (Figure 2). Scenario-1 considers IMERG and MODIS (COT, CER and CWP); scenario-2 considers IMERG, MODIS (COT, CER and CWP) and ALOS-DEM (elevation, slope); and scenario-3 considers IMERG, MODIS (COT, CER and CWP), ALOS-DEM (elevation, slope) and CHIRTS (Tx, Tn). The performance of each scenario is verified consid-ering (i) the 20 pixels containing the rain gauges used for RF model training (i.e., the training dataset), and (ii) the 8 pixels containing the rain gauges not used for RF model training (i.e., the validation dataset).
In addition, IMERG at its native spatial resolution (i.e., 10 km) and at a spatial resolution of 1 km (after the application of a bilinear interpolation) are used as benchmarks to evaluate the different scenarios. Hereafter, these benchmarks will be referred to as IMERG-10 km and IMERG-1 km, respectively.
As it could be argued that using the precipitation time series of neighboring stations could lead to overly optimistic results, as those series are to be correlated, we tested and validated alternative RF models using different time periods rather than different stations. Their inputs were the same datasets applied in scenarios 1, 2 and 3, randomly distributed into a training set consisting of 75% of the observations and a validation set consisting of 25%. These models yielded similar metrics as the proposed methodology. Hence, we chose to keep the proposed method in the next stages of the study.
where and are the values of the scenarios' outcomes or benchmarks and the rain gauge observations, respectively, is the length of the overall data, and, ̅ and ̅ are the average values of the scenarios' outcomes or benchmarks and rain gauge observations.
For each scenario, CC, MAE and RMSE were computed after the application of the RF models at 1 km considering (i) the 20 pixels including the rain gauges used for the RF model training (training dataset), and (ii) the 8 pixels including the rain gauges not used for RF model training (validation dataset).

Categorical Metrics
Categorical metrics were considered to assess the capability of the RF models to detect daily precipitation events. The categorical metrics consider precipitation as a discrete value with four possible situations: both the scenario (or benchmark) and the rain gauge detect a precipitation event (a), only the scenario (or benchmark) reports a precipitation event (b), only the rain gauge detects a precipitation event (c), or neither the scenario (or benchmark) nor the rain gauge report a precipitation event (d). To avoid the inconsistencies that could be produced by light precipitation events, a threshold of 1 mm.day -1 was used in the precipitation estimates to separate rainy from non-rainy days.
The metrics considered in the current work are: (i) the probability of detection (POD); (ii) the false alarm ratio (FAR); and (iii) the critical success ratio (CSI). POD assesses the ability of the scenario (or benchmark) to forecast precipitation events (Equation (5)), whereas FAR evaluates how often the scenario (or benchmark) detects an event when it should not do so (Equation (6)). POD and FAR values vary from 0 to 1, with perfect scores being 1 and 0 for POD and FAR, respectively. CSI computes the ratio between the number of precipitation events correctly detected by the scenario (or benchmark) and the total number of precipitation events (Equation (6)). Values vary from 0 to 1, with a perfect score being 1.
To facilitate the inter-comparison of the scenarios and the benchmarks, the results are presented in the form of a performance diagram [59]. Performance diagrams have previously been used for daily precipitation analysis in Brazil, Bolivia and Pakistan [19,60]. Figure 3 shows the median statistical metrics (i.e., CC, MAE and RMSE) obtained after the application of the RF model (Stage 2), considering both the training and validation datasets. Applying bilinear interpolation to downscale IMERG to 1 km (i.e., IMERG-1 km) led to very slight improvements compared to the use of IMERG data at its native spatial resolution (i.e., IMERG-10 km) (Figure 3). IMERG-10 km and IMERG-1 km CC were higher for the training dataset (0.35, 0.38) than for the validation dataset (0.32, 0.34).

Evaluation of Results from Quantitative Indices
All of the scenarios studied had similar CC magnitudes in both analyses (i.e., considering the training and validation datasets), which indicates that there were no overfitting issues in the trained RF models. It is worth mentioning that the lower statistical scores obtained by each scenario when considering the validation dataset are related to the lower scores overall that were also observed for the benchmarks.
The Taylor diagrams show that the scenarios reached more accurate results than the benchmarks according to the statistical indices (i.e., CC, RMSE and Standard Deviation), as all scenarios are closer to the reference dot. The results of the scenarios are ranked from scenario-3 to scenario-1 (in reverse order). Indeed, scenario-3 achieved the best results for both the training and validation datasets, with the highest CC scores (0.58 and 0.48) and the lowest MAE values (2.1 and 2.48 mm) and RMSE values (3.91 and 4.68 mm). Considering the validation dataset, these correspond to improvements of 50%, 15% and 29% when compared to IMERG-10 km and to improvements of 41%, 12% and 26% when compared to IMERG-1 km.
Scenario-1, which only considers the MODIS and IMERG variables as regressors (COT, CER, CWP, IMERG), also considerably outperformed both benchmarks. For example, when considering the validation dataset, Scenario-1 reached CC and RMSE values of 0.43 and 4.82, whereas IMERG-10 km reached values of 0.32 and 6.60, respectively. This corresponds to improvements of 34% and 27%, respectively. Therefore, even if the temporal mismatch between IMERG and MODIS observations (IMERG relies on 48 half-hour observations, whereas MODIS relies on two five-min observations) can result in a "nonrainy" cloud cover condition, despite detection of a precipitation event by IMERG (before or after the MODIS Terra and/or Aqua observations), the MODIS cloud optical and microphysical properties still constitute very valuable information for downscaling SPEs at the daily time step. All scenarios present higher CC and lower RMSE and MAE than the benchmarks, highlighting the benefits of the proposed downscaling method. Figure 4 shows the statistical metrics (i.e., CC, MAE and RMSE) obtained after RF model application (Stage 2) at the 28 pixels that contain rainfall stations.
When considering the training dataset, the max and min CC for IMERG-10 km were 0.52 and 0.22, and increased to 0.71 and 0.46 in scenario-3 (i.e., the scenario with the best results). This corresponds to improvements of 36% and 109%, respectively. Similarly, the scenarios' RMSE and MAE values were lower than those of the benchmarks (Figure 4). For example, the max and min RMSE (MAE) values decreased from 7.64 mm and 4.09 mm (4.39 mm and 1.71 mm) for IMERG-10 km to 6.65 mm and 2.78 mm (3.86 mm and 1.47 mm) for scenario-3. This corresponds to improvements of 13% and 32% for max and min RMSE, respectively, and to improvements of 12% and 14% for max and min MAE, respectively.
To assess the overall performance of the scenarios, it is convenient to focus on the results obtained with the eight stations not considered in the training procedure (i.e., the validation dataset). The scenario-3 outperforms IMERG-10 km for all the statistical indices. The improvement rates in each maximum/minimum statistical metric were 67%/52% for the CC metric, 24%/16% for the MAE metric, and 31%/33% for the RMSE metric.
It is worth mentioning that scenario-1, which considers MODIS and IMERG variables as regressors (COT, CER, CWP, IMERG) also considerably outperformed the benchmarks. For example, when considering the validation dataset, scenario-1 presented max and min CC values of 0.60 and 0.33, corresponding to improvements of 50% and 32% when compared to IMERG-10 km and to improvements of 39% and 27% when compared to IMERG-1 km. The inclusion of topography features (slope and elevation) in scenario-2 did not improve the precipitation estimates compared to scenario-1. Indeed, the max and min values of the indices for both scenarios were very close. As topography might control the cloud dynamics and accumulation [61][62][63], topography and cloud properties may be redundant, which could explain why scenario-1 and scenario-2 performed similarly.
On the other hand, the consideration of air temperature in scenario-3 led to a significant improvement in precipitation estimates. In accordance with the Clausius-Clapeyron relation, air temperature changed over time according to the climate conditions. Those variations influence the atmospheric humidity and help in the formation of clouds, and thus precipitation events. Therefore, air temperature may be complementary to MODIS cloud optical and microphysical properties [52,53], as it supports the discretization of rainy from non-rainy clouds.  Figure 5 illustrates the precipitation event forecasting capability exhibited by the scenarios and the benchmarks. The POD was higher for the scenarios than for the benchmarks, demonstrating that the scenarios have better forecasting performance. Nevertheless, the FAR values indicate that the scenarios detect some precipitation events when they should not. Despite these shortcomings, the overall performance of the scenarios is comparable to that of the benchmarks, since the critical success indexes (CSI) computed for both groups (i.e., the benchmarks and the scenarios) have similar values (around 0.5).

Evaluation from Qualitative Indices
A recent study in China proposed two-step machine-learning modeling to merge six SPEs with rain gauge information [9]. First, a classification model was trained to identify rainy and non-rainy days, and secondly, a regression model was used to estimate the rainfall amounts of the classified rainy days. The resulting merge datasets outperformed the six studied SPEs in terms of categorical metrics (including POD, FAR and CSI) and substantially alleviated the temporal and spatial biases in the rainfall estimates. This two-step modeling approach could be used as a guideline for further studies aiming to improve the categorical (and quantitative) metrics of the proposed scenarios. Figure 5. Performance diagram used to assess the SPEs' daily prediction performance using the training dataset (a) and the validation dataset (b). Curved lines and straight lines represent the CSI and BIAS values, respectively. The reference black dot represents a perfect statistical score.

Assessment of the Precipitation Spatial Pattern
The climate variability over the study area depends on the easterly flow that transports humid air from the Amazon region (summer) and the westerly flow that comes from the Pacific region (winter). The westerly flow inhibits moisture transport toward Lake Titicaca, and instead provides the study region with dry air [64] The model presented in the current study does not consider the previously mentioned atmospheric circulation, instead focusing on the analysis of precipitation patterns contemplating only the amount of rainfall that fell in 2012. In that sense, the scenarios and the benchmarks allow us to quantify the annual precipitation over the study area ( Figure  6). The spatial pattern of the annual precipitation is more conspicuous in the benchmarks than in the scenarios. Indeed, all scenarios underestimate the amounts of precipitation observed by the benchmarks across the northern Amazonian region, whereas they overestimate the precipitation amounts across the rest of the study area ( Figure 6).
According to the isohyet maps of the Titicaca-Desaguadero-Poopó-Salar de Coipasa system [65], a positive precipitation anomaly can be observed over Lake Titicaca (i.e., there is more precipitation over the lake surface than in the surrounding land). Indeed, Lake Titicaca allows more solar radiation absorption, resulting in higher temperatures on the lake surface than for the surrounding land. The air masses pick up lake moisture while their temperature rises, increasing rainfall above the lake [66]. This strong land-water contrast (in terms of radiation absorption, temperature, and humidity) is a source of error for SPEs, which underestimate precipitation over the Lake Titicaca [6]. This climatic behavior is correctly captured by scenario-3 (Figure 6e), which estimates an annual precipitation of approximately 1000 mm across the lake. The consideration of Tx and Tn in scenario-3 allows a better representation of that behavior and improves the precipitation estimates across Lake Titicaca (i.e., positive rainfall anomaly).

Conclusions
This study assessed the benefits of using cloud optical and microphysical properties derived from MODIS along with temperature estimates derived from CHIRTS as input variables in a random forest (RF) model to downscale daily satellite-based precipitation estimates (SPEs) (i.e., IMERG) from 10 km to 1 km. Unlike others regressors generally used in statistical downscaling approaches (e.g., NDVI datasets), cloud optical and microphysical properties and air temperature have the substantial advantages of providing variables that are closely related with precipitation events at finer temporal scales (e.g., daily).
The results showed that the consideration of MODIS cloud optical thickness (COT), cloud effective radius (CER) and cloud water path (CWP) as auxiliary data provides a suitable alternative for daily SPE downscaling. Indeed, the downscaled precipitation estimates obtained with this dataset improved the reliability of IMERG by 34% and 27% for CC and RMSE, respectively. The additional use of CHIRTS maximum and minimum air temperature estimates (Tx and Tn, respectively) led to additional improvements. The combination of the MODIS and CHIRTS datasets improved the CC and RSME by 50% and 29%, respectively, compared to IMERG at its native spatial resolution (i.e., 10 km). This result suggests that the downscaled precipitation estimates are prone to improvement with the inclusion of a new set of regressors. Along these lines, the consideration of additional variables that are linked to precipitation at the daily time step (e.g., aerosol concentration, total water content, and wind flow patterns) should be considered as additional model inputs for improving the model outcomes. A future challenge will be the use of two RF models, one to detect rainfall events (i.e., a classification model) and another to estimate the amount of precipitation in the detected rainfall events (i.e., a regression model) in order to improve SPEs' daily precipitation forecasting. Funding: This work was funded by the Centre National d'Etudes Spatiales (CNES) in the framework of the IBYS project (Irrigation by satellite). The first author is grateful to the IRD (Institut de Recherche pour le Développement) for its financial support.