Soil moisture (SM) and evapotranspiration (ET) are key variables of the terrestrial water cycle with a strong relationship. This study examines remotely sensed soil moisture and evapotranspiration data assimilation (DA) with the aim of improving drought monitoring. Although numerous efforts have gone into assimilating satellite soil moisture observations into land surface models to improve their predictive skills, little attention has been given to the combined use of soil moisture and evapotranspiration to better characterize hydrologic fluxes. In this study, we assimilate two remotely sensed datasets, namely, Soil Moisture Operational Product System (SMOPS) and MODIS evapotranspiration (MODIS16 ET), at 1-km spatial resolution, into the VIC land surface model by means of an evolutionary particle filter method. To achieve this, a fully parallelized framework based on model and domain decomposition using a parallel divide-and-conquer algorithm was implemented. The findings show improvement in soil moisture predictions by multivariate assimilation of both ET and SM as compared to univariate scenarios. In addition, monthly and weekly drought maps are produced using the updated root-zone soil moisture percentiles over the Apalachicola–Chattahoochee–Flint basin in the southeastern United States. The model-based estimates are then compared against the corresponding U.S. Drought Monitor (USDM) archive maps. The results are consistent with the USDM maps during the winter and spring season considering the drought extents; however, the drought severity was found to be slightly higher according to DA method. Comparing different assimilation scenarios showed that ET assimilation results in wetter conditions comparing to open-loop and univariate SM DA. The multivariate DA then combines the effects of the two variables and provides an in-between condition.
Climate extremes such as droughts and floods are becoming more severe and frequent, causing unprecedented threats to food and water security (Hameed et al. 2019; Alipour et al. 2020; Rammig et al. 2020). Drought is a natural climate extreme occurring in virtually all climatic zones. Drought is considered to be a complex phenomenon classified into four major types including meteorological, agricultural, hydrological, and socioeconomic drought. Among these four types, agricultural drought refers to a period of time with a deficit in soil moisture (SM), which could consequently result in crop failure. With the purpose of monitoring agricultural drought, several indices have been proposed based on a combination of temperature, precipitation, evapotranspiration (ET), and soil moisture (Ahmadalipour et al. 2017; Vicente-Serrano et al. 2010; Mu et al. 2013; Almamalachy et al. 2020).
Satellite data imagery for drought monitoring can be applied for either gathering atmospheric data such as precipitation and relative humidity, or land surface data acquisition such as SM and ET. The latter can be indirectly assimilated into the land surface models to achieve more accurate and reliable predictions of hydrologic fluxes as well as for monitoring purposes (Kumar et al. 2014; Pan and Wood 2006; Pipunic et al. 2008; Reichle et al. 2014; Sawada et al. 2015; Xu et al. 2020). SM prediction using land surface models driven by meteorological forcing carries considerable uncertainty due to spatiotemporal variability of forcing data. This uncertainty is magnified by large spatial variability of the land surface processes such as exchanges of energy, mass, and momentum. To improve the model predictions, an accurate characterization of the uncertainty in the model state, i.e., SM, is critical. This is particularly important for highly dynamic systems and those with sensitivity to initial conditions (Kumar et al. 2014; Mo et al. 2012).
As discussed earlier, land surface and hydrologic properties can be simulated and predicted by land surface models which provide a simplified representation of physical processes. However, an accurate prediction of these components, such as SM, ET, and streamflow, is highly dependent on the quality of model forcing data, the model parameters (measured or estimated through calibration), initial and boundary conditions, and model structure. For land surface and hydrologic models, the integration of data assimilation (DA) techniques with these models has been highly recommended because it improves the accuracy of water and energy balance computations and increases the model’s predictive skills (Reichle et al. 2014; Sawada et al. 2015; Seo et al. 2003). Earth system DA seeks to exploit real-time observations for more accurate hydrologic forecasts (Kumar et al. 2014; Reichle et al. 2014). DA aims at merging current and past observations with a dynamical model, using the model’s prognostic equations to estimate more accurate and reliable model state variables and parameters. Additionally, it provides a mathematical framework through which to optimally combine observations and model simulations (usually considered as prior information) based on their respective uncertainties (Moradkhani et al. 2018). For forecasting applications, it is used to characterize initial condition from available observations, so that more accurate forecasts can be generated (Yan et al. 2017; Davolio et al. 2017; Sahlaoui et al. 2020; Poletti et al. 2019; Abbaszadeh et al. 2020). These improvements are also dependent on the chosen variables to be assimilated into the system and their temporal and spatial relationships (Kumar et al. 2014).
In addition to SM, ET is also a crucial component of the terrestrial water cycle as a considerable amount of solar energy alongside 60% of total precipitation is consumed by ET (Trenberth et al. 2009). ET simulation is a complex task due to its dependence on many land–atmosphere interaction processes (Walker et al. 2019). The literature shows that ET variation is highly dependent on SM (Berg and Sheffield 2018; Jung et al. 2010; Purdy et al. 2018; Walker et al. 2019). Soil moisture controls latent and sensible heat exchange between land and atmosphere which causes feedback mechanisms in land–atmosphere interactions (Brutsaert and Stricker 1979). Several studies have shown that ET significantly contributes to the improvements of SM estimations (Berg and Sheffield 2018; Jung et al. 2010; Purdy et al. 2018; Walker et al. 2019) and plays a crucial role in predicting flash droughts (Chen et al. 2019). This strong SM–ET relation stimulates interest in the joint assimilation of observation of SM and ET into the land surface models for more accurate and reliable prediction of hydrologic fluxes and drought monitoring purposes.
The increasing availability of new types of satellite observations, in particular, SM and ET, provides the opportunity to explore their synergistic use (Su et al. 2014). However, the literature shows that the majority of efforts have gone into individual assimilation of SM or ET into land surface models and few evaluated the merits of dual assimilation of these variables (see, e.g., Hain et al. 2012). Several studies have reported the SM assimilation for improving drought monitoring and forecasting purposes (Bolten et al. 2010; Kumar et al., 2014; Yan et al. 2017; Yan et al. 2018). Kumar et al. (2014) showed that shorter-time-scale drought estimation can be improved by SM data assimilation. For ET assimilation, many efforts have gone into univariate assimilation of this prognostic variable to improve model simulations such as SM, latent heat flux, and ET itself (Pan et al. 2008; Pipunic et al. 2008; Qin et al. 2008; Schuurmans et al. 2003). It is noted that most of these studies provided results at relatively coarse spatial resolutions (12.5–25 km) that do not capture local to finescale spatial variations. Essentially, we need to downscale the coarse-scale microwave SM retrievals because more detailed information is required for improved drought monitoring and assessment at regional or local scales, especially for agricultural applications (Yoon et al. 2012). Considering that little attention has been given to multivariate DA in land surface models in particular at high resolution, this study investigates the benefits of joint assimilation of satellite ET and SM at 1 km and the framework is designed on a fully parallelized system to cope with the computational complexity. At this resolution, large-scale as well as local drought conditions can be monitored, hence a more realistic assessment of actual drought can be made.
Optimal assimilation methods are needed to maximize information content from observations and model simulations (Abbaszadeh et al. 2019a). However, most assimilation algorithms are suboptimal for complex real-world problems to which they are applied. However, the sequential Monte Carlo (SMC) methods based on particle filters (PF) do not rely on some restrictive assumptions such as the form of the probability distributions and Gaussian error assumption in model and observation, rather they propagate the full distribution of variable of interest over time in a nonlinear dynamical system (e.g., Moradkhani et al. 2018; Matgen et al. 2010; Sawada et al. 2015). As a successor version of SMC methods based on particle filter and Markov chain Monte Carlo (PF-MCMC), Moradkhani et al. (2012) and Abbaszadeh et al. (2018) embedded genetic algorithm in the PF-MCMC and developed the evolutionary PF that shows improvement in the state and parameter estimation in high dimensional systems. Therefore, in this study, we use this approach to jointly assimilate soil moisture and evapotranspiration data into a land surface model with the aim of improving drought monitoring.
The remainder of the paper is organized as follows. Sections 2 and 3 describe the study area and datasets. This includes North American Land Data Assimilation System (NLDAS), Soil Moisture Operational Product System (SMOPS), and MODIS evapotranspiration (MOD16A2) as well as datasets used for downscaling SMOPS to 1-km resolution. Sections 4 and 5 briefly explain the land surface model used in this study and the procedure for parallel joint DA. Section 6 assesses the performance of the proposed methodology and section 7 discusses the results of the DA performance and drought monitoring over the study area. Finally, section 8 provides a summary of the findings with the conclusion and potential for future expansion of this work.
2. Study area
The Apalachicola–Chattahoochee–Flint (ACF) basin, located in the southeastern United States, encompasses the three states of Alabama, Georgia, and Florida. The total area of the basin is 50 800 km2 most of which is located in western Georgia (Fig. 1). The Chattahoochee and Flint are the two major rivers in this basin originating from the north of Lake Sidney Lanier and south of Atlanta, respectively. These two rivers merge at Lake Seminole to form the Apalachicola River, which drains into the Apalachicola Bay and the Gulf of Mexico. The region has undergone extensive irrigation expansion between the 1970s and 1980s (Singh et al. 2017). Irrigated farmlands make up 452 000 acres within the region, making the area susceptible to huge economical losses during drought periods. The average annual precipitation is about 1400 mm which results in about 38–1020 mm annual runoff. Despite high average rainfall and runoff, the basin has experienced three extensive multiyear droughts between 2000 and 2015 (2000–01, 2007–08, and 2011–12). In 2011–12, the exceptional drought condition prevailed which led to significant ecosystem and economic losses (Leitman et al. 2016).
a. Forcing data
The meteorological forcing data are acquired from phase 2 of the North American Land Data Assimilation System (NLDAS-2; Xia et al. 2012b). More specifically, this study uses the NLDAS_FORA0125_H_002 dataset, which covers data from 1979 to present with spatial resolution of 1/8° and temporal resolution of 1 h. The land surface model used in this study is Variable Infiltration Capacity (VIC) model (Liang et al. 1994; Wood et al. 1992). VIC is forced with the air temperature at a height of 2 m above ground, total precipitation, incoming shortwave and longwave radiation at the surface, atmospheric pressure, vapor pressure and wind speed. The entire forcing data are available through the North American Land Data Assimilation System (NLDAS; Xia et al. 2012a).
b. Soil Moisture Operational Products System
The SMOPS developed by NOAA NESDIS (https://www.ospo.noaa.gov/Products/land/smops/index.html) is a combination product of multiple satellites and sensors providing global SM maps generated in 6-h and daily intervals with 0.25° × 0.25° spatial resolution (Zhan et al. 2011). SMOPS uses retrievals from six satellites including GPM, SMAP, GCOM-W1, SMOS, MetOp-A, and MetOp-B. The output product includes volumetric SM of the top surface layer (1–5 cm) alongside with quality information and metadata. Some examples of SMOPS applications in SM assimilation into land surface models include Nair and Indu (2016) and Yin et al. (2015, 2019).
The focus of agricultural drought monitoring has been mainly on broadscale conditions whereas more detailed information is required for improved drought monitoring and assessment at the small regional or local scale (Yoon et al. 2012). For this, SM products at finer resolutions (usually kilometer to subkilometer) are needed to better capture the spatial variability of the drought impacts, especially for agricultural applications. Original SMOPS images are available at 0.25° spatial resolution across the conterminous United States (CONUS). This resolution is too coarse to be applied over the ACF basin and to estimate the finescale drought extent over the area. Therefore, spatial downscaling can be applied to provide finescale regional drought information. Using a method proposed by Abbaszadeh et al. (2019b), the SMOPS product is downscaled to 1-km spatial resolution. In this method, high-resolution remotely sensed satellite imagery including MODIS NDVI (MOD13A2), MODIS LST (MOD11A2), Integrated Multisatellite Retrievals for GPM (IMERG) precipitation data, digital elevation maps (DEM), and ground-based soil moisture observation networks [i.e., U.S. Climate Reference Network (USCRN) and Soil Climate Analysis Network (SCAN)] are used to train a random forest (RF) machine learning algorithm (Abbaszadeh et al. 2019b). Out of 313 SCAN and USCRN stations, 66 stations were chosen to be included in the test set to evaluate the performance of the downscaling algorithm. Table 1 shows some statistics of the downscaled and the original SMOPS product over these stations. Figure 2 shows the spatial variability of the proposed product comparing it with the original SMOPS data.
c. MODIS evapotranspiration
The MODIS global evapotranspiration product MOD16 is a 500-m gridded land surface ET dataset for global land areas available at 8-day, monthly, and annual intervals (Mu et al. 2007, 2011). The output variables of the MOD16 product include 8-day, monthly, and annual ET, λE (latent heat flux), PET (potential ET), PλE (potential λE), and ET_QC (quality control). The data used in this study are the MOD16A2 ET product, which is produced at 500-m spatial resolution and 8-day compositing periods in the Sinusoidal projection (Running et al. 2017). Given the temporal and spatial scale mismatch between the MODIS and SMOPS data, different scenarios for DA can be considered. For instance, to improve the SM estimation, the assimilation of SMOPS can be done at a daily time scale, and then the MODIS ET is assimilated every 8 days. Another scenario could be the decomposition of 8-day ET data into daily maps by filtering the maps with the day of the year (DOY) map provided with the product. In the latter scenario, the resulting daily maps contain only a small number of pixels with data, resulting in limited coverage of the area, especially in the winter period. Thus, in this study, the first approach is chosen. To cope with the spatial resolution discrepancy between the observation (MODIS ET) and model prediction, the MODIS ET at 500 m was upscaled to 1-km spatial resolution using the mean aggregation technique. Figure 3 shows the annual average of the datasets used in this study.
4. VIC land surface model
To predict the terrestrial water, energy, and biogeochemical processes, land surface models (LSMs) are used where the governing equations of the soil–vegetation–snowpack medium are solved (Mo et al. 2012). An accurate representation of land surface processes is critical for improving the coupling of land and atmosphere at various spatial and temporal scales and over heterogeneous areas. The driving components of a typical land surface model are the initial conditions (states), boundary conditions including forcings and fluxes or states, and the soil, vegetation, and other surface parameters. In this study, we used Variable Infiltration Capacity version 5 (VIC-5) model (Liang et al. 1994) as the LSM. VIC is a semidistributed hydrologic model used in many applications including dataset construction (Nijssen et al. 2001), historic trend analysis, data evaluation, data assimilation, drought monitoring and forecasting (Nijssen et al. 2014; Shukla et al. 2011), and climate change impact analysis (Hamman et al. 2018). Unlike previous versions, VIC-5 supports parallel processing by utilizing the MPI standard and shared memory threading via OpenMP. This new feature enables us to implement the VIC model on the high-performance computing (HPC) system, which significantly improves the run time of the modeling. For more detailed information about the VIC algorithm, its main governing equations, and model state variables and parameters, we refer the interested readers to Gao et al. (2010).
a. Sequential Bayesian theory
where and ut are the vectors of the uncertain state variables and forcing data at time step t, respectively. The terms and represent the vectors of model parameters and observation data while ωt and υt denote model structural and measurement errors, respectively. These errors are generally assumed to be independent white noises with mean zero and covariance
where p(yt|xt) denotes the likelihood at time step t and p(xt|y1:t−1) is the prior distribution. The term p(yt|y1:t−1) is the normalization factor, and p(y1:t) is the marginal likelihood function, both of which can be determined using Eqs. (5) and (6):
Solving Eq. (3) analytically is only available for special cases such as linear systems with Gaussian assumption of noises in the system (i.e., the Kalman filter). Thus, for practical reasons this equation is generally approximated using a set of random replicates with associated weights:
where wi+ is the posterior weight of the ith particle. Parameters δ and N are the Dirac delta function and the number of particles, respectively. The normalized weights are determined by
In which, wi− is the prior weight of the ith particle. The term can be calculated using the likelihood as follows:
According to the posterior weights, the sequential importance resampling (SIR) is used to resample the state variables and parameters. Afterward, the proposal parameter distribution is generated using Eq. (10):
In this equation, and are the parameters before and after SIR implementation and st is a small tuning factor. To accept or reject the proposal parameter samples , a metropolis acceptance ratio α is calculated:
where is the proposed joint probability distribution:
where is a sample from the state proposal distribution and is the resampled forcing data at time step t. The tuning factor st is a time-variant unknown variable that can be estimated using the variable variance multiplier (VVM) method (Leisenring and Moradkhani 2011).
Evolutionary Particle Filter
The evolutionary PF-MCMC (hereafter EPFM) is a formal Bayesian DA technique that is implemented by a combination of sequential Monte Carlo technique genetic algorithm and MCMC that provides a full probability distribution of state variables and parameters and consequently their predictive uncertainty. The EPFM DA utilizes the MCMC technique twice in a sequential framework, once before resampling, in order to crossover and mutate the particles, and consequently produce a more informative prior distribution for state variables, and a second time after resampling to generate proposal parameter. While improving DA performance, EPFM significantly precludes the particle degeneracy and sample impoverishment problems that had been the main concerns in using the particle filter method. Here, we only summarize the four main steps of EPFM DA approach, and refer the readers to Abbaszadeh et al. (2018) for more detailed information.
Particles are selected for crossover operation from the original ensemble pool. To select these particles, we use the roulette wheel selection approach. Then we assign a fitness value to each ensemble member. The values of weights are representative of the quality of each particle; therefore, they can be directly used as the fitness value.
For crossover operation, we use the arithmetic crossover procedure to linearly combine the pair of selected particles. This process is shown by the following equations:
where and are the parent particles and and are the pair of newly generated offspring particles. Parameter ξ is a uniform random value ranging from 0 to 1.
To further intensify the diversity of the newly generated particles, we use GA mutation operator as follows:
where η is a random sample from the Gaussian distribution with zero mean and variance . The term is the variance of the prior states at the time t, and φ is a small tuning parameter.
Finally, we implement the MCMC approach to accept or reject the new ensemble members generated by GA crossover and mutation operators.
b. Performance measures
In this study, we use two deterministic performance measures including Kling–Gupta efficiency (KGE; Gupta et al. 2009) and unbiased root-mean-square error (ubRMSE; Entekhabi et al. 2010), as well as reliability (Renard et al. 2010), as a probabilistic measure, to assess the performance of the EPFM data assimilation approach. These measures are defined in Eqs. (17)–(19) as follows:
where yt and are observations and model predictions, respectively. denotes the covariance calculated between the observations and model predictions. Parameters σ and σ′ are the standard deviations of observations and model predictions, and μ and μ′ denote the average of observations and model predictions, respectively, and E[.] is the expectation operator. Parameter Zt is the score that represents the quantile of every observation time step, and Ut is the uniform distribution between 0 and 1. For more information, we refer the readers to Renard et al. (2010).
6. Research framework
The schematic view of the proposed framework is presented in Fig. 4. As is shown in this figure, first, the VIC model is driven with an ensemble of NLDAS-2 forcing data. To generate this ensemble, we use lognormal error distributions with a relative error of 25% for precipitation, and for temperature, the error is assumed to be homoscedastic and normally distributed with mean zero and standard deviation of 5°C. It is assumed that uncertainty in other forcing data including incoming shortwave and longwave radiation, specific humidity, surface pressure, near-surface wind in u and υ components are negligible. To characterize uncertainty in the initial condition, we assume a normal error distribution with a relative error of 4% (0.04 m3 m−3) for the entire soil column with three layers. Having the model ensemble for ET and SM, the downscaled SMOPS satellite data and MODIS ET are assimilated into the VIC model using the EPFM algorithm. It is also assumed that the MODIS ET observation and the downscaled SMOPS SM errors follow a normal distribution with a relative error of 10% and 4%, respectively. The choice of the 4% for observed SM is based on the ubRMSE of the downscaled product over the test stations. In real case experiments, the remaining errors arise from the model structural uncertainty, which herein is represented by adding white noise with a relative error of 25% to the SM at three layers and ET model outputs. The error value assumptions for forcings, ET observations, and model outputs are based on multiple trial-and-error experiments and our previous study (Abbaszadeh et al. 2018).
The main challenge in applying the particle filter data assimilation is the computational intensity. A great percentage of these computations are carried out inside the model, and fortunately, this provides a fertile ground for parallelization schemas such as model or domain decomposition (Yan et al. 2018). Although parallelization will significantly improve the framework performance, inefficient implementation of the assimilation could result in the run time being dominated by the assimilation part instead of the model. To overcome this issue, we propose a framework that incorporates both model and domain decomposition in two steps. First, each ensemble member is assigned to one core and the VIC model is run in parallel to provide an ensemble of model results (model decomposition). Then, having the state variables and model outputs of the entire ensemble, a particle tracking system tracks the index of each particle before and after the resampling (earlier discussed in section 5) step and stores the information of the prevailing particles in memory to be used in the next time step. At the same time, it will recursively partition the domain into smaller pieces and performs EPFM calculations (i.e., likelihood function evaluation, crossover, mutation, kernel density estimation, etc.) on each piece and then merge all the sections results using a divide-and-conquer parallel algorithm (domain decomposition) (Bentley 1980; Blanes et al. 2012; Prasad and Bruce 2011). The approach significantly improves model performance in terms of model run time. Figure 5 shows a schematic view of the parallelization framework used in this study.
7. Results and discussion
a. Performance of the data assimilation approaches
The spatial distribution of the performance measures over the study area is shown in Fig. 6. The top panel compares the KGE of predicted versus observed soil moisture for four different cases of open loop (OL), univariate assimilation of ET (ET-DA), univariate assimilation of SM (SM-DA), and multivariate assimilation of ET and soil moisture. In the OL run, the VIC model is run without assimilation whereas univariate and multivariate DA correspond to individual and joint assimilation of ET and SM, respectively. This figure shows that the OL run has an average KGE of 0.5 or less over the study area, and it has a low value of less than 0.10 in the shoreline. In the ET-DA scenario, MODIS ET observations were independently assimilated into the VIC model. Univariate SM-DA corresponds to the scenario in which only SMOPS downscaled data at 1 km are assimilated into the VIC model. As can be seen in Fig. 6, ET assimilation improved the top soil layer SM based on the downscaled SMOPS data as compared to the OL run, especially in the east and northeast of the region. In the SM-DA scenario the results were significantly improved that, on average, more than 78% of the study area has a KGE of 0.70 or higher. Further improvement is achieved by multivariate assimilation of ET and SM observations into the VIC model. The performance measures corresponding to the multivariate assimilation are provided under the multivariate column of Fig. 6. It is seen that the multivariate DA outperforms the OL, ET-DA, and SM-DA according to all three performance measures. Except for a small area near the outlet of the basin, the ubRMSE is less than 0.04 m3 m−3 over the entire study area. While this is true for both SM-DA and multivariate DA, the later performance shows some improvements. Comparing the reliability of the three cases shows the superiority of the multivariate DA as compared to OL and univariate DA assimilations as well.
To provide an independent measure of performance for the root-zone SM predictions we compared the multivariate DA results with the in situ measurements of two stations inside the ACF region. Two of the USCRN and one of the SCAN stations are located inside the ACF basin, although only one of the USCRN stations provided data for the analysis period of the current study. The comparison is shown in Fig. 7. As shown in this figure, SM predictions are consistent with the in situ measurements at both sites. The correlation coefficient for SCAN and USCRN stations are 0.629 and 0.631 and the ubRMSEs are 0.065 and 0.010, respectively, for both sites. Table 2 compares the performance metrics against all four scenarios in these two stations. The table shows that the multivariate DA outperforms all other scenarios. Comparing the ET-DA with the SM-DA, univariate assimilation of SM presumably results in better performance in SM estimation. Both univariate and multivariate scenarios significantly improve the OL results at the top soil layer and root-zone SM. This indicates that the developed framework is successful in representing the SM conditions at different soil layers which will be used later for drought monitoring over the ACF basin.
b. Drought monitoring over the ACF basin
A typical approach to characterize drought is through normalized indices that represent the deficit of hydrologic variables, such as precipitation, soil moisture, or streamflow. A variety of metrics have been developed for quantifying drought (Heim 2002). In this section, we evaluate the impact of multivariate DA using EPFM on categorizing drought through percentile-based drought indices that have been used in the NLDAS drought monitoring system (Sheffield et al. 2012). Using a method described by Kumar et al. (2014), we calculate the daily percentile values of SM using NLDAS data from 1979 to 2018.
The U.S. Drought Monitor (USDM) has been providing weekly estimates of drought conditions since 1999 (Svoboda et al. 2002). USDM depicts the drought intensity by classifying it into five categories: abnormally dry for percentile below 30% (D0), moderate drought, percentile ≤ 20% (D1), severe drought, percentile ≤ 10% (D2), extreme drought, percentile ≤ 5% (D3), and exceptional drought, percentile ≤ 2% (D4). In this section, we compare the drought conditions determined from the SM percentiles of multivariate DA against USDM archives (http://droughtmonitor.unl.edu/MapsAndData/MapArchive.aspx) during the period of 2018–19.
Root-zone SM is one of the leading indicators of agricultural drought, especially in warm seasons and climates (Bolten and Crow 2012). Figure 8 shows the monthly drought extent, based on estimated root-zone SM percentiles derived from multivariate DA results over the ACF basin. The first row shows that the drought is the most severe during the winter of 2018. After the winter, the drought extent shrinks and the severity decreases, from south to north, and then again increases during summer. After that, the conditions become completely normal during fall. One reason to explain this could be Hurricane Michael, which made landfall early and brought rainfall to Florida and Georgia, but not really to Alabama. Hence, the patterns in October show normal conditions for the eastern part of the basin but the western and northwestern parts show small areas with D0 drought categories.
Figure 9 shows a comparison between the USDM and the drought categories derived from root-zone SM percentiles after multivariate DA for four different weeks and scenarios. It is noted that the spatial resolution of multivariate DA is 1 km providing a more detailed depiction of the drought extension over the ACF region. This could help capture local drought conditions, whereas the USDM focuses on broadscale conditions only. The findings of this study are consistent with the USDM maps during the winter and fall season in terms of drought extension, however, DA results indicate a slightly more severe drought. In Fig. 9, the week of 6–13 February shows the most severe drought status of the year 2018. The drought extends almost over the entire domain, however, except for ET-DA, all the other scenarios extend the categories to D4 showing that the conditions are more severe based on these scenarios. The ET-DA results show normal conditions in the southern and northwestern portion of the region while other areas exhibit D0–D1 drought categories. The OL follows a somewhat similar pattern but with more severe conditions. This means that the ET assimilation results in wetter conditions comparing to the OL whereas the SM assimilation shows a dryer soil condition. This may be due to the fact that SMOPS shows lower values for SM comparing to the VIC predictions. Hence, the drought pattern leans toward D3–D4 categories in the SM-DA scenario. The multivariate DA combines the ET-DA and SM-DA conditions and is showing an in-between pattern, although the SM assimilation effect is more dominant. The second row of this figure shows the last week of March in which the conditions were normal in the southern and northern part of the region but D0–D1 categories in the central portion. However, the OL results indicate drought conditions all over the region. The drought extent significantly shrinks after ET assimilation and further expands in the SM-DA scenario. After the joint assimilation of both variables, the patterns are more consistent with the USDM map, although more severe drought is identified. This may be attributed to the coarser resolution information that the USDM is using to characterize drought. Also, USDM does not strictly map SM rather it considers streamflow and precipitation at different time scales, local experts, etc. (Svoboda et al. 2002). However, the findings of this paper rely only upon SM percentiles.
To further investigate the impact of torrential rainfall of Hurricane Michael on drought propagation, Fig. 10 compares the USDM maps and multivariate DA results for five weeks from 18 September to 16 October, four weeks prior, and during the Hurricane. It is shown in the top panel of this figure that a moderate drought is propagating throughout the region even after the occurrence of this hurricane on the 7 October and shrinks a week after with a one week delay. The multivariate DA, on the other hand, immediately detects the termination of drought due to the heavy rainfall of Hurricane Michael as it is shown in the week of 9 October the region is completely drought-free because of the early landfall of this hurricane.
8. Concluding remarks
This study examines the advantage of multivariate assimilation of SM and ET into the VIC land surface model to improve drought monitoring over the ACF region. The proposed framework efficiently performs this task by using a fully parallelized algorithm to account for the computational intensity of the multivariate DA in the high spatial resolution of 1 km. The study is conducted over the ACF region using the NLDAS-2 datasets as the forcing variables. SMOPS and MODIS ET are used as SM and evapotranspiration observations. The original SMOPS data are downscaled to 1-km resolution and then alongside ET, are assimilated into the VIC land surface model using the EPFM data assimilation scheme to improve root-zone SM predictions.
To examine the added value of multivariate assimilation, its performance was compared with that of the univariate assimilation of ET and SM. Both the probabilistic (i.e., reliability) and deterministic (i.e., KGE and ubRMSE) measures confirmed that the multivariate assimilation of ET and SM contributes more to improving drought prediction than any other univariate assimilation configuration. The posterior soil moisture estimates provided by the multivariate DA and NLDAS soil moisture data from 1979 to 2018 were used to calculate the daily soil moisture percentile over the ACF region for the year 2018. These values were transformed into five drought categories defined by the USDM that depicts the drought severity from abnormally dry (D0) to exceptional drought (D4). The generated drought maps were compared with those reported by the USDM. The findings of this study also revealed that the univariate assimilation of ET results in higher SM values, hence the wetter conditions compared to the OL, SM-DA, and multivariate DA. However, the SM-DA results in the driest condition most likely due to the soil moisture values of SMOPS, which are generally less than that of the model predictions. The multivariate DA, which simultaneously incorporates the ET and SM observations into the model, results in a more realistic soil moisture condition that is corroborated with in situ data and USDM maps. Since the USDM uses multiple indicators along with expert knowledge for drought monitoring, it is therefore not expected that our single-variable-based drought estimation perfectly matches with its drought maps.
This study applies a DA framework for drought monitoring over the ACF basin in the southeastern United States. Future studies will extend the geospatial scale of the study region and apply this framework over the southeastern United States and CONUS. Furthermore, in this study, it is assumed the model and observation errors are independent at each grid cell. To further enhance this multivariate assimilation experiment, future studies can include the cross correlations among the neighboring grids. Although this might improve the data assimilation effectiveness, it will dramatically increase the computational complexity making the data assimilation less favorable in large-scale studies. It is also worthy to investigate the impact of ET frequency in improving the performance of the data assimilation.
Partial financial support for this project was provided by the National Oceanic and Atmospheric Administration (NOAA) Modeling, Analysis, Predictions, and Projections (MAPP) (Grant NA18OAR4310319). Also, partial financial support was provided by the National Science Foundation, Grant 1856054.
Denotes content that is immediately available upon publication as open access.