An ensemble of 48 physically perturbed model estimates of the 1∕8° terrestrial water budget over the conterminous United States, 1980–2015

Terrestrial water budget (TWB) data over large domains are of high interest for various hydrological applications. Spatiotemporally continuous and physically consistent estimations of TWB rely on land surface models (LSMs). As an augmentation of the operational North American Land Data Assimilation System Phase 2 (NLDAS-2) four-LSM ensemble, this paper describes a dataset simulated from an ensemble of 48 physics configurations of the Noah LSM with multi-physics options (Noah-MP). The 48 Noah-MP physics configurations are selected to give a representative cross-section of commonly used LSMs for parameterizing runoff, atmospheric surface layer turbulence, soil moisture limitation on photosynthesis, and stomatal conductance.

The dataset spans from 1980 to 2015 over the conterminous United States (CONUS) at a monthly temporal resolution and a $\mathrm{1}/\mathrm{8}$^∘ spatial resolution. The dataset variables include total evapotranspiration and its constituents (canopy evaporation, soil evaporation, and transpiration), runoff (the surface and subsurface components), as well as terrestrial water storage (snow water equivalent, four-layer soil water content from the surface down to 2 m, and the groundwater storage anomaly). The dataset is available at https://doi.org/10.5281/zenodo.7109816 (Zheng et al., 2022). Evaluations carried out in this study and previous investigations show that the ensemble performs well in reproducing the observed terrestrial water storage, snow water equivalent, soil moisture, and runoff. Noah-MP complements the NLDAS models well, and adding Noah-MP consistently improves the NLDAS estimations of the above variables in most areas of CONUS. Besides, the perturbed-physics ensemble facilitates the identification of model deficiencies. The parameterizations of shallow snow, spatially varying groundwater dynamics, and near-surface atmospheric turbulence should be improved in future model versions.

1 Introduction

Estimates of the terrestrial water budgets (TWBs) – evapotranspiration, runoff, terrestrial water storage, and their constituents – over continental domains are of high interest for a broad range of hydrological applications. Publicly available data have been applied to investigate the state of the terrestrial water cycle (Trenberth and Fasullo, 2013a; Rodell et al., 2015; Scanlon et al., 2018; Yin and Roderick, 2020); to understand the interactions among hydrological processes, vegetation, climate, and human activities (Trenberth and Fasullo, 2013b; LaFontaine et al., 2015; Ward et al., 2014; Levia et al., 2020); to examine the availability and variability of water resources and use (Wu et al., 2021; Hejazi et al., 2014; Scanlon et al., 2012; Voss et al., 2013; Lv et al., 2019; Le et al., 2011; Rodell et al., 2009); and to assess the risk of extreme events such as droughts (Peters-Lidard et al., 2021; Prudhomme et al., 2014; Dai, 2013; Su et al., 2021) and floods (Emerton et al., 2017; Lin et al., 2018).

As the applications have expanded, the availability of TWB estimates has increased rapidly (Peters-Lidard et al., 2018; Saxe et al., 2021; Zhang et al., 2018). Commonly used estimation methods include remote sensing, in situ observations, and model simulations (Saxe et al., 2021; McCabe et al., 2017; Pan et al., 2012; Gao et al., 2010; Trenberth et al., 2007). Among these methods, land surface models (LSMs) are apt for continuously producing physically consistent TWBs over a large domain and long period, and their characteristics are particularly favorable for certain circumstances. For instance, LSMs can estimate various TWB components simultaneously; whereas, for some components, such as runoff (Lin et al., 2019; Beck et al., 2017), root-zone soil moisture (Xia et al., 2015b, a), and transpiration (Lian et al., 2018), direct remote sensing is either unavailable or highly uncertain. Additionally, LSMs are valuable in remote or topographically complex regions because of the sparseness of in situ observations (Kim et al., 2021). Estimations based on remote sensing and in situ observations are often impeded by scale mismatches and observation gaps, whereas these issues are rarely an impairment for LSM simulations. Besides, LSM simulations can complement remote sensing and in situ observations well. Combinations of estimates from different techniques can improve the estimation accuracy (Zhang et al., 2018; Pan et al., 2012; Zhao and Yang, 2018), while comparisons between model-simulated estimates and observations can reveal the impacts of human activities (Zaussinger et al., 2019) and underground processes (Zheng et al., 2020).

Several operational LSM simulation systems have been set up over different regions of the globe (Xia et al., 2019; Shi et al., 2011; Carrera et al., 2015; Rodell et al., 2004). The systems combine an ensemble of LSMs to utilize the competitive strengths of different LSMs and eliminate the weakness associated with individual ones. Among them, the North American Land Data Assimilation System (NLDAS) (Xia et al., 2012a, b; Mitchell et al., 2004) stands as a pioneering and successful one. The NLDAS Phase 2 (NLDAS-2) operates over the conterminous United States (CONUS) from 1979 to near real time at a spatial resolution of $\mathrm{1}/\mathrm{8}$^∘. The system generates a set of multi-source synthesized data of surface meteorology, vegetation, and soils and uses them to drive an ensemble of four different LSMs. The four LSMs – namely Noah version 2.8 (Ek et al., 2003; Chen and Dudhia, 2001a, b; Chen et al., 1997), Variable Infiltration Capacity (VIC) version 4.0.3 (Liang et al., 1994), Mosaic (Koster and Suarez, 1992), and Sacramento Soil Moisture Accounting (SAC) model (Burnash et al., 1973), were selected to give a good cross-section of the diverse range of LSMs (Mitchell et al., 2004). The models have varying strengths and weaknesses in process parameterizations and modeling skills (Kumar et al., 2017). A combination of multiple models can produce an aggregated estimate that outperforms most of the individual constituents (Fei et al., 2021; Beck et al., 2017; Guo et al., 2007; Ajami et al., 2007). An ensemble can also quantify the estimation uncertainty resulting from different model formulations (Troin et al., 2021; Cloke and Pappenberger, 2009). Evaluations of the NLDAS-2 four-LSM ensemble estimates have shown satisfactory performance in matching the observed evapotranspiration (ET) (Zhang et al., 2020; Xia et al., 2012b; Kumar et al., 2018), runoff (Xia et al., 2012a), and soil moisture (Xia et al., 2015b, a).

We have enriched the NLDAS-2 four-model ensemble with 48 perturbed-physics configurations of the Noah LSM with multi-physics options (Noah-MP) (Fei et al., 2021; Zheng et al., 2020, 2019). Noah-MP has more physically realistic representations of vertical stratification than the NLDAS-2 models have. A column of land in Noah-MP consists of a vegetation canopy layer, three snowpack layers, four soil layers, and a groundwater component (Niu et al., 2011; Yang et al., 2011). Conceptual (e.g., the five water tanks of SAC) and lumped (e.g., the combined vegetation–soil surface layer of Noah) representations of the stratification of vegetation and soil, as used in the NLDAS-2 models, are minimized. Moreover, Noah-MP has a more comprehensive representation of various land surface processes that are evident at different depths. The modeled processes include snow accumulation and ablation, infiltration, percolation, retention, freeze–thaw of snow or soil water, groundwater recharge/discharge, and energy constraints (Niu et al., 2011). These improvements in vertical stratification and process parameterizations are expected to better estimate TWBs. Indeed, previous comparisons between Noah-MP and the four NLDAS-2 LSMs have shown that Noah-MP is comparable or better when it comes to estimating soil moisture (Cai et al., 2014b), runoff (Cai et al., 2014a; Fei et al., 2021), and ET (Zhang et al., 2020).

Our enrichment also features a single-model perturbed-physics ensemble, which is different from the widely used multi-model ensemble approach. The Noah-MP ensemble is constructed by shuffling the available parameterization options of several selected processes. The ensemble size grows exponentially as a multiplication of the available parameterization options of different processes (Yang et al., 2011; Zhang et al., 2016; Gan et al., 2019). A large ensemble should give a broad cross-section of feasible model formulations to account for the model uncertainty in TWB estimation (Telteu et al., 2021; Mitchell et al., 2004) and is critical for a statistically reliable estimation of the probability of hydrological events such as floods and droughts (Troin et al., 2021). The single-model perturbed-physics ensemble also facilitates uncertainty attribution and reduction. The ensemble consists of pairs that are different in the parameterization of one process and the same for another. Variance analysis of the ensemble can quantify the contribution of the parameterization of a process and compare the relative importance of two processes (Zheng et al., 2019; Clark et al., 2011). Such quantification could inform further model development to reduce the model uncertainty. However, there are also pitfalls unique to the single-model perturbed-physics ensemble. Fei et al. (2021) found that the ensemble members generated by naive perturbation of the Noah-MP physics are not independent enough from each other. The low independence hinders the skill gained from the ensemble method. The finding suggests that advanced techniques of physics perturbations should be developed to maximize the ensemble skill and minimize the ensemble size. An open-accessible dataset should facilitate the research that leverages the advantages and addresses the issues of the perturbed-physics ensemble.

We have previously evaluated the runoff and compared it with NLDAS-2 (Fei et al., 2021). This paper describes the estimation of all the TWB variables, along with the description of the spread among the ensemble members, the difference with the NLDAS models, and the performance with reference to various observations are presented. The paper is organized as follows. Section 2 presents the information necessary for using the dataset, including the dataset variables, file organization, and the source data and models used for data generation. Section 3 describes the intercomparison methods, evaluation metrics, and reference datasets. Section 4 presents the results and discussion. Finally, after stating the data availability in Sect. 5, Sect. 6 draws conclusions.

2 Data description

The dataset contains gridded water budget variables over CONUS. Section 2.1 describes the dataset variables and their physical relationships. The 48 Noah-MP physics configurations used to create the dataset are detailed in Sect. 2.2. Section 2.3 briefly covers the atmospheric forcing, the static parameters of vegetation and soil, and the simulation settings.

