Ensemble Downscaling of Soil Moisture Data Using BMA and ATPRK

Abstract

Soil moisture is essential information for meteorological and hydrological analyses. To date, many efforts have been made to achieve the two goals for soil moisture data, i.e., the improvement of accuracy and resolution, which is very challenging. We presented an ensemble downscaling method for quality improvement of gridded soil moisture data in terms of the accuracy and the spatial resolution by the integration of BMA (Bayesian model averaging) and ATPRK (area-to-point regression kriging). In the experiments, the BMA ensemble showed a 22% better accuracy than the data sets from ESA CCI (European Space Agency-Climate Change Initiative), ERA5 (ECMWF Reanalysis 5), and GLDAS (Global Land Data Assimilation System) in terms of RMSE (root mean square error). Also, the ATPRK downscaling could enhance the spatial resolution from 0.25° to 0.05° while preserving the improved accuracy and the spatial pattern of the BMA ensemble, without under- or over-estimation. The quality-improved data sets can contribute to a variety of local and regional applications related to soil moisture, such as agriculture, forest, hydrology, and meteorology. Because the ensemble downscaling method can be applied to the other land surface variables such as temperature, humidity, precipitation, and evapotranspiration, it can be a viable option to complement the accuracy and the spatial resolution of satellite images and numerical models.

1. Introduction

Recent climate change has increased the magnitude and frequency of drought, flood, and wildfire. For the management of such disasters related to water, it is crucial to understand the characteristics of soil moisture and quantify it precisely. Because soil moisture is essential information for meteorological and hydrological analyses (Western et al., 2002; Ochsner et al., 2013), the measurement of regional soil moisture is also necessary for studies of water resource and disaster (Kim and Lee, 2004), in addition to the remote sensing for extensive areas. In-situ measurement of soil moisture is fast and accurate by the TDR (time domain reflectometry) sensors buried underground. However, the number of measurement points is not enough to ensure the spatial continuity of data.

On the other hand, the PMW (passive microwave) satellite sensor that covers vast areas on the Earth can retrieve the soil moisture using the relationship between soil dielectric constant and water content (Schmugge et al., 2002; Mohanty et al., 2017; Peng et al., 2017). Global data sets of soil moisture have been provided by the ASCAT (Advanced Scatterometer) onboard the MetOp-A (Meteorological Operational Satellite–A) operated by ESA (European Space Agency) at the resolution of 25 to 50 km since 2006 and by the MIRAS (Microwave Imaging Radiometer with Aperture Synthesis) onboard the SMOS (Soil Moisture and Ocean Salinity) operated by ESA at the resolution of 35 to 50 km since 2009. Also, the AMSR2 (Advanced Microwave Scanning Radiometer 2) onboard the GCOM-W1 (Global Change Observation Mission– Water) operated by JAXA (Japanese Aerospace ExplorationAgency) has been producing soil moisture data at the resolution of 10 to 25 km since 2012. In 2015, NASA (National Aeronautics and Space Administration) had launched a next-generation satellite SMAP(Soil Moisture Active Passive)that can produce 3-km and 36-km data using a radar and a PMW sensor, respectively. However, the 3-km soil moisture data was not provided finally because of the malfunction of the radarsensor. Despite the many efforts, the accuracy of the PMW observation for soil moisture is not so satisfactory yet. The PMW detects the land surface or the top 5 cm of soil column while the TDR sensor is buried at a depth of 10 to 30 cm for consideration of the root zone of vegetation. Such difference can lead to the inconsistency between the satellite observation and the in-situ measurement (Kerr, 2007; Collow et al., 2012; Kim et al., 2018).

