Global Downscaling of Remotely-Sensed Soil Moisture using Neural Networks

. Characterizing soil moisture at spatio-temporal scales relevant to land surface processes (i.e. of the order of a kilometer) is necessary in order to quantify its role in regional feedbacks between land surface and the atmospheric boundary layer. Moreover, several applications such as agricultural management can benefit from soil moisture information at fine spatial scales. Soil moisture estimates from current satellite missions have a reasonably good temporal revisit over the globe (2-3 days repeat time); however, 15 their finest spatial resolution is 9km. NASA’s Soil Moisture Active Passive (SMAP) satellite estimates soil moisture at two different spatial scales of 36km and 9km since April 2015. In this study, we develop a neural networks-based downscaling algorithm using SMAP observations and disaggregate soil moisture to 2.25km spatial resolution. Our approach uses mean monthly Normalized Differenced Vegetation Index (NDVI) as an ancillary data to quantify sub-pixel heterogeneity of soil moisture. Evaluation of the downscaled soil moisture estimates against in situ observations shows that their accuracy is better than or equal 20 to the SMAP 9km soil moisture estimates.


Introduction
Soil moisture is a key variable constraining the fluxes between the land surface and atmosphere boundary, and therefore it plays a key role in regulating the feedbacks between the terrestrial water, carbon and energy cycles (Berg et al., 2014;McColl et al., 2017;. Soil moisture partitions the surface energy between latent heat and sensible heat fluxes (Entekhabi et al., 25 1996;Gentine et al., 2007Gentine et al., , 2011Koster et al., 2010). Moreover, plant photosynthesis is regulated by the water available to them through their roots along with atmospheric conditions (Rodríguez-Iturbe and Porporato, 2007;Volk et al., 2000). Finally, soil moisture regulates the surface water fluxes including infiltration and surface runoff generation (Salvucci, 1993;Salvucci and Entekhabi, 1994;Sun et al., 2011). Therefore, there is a need to characterize soil moisture at spatial scales relevant to the representation of land-surface and mesoscale processes in the atmosphere. This can potentially improve representation of 30 evapotranspiration, runoff and precipitation in hydrologic and weather prediction models and result in improved predictive skills (Gedney and Cox, 2003). In addition, knowledge of soil moisture at fine spatial scales is necessary to improve farming practices and optimizing irrigation scheduling.
Soil moisture downscaling can be conducted through different strategies. The first one is to use the synergy of active and passive observations to take advantage of highly accurate but coarse resolution passive observations and active measurements available at higher spatial resolution but not as directly related to soil moisture as passive microwave measurements (Das et al., 2011;Jagdhuber et al., 2015Jagdhuber et al., , 2016Leroux et al., 2016;Montzka et al., 2016;Njoku et al., 2002;Piles et al., 2009;Wu et al., 2017). However, due 10 to the failure of SMAP active instrument in July 2015 these approaches are not directly applicable to SMAP.
Another approach to downscaling is to use ancillary data in combination with coarse scale soil moisture estimates to describe the spatial heterogeneity of soil moisture. Several studies have used this approach to develop soil moisture at fine spatial scales (Chakrabarti et al., 2016(Chakrabarti et al., , 2017Mascaro et al., 2011;Piles et al., 2011Piles et al., , 2014Srivastava et al., 2013;Verhoest et al., 2015).
However, there are limitations in these techniques. Some of them use a linear or non-linear relationship (i.e. a projection) to define 15 the impact of spatial heterogeneity using ancillary data, typically in combination with an energy balance model to relate surface temperature and soil moisture (Colliander et al., 2017a;Merlin et al., 2005Merlin et al., , 2008cMerlin et al., , 2008aMerlin et al., , 2008b. A major issue is that surface temperature at finer spatial scales from satellites cannot be estimated under cloudy conditions. Moreover, surface temperature cannot be used as a signature of soil moisture in energy-limited conditions. The other group of downscaling methods use complex and computationally intensive disaggregation algorithms that are typically unsuitable for global applications. Another group of 20 methods use physical land surface models at fine spatial scales as an information useful for downscaling which adds significant uncertainty to the downscaling due to the errors of the model parameter estimates (Ines et al., 2013;Merlin et al., 2006;Roerink et al., 2000;Shin and Mohanty, 2013). Merlin et al., 2009 andMerlin et al., 2013 also developed a sequential technique to use satellite observations at different scales and sequentially downscale soil moisture observations from a 40km scale to 1km scale and then 100m. Peng et al. (2017) provide a more detailed review of the current status of soil moisture downscaling algorithms and 25 their advantages and disadvantages.
Our approach in this study is to use Neural Networks (NN) to develop a downscaling algorithm for SMAP soil moisture estimates at the global scale. The capability of NNs in learning complex relationships between inputs and target data as well as their quick run-time after training are among the reasons for their popularity in Earth science and remote sensing problems. In recent years, NN has been used to develop soil moisture retrieval algorithms from either passive or active instruments or combination of both 30 (Aires et al., 2012;Alemohammad et al., 2017a;Jiménez et al., 2013;Kolassa et al., 2013Kolassa et al., , 2016Kolassa et al., , 2017Kolassa et al., , 2018Rodríguez-Fernández et al., 2017;Rodríguez-Fernández et al., 2015). They have also been used to retrieve surface turbulent fluxes from remote sensing observations (Alemohammad et al., 2017b;Jiménez et al., 2009).
Our final product is soil moisture estimates at 2.25km spatial resolution with full global coverage every 2-3 days as dictated by the SMAP orbit from April 2015 until end of March 2017 (the first two years of SMAP mission). The 2.25km spatial resolution is 35 chosen since we use SMAP 36km and 9km (¼ of 36km) soil moisture products for training, and developing the scaling relationship.
This relationship is used to estimate soil moisture at the finer 2.25km resolution which is ¼ of 9km. We explicitly emphasize that we assume that the scaling between 36km and 9km is similar to between 9km and 2.25km. Unfortunately, there is no current way to systematically test a potential scale-dependence of such scaling but the systematic comparisons with observations are used as 3 an indirect validation of this assumption. Moreover, we use ancillary data that are available at global scale at all times (monthly NDVI and topographic index) to mitigate issues with current downscaling algorithms.

Datasets
In order to perform the downscaling algorithm, we use three remote sensing based observations: (1) SMAP estimated coarseresolution soil moisture (Section 2.1), (2) Normalized Difference Vegetation Index (NDVI) as a vegetation index that is strongly 5 correlated with soil moisture spatial patterns (Section 2.2) and (3) a topographic index to account for the location of water in the landscape (Section 2.3).
These data are provided on different spatial grids and are re-projected to a common grid for our analysis (Section 3.1). Moreover, in situ soil moisture data that are used for evaluation are introduced in Section 2.4.

SMAP Soil Moisture 10
SMAP satellite estimates surface soil moisture (within the top 5cm of the soil) across the globe with a 2-3 days repeat time (Entekhabi et al., 2010). SMAP uses a passive microwave radiometer in L-band (1.4 GHz) which is known to be sensitive to surface soil moisture and transparent to clouds and atmospheric moisture. Brightness temperature observations from SMAP are used together with ancillary data on vegetation condition and surface temperature to estimate soil dielectric constant from a zerothorder radiative transfer model, commonly referred to as τ-ω model (Jackson, 1993;Jackson and Schmugge, 1991;Kurum et al., 15 2011). Finally, soil dielectric constant estimates are converted to volumetric soil moisture estimates using soil texture data (Chan et al., 2016;O'Neill et al., 2015). Measurements from SMAP are available from April 2015 to present.
In this study, we use two SMAP Level-3 (L3) passive soil moisture products. The first is the version 4.0 SPL3SMP product, which estimates soil moisture on a 36km the Equal-Area Scalable Earth Grid (EASE-grid 2.0) (O'Neill et al., 2016b). This product provides soil moisture estimates for both the morning (6 AM) and evening (6 PM) overpasses of the satellite; however, we only 20 use the morning overpasses in this study since it minimizes observation errors due to Faraday rotation and soil and canopy temperatures are in equilibrium in the morning. The second is the version 1.0 SPL3SMP_E product, which estimates soil moisture with an enhanced spatial resolution posted on 9km EASE-grid 2.0 (O'Neill et al., 2016a). Both of the products use the same brightness temperature observations from SMAP radiometer; however, in the enhanced product spatial resolution is enhanced using the Backus Gilbert interpolation technique Chaubell et al., 2016). While the SPL3SMP_E product is 25 posted on a 9 km EASE 2.0 grid, its native resolution is coarser (~ 33 km). , provides a detailed explanation of how the native resolution and grid spacing of this product are different. Both products are downloaded from National Snow and Ice Data Center (NSIDC) at https://nsidc.org/data/smap/smap-data.html .
To be consistent with the SMAP grid, we re-project our ancillary data to EASE-grid 2.0 at the respective spatial resolutions.
Moreover, our final downscaled soil moisture estimates are at 2.25km spatial resolution on the EASE-grid 2.0 which is nested 30 within the 9km and 36km grids of SMAP L3 soil moisture estimates (Section 3.1).

MODIS NDVI
Plant's photosynthesis and subsequently vegetation greenness are regulated by the moisture available to plants through their roots, as well as atmospheric conditions (Béziat et al., 2013;Kool et al., 2014;Wang et al., 2014). In turn, vegetation modifies its soil moisture environment through changes in the partitioning of evapotranspiration and through precipitation interception (Coenders-35 4 Gerrits et al., 2013;Markewitz et al., 2010). Therefore, surface vegetation cover can be a proxy for the spatial heterogeneity of soil moisture. Vegetation greenness indices such as NDVI are estimated by observations at red and near-infrared frequencies and have relatively high spatial resolution. In this study, we use monthly mean values of NDVI from the MODerate resolution Imaging Spectroradiometer (MODIS) instrument on board Terra satellite. We use version 006 of MYD13A3 product that provides vegetation indices on a 1km spatial resolution at global scale (Didan, 2015). 5 MODIS estimates are provided on a sinusoidal grid, and due to its high resolution, data are divided into 10° × 10° tiles. In order to re-project these data to SMAP EASE-grid 2.0, we used two open source libraries in Python. First, the pyModis library was used to mosaic all the tiles together and generate one global map for each measurement. Next, GDAL library was used to reproject MODIS estimates from its original 1km sinusoidal grid to the EASE-grid 2.0s.

Topographic Index 10
Topographic index (TI) or topographic wetness index is a measure of the soil's saturation tendency given the upstream drainage area and the slope of the local outflow (Marthews et al., 2015). This index is a good proxy for the heterogeneity of soil moisture at the watershed scale and thus could be expected to provide useful information on landscape moisture organization. Here, we use TI at 15 arcsec spatial resolution developed by Marthews et al., 2015. We upscaled the TI values to the EASE-grid 2.0 at 36km, 9km and 2.25km to be used as an ancillary data. 15

International Soil Moisture Network (ISMN)
In situ soil moisture observations from a set of local networks are provided through the ISMN network and are used to independently evaluate the performance of our downscaling algorithm. ISMN collects and standardizes soil moisture observations from several networks around the globe (Dorigo et al., 2011(Dorigo et al., , 2013. The spatial density of probes and their measurement depth varies across different networks. In this study, we only use data from networks that have measurement at a depth of 5cm below  (Albergel et al., 2008;Bell et al., 2013;Calvet et al., 2007;Moghaddam et al., 2010Moghaddam et al., , 201625 Zacharias et al., 2011). Figure S6 shows the spatial distribution of these stations.

Methodology
Our downscaling algorithm is based on an NN approach that uses the two soil moisture estimates from SMAP at 36km and 9km for training, and then retrieves soil moisture at 2.25km spatial resolution using the 9km estimates from SMAP. The NN relates the coarse scale soil moisture and NDVI estimates as well as the fine scale NDVI estimates as input to the fine scale soil moisture 30 estimates as output. To do so, we assume that the scaling relationship between 36km and 9km soil moisture estimates is the same as the scaling relationship between 9km and 2.25km estimates. To the best of our knowledge, this is the first time that the assumption of similar scaling relationship is used to downscale soil moisture, while it has been shown that soil moisture has a fractal scaling property across scales (Famiglietti et al., 2008). In the following the spatial grid setup, data preprocessing and NN training are explained. 35

Spatial Grids
SMAP observations are provided on the EASE-grid 2.0 at different spatial resolutions. The advantage of this grid is the possibility of having nested grids at different spatial resolutions. In this case, SMAP's original radiometer-only soil moisture estimates are at 36km resolution and the enhanced estimates are at 9km resolution which is nested within the 36km one. Since we use the relationship between the 36 and 9km grids (¼ of the original scale) to train our NN algorithm, our final downscaled estimates also 5 have a spatial resolution of ¼ of the 9km product which is used as input in the retrieval step. This results in a spatial resolution of 2.25km on the EASE-grid 2.0. Figure 1 shows the setup of grids used in the training and retrieval steps.

10
To the best of our knowledge, EASE-grid 2.0 at 2.25km has not been used in any other product before. Therefore, we defined the grid using the same parameters as for 36km and 9km, listed in (Brodzik et al., , 2014, other than the grid spacing which is set to 2252.0138025365 m (equals to ¼ of the grid spacing for the 9km grid).
Our NN algorithm estimates soil moisture at a given pixel of the finer resolution grid (the hatched pixel on the black grid in Figure  15 1). To generate the coarse scale soil moisture estimates that are used as input to the NN algorithm, we use a moving window averaging at the coarser resolution grid (the green pixel in Figure 1) that is centered on the target pixel at the finer grid. In the training step, this moving window is 45 km (applied over the 36 km product) and in the retrieval step, the moving window is 11.25km (applied over the 9 km product). All the inputs to the algorithm are averaged at the scale of the moving window from the coarse resolution grid using an area-weighted averaging. The only disadvantage of the moving window technique is that we cannot 20 retrieve soil moisture at the finer scale near coastal regions or places where SMAP does not provide estimates of soil moisture (such as urban areas, ice covered regions and frozen soil conditions). While this is a limitation of our approach, its impact on the spatial coverage of our downscaled soil moisture estimates is very limited geographically. This limitation can potentially be improved by incorporating an extrapolation algorithm that uses NDVI to estimate soil moisture in coastal pixels . Moreover, we used the static water bodies' mask from SMAP 9km soil moisture product and masked out all pixels that 25 have a water fraction of more than 10%. The reason is that high water fraction impacts brightness temperature observations and results in a biased estimate of soil moisture. We also observed higher discrepancies between the SMAP 36km and 9km soil moisture estimates in regions with high water fraction (not shown here). 6

Neural Networks Setup
Our retrieval algorithm is a statistical approach based on NN. NN uses a set of training data to learn the relationship between inputs and outputs without explicitly modelling the physical relationship between them. This is a powerful statistical approach and has been shown to be able to approximate any continuous function (Cybenko, 1989). NN consists of a set of layers and neurons in each layer that mimic the neural system in humans' brain. Each neuron has a weight and bias corresponding to the neurons from its 5 previous layer. The output of a neuron is a linear summation of the weighted inputs plus the bias that goes through an activation function. In each Multi-Layered Perceptron (MLP) NN (Rumelhart et al., 1985), there is an input layer, an output layer and one or more hidden layers in between. Number of neurons in the input layer is equal to the number of inputs provided, and number of neurons in the output layer is equal to the number of outputs given to the NN during training.
The right choice of the number of neurons in the hidden layer and number of hidden layers depend on the complexity of the 10 relationship between inputs and outputs. While having more hidden layers and neurons can potentially increase estimation accuracy, it might result in over-fitting to training data. This means that the NN has lower accuracy when applied to other independent data. To avoid over-fitting, we use the simplest network (i.e., minimum of hidden layers and neurons) that has an acceptable estimation accuracy. Moreover, we use an extensive and representative set of data for training that covers the full range of soil moisture dynamics across different climates and surface conditions. Furthermore, we test the generalization ability of the 15 NN in an independent dataset.
We use the first two years of data from SMAP beginning on April 1 st , 2015 until March 31 st , 2017. For training data, we sample SMAP soil moisture estimates every 10 days at global scale. This results in 46.3 million data points after removing pixels with high percentage of water bodies. For validation, we sample the SMAP data every 5 days with one day lag with respect to the training data at global scale; therefore, training and validation samples are mutually exclusive. This results in 92 million data points 20 after removing pixels with high percentage of water bodies.
During training the weights and biases of each neuron are estimated iteratively by minimizing a mean squared error cost function using a gradient-descent algorithm with back propagation (Hagan and Menhaj, 1994;Rumelhart et al., 1985). For this purpose, training data are divided to three categories of training (60%), validation (20%) and testing (20%). After each iteration on the training category, the validation data are used to check for over-fitting and test data are used to check convergence. When changes 25 in the cost function are smaller than a threshold, training stops, and weights and biases that resulted in the best performance are selected as the parameters of the network. The validation data used here should not be mistaken with the independent validation dataset sampled mutually exclusive to the training dataset. Throughout the paper, we will use "validation data" terminology to refer to the dataset sampled mutually exclusive to the training data (explained in the previous paragraph) and will be used to evaluate the performance of the NN training in terms of over-fitting. 30 We trained a set of networks with number of hidden layers between 1 and 5 and number of neurons in each hidden layer between 1 and 15. After evaluating their performance, we selected the network with one hidden layer and 5 neurons, which has a good accuracy, while adding more hidden layers or neurons did not change the performance. We use the hyperbolic tangent sigmoid function as the activation function for the hidden layer, a standard practice, and a linear function as the activation function for output layer. 35 7

Inputs and Downscaling Schemes
We develop four different schemes with increasing number of ancillary inputs for downscaling SMAP soil moisture estimates. In this section, we introduce each of them, and we compare their performance in Section 4. At the end, only one scheme was selected to provide the final downscaled soil moisture estimates. Scheme R1 has the three following inputs: soil moisture estimates from SMAP on the moving window grid at the coarse scale 5 resolution (45km for training and 11.25km for retrieval), mean NDVI at the coarse scale moving window, and NDVI in the target pixel at the fine scale grid (9km for training and 2.25km for retrieval).
Scheme R2 has all the inputs in scheme R1 plus standard deviation of NDVI at the fine scale pixels within the moving window.
For example, in the training step there are 25 9km pixels within the 45km moving window. We estimate the standard deviation of those, and input that to the network. This estimate provides a proxy of the heterogeneity within the coarse scale grid. 10 Scheme R3 has all the inputs in scheme R1 plus TI at the moving window of the coarse scale and TI in the target pixel at the fine scale grid.
Scheme R4 has all the inputs from schemes R1-R3. Table 1 lists all the four schemes and their respective inputs. For each scheme, the 9km SMAP soil moisture estimates are used as the NN target data.

Results
In this section, we first present the results of NN training and evaluation of the downscaling from 36km to 9km. We only present the results from Scheme R1, since all four schemes have similar performances. Next, we apply the downscaling scheme to the entire 2-year SMAP data record and generate soil moisture at 2.25km spatial resolution. Finally, we evaluate the accuracy of downscaled soil moisture from all the four schemes using in situ soil moisture estimates from the ISMN dataset. 20

Evaluating NN Training
To evaluate the success of the NN training, we compute the correlation coefficient (R 2 ) and the unbiased Root Mean Square Difference (ubRMSD) between the NN estimates and the target data for the training and validation data. We compare the metrics aggregated across all the data and spatially at each pixel. Figure 2 shows the density scatterplot of 9km NN estimates versus target 9km soil moisture data in the training and validation 5 datasets. Both datasets have similar ubRMSD and R 2 which shows that the NN is able to generalize beyond the training data. R 2 of ~0.98 is an almost perfect correlation between the target and estimates, and indicates that the NN setup with the inputs provided is capable of learning the relationship between coarse and fine scale soil moisture. In order to compare the performance of NN with other techniques, we present similar results for two other approaches. The first approach assumes a uniform distribution of soil moisture (i.e. no heterogeneity) within each 36km grid pixel. This means that all the 9km pixels that fall within a 36km grid pixel are attributed the same soil moisture value as the 36km one (i.e. no downscaling). 15 The second approach is to use a linear interpolation (i.e. a crude downscaling). This is used as a reference to a simpler approach than the NN. In this approach, soil moisture values on the 36km grid are assumed to be at the center point of the pixel. Then using a linear interpolation, soil moisture is estimated at the center point of each 9km grid pixel. We implement the interpolation at the global scale for each SMAP observation. Figure 3 shows the density scatterplots for these two approaches. As the scatterplots and performance metrics show, the NN algorithm has higher correlation and lower ubRMSD compared to both no heterogeneity and 20 interpolation approaches. As expected, the homogeneous algorithm has the worst performance. To further evaluate the performance of the NN downscaling algorithm, we evaluate the correlation between target soil moisture and NN estimates at 9km in each pixel (Figure 4). These estimates are based on the validation dataset, and on average 50 data 5 points are used to calculate the correlations (in some snowy regions like Himalaya Mountains as few as 20 data points were used, the minimum number of samples set to estimate the correlation in this evaluation). Correlations range between a minimum of 0.6 in dense tropical forests (where NDVI saturates and limited spatial inhomogeneity can be found) to ~1 in most other parts of the world, except in very northern regions were the correlation is closer to 0.9. The smaller value of correlation in the tropics is a result of the low seasonality and saturation effect of NDVI in those regions, so that NDVI does not provide useful information for 10 downscaling (Morton et al., 2014). These regions have relatively high soil moisture and relatively low variability in NDVI throughout the year. In addition, NDVI and EVI tend to saturate at high vegetation cover so that there is only limited variability in those regions. As a result, disaggregating soil moisture using NDVI has a lower accuracy compared to other regions. Nevertheless, NDVI and EVI are the only global fine scale vegetation indices that can be used as ancillary data to describe the sub-pixel heterogeneity of soil moisture, and both lack seasonality in the dense tropical forests. Moreover, Figure S1 shows the percentage 15 difference between the NN retrieval and SMAP estimates at 9km. Largest differences are observed in mountainous and coastal regions where soil moisture variations are influenced by factors such as the terrain slope. However, inclusion of other ancillary data such as slope in other downscaling schemes did not improve the bias and schemes R2-R4 have similar performance both in terms of aggregated R 2 and ubRMSD and spatial patterns of correlation between NN estimates and target soil moisture.
Moreover, Figures S2-S5 show the spatial correlation maps and percentage difference between SMAP estimates at 9km and the 20 estimates from Interpolation and No Heterogeneity approaches. These figures show that both approaches have significantly higher bias with respect to the SMAP estimates compared to the NN retrieval.

Downscaled Soil Moisture
Using the NN that was trained and evaluated in the previous section, we then estimate soil moisture at 2.25km EASE-grid 2.0 resolution. For this, we use the SMAP 9km observed soil moisture as input together with ancillary data for each scheme as described in Section 3.3.
10 Figure 5 shows the average global soil moisture at three different spatial resolutions (36km and 9km from SMAP observations, and 2.25km downscaled from SMAP observations). Overall, the spatial patterns are as expected, with some heterogeneity being added as the resolution increases. These maps show that our downscaling algorithm preserves the large-scale spatial patterns of soil moisture while disaggregating it at local scales. Moreover, the latitudinal average plots (on the right side of each panel of Figure 5) show that at higher spatial resolutions there is more spatial heterogeneity. The latitudinal average for the 36 km product 5 is much smoother than the 2.25 km one. Figure 6 shows a snapshot of soil moisture estimates from SMAP 36km product, SMAP 9km product and the downscaled product at 2.25km across western Africa from Oct. 18 th , 2015, in order to further demonstrate the added information of the downscaling algorithm. From left to right in Figure 6 with increasing spatial resolution, finer spatial patterns are revealed, with finer structure, while respecting the overall 9km structures. We note that there is no random spatial noise in the higher spatial resolution data 10 emphasizing that the additional NN input (NDVI) is useful to preserve the spatial connectivity of different high-resolution pixels.

Large Scale Evaluation of Downscaled Soil Moisture
In this section, we analyze large-scale patterns of the downscaled soil moisture estimates at 2.25km spatial resolution. Our first 5 analysis is the spatial heterogeneity of the downscaled soil moisture. We calculate coefficient of variation (CV), which is the spatial standard deviation divided by the mean, of the 16 2.25km pixels within each 9km pixel at each time. Figure 7 middle panel shows the mean CV at each pixel for the downscaled soil moisture. For comparison, we also calculate CV for the 9km soil moisture estimates from SMAP at the 36km grid (Figure 7 top panel), and CV for the 2.25km downscaled estimates at the 36km grid ( Figure   7 bottom panel). Panels in Figure 7 have different range of CV which is expected given the difference in their spatial scales, but 10 inset PDFs have similar scales for better comparison. These plots reveal similar spatial patterns across different scales, yet with finer and contiguous structures at higher resolution. This means that in regions that have high spatial heterogeneity in the 9km SMAP estimates, we also observe high heterogeneity in the 2.25km downscaled soil moisture estimates. Meanwhile, the 2.25km estimates have higher variation (i.e. heterogeneity) at the 36km grid scale compared to the 9km product. This indicates that the 2.25km product is not just a smoothed version of the 9km product but can add spatial variability not captured by the coarser 9km 15 pixels. Therefore, our NN algorithm is appropriately explaining the spatial variability of soil moisture using NDVI as ancillary data. The other downscaling schemes (R2-R4) exhibit similar spatial patterns for CV. This is an indication that the other inputs to these schemes do not provide any extra information on the heterogeneity of soil moisture at the finer spatial scale.
Our second analysis evaluates the conservation of the water balance within each coarse scale pixel. Since we are disaggregating soil moisture from 9km grid to 2.25km grid, the downscaling algorithm may not preserve the mean soil moisture within each coarse 20 scale pixel (water balance). Figure 8 shows the spatial map of the percentage average difference (with respect to their mean) between the mean soil moisture from the 16 2.25km pixels within each 9km pixel and soil moisture estimates at 9km from SMAP.
The inset PDF also shows the probability distribution of these differences. On average, there is a positive bias of ~1.5% in the differences, and the downscaled soil moisture estimates at 2.25km are slightly positively biased, i.e. more humid. However, the differences are small (mostly less than 3% in absolute value) and within the SMAP retrieval mission accuracy (<0.04 m 3 m -3 ). This 25 analysis also reveals that our downscaling algorithm is robust in preserving the water balance within each coarse scale pixel and still within the original requirements of the SMAP mission.

Comparison against ISMN Data
In this section, we evaluate downscaled soil moisture estimates against ISMN dataset. We calculate R 2 , anomaly R 2 and ubRMSE for each of the four downscaling schemes as well as SMAP estimates at 36km and 9km. Anomalies are calculated based on a 30 5 day moving average window. We use stations from 10 networks within the ISMN dataset. For networks that have more than one station within each SMAP pixel, we average the stations within the 2.25km pixels and then calculate the metrics for each of the 36km, 9km and 2.25km products. Figure 9 shows the average of each metric across stations for different networks and the average among all of the networks.
In general, all four schemes have quite similar performance within each network and across all networks. During the NN training, 10 all four schemes had similar performance. In some cases, scheme R1 had slightly better performance (i.e. ~2% higher R 2 ). For these reasons, and to reduce complexity of the downscaling algorithm, our final downscaling algorithm is scheme R1.
Our downscaling scheme always has equal or better performance compared to the 9km product across individual networks. Since the 9km product is the input to our downscaling, this shows that in terms of temporal correlation and ubRMSE our downscaling algorithm either improves the 9km product or is similar in terms of accuracy. In comparison to the 36km product, our downscaling 15 algorithm's performance follows the 9km product performance. At stations where the 9km product has better performance compared to the 36km product (e.g. ubRMSE in iRON), the 2.25km product also has better performance compared to the 36km product, and vice versa. The main added value of the 2.25km product is its enhanced spatial resolution that provides a measure of soil moisture heterogeneity at finer spatial resolution using satellite observations while having an acceptable temporal accuracy.

Conclusions 5
In this study, we have developed a statistical disaggregation algorithm using neural networks to downscale soil moisture observations from SMAP to a spatial resolution of 2.25km. We use the two level-3 soil moisture estimates from SMAP at 36km and 9km spatial resolutions to train our downscaling algorithm. Assuming the relationship is consistent between scales, we apply the downscaling algorithm to 9km soil moisture estimates from SMAP and disaggregate those to 2.25km. Our downscaling algorithm uses only mean monthly NDVI estimates as ancillary data at both the coarse and fine spatial scales to estimate the 10 heterogeneity of soil moisture at sub-pixel level. Our investigation shows that topographic index and variability of NDVI within the coarse scale pixel do not provide additional useful information on the spatial heterogeneity of soil moisture for the downscaling using our proposed NN technique and are thus omitted in the final product. This is the first study to estimate soil moisture at global scale at very fine spatial resolution (~2km). This new product can be used in a range of other studies and applications including land surface -atmosphere interaction modeling, evapotranspiration modeling and agricultural management. However, in some 15 parts of the world that farm sizes are smaller than the resolution of this product, further disaggregation of soil moisture estimates is needed. Moreover, we will extend these estimates to near real time as long as SMAP soil moisture estimates are available.
Our evaluation shows that the performance of the downscaled soil moisture estimates is equal or better than that of the SMAP 9km estimates in terms of temporal correlation, anomaly correlation and ubRMSE when compared to in situ soil moisture estimates 16 from ISMN. Moreover, averaging the estimates at 2.25km within each 9km pixel shows that our downscaling algorithm has high accuracy in preserving the water balance (less than 5% error, within the SMAP mission requirements).
In this study, we use the relative value of NDVI (in a given pixel with respect to the neighboring pixels) as an auxiliary information to predict spatial variability of soil moisture in each coarse-scale pixel. While the relationship between soil moisture and NDVI at phenological time scales may not be valid at the temporal scale of SMAP observations (couple of days), our assumption builds on 5 the relative value of NDVI within a small region and not the absolute value. Therefore, it is reasonable to use NDVI as a predictor in this case. Use of NDVI as an ancillary measurement to disaggregate soil moisture builds on the assumption that there is a moderate vegetation cover in the pixel of interest. Therefore, this lowers the quality of the downscaling algorithm in bare soil or sparsely vegetated regions. Figure 8 shows a dry bias in arid regions which is an indication of lack of NDVI -soil moisture relationship. Moreover, NDVI estimates tend to saturate in highly dense vegetated regions such as in tropical forests, which results 10 in limited ability of the downscaling algorithm to resolve the heterogeneity of soil moisture in those regions.
This study shows the potential of using neural networks with a large number of training data to develop a downscaling algorithm for soil moisture. Data generated using this algorithm can be used to provide an improved understanding of the dynamics of soil moisture and land surface atmosphere interactions at global scale.