improve diagnosis of weather and climate model biases and enhance the accuracy of their future projections (Williams, 2019).

Remote-sensing data products of SM and cloud vertical profiles from polar-orbiting Earth-observing satellites have advanced considerably in the past decades and have the potential to substantially advance our understanding of SM-cloud-precipitation relationships. Specifically, the 3-hourly/9 km Soil Moisture Active Passive (SMAP/ L4) and 1.5-hourly/5 km Global Precipitation Measurement Dual-Frequency Precipitation Radar (GPM/DPR/ L2A) provide high-resolution estimates of the topsoil moisture content and the vertical distribution of hydrometeors within and above the PBL, respectively, at a global coverage. Many studies have confirmed the accuracy and reliability of SMAP/L4 (Koster et al., 2018;Reichle et al., 2017;Tavakol et al., 2019;L. Zhang et al., 2017;X. Zhang et al., 2017) and GPM/DPR/L2A (Lasser et al., 2019;Liao & Meneghini, 2022;Pejcic et al., 2020) data products.

In this paper, we use functional decomposition of a large database of collocated SMAP/L4 surface SM and GPM/ DPR/L2A cloud vertical profiles (CVPs) to investigate the functional relationship between antecedent SM and the vertical distribution of cloud water and reflectivity in the lower troposphere. Our method, called high-dimensional model representation or HDMR (Gao et al., 2023;Li & Rabitz, 2010, 2012) expresses all variable interactions in a hierarchical order and uses the superposition principle to disentangle individual and cooperative contributions of SM and other land-surface variables to the CVP. We are particularly interested in the so-called first-order component functions as they quantify each land-surface variable's direct contribution to the CVP. As byproduct, HDMR yields maps for our study area of the SM contribution to cloud reflectivity and rainfall, as a function of cloud height and SM time lag. Note that the CVP state decomposition is necessarily incomplete as (a) CVPs are highly variable due to turbulent local dynamics (Heinze et al., 2017) and (b) we consider only a handful of governing variables. The resulting component functions are not intended for CVP prediction, hence we do not report summary statistics of the quality of fit of the HDMR decomposition in the remainder of this paper. The coefficient of determination is low (R 2 = 0.14 for HDMR and R 2 = 0.05 for linear regression). Rather, we aim to unravel the contribution of each macroscopic land-surface and atmospheric variable to the CVP, which can be robust and statistically significant when analyzing sufficient observational samples. This can help diagnose weather and climate model biases.

This paper is organized as follows. Section 2 discusses the SMAP/L4 SM and GPM/DPR/L2A satellite products and study region. Section 3 describes the data preprocessing steps and briefly reviews the HDMR method. Section 4 presents the results of our analysis and documents the relationship between SM and the CVP as a function of cloud height, time lag, and location in our study region. Section 5 summarizes our main findings and presents suggestions for future work.

We use the publicly available 3-hourly/9 km SMAP/L4 (March 2015-present) and 1.5-hourly/5 km GPM/DPR/ L2A (March 2014-present) data products and select samples from our study region in the warm seasons (April to October) from 2016 to 2019, focusing on convective precipitation in the afternoon hours until midnight (14:00-24:00 CDT). The altitude extends from 1 to 5 km with Findell and Eltahir (2003a) identifying the 1-3 km zone as a critical region for convective triggering. We briefly discuss the SMAP/L4 and GPM/DPR/L2A products and our study region. A more detailed description of the satellite data products is found in Text S1 in Supporting Information S1.

