The study sites lie in the Otway and Murray Basins, in south-east South Australia (Figures 1a,b). The region has a Mediterranean type climate, with annual rainfall of 640–750 mm year⁻¹, and potential ET of 960–1,400 mm year⁻¹. The region is generally flat, with undulations in the south-west forming Mount Gambier and Mount Burr (Doble et al., 2017). Up to 95% of regional groundwater resources are found in the shallow unconfined aquifer, with WT levels being within 2 m from the surface at their highest (Benyon et al., 2006). Spatial variations in the depth to WT are mainly driven by surface topographic variation between flats and interspersed sand ridges (see Figure 1c).

Softwood and hardwood forestry for timber production has been practiced in the area for more than a century. The dominant plantation timber species are Monterey Pine, Pinus radiata and Tasmanian Bluegum, Eucalyptus globulus (Benyon et al., 2006). Where these plantations overlie shallow WTs, which are within 6 m of the surface, it is likely that the trees directly access groundwater (Canadell et al., 1996; Benyon and Doody, 2004, 2015; Benyon et al., 2006, 2008).

Ecohydrological studies in the area have previously characterized vegetation, soil, climate, net irrigation, potential and actual evaporation, and depth to the WT, from silvicultural sites (Benyon et al., 2006). We selected four of these sites to answer the research questions. Two sites (“Woakwine” and “Jack”) had shallow WT depths (2.5 and 3m), Woakwine with a homogeneous soil profile and Jack a heterogeneous profile. One site (“Jill”) has a deep WT (10 m). The fourth site (“Piney”) records an average WT depth of 6 m (see Figure 1c for a map of the depth to WT). In general 6 m is the maximum WT depth at which the local forestry plantations can directly access groundwater (Benyon et al., 2006), thus this location is at the threshold between deep and shallow WT. Of these sites, Jill has a homogeneous soil profile, and Piney a heterogeneous profile. Climate is similar between the sites, and all contain mature trees: E. globulus except Piney, which contains P. radiata. Site details are presented in Table 1. Note that site names correspond to the those in Benyon et al. (2006), with the codenames in that paper being Woakwine: 4.EGl; Jack 5.EGl; Jill 6.EGl; and Piney 11.PR.

Remotely-sensed ET estimates were obtained from CMRSET (Guerschman et al., 2009). This data estimates actual ET based on surface reflectance from MODIS-Terra and interpolated climate data. It uses reflectance to rescale the Priestley-Taylor potential ET and calculates actual ET through the global vegetation moisture index (GVMI) and the enhanced vegetation index (EVI). Originally, the CMRSET algorithm was devised to produce monthly values of actual ET, using 16-day reflectance values, upscaled in time to a monthly value at 1 km resolution. In this version of the dataset, the CMRSET reached an averaged RMSE of 14.9 mm month⁻¹ when compared to five flux towers. A finer temporal and spatial resolution data-set, with a 250 m grid at a 8-day composite time step, is used in this study. Evapotranspiration values are rescaled to the model grid size (1 km²) as a result of averaging sixteen 250 × 250 m tiles.

The polar chart of Figure 2B shows, for Jack and Jill, the field-measured actual ET and the 1 km² averaged actual ET value from CMRSET. The spatial resolution of the resulting CMRSET dataset is well-suited to distinguish differences in actual ET between distant sites, but not the paired Jack and Jill sites as results are located within the same satellite observation grid. However, actual ET dynamics at the two sites differ. Jack has higher actual ET rates in summer, which are conceptually attributed to direct use of groundwater by vegetation, but not seen at Jill due to the deeper WT (Benyon et al., 2006). Jill, conversely, displays higher actual ET than Jack during May to August, presumably due to increased rainfall recorded at that site relative to Jack (Table 1). Given the differences in these two proximal sites, the simulations in this study can provide insight into how CMRSET resolution might impact the value of calibration to an assimilation of remotely sensed actual ET.

Simulation results were evaluated using several different metrics. These included: (i) RMSE (see Equation 2), (ii) the Pearson correlation coefficient (r) as in Pearson (1920), and (iii) the Continuous Ranked Probability Score (CRPS—Hersbach, 2000), a measure of the difference between the predicted and observed cumulative distributions. CRPS avoids some of the limitations of sum of squared error metrics (like RMSE or r) which can be highly sensitive to the largest errors, potentially biasing the model skill evaluation (Schneider et al., 2020).

The CRPS for the probability density function P(x) given by the ensemble simulation for the variable of interest x (i.e., ET and WT levels) and calculated at a specific timestep k is defined by:

where P_o is distribution of observations. As the observation at a given time is usually a single value, P_o(x) is equivalent to H(x−x_a), where H is the Heaviside function (Bracewell, 2000), applied similarly to Gelsinari et al. (2021).