Another effort for more reliable soil moisture data is associated with the meteorological reanalysis system such as ERA5 (ECMWF Reanalysis Version 5) (Hersbach et al., 2018) by ECMWF (European Centre for Medium-range Weather Forecasts) and GLDAS (Global Land Data Assimilation System)(Rodell et al., 2004), NLDAS(North American Land Data Assimilation System) (Mitchell et al., 2004; Xia et al., 2012), and MERRA (Modern-Era Retrospective Analysis for Research and Applications)(Rienecker et al., 2011) by NASA. The reanalysis systems combine a physics based model with in-situ measurements to produce reanalyzed soil moisture data at the resolution of 30 to 50 km. They provide a globally-covered reanalysis every 3 hours with the RMSE (root mean square error) of approx. 8% (Kimet al., 2017). Although the temporal resolution is excellent, the accuracy and the spatial resolution are not sufficient for local-scale applications in small countries like Korea. In this sense, the improvement of both accuracy and spatial resolution is crucial for reliable soil moisture data.

For the accuracy improvement of gridded data, a BMA (Bayesian Model Averaging) method, which produces a weighted ensemble using the posterior probability of individual data members, has been effectively applied in the applications of numerical weather prediction (Raftery et al., 2005; Min et al., 2007; Kim et al., 2016). For the blending of soil moisture data, however, few studies have employed the BMA ensemble (Kim et al., 2018). Moreover, the effects of the accuracy improvement by the BMA were not combined with the enhancement of spatial resolution and vice versa.Although not integrated with accuracy improvement, the studies for the statistical downscaling of soil moisture have been conducted using a linea rmode lwith multiple land surface variables (Chauhan et al., 2003; Ray et al., 2010; Merlin et al., 2013; Song et al., 2014). Also, AI(artificial intelligence) techniques were recently adopted for the downscaling of soil moisture (Im et al., 2016; Liu et al., 2017). However, the statistics-based approaches have the residual problem, so the data patterns before and after downscaling could be different. To overcome this, we need a residual correction by kriging to be integrated into the statistical downscaling, which could maintain the spatial pattern of original data (Jin et al., 2017; Kim et al., 2017), but the accuracy improvement was not accompanied.

Thus, we propose an ensemble downscaling method to achieve the two goalsforsoil moisture data, i.e., the improvement of accuracy and resolution. We aim to improve the accuracy of soil moisture data using the BMA ensemble for the three gridded data sets and enhance the spatial resolution byATPRK (area-to-point regression kriging) that combines the MLR (multiple linear regression) and ATPK (area-to-point kriging)for a data consistency after downscaling (Kyriakidis, 2004; Park, 2013). We first made sure the data accuracy by creating a BMA ensemble from the soil moisture products of the ESA CCI (European Space Agency– Climate Change Initiative), GLDAS, and ERA5 at the resolution of 0.25°. Then we conducted a spatial downscaling byATPRK using the land surface variables such astemperature, precipitation, vegetation greenness, and evapotranspiration, which produced 0.05° soil moisture data from the BMA ensemble through a regressionmodel accompanied by the residual correction. Our experiment focused on South Korea that has periodical TDR measurements of soil moisture. We used a 10-day composite of soil moisture from April to September between 2014 and 2018 by considering data availability and the main seasonsfor hydrological applications.

2. Theoretical background

1) The blending of multi-source data with BMA

BMA is a statistical method to mitigate the uncertainty problem in the blending of multi-source data by setting up the weights for individual members using a posterior probability (Hoeting et al., 1999). The posterior probability is the conditional probability for the value of an individual member to be the true value, i.e., the suitability of an ensemble member (Raftery et al., 2005). BMAis a mixture model using a weighted combination of multiple PDFs (probability density functions) (Claeskens and Hjort, 2008). In actual implementations, a simplified method using a weighted average approximated from the PDFs is preferred for the sake of convenience (Dempster et al., 1977; Kim et al., 2016). The BMA approximation can be expressed as

\(p\left(y \mid f_{1}, \ldots, f_{K}\right)=\sum_{k=1}^{K} w_{k} g_{k}\left(y \mid \tilde{f}_{k}\right)\)       (1)