The SMAP mission Level 4 SM (L4_SM) product gives 3-hourly estimates of surface (0-5 cm) SM (see Figure 1a), root-zone SM, and other land-surface variables at 9-km spatial resolution and global coverage (Reichle et al., 2015). As our preliminary data analysis showed that cloud reflectivity was more responsive to topsoil SM than to root-zone SM, we chose topsoil SM as the key input variable in our CVP functional decomposition. The reasons for this are a stronger diurnal variation and larger root density at the surface (Oerter et al., 2021). As auxiliary land variables, we consider land-surface temperature (LST) and leaf area index (LAI). In addition, hourly estimates of low-level atmospheric temperature (AT) and total precipitable water (TPW) from 0.25 °× 0.25°ERA-5 reanalysis are adopted as precursors to mesoscale convective events (Findell & Eltahir, 2003a;Holloway & Neelin, 2010;Sherwood, 1999). We use the mean AT for the critical region, 1-3 km above the soil surface. Section 3.2 discusses in more detail our selection of auxiliary land-surface and atmospheric variables.

The GPM/DPR/L2A product (GPM_2ADPR) provides a swath of precipitation profiles (see Figure 1b) every 1.5 hr at a horizontal resolution of 5 km and vertical increment of 125 m. The major data fields 'zFactorFinal (dBZ)' and 'typePrecip' provide vertical profiles of Ka-band cloud reflectivity (see Figure 1c) and an 8-digit precipitation type ID. We only select samples classified as convective precipitation and use 250-m vertically averaged cloud reflectivities to minimize the impact of measurement errors. This averaging demonstrated to be rather inconsequential in our analysis.

Our study region is displayed in Figure 1a (95°W-105°W, 32°N-40°N) and is a hot spot for SM-precipitation coupling (Findell & Eltahir, 2003b;Ford et al., 2023;Koster et al., 2004) with large spatial variability in climatological sign and strength of the SMCPF (Findell et al., 2011;Ford et al., 2023;Frye & Mote, 2010;Su & Dickinson, 2017;Yuan et al., 2020). This central US region offers an excellent demonstration of our methodology and allows us to benchmark the inferred patterns of the SMCPF sign and magnitude against literature findings.

We extract the GPM/DPR/L2A swaths that overpass our study region and single out samples classified as convective precipitation in the 'typePrecip' data field. This type classification is crucial to an accurate characterization of the antecedent atmosphere using ERA-5 reanalysis AT and TPW data. To avoid water from interception evaporation, we discard samples which, in the 18 hr leading up to the DPR's scan, received more than 0.5 mm of precipitation according to the Multi-Radar Multi-Sensor (MRMS) Gauge-corrected Quantitative Precipitation Estimates (J. Zhang et al., 2016). This should also reduce the impact of large-scale synoptic systems (Findell et al., 2011). Next, we project SMAP/L4 and ERA-5 data onto GPM/DPR/L2A coordinates using linear interpolation. Time lags Δt = t dprt smap or Δt = t dprt era5 of 7 and 10 hr are chosen, where t dpr is the DPR scan time, t smap denotes the midpoint of the 3-hr averaging interval, and t era5 is the ERA-5 time of hourly snapshots of the atmospheric conditions. In doing so, we allow for a 2-hr grace period so as to maximize the sample size n at each DPR measurement height. For example, SM data with a time lag 6.01 ≤ Δt ≤ 7.99 are pooled together in the 7-hr time lag. The two SM time lags find support by our treatment of the synoptic effect (Δt < 18 hr), and are in agreement with data-driven studies on the diurnal SMCPF (Findell et al., 2011;Welty & Zeng, 2018). Figure 1d displays the number of DPR-measured cloud reflectivities n in the months of April-October (2016-2019) as a function of cloud height. The sample size is not constant with height due to for instance, cloud absence, radar detection limitations and path attenuation (Iguchi et al., 2010). Nevertheless, the pooled reflectivities from April-October guarantee a sufficiently large sample size at each cloud height. Next, we decompose this final collection of SMAP/L4 -GPM/DPR/L2A samples using HDMR and expand the DPR-measured cloud reflectivities at each cloud height as a sum of first-and higher-order structural and correlative contributions of SM and the auxiliary variables.