2.1 Dataset variables

Table 1 lists the dataset variables. The variables are available at each $\mathrm{1}/\mathrm{8}$^∘ grid point in NLDAS-2, indicated by a land–water mask (X). The surface water budgets of each grid cell are represented as follows.

Table 1 The dataset variables.

Neglecting horizontal water exchange between adjacent grids, the water budget closure can be obtained among the precipitation (P; $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), ET (E), runoff (R), and terrestrial water storage change ($W^{\prime }=\frac{\mathrm{d}W}{\mathrm{d}t}$, $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) (Zheng et al., 2020):

$$\begin{array}{}\text{(1)}& P=E+R+W^{\prime },\end{array}$$

where precipitation (P) is from NLDAS-2 (described in Sect. 1) and used as the model input.

Noah-MP resolves the components of the water budget closure Eq. (1). ET (E) consists of canopy evaporation (E_can), ground evaporation (E_gnd), and transpiration (E_tran):

$$\begin{array}{}\text{(2)}& E=E_{\mathrm{can}}+E_{\mathrm{gnd}}+E_{\mathrm{tran}}.\end{array}$$

Runoff (R) has a surface (R_srf) and subsurface (R_sub) component:

$$\begin{array}{}\text{(3)}& R=R_{\mathrm{srf}}+R_{\mathrm{sub}}.\end{array}$$

Terrestrial water storage (TWS; W) is the sum of snow water equivalent (SWE; W_snow), groundwater storage in unconfined aquifers (W_gw), and soil water content in the four model layers (W_soil,i):

$$\begin{array}{}\text{(4)}& W=W_{\mathrm{snow}}+W_{\mathrm{gw}}+\sum _{i=\mathrm{1}}^{N_{\mathrm{soil}}}{W}_{\mathrm{soil},i},\end{array}$$

where N_soil=4 is the number of soil layers. Soil water storage (W_soil,i) is not included in the dataset but can be calculated from the volumetric water content (w_soil,i) as follows:

$$\begin{array}{}\text{(5)}& W_{\mathrm{soil},i}={\mathit{\rho }}_{\mathrm{wat}}{w}_{\mathrm{soil},i}\mathrm{\Delta }{z}_{\mathrm{soil},i}\text{for}i=\mathrm{1},\mathrm{\cdots },N_{\mathrm{soil}},\end{array}$$

where ${\mathit{\rho }}_{\mathrm{wat}}=\mathrm{1000}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ is the water density; and $\mathrm{\Delta }{z}_{\mathrm{soil},\mathrm{1}}=\mathrm{0.1}$ m, $\mathrm{\Delta }{z}_{\mathrm{soil},\mathrm{2}}=\mathrm{0.3}$ m, $\mathrm{\Delta }{z}_{\mathrm{soil},\mathrm{3}}=\mathrm{0.6}$ m, and $\mathrm{\Delta }{z}_{\mathrm{soil},\mathrm{4}}=\mathrm{1}$ m are the thicknesses of the four soil layers.

2.2 The 48 Noah-MP physics configurations

The Noah-MP LSM version 3.6 is used. We perturbed the parameterization of runoff, stomatal conductance, soil moisture stress factor, and near-surface atmospheric turbulence. The processes are selected, as they directly control runoff generation and evapotranspiration. Their importance has been shown in global simulations (Yang et al., 2011). The perturbation creates an ensemble of 48 members (48 = 4 runoff × 2 stomatal conductance × 3 soil moisture stress × 2 turbulence). Limited by computational resources, we did not perturb the parameterization of the cryosphere hydrological processes such as snow albedo (Chen et al., 2014; He et al., 2019) and rain–snow partitioning (Wang et al., 2019), which may limit the usage of the dataset in cryosphere hydrology. Appendix A details the formulation of and parameters in the selected parameterization. The parameters use the Noah-MP default values.

2.3 Domain, temporal span, atmospheric forcings, and static parameters

The simulation domain covers the CONUS (25–53^∘N, 125–67^∘W) at a spatial resolution of 0.125^∘.

The hourly NLDAS-2 atmospheric forcings were used to drive the 48 Noah-MP configurations. This study used seven forcing variables: downward solar radiation, downward longwave radiation, air temperature, surface pressure, specific humidity, wind speed, and precipitation rate. The static datasets, including topography (https://ldas.gsfc.nasa.gov/nldas/elevation, last access: 6 February 2016), predominant vegetation class (https://ldas.gsfc.nasa.gov/nldas/vegetation-class, last access: 6 February 2016), and soil texture type (https://ldas.gsfc.nasa.gov/nldas/soils, last access: 6 February 2016), are also the same as those in NLDAS-2. We used the default Noah-MP lookup tables to convert the input vegetation and soil types to parameter values.

The simulation spans 36 years from January 1980 to December 2015 at a time step of 15 min. The initial states on 1 January 1980 were obtained by cycling the year 1979 one hundred times.

3 Intercomparison and evaluation methods

The evaluations and intercomparisons in this paper are performed for 12 River Forecast Centers (RFCs): Northeast (NE), Mid-Atlantic (MA), Ohio (OH), Lower Mississippi (LM), Southeast (SE), North Central (NC), Northwest (NW), Arkansas (AB), Missouri (MB), West Gulf (WG), California–Nevada (CN), and Colorado (CB). Figure S1 in the Supplement displays the geographical delineation of the RFCs. More details on the RFCs, such as their multi-year average precipitation, potential evaporation, and topography, can be found in Fei et al. (2021, Fig. 1).

The intercomparison and evaluations were conducted at different timescales – the long-term climatological mean, annual cycle, and interannual anomaly. Section 3.1 details how the temporal variations and ensemble spread are derived for the timescales. We utilized the Taylor diagram and Taylor skill score (TSS) to measure the performance of Noah-MP against various reference datasets. The evaluation methods are shown in Sect. 3.2, and the reference datasets are described in Sect. 3.4. In addition to the intercomparisons and evaluations, Sect. 3.3 introduces the Sobol' sensitivity index for the process sensitivity analysis.

3.1 Temporal variability and ensemble spread

We calculated the total temporal variability (σ_total), the variability of the annual cycle (σ_ancy), and the interannual variability (σ_anom) for each ensemble member and the ensemble arithmetic mean following Dirmeyer et al. (2006):

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{total}}& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{12}Y}}\sum _{y,m}(x_{y,m}-x_{\mathrm{clim}}{)}^{\mathrm{2}}},\text{(7)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{ancy}}& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\sum _{\mathrm{m}}(x_{\mathrm{ancy},m}-x_{\mathrm{clim}}{)}^{\mathrm{2}}},\text{(8)}& {\displaystyle }{\mathit{\sigma }}_{\mathrm{anom}}& {\displaystyle }=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{12}Y}}\sum _{y,m}(x_{y,m}-x_{\mathrm{ancy},m}{)}^{\mathrm{2}}},\end{array}$$

where x_clim, x_ancy, and x_anom are the climatology, annual cycle, and interannual anomaly of the monthly time series x_y,m (month m of the year y; $m=\mathrm{1},\mathrm{\cdots },\mathrm{12}$; $y=\mathrm{1},\mathrm{\cdots },Y$).

The ensemble spread for the total time series (${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{total}}$), annual cycle (${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{ancy}}$), and interannual anomaly (${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{anom}}$) are derived also following Dirmeyer et al. (2006):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{total}}={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T\mathit{\sigma }\left(x_{y,m}\right),\text{(10)}& {\displaystyle }{\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{ancy}}={\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}\sum _{m=\mathrm{1}}^{\mathrm{12}}\mathit{\sigma }\left(x_{\mathrm{ancy}}\right),\text{(11)}& {\displaystyle }{\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{anom}}={\displaystyle \frac{\mathrm{1}}{\mathrm{12}Y}}\sum _{y=\mathrm{1}}^Y\sum _{m=\mathrm{1}}^{\mathrm{12}}\mathit{\sigma }\left(x_{\mathrm{anom}}\right),\end{array}$$

where $\mathit{\sigma }\left(x\right)=\sqrt{\left[{\sum }_{n=\mathrm{1}}^N(x-\stackrel{\mathrm{ˇ}}{x}{)}^{\mathrm{2}}\right]/N}$ denotes the standard deviation across the ensemble members at each time step. $\stackrel{\mathrm{ˇ}}{x}$ denotes the arithmetic ensemble mean.

We calculated the ratio ($r=\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}/\mathit{\sigma }$) of ensemble spread $\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}$ to temporal variability σ for the total monthly time series (with a subscript of “total”), annual cycle (“ancy”), and interannual anomaly (“anom”). The ratio enables the intercomparison of the ensemble spread among the selected timescales. A grade is assigned according to the ratio following Dirmeyer et al. (2006): a grade “A” for r<0.316; “B” for $\mathrm{0.316}\le r<\mathrm{1}$; “C” for $\mathrm{1}\le r<\mathrm{3.16}$; “D” for $\mathrm{3.16}\le r<\mathrm{10}$; and “E” for r>10. A lower (higher) r or a higher (lower) grade denotes a lower (higher) ensemble spread.

3.2 Taylor diagram and skill score

The Taylor diagram (Taylor, 2001) is a graphical representation of how closely a model simulation matches observations in terms of correlation coefficient (R), normalized unbiased root mean square error (nuRMSE), and normalized standard deviation (${\widehat{\mathit{\sigma }}}_{\mathrm{f}}$). The TSS is an index that measures the distance between a model simulation and the observations in the Taylor diagram. The TSS is defined as follows:

$$\begin{array}{}\text{(12)}& {\displaystyle }\mathrm{TSS}={\displaystyle \frac{\mathrm{4}(\mathrm{1}+R)}{({\widehat{\mathit{\sigma }}}_{\mathrm{f}}+\frac{\mathrm{1}}{{\widehat{\mathit{\sigma }}}_{\mathrm{f}}}{)}^{\mathrm{2}}(\mathrm{1}+R_{\mathrm{0}})}},\text{(13)}& {\displaystyle }{\widehat{\mathit{\sigma }}}_{\mathrm{f}}={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{f}}}{{\mathit{\sigma }}_{\mathrm{o}}}},\end{array}$$

where σ_f and σ_o are the standard deviations of the model simulation and the observation, and R₀ is the maximum correlation coefficient attainable (in this study, R₀=1). The value range of TSS is from 0 to 1. A higher TSS indicates a higher overall performance of model prediction with reference to the observations.

3.3 Sobol's total sensitivity index

The sensitivity of the Noah-MP ensemble to a physical process can be quantified by the Sobol' total sensitivity index (Sobol', 1993; Zheng et al., 2019). The Sobol' total sensitivity index measures the proportion of the variance of different processes to the total variance, which is defined as follows:

$$\begin{array}{}\text{(14)}& S_{\mathrm{\Lambda }}={\displaystyle \frac{E_{\sim \mathrm{\Lambda }}\left({\mathrm{Var}}_{\mathrm{\Lambda }}\right(Y|\sim \mathrm{\Lambda }))}{\mathrm{Var}\left(Y\right)}},\end{array}$$

where S_Λ is the Sobol' total sensitivity index for one process Λ; ∼Λ represents the other processes except for Λ; Y represents the 48 Noah-MP ensemble members; Var(Y) is the variance of Y; ${\mathrm{Var}}_{\mathrm{\Lambda }}\left(Y\right|\sim \mathrm{\Lambda })$ denotes the variance among different parameterization schemes of the process Λ, and E(∼Λ) denotes the arithmetic average across all combinations of the other processes except for Λ. Detailed calculations and examples can be found in Zheng et al. (2019, Appendix A).

3.4 Reference data

3.4.1 Terrestrial water storage