where y is the target variable, and \(\tilde{f}_k\) is the bias-corrected value of each member, which affects the conditional PDF \(g_{k}\left(y \mid \widetilde{f}_{k}\right)\). w_k is the weight for each member determined by EM(expectation–maximization)algorithm. We can assume that the errors of \(\tilde{f}_k\) follow a normal distribution, and the conditional PDF \(g_{k}\left(y \mid \widetilde{f}_{k}\right)\) can be approximated by a normal distribution centered on \(\tilde{f}_k\) with the variance σ² (Raftery et al., 2005). Then, the BMA approximation is the conditional expectation of y:

\(E\left[y \mid f_{1}, \ldots, f_{K}\right]=\sum_{k=1}^{K} w_{k} \tilde{f}_k\)       (2)

The weight (w_k) and variance (σ²) of the ensemble members can be estimated by maximizing the log-likelihood from the training set:

\(l(\theta)=\ln \left(\sum_{k=1}^{K} w_{k} g_{k}\left(y \mid \widetilde{f}_{k}\right)\right)\)       (3)

where θ is a parameter vector for the weight and variance of the \(\tilde{f}_k\) to maximize the log-likelihood function l. The EM algorithm is an efficient procedure to derive the maximum likelihood estimator (Dellaert, 2002) by the iterations between an expectation (E) step, which creates a function for the expectation of the log-likelihood evaluated using the current estimate for the parameters, and a maximization (M)step, which computes parameters by maximizing the expected log-likelihood found on the next E step.

2) ATPRK downscaling by the integration of MLR and ATPK

Downscaling is a procedure for the creation of fine-resolution data from coarse-resolution data. It can be divided into two broad categories: (1) dynamical downscaling that needs a model-based computation within a fine-resolution grid (Ribalaygua et al., 2013; Peng et al., 2017) and (2) statistical downscaling that facilitates statistical relationships between the original data at a coarse resolution and the multiple explanatory variables at a fine resolution (Park, 2010; Gevaert et al., 2016; Mandal et al., 2016).We built an MLR model using the BMA ensemble at a coarse resolution as a response variable, and the land surface variables aggregated from the fine resolution into the coarse resolution as explanatory variables. However, any statistical model has inevitable residuals, i.e., the differences between original data and estimated values, so the residuals occurring from the MLR model at the fine resolution can also cause the value inconsistency between the original and the downscaled data. To overcome this, we need an accurate residual grid at the fine resolution of which spatial pattern is in accordance with the residual grid at the coarse resolution.Based on the assumption that the residual is an intrinsically stationary process, the ATPK creates an interpolated residual grid at the fine resolution. Residual is a deficit that was not explained by an MLR model. So, the values identical to the original data can be derived by adding the residual to the MLR estimate. Such a residual correction is crucial to the data consistency in downscaling.

The ATPK is a delicate statistical process for the creation of a fine-resolution grid by a full consideration of the spatial association within a coarse-resolution pixel. Fine-resolution residual R(u_j) can be expressed as the multiplication of the coarse-resolution residual R(v_i) and the kriging weight λ(u_j):

\(R\left(u_{j}\right)=\sum_{i=1}^{N} \lambda_{u_{j}} R\left(v_{i}\right), \quad v_{i} \in u_{j} \text { and } \sum_{i=1}^{N} \lambda_{u_{j}}=1\)       (4)

where the kriging weight \(\lambda_{u_{j}}\) is calculated as follows (Kyrialkidis, 2004; Park, 2013; Wang et al., 2015)