SMCPFs are notoriously difficult to observe outside of model environments (Ford et al., 2023), hence innovative analytical approaches are required to study them (Berg et al., 2013;Findell et al., 2011;Guillod et al., 2014;Knist et al., 2017;Koster et al., 2004;Seneviratne et al., 2006). HDMR is particularly appealing in this context. Unlike methods such as multivariate linear regression and correlation analysis (Ford et al., 2023;G. Wang et al., 2024;Welty & Zeng, 2018), HDMR accounts explicitly for correlation among the input variables (which is common) and expresses all variable interactions in a hierarchical order. This allows us to disentangle direct and cooperative contributions of individual and groups of dependent input variables to the CVP.

Suppose we group all land-surface and atmospheric variables that govern cloud reflectivity y = f (x) at a given altitude in a d × 1 vector x = (x 1 , … , x d ) ⊤ . HDMR builds on the multivariable function expansion of Soboľ. (1993) and decomposes the output, y = f (x), of a scalar-valued square-integrable function, (Li & Rabitz, 2012)

where f 0 denotes the mean output and ϵ ∼ N(0, σ 2 ϵ ) is a normally distributed residual with zero-mean and constant variance, σ 2 ϵ . The n 1 = d first-order functions, f i (x i ) , characterize the individual effects of the input variables on cloud reflectivity, y. The

, characterize the cooperative contributions of two, three, up to all land-surface variables combined to y. In physical systems, third-and higher-order effects are usually negligible (Falchi et al., 2018;Gao et al., 2023;Kucherenko et al., 2011;Rabitz & Aliş, 1999;Shereena & Rao, 2019;H. Wang et al., 2017) and, therefore, we consider only the n 12 = n 1 + n 2 first-and second-order component functions in our CVP function expansion

where subscript u stands for component function index rather than its order as in Equation 1. Thus, f 0 signifies the mean reflectivity in units of dBZ, and f 1 , … , f d , and, f d+1 , … , f d + d(d-1)/ 2 , are the first-and second-order component functions, respectively, which quantify the individual and bivariate contributions of the land-surface and atmospheric variables to cloud reflectivity, y.

We follow Gao et al. (2023) and write the first-and second-order component functions of Equation 2as a sum of linear multiples of orthonormalized polynomial functions of degrees 1 to p = 3 (Li & Rabitz, 2012). This formulation with extended bases and expansion coefficients (linear multiples) helps satisfy a so-called relaxed vanishing condition (Hooker, 2007) Geophysical Research Letters

where u is a subset of superset U = {1, … , d}, x u denotes the dimensions u of the input vector and w u (x u ) is the probability density function (pdf) of x u . The vanishing condition (3) dictates that a second-order component function, f ij (x i , x j ) must be orthogonal to its lower order component functions, f i (x i ) and f j (x j ) , and guarantees an exact delineation of the structural and correlative contributions of single and groups of input variables to y (Gao et al., 2023;Li & Rabitz, 2012). D-MORPH regression (Li & Rabitz, 2010) enforces the vanishing condition (Li & Rabitz, 2010) in its search for the optimum expansion coefficients. This method is described in Text S2 in Supporting Information S1.

The statistical significance of a particular component function is readily determined by comparing the performance of the function expansion with and without this component function. Let SSR 1 be the sum of the squared residuals of the function

is the quantile function of the Fisher-Snedecor distribution with ν 1 and ν 2 degrees of freedom at p α = 1α and significance level α ∈ (0, 1). The

10.1029/2024GL110347 magnitude of the F-statistic conveys the importance of f d (x d ) in explaining the CVP and is, thus, a measure of the feedback strength. Here, we are mainly interested in f 1 (x 1 ), … , f d (x d ) , as these first-order component functions quantify the direct contribution of x 1 , … , x d to the CVP (shown in next section). Most of the second-order component functions, except for f 23 (x 2 , x 3 ) and f 34 (x 3 , x 4 ) , contribute substantially less to the CVP, with Fstatistics (see Figure S4 in Supporting Information S1) that are almost an order of magnitude smaller than their first-order counterparts (see Figure 2). Despite this apparent insignificance, the second-order component functions still exert control on f 1 (x 1 ), … , f d (x d ) through their shared use of basis functions. The extended bases guarantee hierarchical orthogonality.