With this definition, the CRPS is interpreted as the area between the cumulative probability function of the ensemble forecast and the observation function (linear). A CRPS of zero implies a perfect deterministic forecast (Hersbach, 2000). The CRPS varies from 0 to 1, and is usually calculated for all timesteps and then averaged over a simulation period.

RMSE values were calculated for all three calibration strategies used in these experiments. When the simulations were based on a parameter ensemble, RMSE values were calculated for the ensemble average predictions. CRPS values were calculated for the ensemble simulations only (i.e., the open loop and data assimilation simulations), following Equation (4). Verification skills were calculated over a validation period as previously described in Section 2.4.2.

Residuals of actual ET were computed as:

where i_k is the satellite observed actual ET and j_k is the modeled actual ET at time k, and L is the total number of observations. Values greater than zero indicate that the modeled actual ET under-estimates observations, while values higher than zero imply that the modeled actual ET is lower than the remote sensing actual ET observations.

Including actual ET in the models through calibration or assimilation improved simulations of actual ET in the sense that RMSEs declined and r increased, for every site other than Jill. The improvements were most evident at Jack, in which the RMSE decreased by ~30% and the r improved by ~100% upon inclusion of actual ET data in the modeling. Improvements at Woakwine and Piney are more modest, but still represent ~20% reductions in RMSE (see Table 3). The modeled actual ET at Jack is shown in the top panel of Figure 4. This panel shows that the Level-0 modeled actual ET (yellow line) greatly underestimates MODIS-CMRSET observations, an underestimation which is overcome by incorporating actual ET information, in the Level-1 and Level-2 models. The cause of the peak in ET of April 2003 is due to the model input, as this is similarly found at all the three levels, meaning that is not an artifact of the Ensemble Kalman filter update or the input perturbation strategy. The model response has been interpreted as a situation where the soil profile is particularly wet combined with high Potential ET (similar to summer conditions), causing the model to over-estimate actual ET.

As described in Section 2.4.3, the deep WT at Jill occurs within the same MODIS-CMRSET pixel as the shallow WT at Jack, but represents a very different water availability scenario for the plantation. With the model for Jack reproducing actual ET fluxes well when constrained with MODIS-CMRSET observations, it is perhaps unsurprising that the opposite occurs in the different environment at Jill. Level-1 and Level-2 models at Jill had an ~10% higher RMSE and an ~15% lower r relative to the Level-0 model.

CRPS values, aggregated over the validation period for the model runs, are reported in Table 4, with the ET results shown on the left hand side. CRPS values calculated at each time represent the actual ET level distribution of the ensemble around the observation, with a low value indicating good agreement and thus reduced uncertainty. The CRPS suggest that assimilation improves the ET ensemble distributions at Jack, Woakwine, and Piney, although the improvements are relatively small in percentage terms, and marginal at Jack. At Jill, however, the assimilation of actual ET increases the CRPS, actually worsening the model performance relative to calibration.

Figure 5 shows the residuals of 6-month aggregated values of actual ET for the single calibration, open-loop and assimilation runs, compared to CMRSET predictions. Values above zero indicate overestimation of modeled ET compared to remotely sensed observations (CMRSET), while values below zero are underestimates (Figure 5). For almost all runs across sites and years, the groundwater models calibrated only to WT levels underestimate actual ET. The right side of Figure 5 shows the average of these residuals for each 6 month period, showing that Level-0 models underestimate actual ET by between 74 and 140 mm year⁻¹. However, these underestimates were reduced at Level-2 by 25 % (at Jill) to 70 % (i.e., Jack and Piney). Examination of the calibrated parameters of the Level-0 models show that they predict shorter root lengths and a more linear root density distribution than do the Level-1 and Level-2 simulations.

The Level-0 simulations generated the best performance in terms of WT predictions as measured by RMSE in all the locations, as shown in Table 3, with the r being similar for most of the simulations. The CRPS metrics in Table 4, however, show that assimilation can improve ensemble WT predictions compared to the open loop, notably at Piney and Jack, in which quite large percentage improvements in CRPS are achieved by assimilation. The r in the Level-2 Jack simulation is reduced; this is a consequence of the EnKF continuously updating the WT level (see Figure 4B). For example, the Level-2 model shows abrupt changes in WT in the period of June–December 2003. These changes are due to the filter causing rapid and non-smooth changes in WT predictions, which lower performance metrics. Not all abrupt changes in the model predictions, however, are associated with the EnKF. For example, the shifts in WT around September 2002 occur in both Level-1 and Level-2 simulations, and are associated with the gap of the potential ET forcing data. As a results of this gap, the period from 2002/04/01 to 2002/10/31 was removed when calculating performance metrics.

Similarly to Jack, the CRPS for Piney improves between Level-1 and Level-2 and the performance of the model is good overall. However, the r-value is not improved. At Jill and Woakwine, CRPS, RMSE, and r worsened or were close to unchanged between Level-1 and Level-2. At Jill, the declining performance may indicate that actual ET data were disinformative of the local water balance (Beven and Westerberg, 2011).