\(\begin{array}{} \left[\begin{array}{cccc} c\left(v_{1}, v_{1}\right) & \cdots & c\left(v_{1}, v_{N}\right) & 1 \\ \vdots & \ddots & \vdots & \vdots \\ c\left(v_{N}, v_{1}\right) & \cdots & c\left(v_{N}, v_{N}\right) & 1 \\ 1 & \cdots & 1 & 0 \end{array}\right]\left[\begin{array}{c} \lambda_{u_{1}} \\ \vdots \\ \lambda_{u_{N}} \\ \mu \end{array}\right]=\left[\begin{array}{c} c\left(u, v_{1}\right) \\ \vdots \\ c\left(u, v_{N}\right) \\ 1 \end{array}\right] \\\sum_{i=1}^{N} \lambda_{u_{j}} \bar{c}\left(v_{i}, v_{\mathrm{j}}\right)=\bar{c}\left(u, v_{j}\right) \end{array} \)       (5)

The term \(\bar{c}\left(v_{i}, v_{j}\right)\) is the coarse-to-coarse pixel covariance, i.e., between coarse pixels centered at v_i and v_j. \(\bar{c}\left(u, v_{j}\right)\) is the fine-to-coarse pixel covariance, i.e., between fine and coarse pixels centered at u and v_j . The term μ is the Lagrange multiplier (Goovaerts, 2008). To derive the two types of covariance (coarse-to-coarse and fine-to-coarse pixel covariance), we first need the fine-to-fine pixel covariance. However, the fine-to-fine pixel covariance cannot be directly created because the residual values exist only at the coarse resolution. For an optimal derivation of the fine-to-fine pixel covariance from the known coarse-to-coarse pixel covariance, a deconvolution method was employed (Goovaerts, 2008). More specifically, after setting up an initialfine-to-fine pixel covariance, the theoretically regularized fine-to-fine pixel covariance is compared with the known coarse-to-coarse pixel covariance. By calculating the difference between the two covariance models, the parameters for the fine-to-fine pixel covariance are adjusted to minimize the difference criteria; this procedure is repeated until the difference criterion reaches a given tolerance (Fig. 1). In this sense, the ATPK can be interpreted as a coarse-to-fine pixel kriging for gridded data. A few recent studies have begun to employ the ATPRK for downscaling of the spectral reflectance of satellite images (Wang et al., 2015;Wang et al., 2016) and the satellite products such as precipitation and temperature (Park, 2013; Xu et al., 2020).

Fig. 1. The procedure of ATPK (area-to-point kriging) to derive the residuals at a fine resolution.

3. Data and methods

1) Soil moisture data

(1) ESA CCI soil moisture