We are left with the selection of auxiliary land-surface and atmospheric variables (x 2 , … , x d ) which complement SM, x 1 , in explaining DPR-measured cloud reflectivities, y. We settled on LST, LAI, AT, and TPW and, hence, x = (SM, LST, LAI, AT, TPW) ⊤ is a 5 × 1 vector. LST and LAI modulate evapotranspiration under SM-limited or energy-limited regimes (Seneviratne et al., 2010) and AT and TPW contain synoptic information for the SMCPF about atmospheric preconditioning (Ford, Quiring, et al., 2015;Tuttle & Salvucci, 2017). Figure S1 in Supporting Information S1 presents a correlogram of the five input variables. The decision which input variables to use in regression analysis is always somewhat arbitrary. Our five input variables measure different and complementary properties of the land-surface and overlying air column. One advantage of using ERA-5 reanalysis data of AT and TPW in the HDMR decomposition is that we account explicitly for the dependence of the CVP on the coarse-scale atmospheric conditions (Tuttle & Salvucci, 2017). The F-statistic in Equation 4 will then determine the additional direct contribution of local land surface conditions to the CVP beyond what can be explained by synoptic conditions. Hence, the use of AT and TPW should enhance our confidence in the inferred causal relationships between local SM and the CVP. We do not consider the latent heat flux (LHF) in our HDMR decomposition. The LHF is a result of model simulation and controlled by several of our HDMR inputs. Any assumed constitutive relationship between SM, LST and LAI and the LHF will trouble inference of the individual and cooperative effects of these land surface variables on the CVP. Indeed, LHF inclusion substantially reduces the magnitude of the direct contribution of SM and LST to the CVP (not shown). The control that land-surface and atmospheric variables exert on the CVP varies with cloud height. In the case of SM in panel (a) this equates to a heightdependent SMCPF with a bottom-heavy relationship between SM and CVP, at about 1-3 km above the surface. Above this level, the impact of SM on the CVP decreases rapidly with altitude. As we will show in Section 4.2, the first-order SM component function f 1 (x 1 ) displays a positive feedback due to a wet soil. A higher SM implies a larger evaporative fraction, promoting moderate PBL growth (see Figure S2 in Supporting Information S1) and moisture accumulation (Yin et al., 2015). The CVP at higher altitudes is less dependent on surface SM and controlled more by the upper atmosphere at levels of 3 km and above (Findell & Eltahir, 2003a). Furthermore, a capping inversion layer can inhibit the upward motion of warm, moist air from the surface to the free atmosphere (Findell & Eltahir, 2003b). Indeed, the HDMR-inferred relationship between SM and CVP as articulated by the F-statistic is corroborated by simulation analyses (Findell & Eltahir, 2003a;Koukoula et al., 2019). The strong agreement in the results of the two time lags is a result of SM autocorrelation. The Δt = 7 hour time lag displays the largest influence on the CVP at all measurement heights.