We used the $\mathrm{1}{}^{\circ }\times \mathrm{1}{}^{\circ }$ monthly Gravity Recovery and Climate Experiment (GRACE) land water-equivalent thickness surface-mass anomaly, level-3, Release 6.0, version 04, as the reference of TWS (W in Table 1). The GRACE products from different processing centers are slightly different. To reduce the noise of the estimates (Sakumura et al., 2014), we used the arithmetic average of the products from three centers: (1) GeoForschungsZentrum Potsdam (or the German Research Center for Geosciences) (https://podaac.jpl.nasa.gov/dataset/TELLUS_GRAC_L3_GFZ_RL06_LND_v04, last access: 10 February 2022), (2) the Center for Space Research at the University of Texas, Austin (https://podaac.jpl.nasa.gov/dataset/TELLUS_GRAC_L3_CSR_RL06_LND_v04, last access: 10 February 2022), and (3) NASA's Jet Propulsion Laboratory (https://podaac.jpl.nasa.gov/dataset/TELLUS_GRAC_L3_JPL_RL06_LND_v04, last access: 10 February 2022). The GRACE satellites began orbiting Earth on 17 March 2002. We selected the period from 2003 to 2015 for the evaluation. There are 14 missing values during the period, which were filled with a simple linear interpolation.

The GRACE products experience signal leakages between land and water (Save et al., 2016). Such signal leakage could impact the estimation in the RFCs that are adjacent to the Great Lakes (Ma et al., 2017) and oceans. To alleviate the impacts of the Great Lakes, we corrected the GRACE TWS estimates in the NC, OH, and NE RFCs using the water level variations observed by the National Oceanic and Atmospheric Administration (NOAA). Appendix B details the algorithm.

3.4.2 Soil moisture

We used the daily North American Soil Moisture Database (NASMD; http://nationalsoilmoisture.com, last access: 6 February 2016) (Quiring et al., 2016) as the reference for the simulated soil moisture (W_soil,i in Table 1), similar to previous NLDAS evaluations (Xia et al., 2015b, a). NASMD assembles the soil moisture time series at multiple depths of more than 2200 stations of 24 networks with quality control. The observation depth varies with the network. We interpolated the observations to the centers of the Noah-MP soil layers, which are 0.05, 0.25, 0.7, and 1.5 m, respectively. The interpolation is performed only when a valid observation is exactly at, or two valid observations exist above and below, the given depth; otherwise, a missing value is given. We excluded the soil layers with more than 50 % missing values to minimize the impacts of missing values on the evaluation, after which 264 w_soil,1, 214 w_soil,2, 95 w_soil,3, and 23 w_soil,4 valid time series remained. Daily data from 1996 to 2013 were then aggregated into monthly values. For any month, if less than 10 d of valid data are available, a missing value is assigned.

3.4.3 Snow water equivalent

We used the daily Snow Data Assimilation System (SNODAS; https://nsidc.org/data/G02158/versions/1, last access: 13 April 2021) as the reference of SWE (W_snow in Table 1). SNODAS is a data assimilation system developed by the NOAA National Weather Service's National Operational Hydrologic Remote Sensing Center. This system aims to provide a physically consistent framework to combine snow modeling and observations from satellites, airborne platforms, and ground stations (National Operational Hydrologic Remote Sensing Center, 2004). The original spatial resolution is 1 km × 1 km, and we aggregated the data to the 0.125^∘ NLDAS grids. SNODAS began on 28 September 2003, and we selected the period from September 2004 to August 2015 that spans 11 whole snow seasons. Clow et al. (2012) showed that the SNODAS SWE performs well in the forest areas of the Colorado Rocky Mountains but performs poorly in the alpine areas.

3.4.4 Evapotranspiration

We used plot-scale AmeriFlux observations and two gridded products as the reference ET (E in Table 1). The two gridded products are derived from different methods, and a common evaluation period from 1982 to 2015 is selected for this study. The gridded products have different spatial resolutions. We downscaled the data to the NLDAS grids and then aggregated them to the 12 RFCs.

We selected 25 AmeriFlux sites (https://ameriflux.lbl.gov, last access: 18 September 2021). The 25 selected sites were selected because they have the longest observation periods for seven major land cover types (i.e., evergreen forest, deciduous forest, mixed forest, shrubland, savanna, grassland, and cropland). Figure S1 and Table S1 detail the selected sites. The data have been widely used in LSM evaluations (Cai et al., 2014a; Zhang et al., 2020). AmeriFlux provides hourly or half-hourly latent heat measurements. We converted the measurement to ET by dividing the latent heat of water vaporization (2.5104×10⁶ J kg⁻¹). The hourly or half-hourly ET is then aggregated into monthly values. In the process of aggregation, if there were less than 8 valid hours in a day, a missing day was marked; if there were fewer than 10 valid days in a month, a missing month was assigned; if there was less than 50 % valid months, the time series was dropped. In this study, the data serve as the ground truth for the gridded ET products and model estimation.

Table 2 The climatological mean, ensemble spread, and temporal variability of different dataset variables. $\overline{\stackrel{\mathrm{ˇ}}{x}}$ denotes the climatological average of the ensemble mean. ${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{ancy}}$, ${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{anom}}$, and ${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{total}}$ denote the spread of the 48 Noah-MP configurations in simulating the multi-year averaged annual cycle, interannual anomaly, and total 36-year monthly time series, respectively. σ_ancy, σ_anom, and σ_total denote the temporal variability of the annual cycle, interannual anomaly, and the total time series, respectively. The rating of the ratio $\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}/\mathit{\sigma }$ is defined in Sect. 3.1. The units of the variables are presented in Table 1. $\stackrel{\mathrm{̃}}{W}$ (${\stackrel{\mathrm{̃}}{W}}_{\mathrm{gw}}$) denotes the terrestrial water storage (groundwater) anomaly (kg m⁻²), whereas W^′ ($W_{\mathrm{gw}}^{\prime }$) denotes the monthly terrestrial water storage (groundwater) change ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$).

The first gridded ET product is FLUXNET multi-tree ensemble (MTE) (Jung et al., 2009) (https://www.bgc-jena.mpg.de/geodb/projects/Home.php, last access: 3 May 2017). FLUXNET MTE ET is a monthly dataset produced from the FLUXNET eddy covariance measurements, remote sensing, and meteorological data using the multi-tree ensemble statistical method (Jung et al., 2009). The product is widely used in LSM evaluations (Cai et al., 2014a; Ma et al., 2017; Xia et al., 2016; Jung et al., 2019; Fang et al., 2020; Zhang et al., 2020; Pan et al., 2020). FLUXNET MTE ET has a spatial resolution of 0.5^∘ × 0.5^∘. We remap the data to the 0.125^∘ × 0.125^∘ NLDAS grids with a first-order conservative method. We use the FLUXNET ET as the reference for the annual cycle, since it replicates the AmeriFlux observations best among the two gridded products (Tables S2–S4 in the Supplement).

The second gridded ET product is the Global Land Evaporation Amsterdam Model (GLEAM), version 3.3a (https://www.gleam.eu, last access: 7 May 2019), which is another widely used ET product (Xu et al., 2019). GLEAM estimates transpiration, canopy evaporation, soil evaporation, open-water evaporation, and sublimation separately and then sums them as ET. The method aims to maximize the utilization of satellite information. The product estimates monthly ET at a spatial resolution of 0.25^∘ × 0.25^∘. We bilinearly interpolated the data to the NLDAS grids. The GLEAM ET is used as the reference for the interannual anomaly, as it better replicates the AmeriFlux-observed anomaly than FLUXNET (Tables S2–S4).

3.5 NLDAS ensemble

We used three NLDAS-2 models – namely Noah-2.8, VIC-4.0.3, and Mosaic – as the benchmark of the Noah-MP ensemble. Their outputs can be publicly downloaded from the NASA Goddard Earth Sciences Data and Information Services Center (https://disc.gsfc.nasa.gov/datasets?keywords=NLDAS, last access: 6 February 2016). More information on the NLDAS-2 models, the forcing and static datasets, and simulation settings can be found in Xia et al. (2012b, a). The NLDAS-2 datasets have been proven to perform soundly for regional hydrological simulations (Xia et al., 2012b, a, 2016, 2015b, a) and are widely used as a benchmark for LSM evaluations (Cai et al., 2014b; Fei et al., 2021).

4 Results and discussions

In this section, we begin by intercomparing the ensemble spread of the dataset variables (Sect. 4.1). Then, Sect. 4.2–4.5 examines the performance and ensemble spread of the TWS anomaly (TWSA), SWE, soil moisture, and ET, in comparison with the NLDAS ensemble. The performance of runoff can be found in our previous papers Fei et al. (2021) and Zheng et al. (2020).

4.1 Comparison of the ensemble spread

We aggregated the dataset variables across CONUS. Table 2 summarizes the ensemble spread and temporal variability of the ensemble mean for all the dataset variables. Among the variables, runoff (including surface and subsurface components) shows the largest spread, accounting for one-fifth to one-third of the climatological mean. The ensemble spread is comparable to or larger than the temporal variability. The magnitude is similar to that estimated in GSWP-2. The spread of ET and TWS is significantly smaller than those observed in GSWP-2. The high consistency among the Noah-MP configurations could be a sign of the limited sampling of available process parameterizations but could also be a result of continuous model improvements. ET might be the former case, since the ensemble does not perturb several important processes such as roughness sublayer (Abolafia-Rosenzweig et al., 2021) and plant hydraulics (Li et al., 2021). TWS is likely the latter case, since the parameters of different groundwater schemes are calibrated using GRACE (Niu et al., 2007). The spread of snow water equivalent depicts the smallest spread, which is 2.5 % of the climatological mean. The small spread reflects a limited sampling of the cryosphere processes as discussed in Sect. 2.2. For soil moisture, the ensemble spread is marginal at the surface and increases significantly in the deep layers. The difference hints that the controlling processes vary with depth. The surface soil moisture is tightly controlled by the atmospheric forcings, whereas the spread of the subsurface soil moisture hints at the complex interplay among various land surface processes (e.g., root water uptake and subsurface runoff) (Koster, 2015).

Table 2 also shows the comparison between the annual cycle and interannual anomaly for different variables. For runoff and soil moisture, the ensemble spread in the annual cycle is larger (i.e., ${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{ancy}}\ge {\mathit{\sigma }}_{\mathrm{anom}}$). Whereas for transpiration (E_tran), groundwater storage (${\widehat{W}}_{\mathrm{gw}}$ and $W_{\mathrm{gw}}^{\prime }$), and snow water equivalent (W_snow), the spread is larger for the interannual anomaly (i.e., ${\stackrel{\mathrm{ˇ}}{\mathit{\sigma }}}_{\mathrm{ancy}}\le {\mathit{\sigma }}_{\mathrm{anom}}$). The dynamics at different times are modulated by different processes (Dickinson et al., 2003). Since the timescale of the interannual anomaly is smaller than that of the annual cycle, the land surface memory and the associated land model parameterizations may contribute a larger part to the annual cycle, whereas the interannual anomaly is more tightly controlled by atmospheric forcing. The smaller ensemble spread in the annual cycle than in the interannual anomaly is a sign of insufficient representation of the physics uncertainty.

4.2 Terrestrial water storage anomaly

Figure 1 shows the annual cycle of the TWSA estimated from GRACE, Noah-MP, and NLDAS. Figure 2 presents the TSS. In the 12 RFCs over CONUS, the TWSA peaks in spring, declines rapidly in summer, reaches a minimum in autumn, and recovers in winter. In terms of the timing of the peak and trough, Noah-MP and the NLDAS models perform similarly. In terms of the amplitude of variation, Noah-MP generally produces higher values in all RFCs. Previous studies have attributed this difference to the inclusion of a bucket groundwater component in Noah-MP (Cai et al., 2014b; Ma et al., 2017). However, we found the Noah-MP configurations without a groundwater component can still produce a higher amplitude, especially considering the structural similarity between Noah-MP and Noah. Further investigation of the model difference is necessary.

Figure 1 Annual cycle of TWSA from the model estimates and GRACE in the 12 RFCs for the period 2003–2015. The TWSA is calculated from the monthly TWS by subtracting the 13-year average (2003 to 2015). Black dots denote the GRACE observations. Error bars show the standard deviation of the year-to-year differences. The shaded areas denote the range between the maxima and minima of the 48 Noah-MP estimates. The solid red line denotes the Noah-MP multi-physics ensemble mean. The three NLDAS models (Noah, Mosaic, VIC) and their ensemble means are denoted by the blue, green, cyan, and dark golden lines, respectively. The 12 RFCs are sorted based on climatic aridity, i.e., the most humid RFC in the top left and the driest RFC in the bottom right.

Figure 2 shows the Taylor diagram for the annual cycle of TWSA. The Noah-MP configurations generally outperform the NLDAS models in most of the RFCs, which results in the superior performance of the ensemble mean (shown in Table S5). Detailed examination of the TSS reveals that Noah-MP and NLDAS have similar correlation coefficients. Their difference is manifested in the modeled standard deviation (i.e., the amplitude of variation). In NE, MA, NW, and CN, Noah-MP underperforms compared with NLDAS, mainly due to overestimating the standard deviation. However, the interpretation of the overestimation is multifaceted. First, Noah-MP could overestimate the variability due to unsuitable parameters. For instance, specific yield is an important parameter for the simulation of groundwater storage and water level (Lv et al., 2021). The parameter is calibrated as 0.2 by global simulations (Niu et al., 2007). The globally calibrated value may be overestimated in the RFCs, leading to an overestimation in TWSA. Second, the GRACE data could underestimate the temporal variability at these coastal RFCs due to signal leakage from the ocean (Cai et al., 2014b). In AB and MB, Noah-MP performs slightly better than the NLDAS models, but both underestimate the standard deviation. The Ogallala Aquifer encompasses the two RFCs. Noah-MP can better present the groundwater changes than the NLDAS models due to an unconfined aquifer module. But Noah-MP does not include the confined aquifers, leading to an underestimation of the groundwater storage variability.

Figure 2 Taylor diagram of TWSA's annual cycle for the 48 Noah-MP configurations (red dots) and three NLDAS models (blue dots for Noah, green stars for VIC, and cyan square for Mosaic) at the 12 RFCs. Black stars denote the observations. The ensemble mean of Noah-MP and NLDAS is presented by purple and dark golden diamonds, respectively. The distance between the model and observation presents the nuRMSE. The radial lines show the correlation coefficient, while the distance to the origin along the line denotes the normalized variability.

Figure 3 As in Fig. 2, but for the interannual anomaly of TWSA.

Figure 4 Monthly surface (0–0.1 m) and root-zone (0–1 m) soil moisture (left column) and the annual cycle (right column) from the Noah-MP ensemble, the NLDAS models, and the NASMD observations for the period 1996–2013. Only the RFCs with more than 10 observational sites are considered. Black dots denote the arithmetic average of the valid NASMD observations. Error bars in the right columns denote the standard deviation of the year-to-year differences. The shaded areas denote the range between the maxima and minima of the 48 Noah-MP estimates. The solid red line denotes the Noah-MP multi-physics ensemble mean. The three NLDAS models (Noah, Mosaic, VIC) and their ensemble mean are denoted by the blue, green, cyan, and dark golden lines, respectively. The five RFCs are sorted based on climatic aridity.

Figure 5 The first column shows the TSS of the NLDAS ensemble mean in simulating the annual cycle of soil moisture at different depths (0–0.1, 0.1–0.4, 0.4–1, 1–2 m). The second column presents the difference in TSS between the Noah-MP ensemble mean and the NLDAS ensemble mean. The third column depicts the TSS difference between the arithmetic average of the Noah-MP ensemble mean and the NLDAS models to the NLDAS three-model mean. P is the percentage of sites at which a higher TSS appears relative to the NLDAS ensemble mean. P1 (P2) is the percentage of sites at which a higher TSS appear given that the Noah-MP ensemble mean outperforms (underperforms) the NLDAS ensemble mean. The evaluation period is 1996–2013.

Figure 6 As in Fig. 5, but for interannual anomaly.

Figure 3 shows the Taylor diagram for the interannual anomaly of TWSA. Compared with the annual cycle (Fig. 2), both the Noah-MP configurations and the NLDAS models degrade significantly. Noah-MP still performs better than NLDAS in most RFCs. However, the superiority is marginal. In three RFCs – namely NE, MA, and OH – Noah-MP notably underperforms NLDAS. The underperformance is mainly due to higher variability than GRACE. Similar to the annual cycle, possible reasons include: (1) Noah-MP overestimated the variability due to unsuitable parameter values, and (2) GRACE underestimated the variability due to the signal leaked from oceans.

4.3 Soil moisture

Figure 4 presents the time series of the surface (0–0.1 m) and root-zone (0–1.0 m) soil moisture in NC, NW, AB, WG, and CB. These RFCs were selected as they have more than 10 valid sites. Table S6 presents the corresponding TSSs. The Noah-MP configurations are consistent in estimating the surface soil moisture, having a spread remarkably smaller than that among the three NLDAS models. The NLDAS models have more diverse representations of soil moisture. Noah is the same as Noah-MP. Both solve Richards' equation to present the soil moisture dynamics in four layers (Niu et al., 2011). Mosaic also solves Richards' equation but at three soil layers. The top layer is further divided into tiles to better represent spatial heterogeneity (Koster and Suarez, 1992). VIC is different from them, utilizing a conceptual soil water tank (Liang et al., 1994). It could be that Noah-MP underestimated the ensemble spread, especially when considering its inability to represent the subgrid heterogeneity. The NLDAS models could also overestimate the spread when considering the conceptual representation of VIC. The spread among the Noah-MP configurations increases significantly from the surface (Fig. 4e) to the root zone (Fig. 4k). The ensemble spread in the root-zone soil moisture reflects the difference in modeling root-water uptake for plant transpiration and soil-bottom drainage as described in Appendix A. Further investigation (Figs. S2–S5) shows that, in the deep layers (the third and fourth layers), Noah-MP has a comparable or greater spread than NLDAS.

Comparison between Noah-MP and Noah shows that they perform similarly in AB (Fig. 4e and f) and WG (Fig. 4g and h) but are different in NC (Fig. 4a and b), NW (Fig. 4c and d), and CB (Fig. 4i and j). The similarity in AB and WG is reasonable, since the two models have similar soil layer structures and parameterizations. The dissimilarity in NC, NW, and CB is most pronounced in winter. It could result from the different snow parameterizations in Noah-MP and Noah, which is investigated in Sect. 4.4.

In the RFCs and soil layers examined in Fig. 4, the Noah-MP ensemble mean performs similarly or better than the NLDAS ensemble mean except in AB and CB. The superiority of the Noah-MP in simulating soil moisture is also reported in previous evaluations (Cai et al., 2014b). In AB and CB, individual NLDAS models do not show a consistent superiority over Noah-MP. In AB, the best NLDAS model is VIC in winter and Noah in summer. The performance of VIC in winter corresponds to the best-performing snow estimation (Fig. 8h). In CB, the best NLDAS model is Noah in winter and Mosaic in summer. Mosaic carefully considers the subgrid variability of soil moisture, which could lead to better skill in RFCs with complex topography such as CB and NW. The NLDAS ensemble mean takes the advantage of wintertime soil moisture from VIC in AB and summertime soil moisture from Mosaic in CB. Both the advantage of VIC and Mosaic come from the representation of subgrid heterogeneity.

Figures 5 and 6 compare the TSS between the NLDAS and Noah-MP ensemble mean at each NAMSD site for the annual cycle and interannual anomaly, respectively. The comparison varies significantly with site and soil layer depth, revealing two major patterns. First, similar to Fig. 4, NLDAS tends to outperform Noah-MP in AB and CB. The high skill of the NLDAS ensemble mean is likely a result of a high ensemble skill gain (Fei et al., 2021) related to the diversity among the NLDAS models. Noah-MP has both a low ensemble spread and an inadequate representation of subgrid heterogeneity in the two RFCs. Second, Noah-MP outperforms NLDAS in other RFCs. In NC, OH, and MA, the low performance of NLDAS is related to the anomaly in wintertime soil moisture (Fig. S2). The anomaly suggests that the NLDAS models generally have difficulty in modeling snow and snow–soil moisture interactions (refer to Sect. 4.4 for more information). On the other hand, Noah-MP has a better snow module, leading to a higher soil moisture estimation skill.

To maximize the utilization of the NLDAS model diversity and Noah-MP physics improvements, we combine the Noah-MP ensemble mean and the three NLDAS models. The right-hand columns of Figs. 5 and 6 show that the four estimates' arithmetic average outperforms the three-model NLDAS ensemble mean at most NASMD sites. The outperformance suggests an added value of the Noah-MP data. If the Noah-MP ensemble mean already outperforms the NLDAS ensemble mean, the added value appears in almost every site. If the Noah-MP ensemble mean underperforms the NLDAS ensemble mean, the added value can still show up at approximately one-third (one-fourth) of the sites for the annual cycle (interannual anomaly).

4.4 Snow water equivalent

Figure 7a presents the spatial patterns of the multi-year averaged SWE (W_snow) from SNODAS. Snow is mainly distributed in the northeast (NE, NC, OH, and MA) and in the mountains of the west (the Cascade Mountains, Rocky Mountains, and Sierra Nevada in NW, AB, MB, WG, CN, and CB). Figure 7b (c) shows the geographical difference between the Noah-MP (NLDAS) ensemble mean and SNODAS. Both Noah-MP and NLDAS exhibit a considerable underestimation in most areas of CONUS. However, the underestimation of Noah-MP is generally smaller. Figure 7d confirms that both Noah-MP and the NLDAS models tend to underestimate SWE in most areas but exhibit an overestimation when snow is extremely thick (SWE is greater than 400 mm). Noah-MP performs better than the NLDAS models in most cases. The superiority of Noah-MP is likely attributable to the three-layer snowpack module, which can represent snow dynamics better in a wide range of snow depth than the single-layer Noah and Mosaic snow module and the quasi-two-layer VIC snow module. A careful consideration of the surface energy balance can also contribute to Noah-MP's superiority over VIC and Mosaic. Consequently, Noah-MP captures the spatial patterns better than NLDAS, with a spatial correlation of 0.87 versus 0.43. Further examination reveals that the superiority of Noah-MP appears in all elevation bands and is the most significant between 1000 to 2000 m with a spatial correlation of 0.85 versus 0.38 (0.77 versus 0.76 below 1000 m, and 0.89 versus 0.75 above 2000 m).

Figure 7 (a) Spatial distribution of the 11-year averaged (September 2004 to August 2015) SWE from SNODAS. (b) Difference between the NLDAS ensemble mean and SNODAS. (c) Difference between the Noah-MP ensemble mean and SNODAS. The spatial correlation coefficients between the NLDAS–Noah-MP ensemble mean and SNODAS are also presented.

Figure 8 As in Fig. 1, but for SWE between September 2004 and August 2015.

Figure 9 TSS of the ensemble means in simulating the annual cycle (left column) and interannual anomaly (right column) of SWE. (a, b) TSS of the NLDAS ensemble mean. (c, d) The difference in TSS between the Noah-MP ensemble mean and the NLDAS ensemble mean. (e, f) The difference between the arithmetic average of the Noah-MP ensemble mean and three NLDAS models to the NLDAS ensemble mean. The evaluation period is 2004–2015.

Figure 8 compares the annual cycle estimated from SNODAS, NLDAS, and Noah-MP. The annual cycle in the 12 RFCs exhibits a similar pattern: it accumulates in winter, peaks in spring, and melts from late spring to summer. The snow season in the northeastern RFCs (i.e., NE, MA, OH, and NC) spans from October to May, whereas the snow season is longer in the mountainous RFCs of the west (i.e., NW, AB, MB, WG, CN, and CB), lasting to June. From the comparison between Noah-MP and NLDAS, we make three observations. First, consistent with Fig. 7, the NLDAS models underestimate the SWE in all RFCs. Among the three NLDAS models, Noah performs the best in the northeast (e.g., NE, MA, and OH), whereas VIC shows some advantages in the western mountainous RFCs (e.g., NW, CN, and CB). Possible reasons for the Noah superiority in the northeast include the careful consideration of the surface energy balance, whereas, in the western mountains, the elevation bands of VIC can better capture the spatial heterogeneity. Second, Noah-MP outperforms the NLDAS models in almost all RFCs, except for showing a marginal degradation in AB and MB. In these two RFCs, VIC outperforms not only the NLDAS models but also Noah-MP. The terrain is hilly, and the elevation-banded parameterization enables a better representation of subgrid snow variability. Another reason is related to the shallow snow in these two RFCs. Noah-MP has known been experiencing negative biases with shallow snow (i.e., AB, MB, and WG). The bias is attributable to the shallow snow albedo as discussed in Dang et al. (2019) and Wang et al. (2020). Third, the ensemble spread among the Noah-MP configurations is small. The estimates of SWE and its uncertainty should be improved in the future by considering processes such as rain–snow partitioning Wang et al. (2019), snow albedo (Wang et al., 2020; Dang et al., 2019), and roughness length (He et al., 2019).

Figure 9 shows the TSS of the NLDAS and Noah-MP ensembles in estimating the annual cycle and interannual anomaly of SWE. Table S7 summarizes the skill scores for the 12 RFCs. The annual cycle and interannual anomaly exhibit similar spatial patterns. NLDAS performs well in most parts of the CONUS, with TSSs higher than 0.75. In MB, TSS reaches a minimum of 0.60. In comparison with NLDAS, Noah-MP performs notably better in the east and west but marginally worse in parts of the central CONUS. The location of Noah-MP's underperformance coincides well with shallow snow and relatively flat terrain. The coincidence hints at the weakness of Noah-MP in simulating shallow snow and the advantage of VIC in representing subgrid snow variability as discussed previously.

We further averaged the 48 Noah-MP configurations and added their average to the three-model NLDAS ensemble. Figure 9e and f show that the four-estimate ensemble mean outperforms the three-model NLDAS ensemble mean in nearly all areas of CONUS, again proving the added value of the data provided in this paper.

Figure 10 As in Fig. 1, but for ET in the period of 1982–2011.

Figure 11 Differences between the Noah-MP ensemble mean and GLEAM in total ET (E), canopy evaporation (E_can), ground evaporation (E_gnd), and transpiration (E_tran). The units are millimeters per day (mm d⁻¹).

Figure 12 Probability distribution of ET's TSS for the annual cycle (a, b) and interannual anomaly (c, d). The left column compares the Noah-MP (orange) and NLDAS (cyan) ensemble means. The right column reveals the Noah-MP (magenta) and NLDAS (dark blue) ensembles. The upper, middle, and lower quantile lines of the magenta boxes show the 75th, 50th, and 25th percentile values of the Noah-MP ensemble. The upper, middle, and lower lines of the dark blue boxes show the three NLDAS models. The blue, green, and cyan dots denote Noah, VIC, and Mosaic, respectively. The evaluation period can be found in Table S1.

4.5 Evapotranspiration

Figure 10 compares the annual cycle estimated from FLUXNET MTE, NLDAS, and Noah-MP in the 12 RFCs. We choose FLUXNET MTE as the reference here, since its performance is superior when compared to AmeriFlux (Tables S2–S4). In all the 12 RFCs, ET peaks during summer and is lowest during winter. Noah-MP successfully captures the timing of the peak in humid RFCs (i.e., NE, MA, OH, LM, SE, NC, and NW) but shows a 1-month lead in a few semi-arid and arid RFCs (i.e., MB, WG, and CN). The average of the three NLDAS models better reproduces the timing of the peak, but the models differ from each other significantly. Among the Noah-MP and NLDAS ensembles, VIC and Mosaic are notably different. VIC exhibits a systematic underestimation, while Mosaic shows an overall overestimation. The 48 Noah-MP configurations and Noah perform closely during autumn and winter, whereas their differences are pronounced during spring and summer. During spring and summer, Noah is the closest to FLUXNET MTE in most RFCs except NE and MA, whereas Noah-MP constantly overestimates the ET in all RFCs.

The overestimation of Noah-MP was investigated by separately comparing the three components of ET (i.e., transpiration, canopy evaporation, and ground evaporation) with GLEAM (Fig. S6). Figure 11 shows the overestimation of total ET is closely linked to the overestimation of ground evaporation, which could be partially attributable to the overly high roughness length for heat and water, as described in Appendix A8 and A9. Besides, the lack of a litter layer (Decker et al., 2017) in Noah-MP could also play a part.

Figure 13 Normalized ensemble spread of the multi-year averaged annual (first row) and seasonal (second–fifth rows) ET from NLDAS (first column) and Noah-MP (second column). The ensemble spread is normalized by the temporal variability of the FLUXNET MTE ET calculated using Eq. (6).

Figure 14 The Sobol' total sensitivity of the multi-year averaged and seasonal (spring – MAM, summer – JJA, autumn – SON, winter – DJF) ET to the four parameterizations: stomatal conductance, soil moisture limitation to transpiration (β factor), turbulence, and runoff. Higher values indicate higher sensitivities.

Figure 12 evaluates Noah-MP and NLDAS using the 25 AmeriFlux sites. The NLDAS ensemble mean outperforms the Noah-MP ensemble mean for the annual cycle, and this outperformance results from two causes. First, an NLDAS member, Noah, performs the closest to the observations, as shown in Figs. 12b and 10. Second, the three NLDAS models are remarkably different from each other. The diversity of the ensemble gives a higher skill gain by combining them, as shown by the difference between the ensemble mean skill (Fig. 12a) and the median TSS (Fig. 12b). On the other hand, the Noah-MP configurations are too similar to each other, and all have a positive bias (Fig. 10). However, for the interannual anomaly, the Noah-MP ensemble mean slightly outperforms the NLDAS ensemble mean (Fig. 12c). Figure 12d shows that the Noah-MP configurations marginally outperform the NLDAS models. Among the NLDAS models, VIC performs the best, and Noah does not exhibit the same superiority shown in the annual cycle. The difference between the NLDAS ensemble mean performance (Fig. 12c) and median performance (Fig. 12d) is marginal, suggesting that the NLDAS ensemble skill gains are not notable for the interannual anomaly.

Figure 13 examines the ensemble spread of Noah-MP and NLDAS. The ensemble spread is normalized by the temporal variability calculated using the FLUXNET MTE ET. NLDAS has a significant spread in the southeast and west in all seasons, while spring shows the largest value. As seen in Fig. 10, the NLDAS ensemble spread mainly reflects the differences between VIC and Mosaic. The Noah-MP ensemble has a notably smaller spread than NLDAS. The Noah-MP ensemble spread is manifested in spring and summer in the southeastern (SE, LM, and WG) and western (CN, CB, and NW) RFCs.

We can decompose the Noah-MP ensemble spread and pinpoint the dominant process using Sobol' sensitivity analysis (Zheng et al., 2019). Figure 14 delineates the Sobol' total sensitivity index of total ET to the four processes described in Appendix A. In spring and summer, for the regions where the Noah-MP configurations show significant spread (SE, LM, WG, CB, CN, and NW) (Fig. 13), ET is most sensitive to the parameterization of stomatal conductance (Fig. 14e and i) and then to the β factor (Fig. 14f and j). However, for regions with positive biases (NC, OH, and LM, as shown in Fig. 11b and c), the Noah-MP estimation is more sensitive to the turbulence parameterizations (Fig. 14c and g). During autumn and winter, the parameterizations of stomatal conductance (Fig. 14m and q) and β factor (Fig. 14n and r) still have significant impacts on the estimation of ET, and these impacts could be a result of the “memory” of TWS (Zheng et al., 2019). Besides these two processes, the runoff parameterization is dominant during autumn in the east (Fig. 14p), and the turbulence parameterization is dominant during winter (Fig. 14s).

5 Code and data availability

The dataset is freely available for download from the Zenodo online repository at https://doi.org/10.5281/zenodo.7109816 (Zheng et al., 2022). The dataset (along with datasets on which it is based) is subject to a Creative Commons BY (attribution) license agreement (https://creativecommons.org/licenses, last access: 16 August 2021).

6 Conclusions

This paper describes a $\mathrm{1}/\mathrm{8}$^∘ dataset of the TWB over the CONUS from 1980 to 2015 simulated from an ensemble of 48 perturbed-physics configurations of Noah-MP. This Noah-MP multi-physics ensemble features an enrichment of the NLDAS-2 four-model ensemble and brings convenience for multi-model comparison. The dataset has already been used in the monitoring of groundwater storage change (Rateb et al., 2020), the analysis of LSM parameterization sensitivity (Zheng et al., 2019), the development of model evaluation method (Zheng et al., 2020), and hydrological ensemble simulations (Fei et al., 2021). This paper details the Noah-MP parameterizations employed and evaluates the estimated TWSA, soil moisture, SWE, and ET in comparison with the NLDAS ensemble.

The spread of the ensemble estimation is the largest for the runoff. The spread in surface runoff accounts for 34 % of its climatological mean. The spread is comparable to the previous estimates for multi-model ensembles (Dirmeyer et al., 2006). The ensemble spread in snow water equivalent is the smallest, 2.5 % of its climatological mean. The ensemble has not included different parameterizations of several snow processes such as rain–snow partitioning, snow albedo, and roughness length, which could lead to an underestimation of the ensemble spread. The underestimation of the ensemble spread becomes more apparent when the Noah-MP ensemble mean is biased (Fig. 8). The bias is more pronounced relative to the NLDAS models in parts of AB and MB where snow is shallow and the terrain is hilly (Fig. 9). VIC performs better there, suggesting the importance of considering subgrid variability.

Evaluation against various reference data shows that Noah-MP generally performs better than the NLDAS models. The augmented three-layer snow module of Noah-MP significantly improves the estimation of snow and also wintertime soil moisture. On the other hand, the NLDAS models outperform Noah-MP in AB and CB for surface soil moisture (Figs. 4 to 6), in AB and MB for snow (Figs. 8 and 9), and in OH, LM, and NC for the annual cycle of ET (Fig. 10). The outperformance of NLDAS is likely attributable to the consideration of subgrid variability in soil moisture with Mosaic and in snow with VIC. The Noah-MP ensemble could be improved by increasing the spatial resolution or developing parameterizations of subgrid heterogeneity. Noah-MP also underestimates the temporal variability of TWS in coastal RFCs (Figs. 2 and 3). Correction of the ocean signal leakage in the GRACE data and representation of spatially varying unconfined aquifers' parameters should be beneficial. For the annual cycle of ET, there is a systematic overestimation in spring and summer. Sobol' sensitivity analysis of the Noah-MP ensemble reveals that the bias is mainly related to the parameterization of turbulence. We have examined the code and found that the implementation is inconsistent with the literature. The parameterization of the roughness length of heat and water vapor likely contributes to the ET overestimation.

The Noah-MP ensemble shares the same atmospheric forcing and static parameters with the NLDAS models. The similarity enables the comparison between the multi-model and perturbed-physics ensemble methods as shown in Fei et al. (2021). Besides, Noah-MP complements the NLDAS models well. Adding Noah-MP to the NLDAS ensemble can consistently improve the TWB variables in most areas of CONUS.

The Noah-MP model has been undergoing rapid development. New components such as plant hydraulics (Li et al., 2021), roughness sublayer (Abolafia-Rosenzweig et al., 2021), crops (Liu et al., 2016), and dynamic rooting depth (Liu et al., 2020) have been added. New schemes for the processes such as rain–snow partitioning (Wang et al., 2019) have been included. Parameterizations of surface roughness length (He et al., 2019; Zhang et al., 2021), snow albedo (Wang et al., 2020), and vertical soil layers (Zhao et al., 2022; Shellito et al., 2020) have been refined. The dataset can be improved by using the updated model and including more perturbations.

Appendix A: Formulation of the used Noah-MP parameterization schemes

A1 SIMGM runoff parameterization scheme

SIMGM is a TOPMODEL-based runoff model (Niu et al., 2007). The scheme parameterizes runoff (R_srf and R_sub) as an exponential function of groundwater table depth (z_wt, m, positive down) as follows.

$$\begin{array}{}\text{(A1)}& {\displaystyle }{R}_{\mathrm{srf}}=Q_{\mathrm{soil},\mathrm{srf}}\left[\right(\mathrm{1}-f_{\mathrm{frz},\mathrm{1}})f_{\mathrm{sat}}+f_{\mathrm{frz},\mathrm{1}}],\text{(A2)}& {\displaystyle }{f}_{\mathrm{sat}}=f_{\mathrm{sat},\max }\exp [-\mathrm{0.5}f(z_{\mathrm{wt}}-z_{\mathrm{bot}}\left)\right],\text{(A3)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{sub}}& =[\mathrm{1}-\underset{i=\mathrm{1},\mathrm{\dots },N_{\mathrm{soil}}}{max}(f_{\mathrm{frz},i}\left)\right]R_{\mathrm{sub},\max }\\ & \exp [-\mathrm{\Lambda }-f(z_{\mathrm{wt}}-z_{\mathrm{bot}}\left)\right],\end{array}\end{array}$$

where Q_soil,srf is the water incident on the soil surface (the sum of precipitation throughfall, snowmelt, and dewfall; $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$); f_frz,i is the fractional frozen area of the ith soil layer (m² m⁻²), which is parameterized using the frozen water content of the soil layer following Niu and Yang (2006); f_sat is the saturation fraction of the grid cell (m² m⁻²); and z_bot is the depth of the soil column bottom (2 m in this study), and z_wt is the groundwater table depth (m), which is converted from the groundwater storage by a specific-yield parameter. The groundwater storage is predicted using a dynamic groundwater model interacting with the soil column bottom (Niu et al., 2007).

The scheme has four calibratable parameters: (1) f_sat,max, the maximum saturation area fraction (m² m⁻²), which is defined as the cumulative distribution function of the topographic index when the grid-cell-mean water table depth is zero; (2) f, a runoff decay factor (unitless); (3) R_sub,max, the maximum subsurface runoff when the grid-cell-mean water table depth is zero ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$); and (4) Λ, the grid-cell-mean topographic index (unitless). In this study, the parameters have the following values: $f_{\mathrm{sat},\max }=\mathrm{0.38}$ m² m⁻², f=6, $R_{\mathrm{sub},\max }=\mathrm{5}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and Λ=10.5.

A2 SIMTOP runoff parameterization scheme

SIMTOP is also a TOPMODEL-based runoff parameterization scheme, the same as SIMGM (Eqs. A1–A3). The major difference between SIMTOP and SIMGM is that SIMTOP parameterizes the groundwater table depth (z_wt) using the soil liquid water content by assuming the water head is at equilibrium throughout the soil column down to the water table (Niu et al., 2005). Although SIMTOP and SIMGM share the same conceptual model of runoff generation, implementation differences exist. First, in contrast to Eqs. (A2) and (A3), SIMTOP does not use the soil column bottom depth (z_bot) in calculating the saturation area fraction (f_sat) and subsurface runoff:

$$\begin{array}{}\text{(A4)}& {\displaystyle }{f}_{\mathrm{sat}}=f_{\mathrm{sat},\max }\exp (-\mathrm{0.5}f{z}_{\mathrm{wt}}),\text{(A5)}& {\displaystyle }\begin{array}{rl}{R}_{\mathrm{sub}}& =[\mathrm{1}-\underset{i=\mathrm{1},\mathrm{\cdots },N_{\mathrm{soil}}}{max}(f_{\mathrm{frz},i}\left)\right]R_{\mathrm{sub},\max }\\ & \exp (-\mathrm{\Lambda }-f{z}_{\mathrm{wt}}).\end{array}\end{array}$$

Second, parameter values are slightly different for the runoff decay factor and maximum subsurface runoff: f=2, and $R_{\mathrm{sub},\max }=\mathrm{4}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$.

A3 NOAHR runoff parameterization scheme

The Noah runoff (NOAHR) scheme parameterizes surface runoff (R_srf) as infiltration excess:

$$\begin{array}{}\text{(A6)}& R_{\mathrm{srf}}=Q_{\mathrm{soil},\mathrm{srf}}-Q_{\mathrm{soil},\mathrm{in}},\end{array}$$

where Q_soil,in is the infiltration into the soil ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). The infiltration is derived from the approximate solution to the Richards equation following Philip (1969) with additional considerations of the spatial variability of precipitation and infiltration capacity. By assuming exponential and independent distributions of precipitation and infiltration capacity within a model grid cell, NOAHR formulates the soil infiltration as follows:

$$\begin{array}{}\text{(A7)}& {\displaystyle }{Q}_{\mathrm{soil},\mathrm{in}}=Q_{\mathrm{soil},\mathrm{srf}}{\displaystyle \frac{I_c}{Q_{\mathrm{soil},\mathrm{srf}}\mathrm{\Delta }t+I_c}},\text{(A8)}& {\displaystyle }{I}_{\mathrm{c}}=w_{\mathrm{d}}[\mathrm{1}-\exp (-K_{\mathrm{\Delta }t}\mathrm{\Delta }t\left)\right],\text{(A9)}& {\displaystyle }{w}_{\mathrm{d}}=\sum _{i=\mathrm{1}}^{N_{\mathrm{soil}}}{\mathit{\rho }}_{\mathrm{wat}}(w_{\mathrm{sat},i}-w_{\mathrm{soil},i})\mathrm{\Delta }{z}_{\mathrm{soil},i},\end{array}$$

where I_c is the soil infiltration capacity of the model grid cell (kg m⁻²), w_d is the water deficit of the soil column (kg m⁻²), and Δt is the model time step (s). Following Chen and Dudhia (2001a), the parameter is assumed as propositional to the saturated hydraulic conductivity of the first soil layer (K_sat,1; $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$):

$$\begin{array}{}\text{(A10)}& K_{\mathrm{\Delta }t}={\displaystyle \frac{K_{\mathrm{\Delta }t,\mathrm{ref}}}{k_{\mathrm{ref}}}}{K}_{\mathrm{sat},\mathrm{1}},\end{array}$$

where K_Δt,ref and k_ref are two parameters. In Noah-MP (and Noah), $K_{\mathrm{\Delta }t,\mathrm{ref}}=\frac{\mathrm{3}}{\mathrm{86}\mathrm{400}}{\mathrm{s}}^{-\mathrm{1}}$, and $k_{\mathrm{ref}}=\mathrm{2}\times {\mathrm{10}}^{-\mathrm{3}}$ kg m⁻² s⁻¹. K_sat,1 is assigned using a soil parameter lookup table according to the soil texture type.

NOAHR assumes free drainage at the soil column bottom. The subsurface runoff is calculated as

$$\begin{array}{}\text{(A11)}& R_{\mathrm{sub}}={\mathit{\alpha }}_{\mathrm{slope}}{K}_{\mathrm{soil},N_{\mathrm{soil}}},\end{array}$$

where α_slope is the terrain slope index, which is arbitrarily given as 0.1 in the adopted version of Noah-MP. $K_{\mathrm{soil},N_{\mathrm{soil}}}$ is the hydraulic conductivity of the bottom soil layer, which is parameterized following Clapp and Hornberger (1978).

A4 BATS runoff parameterization scheme

The BATS scheme parameterizes surface runoff (R_srf) as a function of soil wetness (Yang and Dickinson, 1996):

$$\begin{array}{}\text{(A12)}& {\displaystyle }{R}_{\mathrm{srf}}=Q_{\mathrm{soil},\mathrm{srf}}\left[(\mathrm{1}-f_{\mathrm{frz},\mathrm{1}})f_{\mathrm{sat}}+f_{\mathrm{frz},\mathrm{1}}\right],\text{(A13)}& {\displaystyle }{f}_{\mathrm{sat}}={\mathit{\theta }}^{\mathrm{4}},\text{(A14)}& {\displaystyle }\mathit{\theta }={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{soil}}}\frac{w_{\mathrm{soil},i}}{w_{\mathrm{sat},i}}\mathrm{\Delta }{z}_{\mathrm{soil},i}}{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{soil}}}\mathrm{\Delta }{z}_{\mathrm{soil},i}}},\end{array}$$