The ESACCI SM v04.7, which is created using the data from three AMW (active microwave) sensors (AMI-WS, MetOp-AASCAT, and MetOp-BASCAT) and seven PMW sensors (SMMR, SSM/I, TMI, WindSat, AMSR-E, AMSR2, and SMOS), has three types of soil moisture products (ACTIVE, PASSIVE, and COMBINED) (Hollmann et al., 2013). The depth of the soil layer can be different according to the sensors, but the daily information of the topsoil(approx. 0 to 5 cm)is provided at the resolution of 0.25° (https:// www.esa-soilmoisture-cci.org/). We obtained level 3 COMBINED daily soil moisture and aggregated it into a 10-day composite.

(2) GLDAS soil moisture

GLDASis a reliable assimilation system that integrates and a vast amount ofin-situ observations with multiple LSM (land surface model) such as Mosaic (Koster and Suarez, 1996), VIC (Variable Infiltration Capacity) (Liang et al., 1994), Noah (Ek et al., 2003), and CLM (Community Land Model)(Dai et al., 2003)(Rodell et al., 2004; de Jeu et al., 2012). We used GLDAS Noah V2.1 and extracted the soil moisture variable at the topsoil layer(0 to 10 cm) at intervals of 3 hours and the resolution of 0.25° (https://disc.gsfc.nasa.gov/datasets/ GLDAS_NOAH025_3H_2.1/summary), and aggregated the data into a 10-day composite.

(3) ERA5 soil moisture

ERA5 is a state-of-the-art reanalysis produced by ECMWF, following up on its predecessor ERA-Interim (Hersbach et al., 2020). ERA5 soil moisture data is obtained by running the revised TESSEL(Tiled ECMWF Scheme for Surface Exchange over Land) model (Balsamo et al., 2015). The ERA5 hourly reanalysis on the surface level includes the volumetric soil water at the resolution of 0.25°. We extracted the topsoil layer (0 to 7 cm) variable from the CDS (Climate Data Store) of Copernicus Climate Change Service (https://cds. climate.copernicus.eu/) and aggregated the data into a 10-day composite. ERA5 soil moisture data sometimes has a noisy value because of the mixed-pixel problem along the coastline, unlike GLDAS (Kim et al., 2018). So, we compared the ERA5 and GLDAS reanalysis pixel by pixel and removed the ERA5 pixel that showed a difference over 50% to the GLDAS pixel.

\(\frac{\left|G L D A S_{(i, j)}-E R A 5_{(i, j)}\right|}{G L D A S_{(i, j)}}>0.5\)       (6)

2) Explanatory variables

(1) LST

LST (land surface temperature) is the radiative temperature on the ground obtained from thermal infrared data. LST is commonly used as an explanatory variable for the estimation of soil moisture because it is closely associated with surface energy balance and water content(Wan et al., 2004; Park, 2005;Jang et al., 2019). We used “Level 3 MODIS Land Surface Temperature and Emissivity Product” of the Aqua satellite (MYD11C3) obtained from LP DAAC (Land Processes Distributed Active Archive Center) (https:// lpdaac.usgs.gov/). The generalized split-window algorithm was employed in the retrieval of MODIS LST. T31 and T32 are the brightness temperature of band 31 (approx. 11 μm) and 32 (approx. 12 μm), respectively. The constant C, A₁, A₂, A₃, B₁, B₂, B₃ are obtained from an RTM (radiative transfer model) and the look-up tables created by linear equations(Wan and Dozier, 1996; Wan, 2002).

\(\begin{array}{l l} T_{S}=C+ \left(A_{1}+A_{2} \frac{1-\varepsilon}{\varepsilon}+A_{3} \frac{\Delta \varepsilon}{\varepsilon^{2}}\right) \frac{T_{31}+T_{32}}{2}+\left(B_{1}+B_{2} \frac{1-\varepsilon}{\varepsilon}+B_{3} \frac{\Delta \varepsilon}{\varepsilon^{2}}\right) \frac{T_{31}-T_{32}}{2} \end{array}\)       (7)

We extracted the monthly variable of daytime LST at the resolution of 0.05°, which is very useful because the gap-filling has been completed for the entire area. Particularly in summer, the MODIS products for 8-day LST include many missing pixels because of the influence of the cloud. In worst cases, more than half of the 8-day LST pixels in Korea may be missed in summer. An appropriate database for the explanatory variables of the ATPRK requires the minimized missing pixels. Hence, we alternatively employed the cubic spline interpolation (Aires et al., 2004; Chen et al., 2006) for converting the monthly LST into a 10-day composite. Because the 10-day LST shows a gradual change during a year, the cubic spline can be a reasonable method to interpolate the LST time series. In this way, we could obtain the 10-day LST almost without missing pixels.

(2) NDVI

NDVI (normalized difference vegetation index) is a vegetation index that represents the greenness of vegetation.It is based on the spectral characteristics that vital green vegetation has a higher reflectance at the NIR(near-infrared) band and a lower reflectance at the red band (Huete et al., 1999; Wan et al., 2004). So, the NDVI is the ratio of the difference and the sum of the NIR and red reflectance values. Soil moisture is closely associated with the NDVI, particularly in the growing season of vegetation (Farrar, 1994; Narasimhan et al., 2005). We extracted the 0.05° monthly NDVI variable from the “Vegetation Indices Monthly L3 Global 0.05 Degree Product” of the Aqua satellite (MYD13C2). In the same way as the 10-day LST, we carried out the cubic spline interpolation to convert from the monthly values to the 10-day composite

(3) PET

Part of soil moisture is evaporated into the air and is also absorbed by plant roots. Then the stomatal transpiration of the plant makes the water content emitted to the atmosphere. So, the soil moisture and evapotranspiration are closely related to each other within the hydrological cycle, which is also affected by the conditions of the meteorological factors such as temperature, humidity, wind, and precipitation (Wang and Liang, 2008; Hong et al., 2009; Jung et al., 2010; Ahn et al., 2013, Peng et al., 2017). We used the PET (potential evapotranspiration), a maximum potential of evapotranspiration under a given weather condition, which is widely used for the analyses of hydrology, meteorology, and drought (Lim et al., 2015). We extracted the 500-m PET variable from the “Net Evapotranspiration Gap-Filled 8-Day L4 Global 500 m Product” of the Aqua satellite (MYD16A2GF). The MODIS PET was retrieved by the integration of the PM (Penman-Monteith) equation (Monteith, 1965) and the revised RS-PM (Remote Sensing-based Penman-Monteith) algorithm (Mu et al., 2011). The MODIS data for land cover, surface albedo, and LAI (leaf area index) and the GMAO (Global Modeling and Assimilation Office) climate data were used in the PET retrieval.

(4) SPI3

Because the 10-day soil moisture can be associated with accumulated precipitation for more than ten days, we used the accumulated precipitation in the form of SPI(standardized precipitation index). The SPI represents the degree of dryness or wetness in terms of the accumulated precipitation compared with the usual amount(Mckee et al., 1993; Edwards, 1997). We calculated SPI3 from 3-month accumulated precipitation because the SPI3 is commonly used in the analysis of hydrological characteristics in Korea. We extracted daily precipitation from the TRMM (Tropical Rainfall Measuring Mission) 3B42 product developed by the joint program of NASA and JAXA. The data was downloaded from the website of GES DISC (Goddard Earth Sciences Data and Information Service Center)(https://disc.gsfc.nasa.gov/). We built an ECDF (empirical cumulative density function) from the 3- month accumulated precipitation and derived the SPI3 by converting the ECDF value to the z-score.

\(z_{i}=\frac{x_{i}-\mu}{\sigma}\)       (8)

3) In-situ measurements

As the reference forthe calibration and validation of the BMA ensemble, we gathered the TDR measurement provided by RDA(Rural Development Administration). It is the hourly soil moisture at a depth of 10 to 30 cm for the 63 stations (Table 1 and Fig. 2), and we aggregated the data into a 10-day composite.

Table 1. TDR measurement points operated by RDA (Rural Development Administration)

Fig. 2. Study area with the RDA (Rural Development Administration) stations for soil moisture measurement.

4) The procedure of ensemble downscaling

For the improvement of both accuracy and spatial resolution of soil moisture data, we conducted an experimentfor South Korea (33.75-38.75ºN and 125.5- 129.75ºE) during 2014-2018. We focused on April to September because soil moisture is essential for the vegetation growth in these months. First, the ESACCI daily, GLDAS 3-hourly, and ERA5 hourly data in 0.25° were aggregated into a 10-day composite to produce 4,832 matchups. We divided the 10-day matchups into six groups according to the months (April to September) to take account of the seasonality for the BMA training. Then we set up the weights for each member for each month to produce a BMA ensemble every ten days. Fig. 3 shows the procedure of creating an ensemble grid data for the 10-day soil moisture. If a pixel has three members, the BMA weighting scheme for the three members is applied for the weighted mean. If there are two members for a pixel, the BMA weighting scheme for the two members (e.g.,A:B orA:C or B:C) will be used. If only one memberis available, the pixel is filled with the member value. If no member value is available, a missing flag is recorded for the pixel. We used “EBMAforecast” library in R to create the BMA ensemble.

Fig. 3. Workflow for the map creation with a BMA ensemble.

Then we conducted the ATPRK downscaling from the BMA-ensembled soil moisture. Monthly LST and NDVI at a 0.05° resolution and the 8-day PET on a 500-m resolution were converted to the 10-day interval using the cubic spline interpolation. Daily precipitation at a 0.25° resolution was aggregated into a 10-day interval. The different spatial resolution of the four variables was rearranged on a 0.05° grid.

We upscaled the explanatory variables such as LST, NDVI, PET, and SPI3 from the fine resolution (0.05°) to the coarse resolution (0.25°) to build and MLR model at the coarse level. The residual at the coarse resolution, which was created by subtracting the MLR estimates from the coarse-resolution data, were downscaled by ATPK to produce the residual at the fine resolution. Also, the MLR estimates at the fine resolution were created by inputting the explanatory variables at the fine resolution to the MLR model. Finally, the estimates and residuals at the fine resolution were summed up for the residual correction. Fig. 4 shows the process of the ATPRK that integrates MLR and ATPK. We used “gstat” library in R to conduct the ATPRK downscaling.

Fig. 4. The workflow of the ensemble downscaling of soil moisture data.

4. Results and Discussions

1) BMA ensemble

Because the hydrological characteristics can be different by seasons in Korea, the BMA training was carried out by months. The weight for each member for each month was applied to the creation of the BMA ensembles for the early, middle, and late April to September. To compare the accuracies of the GLDAS, ERA5, ESA CCI, and the BMA ensemble against the TDR measurements, the validation statistics such as MBE (mean bias error), MAE (mean absolute error), RMSE (root mean square error), and CC (correlation coefficient) were presented in Table 2. Also, Table 3 showsthe validation statistics of August when the BMA ensemble notably out performed the other months.

Table 2. Validation statistics for the three members and the BMA ensemble (2014–2018)

Table 3. Validation statistics for the three members and the BMA ensemble (August 2014–2018)

In Table 2, the RMSE of the three members was approx. 7.3 to 7.5%, whereas that of the BMA ensemble was 5.806%. Each member had a certain amount of uncertainty, and any one of themembers was not superior to the others. So, the blending of the members based on the posterior probability brought about an overall decrease in errors. The mean RMSE of the three members was 7.412%, and the BMA ensemble showed an approx. 22% improvement (7.412-5.806)/7.412). The CC of GLDAS and ESA CCI wasrelatively low (0.341 and 0.323,respectively), and that of ERA5 was also not so good (0.500). However, the BMA ensemble had a CC of 0.659, which correspondsto a great improvement compared with the three members. In the case of August in Table 3, the mean RMSE of the three members was 8.410%, and the RMSE of the BMAwas 5.665%, which indicates a 33% improvement((8.410-5.665)/8.410). It means that the BMA system assigned the weights reasonably to each member according to months. Fig. 5 shows that the BMA ensemble maintained a very stable RMSE for the entire period of the experiment, unlike the three members. Fig. 6 illustrates the example maps for the created BMA ensemble. The distributions of values were natural as the result of the appropriate blending of the three different data sets.