While the impact of SM on the CVP is most pronounced near the surface, LST affects the CVP over a much larger range of altitudes (see Figure 2b). Its F-statistic displays a bimodal relationship with height, peaking close to the surface for Δt = 7 h and at much higher altitudes of about 3.5-4.0 km for Δt = 10 h. As will be further discussed in the next section, f 2 (x 2 ) exhibits a positive correlation with LST, suggesting that positive LST anomalies (dry soils) play an important role in shaping the CVP. The finding that the lower level CVP (1.0-2.5 km) is controlled by LST is in qualitative agreement with the pathway of negative SMCPFs. Positive LST anomalies enhance the sensible heat flux and convective triggering potential (CTP), thereby promoting rapid PBL growth. Figure S2 in Supporting Information S1 testifies to this conjecture and displays PBL heights from ERA-5 reanalysis data. The response of the PBL height to wet and dry surfaces is in agreement with results from simulation-based and observational studies (Findell & Eltahir, 2003a;Ford et al., 2023;Xu et al., 2021), which suggest two mechanisms Geophysical Research Letters 10.1029/2024GL110347 for convective initiation: significant moistening of the PBL (over wet soil) and rapid growth of the PBL (over dry soil). Note that the LST contribution to the CVP decreases at 3.0 km height to increase again between 3.5 and 4.0 km. The reason why LST is influential at higher altitudes may be twofold. On the one hand, LST anomalies favor a strong CTP, increasing the capability of air parcels to overcome convective inhibition and reach the level of free convection (Taylor et al., 2012). If f 2 (x 2 ) measures the contribution of near-surface air to cloud reflectivity for a specific height and time lag, then its F-statistic (in Figure 2b) somehow approximates the dynamics of convective updrafts such that the largest F-statistic value shifts from Δt = 7 hours to Δt = 10 hours with height changing from 1.0 to 5.0 km. Local LST may also influence higher-altitude CVP by elevating the altitude of the melting layer, which is characterized by a sharp increase in radar reflectivity and resides between 3.0 and 5.0 km above the surface during the pre-monsoon and monsoon seasons in the central United States (Song et al., 2021).

The F-statistic of the LAI component, f 3 (x 3 ) , has a negligible impact on the CVP in the lower atmosphere, due to the governing effects of SM, LST, and AT (see Figure 2d) on initiating convection and the subsequent formation of clouds and precipitation. LAI has a somehow stronger impact on the higher-level CVP. The f 3 (x 3 ) component function in Figure S3a in Supporting Information S1 suggests that dense vegetation tends to suppress cloud formation in the upper atmosphere at 3-5 km height. A possible explanation is the role of LAI in modulating surface energy partitioning (Lauwaet et al., 2009), with denser vegetation potentially reducing surface heating and constraining the development of deep convection. Text S3 in Supporting Information S1 presents an analysis of these atmospheric controls on the CVP and physical underpinning comparable to that of the land-surface variables.

The second-order component functions have only a small contribution to the CVP with exceptions of the LST-LAI component, f 23 (x 2 , x 3 ) , and the LAI-AT component, f 34 (x 3 , x 4 ) . The joint contributions of LST-LAI and LAI-AT to the CVP at 1-3 km height may be linked to processes such as plant transpiration and boundary layer moistening. This is further illustrated in Figure S4 in Supporting Information S1 and discussed in Text S4 in Supporting Information S1. The comparison between HDMR and traditional data-driven methods (see Figures S5-S6 and Text S5 in Supporting Information S1) reveals that all three approaches offer valuable insights into the diurnal SMCPF. However, HDMR is able to discern time lags in this feedback and captures the evolving physical dynamics of thermal updrafts.

In this section, we focus on the spatial pattern of the SMCPF within the study region. We reiterate that we conduct functional decomposition of the cloud reflectivity using all the samples of April-October (2016-2019) for a specific time lag and cloud height, to guarantee an adequate number of samples and storm events. Our goal here is to present a 4-year averaged spatial distribution of the derived component functions and identify locations of positive and negative SMCPF rather than focusing on interannual and/or cross-season variations.

Figures 3a and 3b present the spatial distribution of the antecedent 7-hr SMAP/L4 soil wetness at the top layer (0-5 cm), collocated at the coordinates of the GPM/DPR/L2A samples of April-October (2016-2019), alongside the corresponding first-order component function, f 1 (x 1 ) (dBZ), evaluated at 2.0 km. The GPM/DPR/L2A samples are spatio-temporally scattered, directly corresponding to the satellite's path and observation times. This results in voided areas on the map. Our examination of SM's feedback strength, conditioned on an altitude of 2.0 km and a 7-hr time lag, is of particular interest upon our prior analysis of the F-statistic in Figure 2a. Panels (b-c) reveal the positive feedback from SM represented by f 1 (x 1 ) . With a degree of saturation exceeding 0.4, wet soil could increase cloud reflectivity by up to 4 dBZ. The fact that the absolute value of f 1 (x 1 ) decreases with height in Panel (c) again lends support to our inferred height-dependent SMCPF in Section 4.1, underscoring the stronger coupling between SM and CVP in the low-level atmosphere. As a byproduct, we demonstrate in Text S6 and Figure S7 in Supporting Information S1 the application of the Marshall and Palmer (1948) formula to transform f 1 (x 1 ) (dBZ) into estimates of rainfall rate.