where θ is the averaged wetness throughout the soil column (m³ m⁻³).

Similar to NOAHR, the BATS scheme also assumes a free drainage boundary condition at the soil column bottom. Subsurface runoff (R_sub) is parameterized as follows:

$$\begin{array}{}\text{(A15)}& R_{\mathrm{sub}}=\left(\mathrm{1}-\underset{i=\mathrm{1},\mathrm{\cdots },N_{\mathrm{soil}}}{max}\left(f_{\mathrm{frz},i}\right)\right)K_{\mathrm{soil},N_{\mathrm{soil}}}.\end{array}$$

A5 Ball–Berry scheme of stomatal resistance

Leaf stomata are the small pores typically found on the underside of leaves. They control the gas exchange of CO₂, H₂O, and O₂ between the internal leaf structure and the external atmosphere. In LSMs, the opening and closing of the stomata are characterized by stomatal conductance.

The Ball–Berry scheme for parameterizing the stomatal conductance (g_s) for H₂O is as follows:

$$\begin{array}{}\text{(A16)}& g_{\mathrm{s}}=m{\displaystyle \frac{A}{c_{\mathrm{s}}}}{\displaystyle \frac{e_{\mathrm{s}}}{e_i}}{P}_{\mathrm{atm}}+b,\end{array}$$