Fig. 5. Monthly RMSE of soil moisture for the period between Apr 2014 and Sep 2018.

Fig. 6. Blended soil moisture with BMA: Case of 2015.

2) ATPRK downscaling for the BMA ensemble

The BMA ensemble at the resolution of 0.25º was downscaled to the resolution of 0.05º by the ATPRK using the four explanatory variables(LST, NDVI, PET, and SPI3). In the case of August, the MLR model had the explanatory power (R2) of 65.4% for the response variable (BMA ensemble). Fig. 7 shows the soilmoisture values for the BMA ensemble and the MLR estimates atthe coarse and fine resolutions, respectively, formiddle August 2015. In the areas with a significant absolute residual such as Gangwon-do, some inconsistency between the BMA soil moisture and the MLR estimates was found (Fig. 8). So, we created the residual grid at the fine resolution using ATPK to avoid the inconsistency and conduct the residual correction.

Fig. 7. Soil moisture values for BMA (left), MLR estimates at coarse resolution (middle), and MLR estimates at fine resolution (right). Case of the middle August 2015.

Fig. 8. Soil moisture values for BMA (left), MLR residuals at coarse resolution (middle), and ATPK residuals at fine resolution (right). Case of the middle August 2015.

The ATPRK downscaling (i.e.,the integration of MLR and ATPK) is accomplished when the fine-resolution residual(Fig. 8) is added to the fine-resolution estimate (Fig. 7). It could ensure data consistency by improving spatial resolution and preserving the original pattern at the same time. Fig. 9 and 10 show the examples of the ATPRK downscaling for the middle August 2015 and late August 2017. In August 2015, the amount of rainfall was only 110 mm, which ranked in the top five for recent 50 years in terms of the small rainfall in August. The soil moisture was relatively low forsummer (approx.24% to 34%)(Fig.9), and the MLR downscaling resulted in the underestimation in Gyeonggi-do and the over-estimation in Gangwon-do. Meanwhile, the late August 2017 was a usual summer with the almost average level of rainfall and temperature. The soil moisture was an ordinary level for summer (approx. 30% to 38%) (Fig. 10), and the MLR produced a wide underestimation for the south part of Korea. Unlike the MLR, however, the residual correction of the ATPRK brought about a consistent spatial pattern, and the coarse-resolution data were well preserved in the fine-resolution data without under- or over-estimation, irrespective of the amount of soil moisture. The data consistency was demonstrated quantitatively by the scatter plots between the 0.25° pixels for the BMA ensemble and the 0.25° pixels upscaled from the 0.05° pixels of the ATPRK result. The scatter plots agreed with the 1:1 line, and the CC was 0.998. Indeed, the MLR residual on the 0.25° grid was accurately transformed into the ATPK residual on the 0.05° grid, with the CC of 1 (Fig. 11). However, one coarse pixel did not always correspond to 25 fine pixels because of the missing grid at the 0.05° resolution. With no missing grid, the CC between the BMA ensemble and the ATPRK result should be 1 theoretically. In actual circumstances, however, a few missing values can exist in the fine resolution pixels, particularly near the coastline. The dots that were not perfectly on the 1:1 line corresponded to the cases that include lessthan 25 pixels. In our experiment, the missing value effect was negligible, and we made sure that the residual-corrected ATPRK could preserve the spatial patterns of the BMA ensemble, despite the various range of soil moisture.