Significant positive feedback of SM is evident in regions such as northern Texas, central Oklahoma, northwestern and southeastern Kansas, and northeastern New Mexico. All these areas, with the exception of northeastern New Mexico, are located inside or close to the 'transitional regions' as delineated by the dashed gray lines as categorized by Findell and Eltahir (2003b). The middle transitional region, spanning from the semi-arid southwestern to the humid southeastern parts of the central United States, is influenced by both the dry and wet soil advantage regimes (Findell & Eltahir, 2003b). Hence, this dual influence explicates the observable positive feedback in the central and eastern sections of the transitional region and negative feedback in the southwestern part (detailed below). These local wet soil anomalies can be attributed to early warm-season mesoscale convective systems (MCSs) (Hu et al., 2021).

Since SM can indirectly exert feedback on cloud and precipitation through heating or cooling the surface (Duerinck et al., 2016), we further delve into examining spatially the samples of antecedent 7-hr LST (K) and their contribution to cloud, f 2 (x 2 ) (dBZ), and rainfall, ΔR (mm/hour), in Figures 3d-3f and Figure S7 in Supporting Information S1, respectively. f 2 (x 2 ) exhibits a non-linear dependence on LST where LST anomalies exert the most significant influence. From Figures 3d and 3e, it is suggested that LST above 305 K accounts for an increase of up to 4.0 dBZ in the cloud reflectivity and 2.0 mm/hr in rainfall rate (see Figure S7 in Supporting Information S1) at both 2.0 and 3.5 km. On the contrary, the samples with a cooler surface (LST<290 K) seem to foster a more stable atmospheric state, thereby reducing cloud reflectivity, especially in the lower troposphere (h ≈ 2.0 km).

Geographically, the most significant effects of these anomalies are evident and clustered in the southwest of the study region (101°W-105°W and 32°N-36°N). Within this area, we find a significant negative correlation (R = -0.41, p < 0.0001, shown in Figure S8a in Supporting Information S1) between surface SM and the LST component function, f 2 (x 2 ) . Moreover, we illustrate in Figure S8b in Supporting Information S1 that LST Findell and Eltahir (2003b). Panel (c) displays f 1 (x 1 ) (dBZ) as a function of antecedent 7-hr SM. Values of f 1 (x 1 ) are evaluated at three separate heights: 2.0 km (red), 3.5 km (yellow), and 5.0 km (blue). The light-colored regions correspond to the 95% bootstrap confidence intervals. The bottom row of panels presents the same content as panels (a-c) but for panel (d) SMAP/L4 LST and (e, f) its associated component function, f 2 (x 2 ).

10.1029/2024GL110347 contributes to CVP preferentially over dry soil with saturation between 0.1 and 0.4. These findings underscore the presence of the intrinsic SM-LST coupling within SMCPF pathways (Seneviratne et al., 2010), and we can conveniently interpret f 2 (x 2 ) as a proxy for the indirect and negative feedback of SM on CVP. Notably, our identified negative feedback region (101°W-105°W, 32°N-36°N) is consistent with the one proposed by Findell and Eltahir (2003b) (represented by the black dashed line in Figures 3d and 3e). Several factors can play a role when it comes to the sources of convective clouds and precipitation over the dry soil. For instance, the monsoonal moisture incursion into New Mexico can increase local humidity and offset the reduced evapotranspiration from the local dry soils (Klein & Taylor, 2020;Wallace et al., 1999). Besides, the Great Plains Low-Level Jet (GPLLJ) can transport abundant moisture southerly from the Gulf of Mexico into the central United States (Feng et al., 2016;Ford, Rapp, & Quiring, 2015).