where g_s is the leaf stomatal conductance ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$), m is a vegetation-type dependent parameter (unitless), A is the leaf photosynthesis rate, c_s is the CO₂ partial pressure at the leaf surface (Pa), e_s is the water vapor pressure at the leaf surface (Pa), e_i is the saturated water vapor at the stomata (Pa), P_atm is the ambient air pressure (Pa), and b is the stomatal conductance at zero photosynthesis ($\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$). The parameters m and b are assigned from a lookup table using the vegetation type.

A6 Jarvis scheme of stomatal resistance

The Jarvis scheme for parameterizing the canopy resistance (R_c) based on the product of four stress factors (s m⁻¹) is calculated as follows (Chen et al., 1996; Sellers et al., 1996; Jacquemin and Noilhan, 1990; Jarvis, 1976):

$$\begin{array}{}\text{(A17)}& {\displaystyle }{R}_{\mathrm{c}}=& {\displaystyle }{R}_{\mathrm{c},\min }{\displaystyle \frac{\mathrm{1}}{f_{\mathrm{1}}{f}_{\mathrm{2}}{f}_{\mathrm{3}}\mathit{\beta }}},\text{(A18)}& {\displaystyle }{f}_{\mathrm{1}}=& {\displaystyle \frac{\frac{R_{\mathrm{c},\min }}{R_{\mathrm{c},\max }}+f}{\mathrm{1}+f}},\text{(A19)}& {\displaystyle }f=& {\displaystyle }\mathrm{0.55}{\displaystyle \frac{\mathrm{2}{R}_{\mathrm{g}}}{R_{\mathrm{gl}}}},\text{(A20)}& {\displaystyle }{f}_{\mathrm{2}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{1}+h_{\mathrm{s}}\left[q_{\mathrm{sat}}\right(T_{\mathrm{l}})-q_{\mathrm{a}}]}},\text{(A21)}& {\displaystyle }{f}_{\mathrm{3}}=& {\displaystyle }\mathrm{1}-\mathrm{0.0016}(T_{\mathrm{ref}}-T_{\mathrm{l}}{)}^{\mathrm{2}},\end{array}$$

