1. Introduction
Plant growth is tightly coupled with soil nutrients, particularly nitrogen (N) and phosphorus (P) [1,2,3,4], and limitation of these nutrients is assumed to be one of the drivers of observed spatial variability in ecosystem carbon (C) fluxes [5,6,7], particularly terrestrial ecosystem productivity [8,9,10,11]. Soil P is also a limiting factor of global terrestrial ecosystem productivity [12,13], and studies have shown that P limitation on plant growth occurs more frequently than previously thought [3]. Because soil nutrients modify plant growth [14], a change in the soil-nutrient composition can alter plant responses to elevated atmospheric carbon dioxide (CO2) concentration and climate change [15,16,17,18,19,20].
Physical and chemical soil properties, such as texture, soil moisture content, and pH can influence soil nutrient concentration [21,22]. For instance, high soil pH can reduce the availability of soil N and P resulting in nutrient deficiencies [23,24]. Soil age is also a key factor determining nutrient content and the type of nutrients that limit vegetation productivity [25,26]. For example, older and more developed soils can be leached of important rock-derived nutrients such as P [27]. Conversely, soil N may be derived from external processes, such as atmospheric N disposition and biological N fixation. The concentrations of soil nutrients also vary with soil depth [28]. For instance, microbial activity can move nutrients upward into the root zone [29,30], while translocation and infiltration may increase the nutrient concentration with depth by moving nutrients downward. Thus, root systems and plant uptake along with hydrology and microbial communities have strong effects on the vertical distribution of soil nutrients [28,31].
Previous studies showed strong relationships between soil nutrient and vegetation carbon use efficiency (CUE) and water use efficiency (WUE) for different ecosystems [32,33,34]. However, the combined effects of multiple nutrients and soil characteristics on CUE and WUE have been rarely examined. CUE, defined as the ratio of net primary production (NPP) to gross primary production (GPP), measures the ecosystem proficiency to store absorbed atmospheric carbon (C) [9]. CUE plays an important role in regulating the terrestrial ecosystem’s carbon balance and determines the amount of carbon allocated to biomass [35]. An additional metric of plant resource economy quantifies the rate of carbon uptake per unit of water use. WUE measures the effect of water on the ecosystem to absorb atmospheric C [36] and varies with soil nutrients content [37]. WUE, a key physiological parameter linking carbon and water cycles, quantifies the amount of water that terrestrial ecosystems use relative to carbon gained [38,39]. CUE and WUE are important indicators of plants’ ability to adapt to changing environmental conditions, such as atmospheric CO2, precipitation, and temperature.
Satellite-based remote sensing of vegetation can be used to estimate global WUE and CUE from NASA’s Moderate Resolution Imaging Spectroradiometer (MODIS) GPP, evapotranspiration (ET), and NPP products starting from year 2000. The MODIS NPP and GPP algorithms were also modified to generate GPP and NPP at 30 m resolution from Landsat and meteorological data for the conterminous United States [40]. Estimates of global GPP are also available from the Max Planck Institute for Biogeochemistry (MPI-BGC) based on upscaling of eddy covariance-based measurements [41]. Currently, NASA ECOSTRESS onboard the space station is producing ET and WUE starting from July 2018 [42] and can also lead to instantaneous GPP estimates [43]. In addition, satellite-based solar-induced fluorescence (SIF) has been related to GPP [44,45,46]. All these satellite-derived data can be used to quantify WUE and CUE. For example, MODIS data have been applied to quantify the spatial variability and drivers of global WUE and CUE [39,47,48,49] and drought impacts on vegetation WUE [50,51]. Nevertheless, the linkage between satellite-derived CUE and WUE, and soil properties have been rarely examined.
This work builds on previous studies of plant-nutrient limitations by exploring multiple nutrient effects, variable soil profiles, and interactions between these factors. The larger effects of multiple nutrient additions (N, P, calcium, or Ca) on forest productivity than a single element was, revealed in a recent study of nutrient availability in northeastern US hardwood forests [27]. Additionally, a strong relationship between satellite-derived canopy greenness and soil N, aluminum (Al), and P was observed in tropical forests [52]. Moreover, many modeling experiments considered soil as a mixed column and thus relationships between soil nutrients and vegetation growth cannot be inferred for soil depth. Here, we examine how relationships of soil variables with CUE, WUE, and greenness indices vary with soil depth and with satellite sensors’ spatial resolutions. We present a novel approach combining CUE and WUE estimates from satellite GPP, NPP, ET, and SIF with soil properties to understand the regulatory mechanisms and effects of the nutrient cycles. Previous studies have shown that WUE varies with land use type, with forest having the highest WUE and croplands having the lowest WUE [39,53,54]. The relationship between CUE and land use type is more complex and is driven by climate controls, such as temperature and precipitation; CUE is higher in a drier and cooler climate than in a wet and warm climate [55]. In general, forest has lower CUE than grasslands, and crops have medium CUE values [55,56]. Based on previous studies, we hypothesize that the CUE and WUE depend on land use and depth-related soil properties. Moreover, this work elucidates how ecosystems are structured and function based on connections between aboveground vegetation productivity and soil properties.
2. Materials and Methods
2.1. Site Selection
We collected soil samples from 19 sites spatially distributed across a warm temperate and fully humid climate [57] in the Commonwealth of Kentucky, USA. The region is characterized by a humid subtropical climate, with cool to moderately cold winters and warm to hot summers. The mean annual temperature ranges from 11.7 °C in the northeast to 15 °C in the southwest. The warmest month is August, while the coldest month is January. Mean annual precipitation is ~1300 mm in the south and 965 mm in the northern areas, with August the driest month and May the wettest month.
The site selection criteria consisted of two stages. Stage 1 involved the selection of potential sites that represented the major land cover types of Kentucky (i.e., agricultural, grassland, and forest). Stage 2 assessed differences in lithologic and topographic effect (slope) to assure analyses of vegetation–soil interaction were not confounded by parent material and relief. Following our site selection criteria, we identified 19 sites for soil sampling (Table 1). We limited our site selection to two soil orders: Alfisol and Ultisol, which minimized differences due to soil physical, chemical, and mineralogical properties (Figure 1). These two soil orders at the 19 sites consist of fine-grained parent material that is slight to strongly leached of nutrients in the topsoil. These soils also have clay-enriched B horizons. Unlike surrounding Inceptisols and Entisols, these two soil orders more likely reflect the surrounding long-term climate and ecological conditions. These two soil orders cover approximately 18% of the global glacial free land area according to the Soil Science Society of America (
2.2. Soil Physical and Chemical Analysis
Continuous soil cores were extracted using a hammer-percussion coring device and the soil horizons were described [58]. The cores were cut in half in the laboratory and an interval-sampling approach was used, where soil samples were extracted from the core every 8 ± 2 cm, to ensure that each horizon was sampled at least once. The soil physical and chemical properties were measured on samples collected to a depth of 60 cm. Particle size, oven-dry bulk density, gravimetric water content, pH, and loss on ignition (LOI) were measured for the 19 sites at Murray State University using previously established methods [59] with modifications [60]. LOI describes the variance in soil organic matter content (SOM) and is herein referred to as OM%. Soil samples were sent to the University of Kentucky Soils Characterization Lab and extractable P, C, Ca, Mg, Zn, and Fe were measured using the Mehlich-3 extraction method [61,62].
We performed soil stable isotopes analysis because they are indicators of soil C and N dynamics in an ecosystem and can be used to determine the sources of SOM and soil C turnover rates [63,64]. Soil samples were sent to the Stable Isotope Ecosystem Lab of the University of California, Merced (SIELO) for isotope analysis and elemental content of organic C and N. Splits of the soil sample (<2 mm size fraction) were prepared for stable isotope and elemental composition analysis with an initial homogenization step where samples were powdered using mortar and pestle for three to five minutes. One portion of this powdered subsample was placed in centrifuge vials and treated with 1N HCl to remove any inorganic C (calcite), which confounds the organic C signal for isotopic and elemental compositional analysis. After acidification, these subsamples were triple-rinsed in deionized water, then dried overnight at 50 °C. The other split was weighed without acid treatment to preserve δ15N values. All powdered samples were weighed to 20–40 mg (based on organic C and N content) into 5 × 9 tin capsules for stable isotope (13C/12C and 15N/14N) and organic elemental (C and N) analysis. The SIELO analyses content and isotopes of C and N with an Elemental Analyzer (Costech 4010) coupled to a Thermo Scientific Delta V Plus isotope ratio mass spectrometer via conflo IV. Reference materials for elemental and isotopic composition were run with all samples and included at least n = 3 of each Costech acetanilide, Peach Leaf (NIST 1547), and USGS 40. Reproducibility for these reference materials was <0.3‰ for both δ13C and δ15N values as well as <3% for C and N content (Table 2). Stable isotope compositions are reported with VPDB and AIR as references for δ13C and δ15N values, respectively, and in per mille (‰). All C:N ratios are reported as weight percent.
2.3. Satellite Data
The MODIS 250 m resolution EVI (MOD13Q1) and NDVI (MOD13Q1) data and 500 m resolution ET data (MOD16A2; hereafter ET) were obtained from the USGS Earth Explorer website (
We also used the GPP product (GOSIF GPP) based on SIF derived from the Orbiting Carbon Observatory-2 (OCO-2). GOSIF is a global, OCO-2 based SIF product, while GOSIF GPP is based on GOSIF and linear relationships between GPP and SIF [65,66]. Both products are available at
CUE based on each of the data inputs was calculated as follows:
(1)
(2)
(3)
where CUE30 is CUE estimated from Landsat 30 m GPP and NPP, CUE250 is CUE derived from MODIS 250 m GPP and NPP, and CUESIF is CUE based on MODIS 250 m NPP and GOSIF GPP. Similarly, WUE for each of the data sets was calculated as:(4)
(5)
We used GOSIF GPP along with MODIS NPP to calculate CUE and WUE to explore how SIF-based GPP can be used to examine the relationships of WUE and CUE with soil properties.
2.4. Statistical Analysis
To identify and rank the soil chemical and physical properties that are significant in predicting CUE and WUE, we used the classification and regression trees (CART) method (the rpart package in R). The advantage of the regression tree is that it uses recursive partitioning to split the data into nodes or branches. This allows the training data to be split into groups that behave independently without the need for data stratification before analysis. For soil data analysis, whereas soil variables are the results of interdependent soil processes, the CART model is very helpful because it reveals variable interactions [67,68]. Thus, CART can help address the interdependence between soil data with depth and CUE and WUE. The predictors used to estimate WUE and CUE were land cover type, soil order, soil horizon, soil depth, layer thickness, bulk density, silt, sand, clay, OM, EC, pH-CaCl2, pH-H2O (1:1 H2O), soil organic carbon (SOC), total C, total N, C:N ratio, δ13C, δ15N, and extractable soil P, K, Ca, Mg, Zn, and Fe. Prediction errors were estimated by randomly dividing the data to 70% training (n = 80 for WUE and WUESIF; n = 72 for CUE30, n = 63 for CUE250, and n = 76 for CUESIF) and 30% testing samples (n = 35 for WUE and WUESIF; n = 31 for CUE30, n = 28 for CUE250, and n = 33 for CUESIF). The differences in the sampling sizes for the CUE data were related to the availability of satellite data for our study sites. We built the regression tree using the training sample and then pruned the tree to minimize overfitting by selecting the number of splits in the CART that is associated with a minimized cross-validated error estimate. The CART results were used to identify the soil chemical and physical variables that were important predictors for CUE and WUE. Then, the mean absolute error (MAE) was calculated by using the testing sample with the pruned CART model. We applied the CART for the 30 m and 250 m CUE and WUE to investigate the impact of spatial resolution on the results of CART. Additionally, CART was developed for CUESIF and WUESIF. Finally, we analyzed the CART results based on GPP derived from MODIS and SIF data to detect any difference in soil chemical and physical variables identified as important predictors for CUE and WUE. The reason for such comparison is that GOSIF GPP is directly estimated from SIF that is a direct measure of energy emitted by plants during photosynthesis, while MODIS GPP and NPP are based on the light use efficiency model and process-based model, respectively. Thus, they represent two different approaches to estimate CUE and WUE. To test for the effects of land cover type, soil order, soil horizon, and soil depth, we performed four additional CART analyses.: (1) excluding only land cover from the CART model (CARTnLC); (2) excluding only soil order from the CART model (CARTnSO); (3) excluding only soil horizon from the CART model (CARTnSH); and (4) excluding only soil depth from the CART model (CARTnSD). This was performed for each of the WUE and CUE data and the model with the highest r2 and lowest MAE was selected as the best model.
Soil physical and chemical properties were compared to satellite CUE and WUE. To achieve this, we summed the 8-day MODIS GPP and ET and 16-day Landsat GPP to produce annual GPP for the study sites. Both MODIS and Landsat-based NPP are directly available as annual products. Then, we used the annual data to calculate WUE and CUE at the annual timescale.
3. Results
3.1. Soil Data
The study sites sampled for this study (7 grasslands, 6 forestlands, and 6 agricultural sites) are dominated by silty clay and silty clay loam soils (Figure 2). Forest sites showed the highest soil SOC and the lowest δ15N values in the top 10 cm of the soil, while agricultural soils had the highest P concentrations, which could be due to fertilization (Figure 3). The low forest δ15N in the top 10 cm of the soil can be related to litter inputs and low microbial processing [69] as evident in higher C:N in the soil A horizon (Table 2). In general, soil physical and chemical properties were not significantly different based on Tukey’s HSD test for land use types (i.e., forest vs. grasslands vs. agriculture) (see Supplementary Data). We note that one of our sites was formerly agricultural land (currently converted to grassland) with soil extractable P values higher than other grassland sites that were not cultivated (Figure 3C).
3.2. Soil Variables Versus CUE and WUE
The pruned CART model for MODIS WUE showed that soil δ13C was the primary splitting node indicating the role of SOC in regulating WUE since soil δ13C is a facet of soil carbon content. MODIS WUE decreased with greater soil δ13C values (less negative) and increased with lower soil δ13C values (Figure 4A). The second split was soil order Alfisol that decreased MODIS WUE (Figure 4A), indicating lower WUE in Alfisol soils. Soil OM was the only split for WUESIF, showing a negative relationship between soil OM and WUESIF (Figure 5A). Table 3 shows the r2 and MAE values for each of the CART models developed to address the impacts of land cover, soil order, soil horizon, and soil depth. The model with the best model was the one with land cover excluded (CARTnLC) from the analysis for MODIS WUE with r2 of 0.39 and MAE of 0.17 gC/KgH2O (Table 3). The WUESIF models showed the highest MAE when all the soil variables were included (Table 3), and thus we selected the model with land cover excluded (CARTnLC) (r2 = 0.19, MAE = 0.64 gC/KgH2O) (Table 3). We found that MODIS WUE outperformed WUESIF likely because MODIS had finer spatial resolution (250 m) than the GOSIF GPP (0.05°).
The pruned CART model for CUE30 showed that land cover was the first split, grassland being the main land cover type in the first split (Figure 6A). The second split was soil extractable Ca that showed a negative relationship with CUE30, with CUE30 increasing at lower Ca values (Figure 6A), and the CART model had r2 of 0.64 and MAE of 0.04 (Table 3). Land cover was the first split and the most important predictor for CUE250, but with the forest as the land cover predictor rather than a grassland as in CUE30 (Figure 7A). CUE250 decreased with forest with higher soil extractable Ca (Figure 7A), and the model had r2 of 0.78 and MAE of 0.06 (Table 3). Since the difference between the original CART and CARTnSH for CUE250 was almost identical, we decided to use the model that included all the soil variables (Table 3). The pruned CART model for CUESIF showed that soil δ13C was the first split showing a positive relationship with CUESIF, that the second split was soil silt content with lower CUESIF at lower soil silt content, and that the third split was soil extractable Mg that a positive relation with CUESIF (Figure 8A); the model had r2 of 0.45 and MAE of 0.1 (Table 3). The variable importance plots for the CART models are shown in the Supplementary Figures S1–S5. As expected, soil order and soil δ13C were the most important predictors of WUE (Figure S1), while OM was the most powerful predictor of WUESIF (Figure S2). Land cover was the most important predictor of CUE30 and CUE250 (Figures S3 and S4). Soil δ13C was the most important predictor of CUESIF (Figure S5). It is important to note that all other soil variables had some amount of predictive capacity for both WUE and CUE (Figures S1–S5).
3.3. Impacts of Different Productivity Measures and Spatial Resolution on the Relationship Satellite-Derived CUE and, WUE Withsoil Properties
We applied the CART for MODIS WUE and WUESIF to investigate the impacts of different productivity measures (based on the MODIS light use efficiency model vs. GPP based on OCO-2 SIF) on the CART results. There were no common variables detected between CART-MODIS and CART-SIF (Figure 4 and Figure 5). We compared GEOSIF GPP and MODIS GPP and found a low correlation (r2 = 0.19) between these two products (Figure S7). In general, GOSIF GPP values were much higher than MODIS GPP which was reported to underestimate site measured GPP, particularly for cropland [40,70]. Besides, we applied CART for Landsat CUE, MODIS CUE, and SIF-based CUE to investigate the impact of spatial resolution on the results of CART. The land cover type and soil extractable Ca were the common variables between MODIS CUE and Landsat CUE (Figure 7A and Figure 8A) highlighting the importance of these soil variables for predicting CUE. No common soil variables were detected among CUESIF, CUE30, and CUE250 (Figure 6, Figure 7 and Figure 8). This could be due to different algorithms used to estimate MODIS and Landsat GPP compared to the GPP derived from SIF and the different spatial resolutions of these data products (Table 4).
4. Discussion
Overall, our approach is robust and applicable to our study sites given the significant relationships of soil chemical and physical properties with CUE and WUE. The recursive partitioning that CART is built upon can model WUE and CUE from soil physical and chemical properties by providing a simplified, but yet powerful perspective on the drivers of WUE and CUE. The CART models showed higher predictability for CUE than for WUE likely due to the uncertainties in MODIS ET or other environmental controls on WUE (e.g., vapor pressure deficit). Our results did not reject our hypothesis that CUE depends on the land cover type with grasslands and forest being the most important predictors of CUE30 and CUE250, respectively. Conversely, our results rejected the hypothesis that WUE is dependent on land cover type.
The CART WUE and CUE results revealed indirectly the correlations of soil δ13C and OM with soil depth. In general, OM was higher and soil δ13C values were lower in the topsoil layers and OM becomes lower and δ13C values increase with an increase in soil depth [71], which was also the case for our study sites (Figure S6). Although CART WUE and CART CUESIF models did not show that soil depth was an important predictor of WUE and CUESIF, they did reveal an indirect indication of the soil depth-related properties in regulating WUE and CUE as evident in Figure 4, Figure 5 and Figure 8. Therefore, our results failed to reject our hypothesis that WUE and CUE are independent of soil depth-related properties.
Our analysis revealed the importance of soil extractable Ca and Mg in regulating CUE and WUE and calls for more data collections and synthesis to better understand the potential role that these soil properties play in moderating the effects of climate change on CUE and WUE. Our methodology demonstrated the ability of remotely sensed data to capture the observed relationships of CUE, WUE, NDVI, and EVI with soil and provides a platform for future exploration of the soil–vegetation interactions using satellite vegetation data and in situ soil data.
The regression tree trained using MODIS data relied on soil δ13C and soil order to predict WUE (Figure 4). Soil δ13C values showed a positive correlation with WUE (Figure S4A) similar to the observed positive correlations between Leaf δ13C and intrinsic WUE [72,73]. This relationship can be attributed to the close relationship between SOC and δ13C values [71,74] and the role that SOC plays in controlling soil water availability. Regression tree analysis showed that soil order Alfisol was important for predicting WUE. This could be driven by our three sites with soil order Ultisol that on average had higher WUE than the other sites. In general, the three forest sites with soil order Ultisol had lower soil extractable Ca and K than the other sites. An increase in the uptake of soil extractable Ca and K can enhance WUE by decreasing ET [37]. We hypothesized that higher WUE in these three sites is due to an increase in the uptake of soil extractable Ca and K. We plan to measure leaf nutrient concentration data from our study sites to test this hypothesis in future studies. Soil OM was the only statistically significant predictor of WUESIF (Figure 5). Additionally, OM showed a negative correlation with WUE and can indicate structural water in clay minerals and/or SOC and can affect soil water holding capacity. This is because an increase in soil OM is critical for forming soil aggregate that can increase the soil-plant available water and the latter controls the ET rates [75]. These findings indicate that an increase in OM could lead to a decrease in WUE.
The regression trees trained using Landsat and MODIS data were very similar, with both relying on land cover and soil extractable Ca to predict CUE. Land cover was the most important predictor for CUE highlighting the difference between forest and grassland ecosystems, which was similar to the observed pattern in MODIS-derived CUE [33,55]. The importance of soil extractable Ca for predicting CUE highlighted the close interactions between plant N uptake, soil N pools, and soil extractable Ca [76,77], and an increase in soil extractable Ca increases GPP and eventually decreased CUE. The CUE increased with a decrease in soil extractable Mg because of the role that Ca plays in regulating physiological responses related to plant stress and growth [37,78]. For instance, the deficiency in soil Ca was related to an increase in root biomass and a decrease in respiration [78], which can lead to an increase in NPP and potentially CUE. The importance of the soil δ13C as a predictor of the CUESIF has not been reported before mainly because of limited site CUE data. As more soil δ13C data becomes available, more in-depth studies can potentially further probe into the linkages between soil δ13C values and vegetation WUE and CUE. The CART CUESIF model showed a positive relationship between CUESIF and soil extractable Mg. Soil Mg is involved in various processes that affect photosynthesis and plant growth, such as chlorophyll a and b production and root and shoot biomass [37,79,80]. A recent meta-analysis study found that an increase in soil extractable Mg leads to an increase in NPP and biomass [81], which can eventually increase CUE as predicted by our CART results.
The regression tree analysis revealed relationships of soil chemical and physical properties with WUE and CUE. The CART analysis allowed us to: unravel the interactions among various soil properties and their and combined effects on WUE and CUE; and to identify soil variables that are most important for predicting WUE or CUE. Our approach provided an understanding of the interrelationships within the datasets (e.g., soil depth) more easily than using, for example, simple or multiple regression analysis. It should be noted that the results shown in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 did not mean that these variables were the only soil variables important for predicting WUE and CUE, but rather that the CART method found these variables to be able to minimize the variability in WUE and CUE between all the tree nodes. Figures S1–S5 show the ranking of the importance of soil variables in predicting WUE and CUE. Thus, all the other soil variables (Figures S1–S5) had certain predictive power for WUE and CUE suggesting the need for more comprehensive soil–vegetation analysis to better untangle the complex soil–vegetation processes that underlie WUE and CUE.
All soil chemical and physical variables included in this study have direct or indirect effects on vegetation processes, including plant growth, WUE, GPP, and ET [37,79,80,82]. Soil extractable Zn and C:N ratio were the only soil variables that did not show any predictive power for either WUE or CUE. (Figures S1–S5). Soil δ15N values appeared once as a variable that had predictive power for WUE (Figure S1). Thus, for our study sites, these soil variables were not important in minimizing the errors in CUE or WUE between the CART tree nodes.
We collected soil samples from two soil orders to minimize the effects of soil order on soil chemical properties. We expect that our soil chemical and physical data to have little variability as they were collected from the center of Landsat 30 m pixels with a homogenous land cover type. This is less likely to be the case for the coarse-resolution MODIS (250 m) and GOSIF GPP (0.05°) data, as individual grid cells could contain multiple land cover types. With time, vegetation type and land use could alter soil properties and result in different soil chemical concentrations, leading to heterogeneous soil properties in a pixel with multiple land cover types. Moreover, soil chemical properties vary significantly as a function of topographic position. Although we controlled for land cover type (for Landsat data), we did not control for topographic variability. In addition, the scale mismatch between Landsat and MODIS can result in differences in the land cover type assigned for calculating GPP and NPP. Thus, MODIS GPP and NPP values might not represent the same land cover type as Landsat derived GPP and NPP. The scale mismatch among site, MODIS, and GOSIF GPP could influence the results since no common soil variables were detected between MODIS and GOSIF GPP CART results for CUE and WUE. Moreover, the difference in algorithms MODIS and GOSIF GPP products (Table 4) could have contributed to the differences in our results as evident in the low correlation between MODIS GPP and GOSIF GPP.
Finally, in implementing our approach, we limited the soil sampling to two soil orders to control for the effect of the variability in soil type on the results. Limiting our sampling to two soil orders could affect the broad application of our results. We anticipated that the observed relationships of CUE and WUE with soil chemical properties could apply to these two soil orders in different climatic regions, but more research is needed to verify the correlations that we observed. The robustness of our results is dependent on the accuracy of satellite data, used to estimate CUE and WUE. Newer, higher spatial resolution data on WUE from ECOSTRESS have the potential to significantly refine these results further. ECOSTRESS produces an operational WUE product at 70 × 70 m spatial resolution every 1–5 days from the International Space Station [42,83]. ECOSTRESS measures thermal infrared radiance, which is converted into a land surface temperature and emissivity (L2_LSTE) product; L2_LSTE is then converted into ET (L3_ET_PT-JPL) using the Priestley–Taylor Jet Propulsion Laboratory (PT-JPL) method [84]. ET is then combined with GPP from MOD17 to create the WUE product (L4_WUE); newer versions of ECOSTRESS replace MOD17 with a native 70 m GPP based on the Breathing Earth System Simulator [85]. ECOSTRESS can detect spatial variability in WUE never seen before (Figure S8). These differences are evident throughout the study region landscape, especially in relation to topography, hydrology, and land cover, and land use. Nevertheless, to our knowledge, our effort is one of the few that relate satellite data to soil properties for investigating the soil–vegetation interactions. CUE and WUE terms are important indicators of plants’ adaptability to changing environmental conditions, such as precipitation and temperature. Understanding changes in CUE and WUE is thus critical to quantify the response of the terrestrial ecosystem to climate change.
5. Conclusions
CUE and WUE are key physiological parameters linking carbon and water cycles. Understanding changes in these efficiencies is critical to quantify future terrestrial ecosystem responses to climate change. Here, we presented a different approach to understand the stand-level soil–vegetation interactions by using remotely sensed data. Data from the 19 sites representative of the major land cover types occurring in a warm-temperate and humid climate were used to investigate the capability of satellite data to derive meaningful information about the soil–vegetation interactions. The results revealed significant relationships between satellite-derived CUE and WUE and measured soil chemical properties. The CART model identified soil variables that were important to predict CUE and WUE, implying that satellite data have the potential for elucidating soil–vegetation interactions.
We have learned from this study that satellite spatial resolution was not important for more accurate detection of relationships between satellite vegetation products and soil chemical properties. Interestingly, coarser-resolution MODIS data produce the best estimates for CUE than high-resolution remote sensing products (e.g., Landsat CUE). Our study is limited by the number of sites used and the reported relationships may differ with more data or additional biomes or soil types. Nevertheless, our study represents a new approach that can be used to better understand the observed spatial variability in satellite-derived CUE and WUE data and can be used to inform land surface models about missing processes that are essential for more accurate predictions of vegetation CUE and WUE [86].
Supplementary Materials
The following are available online at
Author Contributions
B.E.M. and G.E.S. conceived and designed the research; B.E.M., G.E.S. and S.L.K. performed the research and analyzed and interpreted the data; B.F. collected the data and performed lab analysis; B.E.M. was responsible for the research analysis; J.B.F. produced the ECOSTRESS data; J.X. contributed the GOSIF GPP data; B.E.M., G.E.S. and S.L.K. contributed to the writing of the manuscript; H.C., J.B.F. and J.X. reviewed and revised the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding
The material is based upon work supported by NASA Kentucky under NASA award No: NNX15AK28A. J.X. was supported by NASA’s Climate Indicators and Data Products for Future National Climate Assessments (award No: NNX16AG61G).
Data Availability Statement
The data presented in this study is available upon request from the corresponding author.
Acknowledgments
We thank Lance Stewart, Mallory Gerzan, Molly Karnes, Austin Johnessee, and Joalda Morancy for the field, lab, and data assistance. J.B.F. contributed to this research at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Government sponsorship acknowledged. J.B.F. was supported in part by NASA programs: IDS and ECOSTRESS. Copyright 2021. All rights reserved.
Conflicts of Interest
The authors declare no conflict of interest.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figures and Tables
Enlarge this image.Figure 1. Location of the 19 study sites and their land cover types in our study region, (Kentucky, USA).
Enlarge this image.Figure 2. Soil texture pyramid for the study sites showing the dominance of silty clay (SIC) and silty clay loam (SICL) soils. Data represents soil texture for every soil depth sampled for each site. C = clay, SC = sandy clay, SCL = sandy clay loam, SL = sandy loam, LS = loamy sand, S = sand, L = loam, CL = clay loam, SIL = silt loam, and SI = silt.
Enlarge this image.Figure 3. Variability of soil chemical properties with soil depth for the study sites: (A) total soil carbon (%); (B) δ15N; (C) extractable phosphorus (P); (D) extractable potassium (K).
Enlarge this image.Figure 4. Pruned regression tree for (A) MODIS-based WUE predicted as a function of two significant variables. The number in boxes 1 and 2 represent the mean CUESIF of the node, n values represent the terminal node size, soil variables names under the boxes represent the splitting criterion, and the numbers in boxes 3–5 represent the mean estimated CUE. Note that the position of “yes” and “no” at each split represents the decision split at each node. (B) Boxplots showing the distribution of WUE in each terminal node in the regression tree. The model has r2 of 0.39 and a MAE of 0.175 gC/KgH2O for WUE.
Enlarge this image.Figure 5. Pruned regression tree for (A) WUESIF predicted as a function of one significant variable. The number in box 1 represents the mean WUESIF of the node, n value represents the terminal node size, soil variables names under the boxes represent the splitting criterion, and the numbers in boxes 2 and 3 represent the mean estimated CUE. Note that the position of “yes” and “no” at each split represents the decision split at each node. (B) Boxplots showing the distribution of WUE in each terminal branch in the regression tree. The model has r2 of 0.19 and MAE of 0.64 gC/KgH2O for WUE.
Enlarge this image.Figure 6. Pruned regression tree for (A) Landsat-based CUE predicted as a function of two significant variables. The number in boxes 1 and 3 represent the mean CUE of the node, the n value represents the terminal node size, soil variables names under the boxes represent the splitting criterion, and the numbers in boxes 2, 6, and 7 represent the mean estimated CUE. Note that the position of “yes” and “no” at each split represents the decision split at each node. (B) Boxplots showing the distribution of CUE in each terminal node in the regression tree. The model has r2 of 0.64 and MAE of 0.04 for CUE.
Enlarge this image.Figure 7. Pruned regression tree for (A) MODIS CUE predicted as a function of two significant variables The number in boxes 1 and 3 represent the mean CUE of the node, n value represents the terminal node size, soil variables names under the boxes represent the splitting criterion, and the numbers in boxes 2, 6, and 7 represent the mean estimated CUE. Note that the position of “yes” and “no” at each split represents the decision split at each node. (B) Boxplots showing the distribution of CUE in each terminal branch in the regression tree. The model has r2 of 0.78 and MAE of 0.06.
Enlarge this image.Figure 8. Pruned regression tree for (A) CUESIF predicted as a function of three significant variables. The number in boxes 1–4 represents the mean CUESIF of the node, n value represents the terminal node size, soil variables names under the boxes represent the splitting criterion, and the numbers in boxes 5–9 represent the mean estimated CUE. Note that the position of “yes” and “no” at each split represents the decision split at each node. (B) Boxplots showing the distribution of CUE in each terminal branch in the regression tree. The model has r2 of 0.45 and MAE of 0.1.
Latitude, longitude, and land cover type for the 19 sites used in this study.
| Site ID | Latitude | Longitude | Land Cover |
|---|---|---|---|
| DRJ | 37.40N | 88.34W | Agriculture |
| PTN | 37.40N | 87.70W | Agriculture |
| HWY | 37.68N | 86.98W | Agriculture |
| COR | 37.22N | 83.89W | Agriculture |
| LWF | 37.20N | 83.92W | Agriculture |
| UKS | 37.97N | 84.53W | Agriculture |
| ARB | 36.60N | 88.33W | Grassland |
| NFM | 36.63N | 88.33W | Grassland |
| FFG | 37.75N | 86.62W | Grassland |
| HLF | 37.88N | 84.46W | Grassland |
| HLP | 38.13N | 86.71W | Grassland |
| SPF | 38.13N | 84.50W | Grassland |
| BCF | 36.60N | 88.23W | Grassland |
| HBS | 36.76N | 88.12W | Deciduous Forest |
| CR | 36.93N | 88.46W | Deciduous Forest |
| HEA | 37.41N | 88.19W | Deciduous Forest |
| FFF | 37.75N | 86.63W | Deciduous Forest |
| RFS | 37.47N | 83.17W | Deciduous Forest |
| MMF | 36.69N | 88.15W | Deciduous Forest |
Summary of the soil chemical and physical properties for the study sites. Values are shown as the average for the A and B soil horizons. GR = grasslands, FR = forest, AG = agriculture.
| Site ID | Horizon | Top of Subsoil (cm) | Soil Order | BD | LOI | pH | N |
δ13C |
δ15N |
C(%) | C:N | P (g/m2) | K (g/m2) | Ca (g/m2) | Mg (g/m2) | Zn (g/m2) | Fe (mg/kg) | Land Cover |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| DRJ | A | 25 | Alfisol | 1.26 | 5.78 | 6.00 | 0.12 | −21.8 | 7.5 | 1.46 | 11.86 | 31.22 | 36.04 | 474.62 | 16.03 | 2.20 | 196.50 | AG |
| B | 1.34 | 4.75 | 5.75 | 0.05 | −24.1 | 7.5 | 0.33 | 6.68 | 5.21 | 21.88 | 437.97 | 32.62 | 0.68 | 135.25 | ||||
| PTN | A | 20 | Alfisol | 1.24 | 5.16 | 6.70 | 0.10 | −27.1 | 4.7 | 1.21 | 12.34 | 3.36 | 11.32 | 389.78 | 17.88 | 0.47 | 167.50 | AG |
| B | 1.40 | 3.24 | 5.27 | 0.04 | −22.9 | 8.1 | 0.33 | 8.29 | 1.20 | 17.90 | 197.72 | 28.13 | 0.03 | 110.33 | ||||
| HWY | A | 12 | Alfisol | 1.06 | 5.54 | - | 0.15 | −22.6 | 6.0 | 1.53 | 10.00 | 1.79 | 8.74 | 284.81 | 31.94 | 0.89 | 120.00 | AG |
| B | 1.26 | 3.53 | 5.30 | 0.04 | −23.0 | 7.5 | 0.34 | 7.85 | 0.31 | 7.96 | 171.36 | 21.88 | 0.10 | 58.00 | ||||
| COR | A | 30 | Alfisol | 1.15 | 3.75 | 5.77 | 0.07 | −24.8 | 7.7 | 0.67 | 9.84 | 93.70 | 26.60 | 180.08 | 17.86 | 0.35 | 275.33 | AG |
| B | 1.00 | 8.61 | 6.30 | 0.24 | −27.3 | 3.3 | 2.66 | 11.03 | 2.80 | 48.98 | 278.19 | 39.68 | 0.45 | 104.00 | ||||
| LWF | A | 26 | Alfisol | 1.43 | 5.17 | 5.13 | 0.09 | −25.8 | 7.1 | 0.84 | 9.86 | 128.37 | 18.49 | 161.81 | 14.05 | 0.32 | 317.00 | AG |
| B | 1.57 | 4.81 | 5.00 | 0.03 | −23.9 | 7.4 | 0.18 | 6.35 | 2.73 | 12.96 | 138.31 | 7.62 | 0.09 | 50.67 | ||||
| UKS | A | 22 | Alfisol | 1.16 | 6.92 | 5.55 | 0.11 | −23.7 | 7.7 | 1.03 | 8.02 | 21.13 | 79.52 | 311.54 | 32.11 | 0.65 | 151.50 | AG |
| B | 1.22 | 6.18 | 5.63 | 0.08 | −23.5 | 7.8 | 0.56 | 5.83 | 11.46 | 48.95 | 315.21 | 32.36 | 0.33 | 133.00 | ||||
| ARB | A | 24 | Alfisol | 1.23 | 4.97 | 5.70 | 0.13 | −21.3 | 7.2 | 1.19 | 9.05 | 2.35 | 10.54 | 194.84 | 32.09 | 0.70 | 149.33 | GR |
| B | 1.27 | 3.73 | 4.77 | 0.04 | −22.5 | 7.2 | 0.30 | 7.40 | 0.56 | 13.90 | 184.42 | 67.96 | 0.21 | 74.67 | ||||
| NFM | A | 4 | Alfisol | 0.83 | 7.74 | 5.30 | 0.29 | −20.7 | 6.3 | 2.57 | 8.86 | 10.54 | 70.50 | 310.14 | 45.62 | 1.97 | 203.00 | GR |
| B | 1.32 | 3.45 | 5.40 | 0.06 | −21.6 | 8.1 | 0.54 | 8.61 | 1.77 | 34.91 | 204.36 | 40.85 | 0.68 | 70.25 | ||||
| FFG | A | 23 | Alfisol | 1.21 | 4.63 | 5.75 | 0.08 | −25.0 | 5.6 | 0.76 | 9.88 | 2.19 | 9.30 | 233.70 | 12.75 | 0.11 | 107.75 | GR |
| B | 1.30 | 4.42 | 5.60 | 0.06 | −24.5 | 8.8 | 0.64 | 9.43 | 0.84 | 9.75 | 263.29 | 12.55 | 0.08 | 101.00 | ||||
| HLF | A | 23 | Alfisol | 1.09 | 9.42 | 5.85 | 0.26 | −25.1 | 6.3 | 2.48 | 9.58 | 48.53 | 75.21 | 357.72 | 41.92 | 0.85 | 203.00 | GR |
| B | 1.33 | 6.68 | 5.28 | 0.06 | −22.8 | 7.9 | 0.41 | 6.52 | 19.17 | 44.27 | 289.10 | 23.73 | 0.10 | 168.25 | ||||
| HLP | A | 22 | Alfisol | 1.24 | 5.93 | 5.55 | 0.13 | −23.2 | 7.2 | 1.03 | 6.79 | 7.96 | 19.67 | 301.06 | 18.05 | 0.47 | 135.00 | GR |
| B | 1.41 | 4.34 | 5.95 | 0.08 | −23 | 8.0 | 0.56 | 6.13 | 1.20 | 17.43 | 382.13 | 9.89 | 0.10 | 94.75 | ||||
| SPF | A | 29 | Alfisol | 1.41 | 8.73 | 6.90 | 0.20 | −24.1 | 6.5 | 1.76 | 8.93 | 65.27 | 18.01 | 726.50 | 33.66 | 19.00 | 221.33 | GR |
| B | 1.42 | 7.18 | 6.33 | 0.10 | −22.9 | 7.9 | 0.76 | 7.27 | 74.09 | 30.97 | 596.44 | 33.59 | 1.55 | 284.00 | ||||
| BCF | A | 14 | Alfisol | 1.43 | 4.97 | 5.90 | 0.11 | −25.8 | 6.8 | 1.20 | 11.15 | 30.66 | 64.56 | 268.39 | 36.09 | 1.00 | 270.50 | GR |
| B | 1.48 | 4.56 | 5.08 | 0.04 | −24.9 | 7.6 | 0.39 | 8.54 | 2.13 | 43.71 | 169.00 | 80.84 | 0.08 | 113.75 | ||||
| HBS | A | 7 | fUltisol | 0.84 | 25.56 | - | 0.35 | −27.9 | 0.1 | 5.90 | 16.57 | 8.46 | 16.36 | 124.97 | 22.30 | 0.73 | 173.00 | FR |
| B | 1.31 | 3.33 | 3.93 | 0.04 | −25.9 | 7.6 | 0.50 | 12.24 | 1.65 | 16.95 | 42.90 | 42.87 | 0.13 | 94.75 | ||||
| CR | A | 21 | Ultisol | 0.86 | 7.65 | 4.30 | 0.09 | −26.7 | 8.0 | 0.94 | 10.82 | - | - | - | - | - | - | FR |
| B | 1.19 | 4.10 | 4.20 | 0.03 | −25.5 | 6.4 | 0.26 | 9.03 | 3.11 | 14.74 | 22.22 | 63.69 | 0.69 | 214.00 | ||||
| HEA | A | 16 | Alfisol | 1.35 | 2.52 | 4.40 | 0.04 | −26.5 | 5.9 | 0.62 | 13.95 | 3.47 | 13.90 | 37.32 | 17.49 | 0.16 | 91.00 | FR |
| B | 1.24 | 3.47 | 3.98 | 0.02 | −25.4 | 8.3 | 0.27 | 10.57 | 0.78 | 12.69 | 42.56 | 37.07 | 0.15 | 67.75 | ||||
| FFF | A | 7 | Alfisol | - | - | - | - | - | - | - | - | - | - | - | - | - | - | FR |
| B | 1.30 | 4.64 | 7.13 | 0.08 | −25.5 | 3.4 | 0.80 | 9.71 | 0.95 | 41.58 | 805.61 | 33.82 | 0.21 | 73.25 | ||||
| RFS | A | 17 | Ultisol | 1.33 | 6.45 | 5.73 | 0.14 | −27.0 | 6.2 | 1.65 | 12.52 | 1.53 | 23.13 | 231.98 | 26.83 | 0.29 | 52.33 | FR |
| B | 1.24 | 5.09 | 5.68 | 0.08 | −26.4 | 6.3 | 0.84 | 10.53 | 0.92 | 23.09 | 107.91 | 24.71 | 0.09 | 104.00 | ||||
| MMF | A | 11 | Alfisol | 1.20 | 3.31 | 4.60 | 0.06 | −26.1 | - | 0.71 | 12.30 | 2.58 | 18.38 | 27.85 | 14.18 | 0.22 | 116.50 | FR |
| B | 1.40 | 4.89 | 4.00 | 0.04 | −24.3 | - | 0.34 | 8.41 | 0.50 | 35.47 | 89.42 | 108.67 | 0.11 | 102.25 |
Summary of CART r2 and MAE (gC/KgH2O for WUE) for the CART models. The models with the highest r2 and lowest MAE are highlighted in bold. CART = all soil variables used; CARTnLC = land cover excluded from the analysis; CARTnSO = soil order excluded from the analysis; CARTnSH = soil horizon excluded from the analysis; CARTnSD = soil depth excluded from the analysis.
| CART | CARTnLC | CARTnSO | CARTnSH | CARTnSD | ||
|---|---|---|---|---|---|---|
| WUE | r2 | 0.31 | 0.39 | 0.15 | 0.31 | 0.31 |
| MAE | 0.17 | 0.17 | 0.17 | 0.17 | 0.17 | |
| WUESIF | r2 | 0.19 | 0.19 | 0.19 | 0.19 | 0.19 |
| MAE | 1.4 | 0.64 | 0.64 | 0.64 | 0.64 | |
| CUE30 | r2 | 0.62 | 0.52 | 0.62 | 0.64 | 0.59 |
| MAE | 0.04 | 0.03 | 0.04 | 0.04 | 0.04 | |
| CUE250 | r2 | 0.78 | 0.61 | 0.78 | 0.79 | 0.78 |
| MAE | 0.06 | 0.07 | 0.06 | 0.06 | 0.06 | |
| CUESIF | r2 | 0.45 | 0.39 | 0.4 | 0.4 | 0.4 |
| MAE | 0.1 | 0.11 | 0.11 | 0.1 | 0.1 |
Summary of the GPP algorithm for the remote sensing data set used in this study.
| Data | GPP Algorithm | References |
|---|---|---|
| MODIS 250 m | Light use efficiency model using MODIS 250 m derived leaf area index and the fraction of photosynthetically active radiation | [40] |
| Landsat 30 m | Light use efficiency model using Landsat 30 m derived leaf area index and the fraction of photosynthetically active radiation | [40] |
| GOSIF GPP | OCO-2 SIF vs. GPP universal relationships | [66,67] |
© 2021 by the authors.
Details
; Stinchcomb, Gary E 2 ; Cetin, Haluk 1 ; Ferguson, Benedict 3 ; Kim, Sora L 4 ; Xiao, Jingfeng 5
; Fisher, Joshua B 6 1 Department of Earth and Environmental Science, Murray State University, Murray, KY 42071, USA;
2 Department of Earth and Environmental Science, Murray State University, Murray, KY 42071, USA;
3 Department of Earth and Environmental Science, Murray State University, Murray, KY 42071, USA;
4 Department of Life and Environmental Sciences, University of California Merced, Merced, CA 95343, USA;
5 Earth Systems Research Center, Institute for the Study of Earth, Oceans, and Space, University of New Hampshire, Durham, NH 03824, USA;
6 Jet Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91125, USA;