Dual state/rainfall correction via soil moisture assimilation for improved streamflow simulation: Evaluation of a large-scale implementation with SMAP satellite data

The Soil Moisture Analysis Rainfall Tool (SMART) is a rainfall correction scheme developed and updated by Crow et al. (2009; 2011) and Chen et al. (2012). It is based on sequential assimilation of soil moisture (SM) measurements into a simple Antecedent Precipitation Index (API) model to obtain SM increments, and then linearly relates these increments to rainfall accumulation errors. In the study we extended the ensemble Kalman filter (EnKF) version of SMART developed by Crow et al. (2011) to an ensemble Kalman smoother (EnKS) version with probabilistic rainfall estimates.


Tool (SMART)
The Soil Moisture Analysis Rainfall Tool (SMART) is a rainfall correction scheme developed and updated by 2011) and Chen et al. (2012). It is based on sequential assimilation of soil moisture (SM) measurements into a simple Antecedent Precipitation Index (API) model to obtain SM increments, and then linearly relates these increments to rainfall accumulation errors. In the study we extended the ensemble Kalman filter (EnKF) version of SMART developed by Crow et al. (2011) to an ensemble Kalman smoother (EnKS) version with probabilistic rainfall estimates.
Following 2011), the API model is used to capture the response of moisture storage (represented by the API state) to rainfall input: where t is a timestep index; P is the original uncorrected precipitation observation and γ is a loss coefficient (dimensionless) that accounts for storage loss through evaporation, drainage, etc. In the ensemble version of SMART (Crow et al., 2011), Eq. (S1) is converted to: where the superscript (j) denotes the jth ensemble member; η is multiplicative noise with mean 1 added to the observed precipitation to represent random precipitation forcing error; and ω is zero-mean Gaussian noise to represent random API model structure and parameterization error.
The API state can be related directly to SM content via rescaling where K is the Kalman gain; Tim is the covariance matrix between API states at time i and m; R is the measurement error variance for the rescaled SM measurements; the superscript (j) denotes the jth ensemble member; the superscripts "-" and "+" denote API states before and after an update, respectively; к is zero-mean Gaussian noise added to represent the random SM measurement error. Tim is calculated as: where M is the ensemble size; t API − is the ensemble-mean API states before update.
The SMART algorithm then uses ensemble-mean API increment δ to estimate the rainfall correction amount via a simple linear relation. We extended this relation to produce an ensemble of corrected rainfall time series (instead of the single rainfall estimates in past studies) where each ensemble member of the perturbed rainfall time series is corrected by the corresponding member of δ. : where "[ ]" denotes temporally aggregated P or δ (in the SMART study in this paper, this window was set to the 3-hour native SMART timestep without aggregation); l is the new time index for the aggregated windows; is a scaling factor that can either be calibrated or set to a prescribed constant. Finally, negative Pcorr resulted from Eq. (S6) are reset to zero, and the final corrected precipitation time series is (multiplicatively) rescaled to be unbiased over the entire simulation period toward the original precipitation observation time series.

S2. The impact of the λ parameter in the SMART rainfall correction scheme
In the SMART rainfall correction scheme, λ is a scaling factor that linearly relates the API state increment to rainfall correction amount. It can either be calibrated or set to a prescribed constant. We experimented with two strategies of determining λ in this study: 1) calibrating a temporally constant λ at each SMAP pixel separately to optimize the rainfall correlation with respect to the NLDAS-2 benchmark rainfall, and 2) setting λ to a spatial constant of 0.1, which is applicable for any region that may not have a good rain gauge coverage.
The rainfall correction results from the two strategies are shown in Fig. S1, in which Column 1 shows the improvement of correlation coefficient r after SMART correction with λ tuned at each pixel to maximize r (with respect to the NLDAS-2 benchmark), and Column 2 shows results obtained using a domain-constant value of λ = 0.1. Simply setting λ = 0.1 results in slightly smaller correlation improvement compared to the optimal λ case for all temporal accumulation periods (3-hour, 1-day and 3-day), especially for locations in the eastern and western ends of the domain. In general, these reductions are small, and since constant-λ is a more generally applicable case, we decided to use the λ = 0.1 strategies for all the results presented in the main manuscript. Figure S1. Maps of correlation coefficient improvement after SMART EnKS rainfall correction.
The left column shows the results with λ tuned at each pixel to optimize the correlation coefficient of corrected rainfall relative to the NLDAS-2 benchmark, and the right column shows the results with domain-constant λ = 0.1 [-] (this column is identical to the left column in Fig. 3 in the main manuscript). Each row shows results based on different temporal accumulation period: 3-hourly, 1-day and 3-day aggregation, respectively. The number on the lower left corner of each subplot shows the domain-median correlation improvement.