In this data-driven study, we used functional decomposition of a large database of satellite-measured soil moisture (SMAP/L4) and cloud vertical profiles (GPM/DPR/L2A) to quantify the relationship between topsoil moisture content (SM), land-surface temperature (LST), leaf-area index (LAI), atmospheric temperature (AT), total precipitable water (TPW) and cloud reflectivity in the central United States. Results show that the sign and strength of this soil moisture-cloud-precipitation feedback (SMCPF) differs substantially between cloud heights and geographical locations. A significant positive feedback is observed in the lower atmosphere, particularly between 1.0 and 3.0 km and SM time lag of 7 hr. Wet soils (saturation degree of 0.4 and larger) can increase the cloud reflectivity (rainfall rate) by up to 4.0 dBZ (2.0 mm/hr) at a height of about 2.0 km. This is most evident in northern Texas, central Oklahoma, and northwestern and southeastern Kansas. Negative SMCPF feedbacks, are a result of LST anomalies and extend from 1.0 to 4.0 km at a time lag of 7-10 hr. These LST anomalies display comparable increments in cloud reflectivity and rainfall rates to SM, except for northwestern Texas and southeastern and eastern New Mexico. The patterns of SMCPF identified by our satellite data decomposition method are in qualitative agreement with the findings of previous studies (Findell & Eltahir, 2003a, 2003b;Ford et al., 2023;Hu et al., 2021;Koukoula et al., 2019;Qian et al., 2013;Sathyanadh et al., 2017;Su & Dickinson, 2017).

We could have resorted to multivariate linear regression and written the DPR-measured cloud reflectivities as a superposition of linear multiples of the input variables. This approach has elements in common with the multivariable function expansion of Soboľ. (1993) but by pooling together structural, cooperative and/or interaction contributions, linear regression may not accurately portray the direct effects of the input variables on the CVP. This should not discourage anyone from using traditional regression methods in hydroclimatological and hydrometeorological research. For example, G. Wang et al. (2024) uses linear regression to single out the influence of lower-tropospheric moisture on the local SMCPF.

The ANOVA-inspired functional decomposition of HDMR accurately disentangles the individual and correlative contributions of single and groups of correlated input variables to DPR-measured cloud reflectivities. HDMR is a powerful addition to the toolbox of regression methods used by hydrometeorologists and hydroclimatologists for analyzing the complex and intricate relationships between land-surface and atmospheric variables. HDMR results can help diagnose biases in the current generation of weather and climate models and help detect and quantify changes in SMCPFs at regional and global scales as a result of climate change and hydroclimatic extremes.

HDMR functional decomposition, however, requires a relatively large number of samples at each cloud height. This complicates seasonal, interannual, or localized analysis of the SMCPF. Another limitation of HDMR is that the number of input variables should not exceed, say, d = 10, otherwise the matrices involved may become too large-sized complicating the estimation of the component functions' expansion coefficients with D-MORPH regression. Our current set of d = 5 input variables serves as an illustration of the HDMR method but may ignore other variables that control the magnitude and sign of SMCPFs. collocated SMAP and DPR samples, have been archived on Zenodo (Gao et al., 2024) with Creative Commons Attribution 4.0 International license. Precipitation Measurement Mission program (Grant 80NSSC22K0597) and the Weather and Atmospheric Dynamics program (Grant 80NSSC23K1304), as well as the support of the NSF Division of Information and Intelligent Systems (ExpandAI2ES Grant IIS2324008). Yifu Gao wishes to thank the support of the UCI Engineering-Los Alamos National Laboratory (LANL) Fellowship, provided by the UCI Samueli School of Engineering and LANL. Chonggang Xu acknowledges the support through RUBISCO Science Focus Area (SFA) by DOE Office of Science, Biological and Environmental Research, Regional & Global Model Analysis program. The constructive comments from the two anonymous reviewers are highly appreciated and have substantially improved the manuscript.