Fig. 9. Consistency of data values before and after downscaling (middle August 2015)

Fig. 10. Consistency of data values before and after downscaling (late August 2017).

Fig. 11. Consistency between the MLR residual of a coarse pixel and the ATPK residuals of fine pixels.

5. Conclusions

We presented an ensemble downscaling method for quality improvement of gridded soil moisture data in terms of the accuracy and the spatial resolution. Unlike previous studies, we could achieve both goals, that is, the accuracy improvement by BMA and the resolution improvement by ATPRK. In the experiments, the BMA ensemble showed a 22% better accuracy than the data sets from ESACCI, ERA5, and GLDAS in terms of RMSE. Also, the ATPRK downscaling could enhance the spatial resolution from 0.25º to 0.05º while preserving the improved accuracy and the spatial pattern of the BMA ensemble, without under- or over-estimation. In future work for the ensemble downscaling, a longer period data will be required for the BMA ensemble, and the addition of explanatory variables such as land cover, solar radiation, surface albedo, and surface energy balance will be necessary for the ATPRK downscaling. This will contribute to a more stably blended data and a more robust MLR model to produce a better-quality soil moisture. Such data sets can be used in a variety of local and regional applicationsrelated to soil moisture, such as agriculture, forest, hydrology, and meteorology. Because the ensemble downscaling method can be applied to the other land surface variables such as temperature, humidity, precipitation, and evapotranspiration, it can be a viable option to complement the accuracy and the spatial resolution of satellite images and numerical models.