Across both actual ET and WT, and for all levels modeled, the model predictions were most reliable at Piney. This site exhibited the highest stem density of all the plantations, which is likely to have produced near homogeneous conditions of water use, and may mean that satellite observations of actual ET were more reliable than when made over the more variable land surface conditions at the other sites.

The benefits of informing model calibration with actual ET were greatest at the Jack, Piney, and Woakwine sites. These sites had relatively shallow WT that were at the limit of or within the rooting depth of the plantations, suggesting that WT levels might be well coupled to evaporation at these locations. Actual ET assimilation has previously been shown to perform well in areas with significant groundwater-vegetation interactions (Gelsinari et al., 2020).

Information provided from actual ET, reduced model performance only at the Jill site. Presumably this is due to the localized deep WT at Jill, which is not representative of the conditions across the MODIS-CMRSET pixels used for the analysis. Under such a condition, the ET data is disinformative rather than informative, and worsens the model predictions. There is less differentiation between calibration and assimilation under these circumstances, because the ET in all cases does not add a meaningful constraint to the groundwater model. Preliminary analyses of digital elevation models could be useful to identify such locations that might exist in one remote sensing pixel.

Across all sites, incorporating actual ET information into model parameterization reduced the quality of predictions of the WT dynamics while improving the prediction of ET fluxes. This suggests that an improvement in the fidelity of modeling the local water balance might require reduced accuracy in modeling the specific signature of WT levels. It also indicates the risks inherent with using only a single metric to evaluate a model performance: overfitted models based on a single measurement can overestimate a model's ability to replicate all hydrological processes. If the models are only evaluated against WT level results, the single calibration (e.g., WT levels) seems to produce good results.

However, the quantitative management of groundwater resources requires having good estimates of fluxes, such as recharge and ET, and not only of groundwater levels. If the volume of water transpired by plantations is underestimated in models (as it was in all Level-0 models) then management based on model interpretations could be incorrect. Underestimation of water used by forest plantations could contribute to unsustainable management of groundwater resources, which itself can trigger non-linear and hard-to-reverse changes in water resource availability. By accounting for satellite based remote sensing ET, the models better replicated the fluxes of plant transpiration, especially in areas where plant-groundwater interactions were strong. The benefits of informing the model with actual ET greatly improved the estimates of ET volumes, at the cost of losing some precision in the estimates of WT levels. However, this loss of precision is not problematic with the benefit of obtaining better estimates of important hydrologic quantities such as ET or recharge.

Consequently, using actual ET observations from satellites is a valuable approach to improving the fidelity of groundwater models, especially if remote sensing outputs can be bias corrected using field collected data. This approach may be particularly beneficial in sites exhibiting vegetation-groundwater interactions, which closely couple WT levels to ET.

The difference between Level-1 and Level-2 is more subtle than the improvements relative to Level-0 models. There is a substantial increase in effort and skill needed to set up an EnKF scheme, and the benefit of doing so requires careful evaluation. In densely vegetated areas where the WT is within the reach of the plant roots, such as Jack, Woakwine and Piney, data assimilation improved RMSE and CRPS values for WT, actual ET and thus volumes of transpired water. Overall, however, the improvements in simulations between simple calibration and assimilation were relatively modest, and data assimilation may not be necessary for successful use of ET in groundwater models in general.

Although one of the key advantages of remote sensing data is to provide spatially distributed data, in this experiment that piece of information was not used as the focus was on the calibration and data assimilation of models using actual ET data. Regardless of the calibration method, the potential for rapid spatial changes in depth to groundwater and water balance within a satellite pixel is the main constraint on using ET for model calibration. Modelers should determine whether local groundwater dynamics can be reasonably represented at the scale of the satellite ET pixel before proceeding.

Future work in model evaluation could consider longer time periods, more sites, the use of drones for the estimation of actual ET at much higher spatial resolutions (Marzahn et al., 2020), and further refining metrics. For example, representing the CRPS as a time dependent metric allows for more in-depth analyses to identify the specific conditions for which an ensemble is or is not improved (Gelsinari et al., 2021). Similarly, it could further explore the role of the soil profile in moderating model performance. During this experiment, identifying and isolating the effects of integrating ET data over different soil profiles (i.e., Homogeneous vs. Heterogenous) was not conclusive. With the results obtained at Piney and Jack, it was shown that heterogenous profiles (usually considered more difficult to model and calibrate) do not represent an obstacle per se to the application of this methodology.

Here, we focused on measured actual ET fluxes and WT levels, and were not able to examine how ET information altered prediction of the distribution of soil moisture in the profile. A specific assessment of the soil moisture content at all the locations could provide more insights into the effect of using actual ET retrievals on model fidelity in the unsaturated zone. In this study, water content distributions accounted for by using ET fluxes as a proxy for the water content along the soil column.