S3.1. Background and methods
It is well known that correlated errors in different parts of a Kalman filter result in suboptimal filter outputs. Therefore, in the original paper detailing the dual state/rainfall correction system, Crow and Ryu (2009) advised that the corrected rainfall (informed by the SM measurements) should not be fed back into the state EnKF correction scheme into which the 5 same SM measurements are assimilated. Instead, corrected rainfall and states should be combined via an offline model simulation (see Fig. 1

and Sect. 2.4.3 in the main manuscript).
Later studies that applied the dual correction system all followed this general guideline (e.g., Chen et al., 2014;Alvarez-Garreton et al., 2016). However, although this guideline helps avoid first-order error correlation in the system, it does not completely eliminate the possibility of error cross-correlation. Specifically, the corrected rainfall and the updated states are informed by the same SM measurement, thus they potentially inherit the same error from the SM measurement.
When fusing the two schemes together, such inherited error could potentially be amplified, degrading streamflow performance or cause a probabilistic estimate (based on an implicit assumption of independent errors) to be biased or have inaccurate uncertainty spread. In other words, it is possible that the current system still suffers from some second-order issue of overusing the information of SM measurements. Massari et al. (2018) intentionally avoided combining the state update scheme and the rainfall correction scheme in their study due to this legitimate concern.
To further investigate this issue, we designed a set of synthetic experiments and applied in an arbitrary small domain within the Arkansas-Red (a box around the Little Arkansas subbasin, see Table 1 and Fig. 2 in the main manuscript for its location). Synthetic measurements, instead of the real SMAP measurements, were generated and assimilated into the dual correction system so that we have complete control over all the error statistics and correlation, which is impossible in a real-data case. Specifically, a single perturbed VIC realization (with perturbed forcing and states) was treated as the synthetic "truth". Synthetic measurement can then be generated at daily interval by degrading the true surface-layer SM by adding random measurement errors. Precipitation perturbation was assumed to be temporally auto-correlated (first-order autoregressive noise with parameter ϕ = 0.9), and all the other error assumptions and dual correction setup were consistent with those described in Sect. 2.4 in the main manuscript.
We generated two sets of synthetic measurements based on the same truth with the same measurement error statistics but mutually independent realizations of errors. Then, two scenarios of dual correction were designed and carried out (see Fig. S2 for illustration): 6 Scenario 1: the same set of synthetic SM measurement were assimilated into both the state update and the rainfall correction schemes. This scenario mimics the issue in the real-data dual system with error cross-correlation in the two schemes and potentially degraded streamflow; Scenario 2: two sets of synthetic SM measurements (with mutually independent errors) were assimilated into the two schemes separately. This scenario completely avoids the issue of error cross-correlation.
The final runoff performance from the dual correction system were evaluated toward the truth, and the runoff performance from the two scenarios was compared. Differences in the performance of the two scenarios would indicate degradation caused by error cross-correlation.
For these synthetic experiments, runoff was evaluated locally at each grid cell without routing, since we know the true condition locally.

S3.2. Results
Deterministic and probabilistic results from the two scenarios were compared in Fig. S3 and Fig. S4. Clearly, runoff results show only very little difference between the two scenarios in terms of both PER and NENSK (see Sect. 2.5 in the main manuscript for details of the two metrics). This is true for both the total runoff and the fast-and slow-response runoff components separately. This suggests that the streamflow performance is not noticeably degraded by assimilating the same SM retrievals to both the state update and rainfall correction schemes.
Although the cross-correlated error theoretically exists in the system, they are not big enough to cause problematic streamflow results. In other words, we are not over-using the information contained in SM retrievals in the system. This is true both from a deterministic sense and in terms of probabilistic representation. We also experimented the case where the synthetic measurements were assumed to have temporally auto-correlated errors instead of white errors, which in theory creates bigger risk of degradation in the subsequent streamflow, but drew similar conclusions as above (results not shown).
The synthetic results in this section validates that we can safely assimilate the SMAP retrievals into both schemes of the dual correction system without significantly degrading the final streamflow performance.
7 Figure S2. Illustration of the synthetic experiments for investigating error cross-correlation.