where f₁, f₂ and f₃ are the stress factors of solar radiation, vapor pressure deficit, and air temperature, respectively (unitless), which are unitless and range from 0 to 1; β is the soil moisture stress factor, which is detailed in Sect. A7; R_g is the incoming solar radiation (W m⁻²) for unit leaf area index; T_l is leaf temperature (K). q_sat(T_l) (kg kg⁻¹) and q_a (kg kg⁻¹) are the saturated humidity at the temperature of T_l and ambient humidity, respectively. In literature (Chen et al., 1996; Jacquemin and Noilhan, 1990), q_sat(T_l) and q_a are specific humidity, whereas in Noah-MP v3.6, they are implemented as mixing ratio.

The scheme has five parameters: R_c,min, the minimum stomatal resistance (s m⁻¹) per unit leaf area index; R_c,max, the maximum resistance; R_gl, a radiation scaling factor (unitless); h_s, a humidity scaling factor (unitless); T_ref, the optimum temperature (K). Among these parameters, R_c,min, R_gl, and h_s are assigned using a vegetation-parameter lookup table, while R_c,max and T_ref are assumedly assigned to 5000 s m⁻¹ and 298 K, respectively.

A7 Three soil moisture stress factor schemes

The Noah β-factor (NOAHB) scheme parameterizes the soil moisture stress factor controlling transpiration (β factor) as a function of soil moisture, which is calculated as follows:

$$\begin{array}{}\text{(A22)}& \mathit{\beta }=\sum _{i=\mathrm{1}}^{N_{\mathrm{root}}}{\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{soil},i}}{z_{\mathrm{root}}}}min\left(\mathrm{1},{\displaystyle \frac{{\mathit{\theta }}_i-{\mathit{\theta }}_{\mathrm{wilt}}}{{\mathit{\theta }}_{\mathrm{ref}}-{\mathit{\theta }}_{\mathrm{wilt}}}}\right),\end{array}$$

where N_root is the total number of soil layers that contain roots, z_root is the total depth of the root-zone layer (m), and θ_i is the volumetric soil moisture of the ith soil layer (m³ m⁻³). NOAHB has two parameters: θ_ref, the field capacity (m³ m⁻³); and θ_wilt, the wilting volumetric soil moisture (m³ m⁻³).

The Community Land Model (CLM) scheme (Oleson et al., 2004) parameterizes β as a function of soil matric potential, which is calculated as follows:

$$\begin{array}{}\text{(A23)}& \mathit{\beta }=\sum _{i=\mathrm{1}}^{N_{\mathrm{root}}}{\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{soil},i}}{z_{\mathrm{root}}}}min\left(\mathrm{1},{\displaystyle \frac{{\mathit{\psi }}_{\mathrm{wilt}}-{\mathit{\psi }}_i}{{\mathit{\psi }}_{\mathrm{wilt}}-{\mathit{\psi }}_{\mathrm{sat}}}}\right),\end{array}$$

where ψ_i is the water pressure head of the ith soil layer (m), and ψ_i is converted from θ_i using the formula of Clapp and Hornberger (1978). CLM has two parameters: ψ_sat, the saturated water pressure head (m); and ψ_wilt, the wilting pressure head (m).

The Simplified Simple Biosphere (SSiB) scheme (Xue et al., 1991) also parameterizes the β factor as a function of the soil pressure head, similar to CLM. However, the formula is different, as follows:

$$\begin{array}{}\text{(A24)}& \mathit{\beta }=\sum _{i=\mathrm{1}}^{N_{\mathrm{root}}}{\displaystyle \frac{\mathrm{\Delta }{z}_{\mathrm{soil},i}}{z_{\mathrm{root}}}}min\left[\mathrm{1},\mathrm{1}-\exp \left(-c_{\mathrm{2}}\ln \left({\displaystyle \frac{{\mathit{\psi }}_{\mathrm{wilt}}}{{\mathit{\psi }}_i}}\right)\right)\right].\end{array}$$

SSiB has two parameters: ψ_wilt, the wilting pressure head (m); and c₂, a unitless coefficient.

In Noah-MP version 3.6, the parameters θ_sat, θ_wilt, and ψ_sat are assigned using a soil parameter lookup table (Table 2 of Chen and Dudhia, 2001a); ψ_wilt is −10 m, independent of vegetation and soil types (Niu et al., 2011); c₂ is assumed constant at 5.8, whereas in SSiB, this parameter varies with vegetation type (Xue et al., 1991, Table 2).

A8 Chen97 near-surface turbulence scheme

The Chen97 scheme (Chen et al., 1997) parameterizes the surface exchange coefficient for heat (C_h) as follows:

$$\begin{array}{}\text{(A25)}& \begin{array}{rl}{C}_{\mathrm{h}}& ={\mathit{\kappa }}^{\mathrm{2}}{\left[\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}\mathrm{m}}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z}{L}}\right)+{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{m}}}{L}}\right)\right]}^{-\mathrm{1}}\\ & {\left[\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}\mathrm{h}}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z}{L}}\right)+{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{h}}}{L}}\right)\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where κ=0.4 is the von Kármán constant; L is the Monin–Obukhov (M–O) length (m); z is the reference height (m); Ψ_m and Ψ_h are the similarity theory-based stability functions for momentum and heat, respectively; z_0m is the roughness length for momentum (m) and depends on the land cover/land-use type; and z_0h is the roughness length for heat (m). Niu et al. (2011) parameterized $z_{\mathrm{0}\mathrm{h}}=z_{\mathrm{0}\mathrm{m}}\exp (-\mathit{\kappa }C\sqrt{R{e}^{*}})$, where C=0.1 and Re^* is the roughness Reynolds number. However, in the code of Noah-MP version 3.6, z_0h=z_0m.

A9 M–O near-surface turbulence scheme

The M–O scheme is based on the M–O similarity theory (Brutsaert, 1982), which parameterizes C_h as follows:

$$\begin{array}{}\text{(A26)}& \begin{array}{rl}{C}_h& ={\mathit{\kappa }}^{\mathrm{2}}{\left[\ln \left({\displaystyle \frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{m}}}}\right)-{\mathrm{\Psi }}_{\mathrm{m}}\left({\displaystyle \frac{z-d_{\mathrm{0}}}{L}}\right)\right]}^{-\mathrm{1}}\\ & {\left[\ln \left({\displaystyle \frac{z-d_{\mathrm{0}}}{z_{\mathrm{0}\mathrm{h}}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left({\displaystyle \frac{z-d_{\mathrm{0}}}{L}}\right)\right]}^{-\mathrm{1}},\end{array}\end{array}$$

where z_0h=z_0m. d₀ is the zero-displacement height (m),

$$\begin{array}{}\text{(A27)}& d_{\mathrm{0}}=\mathrm{0.64}{z}_{\mathrm{ct}},\end{array}$$

where z_ct is the canopy top height (m).

Appendix B: Estimation of the terrestrial water storage for the RFCs neighboring the Great Lakes

The TWS estimation for the NC, OH, and NE RFCs is performed in two steps: (1) aggregate the GRACE TWS over both the RFC land area and neighboring lakes (lakes Superior, Michigan, and Huron for NC; Erie for OH; and Ontario for NE); and (2) subtract the lake water storage anomaly from the aggregated TWS. The lake water storage in the second step is calculated as the product of the observed water level and the lake area.

The lake water level is an arithmetic average of selected NOAA in situ observations (https://tidesandcurrents.noaa.gov/stations.html?type=Water+Levels, last access: 15 February 2022). For Lake Superior, five observation stations were selected: Point Iroquois, Marquette C.G., Ontonagon, Duluth, and Grand Marais. For Lake Michigan, seven stations were selected: Ludington, Holland, Calumet Harbor, Milwaukee, Kewaunee, Sturgeon Bay Canal, and Port Inland. For Lake Huron, five stations were selected: Lakeport, Harbor Beach, Essexville, Mackinaw City, and De Tour Village. For Lake Erie, eight stations were selected: Buffalo, Sturgeon Point, Erie, Fairport, Cleveland, Marblehead, Toledo, and Fermi Power Plant. And for Lake Ontario, four stations were selected: Cape Vincent, Oswego, Rochester, and Olcott.

The lake area is estimated from the lake boundary data provided by the United States Geological Survey (https://www.sciencebase.gov/catalog/item/530f8a0ee4b0e7e46bd300dd, last access: 15 February 2022). Only the area within the United States is considered, which is within a 150 km radius from the studied RFCs. The lake areas are calculated as follows: 52 441 km² for Lake Superior within the USA, 57 509 km² for Lake Michigan, 23 185 km² for Lake Huron within the USA, 25 494 km² for Lake Erie, and 18 871 km² for Lake Ontario. Month-to-month variations in lake area are neglected in this study for simplicity.