This study aims to identify the impact of vertical support on the information content of soil moisture (SM) for latent heat flux estimation. This objective is achieved via calculation of the mutual information (MI) content between multiple soil moisture variables (with different vertical supports) and current/future evaporative fraction (EF) using ground-based soil moisture and latent/sensible heat flux observations acquired from the AmeriFlux network within the contiguous United States. Through the intercomparison of MI results from different SM–EF pairs, the general value (for latent heat flux estimation) of superficial soil moisture observations , vertically integrated soil moisture observations , and vertically extrapolated soil moisture time series [soil wetness index (SWI) from a simple low-pass transformation of ] are examined. Results suggest that, contrary to expectations, 2-day averages of and have comparable mutual information with regards to EF. That is, there is no clear evidence that the information content for flux estimation is enhanced via deepening the vertical support of superficial soil moisture observations. In addition, the utility of SWI in monitoring and forecasting EF is partially dependent on the adopted parameterization of time-scale parameter T in the exponential filter. Similar results are obtained when analyses are conducted at the monthly time scale, only with larger error bars. The contrast between the results of this paper and past work focusing on utilizing soil moisture to predict vegetation condition demonstrates that the particular application should be considered when characterizing the information content of soil moisture time series measurements.
As most of the impacts of soil moisture for the climate system are induced by its impact on evapotranspiration in water-limited regimes (Seneviratne et al. 2010), soil moisture information is important for the investigation of land–atmosphere coupling. However, it is typically assumed that the shallow vertical penetration depth of microwave-based surface soil moisture retrievals (generally considered to be constrained to the top 5 cm of the soil column) limits their relevance for such investigations. In response, various methods have been developed to vertically extrapolate information contained in surface soil moisture retrievals to recover deeper soil moisture information. The most popular of these approaches include both sophisticated data assimilation strategies (Reichle et al. 2007; Sabater et al. 2007; Reichle et al. 2008; Kumar et al. 2009; Bolten and Crow 2012) and the application of semiempirical low-pass filtering to historical surface soil moisture time series data (Wagner et al. 1999; Ceballos et al. 2005; Albergel et al. 2008; Brocca et al. 2011; Mo et al. 2011; Manfreda et al. 2014; Crow et al. 2015).
A key unanswered question, however, is exactly how the information content of soil moisture, with regards to land–atmosphere interaction studies (especially surface flux prediction) and ecosystem productivity forecasting, increases as soil moisture time series observations are integrated over a deeper vertical depth. Qiu et al. (2014) looked at whether the vertical integration of soil moisture information improves the skills to forecast vegetation states (with lead time up to 40 days). In particular, using long-term soil moisture observations from the U.S. Department of Agriculture’s Soil Climate Analysis Network (SCAN), they calculated mutual information (MI) between multiple soil moisture variables (with various vertical supports) and near-future vegetation conditions to examine the existence of drought information. Results suggest that, relative to superficial soil moisture observations , the vertical integration of soil moisture (i.e., ) does indeed add value to the forecasting of near-future vegetation health anomalies. However, the enhanced information content in can be effectively matched through the low-pass transformation of using an intentionally simplistic global parameterization. This implies that, for the specific target of agricultural drought forecasting, the shallow vertical penetration depth (<5 cm) of microwave-based retrievals does not pose as large a practical limitation as is commonly perceived.
However, it is unclear if the results of Qiu et al. (2014) can be extended to other potential soil moisture applications. Climatologists (e.g., Budyko 1974; Eagleson 1978) have long developed different assumptions about the relationship between soil moisture (SM) and evaporative fraction (EF), that is, the ratio between latent heat flux (LE) and the sum of latent heat flux and sensible heat flux (SH). Among these assumptions, some research has interpreted the relationship as an exponential curve, whereas others have identified this relationship as piecewise linear function (or two-regime section), which transitions from a linear relationship under water-limited condition to an independent relationship under energy-limited condition (e.g., Komatsu 2003; Ford et al. 2014b). For instance, Koster et al. (2009) exploited the two-section evaporative regime to investigate the connection between seasonal meteorological drought and an increase in seasonal air temperature. Their results showed strong agreement between atmospheric general circulation models (AGCMs) and observation-based geographical patterns of drought-induced warming, which supports the idea that the same evaporative controls are also present in nature. Likewise, Ford et al. (2014b) examined seasonal variations on SM–EF interactions over the U.S. Southern Great Plains using both in situ observations and simulations from the Variable Infiltration Capacity model, and confirmed the well-identified two-section evaporative regime. More mechanistic studies include the proposal of a classical resistance-based method that links the soil evaporative efficiency with aerodynamic resistance and soil resistance, the latter of which can be expressed as a function of 0–5 cm soil moisture (Sellers et al. 1992). However, the impact of soil moisture vertical support on the SM–EF interactions over various land covers has received relatively little attention in the existing literature. In addition, the strength of the SM–EF relationship has generally been evaluated using a conventional Pearson product-moment correlation. However, since the relationship between SM and EF is certainly nonlinear, a linear correlation coefficient cannot fully capture the coupling strength between two variables. This issue motivates the use of a more integrative method (i.e., MI) to measure the SM–EF coupling strength in this study.
Here, we apply the methodology of Qiu et al. (2014) to examine the relationship between the vertical support of SM observations and their mutual information with respect to EF. In particular, we attempt to quantify the potential increase in information content, with respect to latent heat flux estimation, of surface soil moisture induced by integration over deeper vertical depths. This is based on analyzing the MI content between current/near-future surface energy flux time series and time series of three soil moisture datasets with various vertical supports: surface soil moisture, vertically integrated soil moisture, and a low-pass transformation of obtained at AmeriFlux sites within the contiguous United States. The extent to which the increased vertical measurement support increases the value of surface soil moisture observations for surface energy estimation is revealed by the differences in MI content (with respect to both current and near-future EF observations) among these three products. In addition, we also examine the impact of potential nonlinearity in the SM–EF relationship by comparing results obtained from the entropy-based metric of MI with comparable inferences based on a conventional Pearson correlation calculation.
2. Data and methods
AmeriFlux network provides continuous measurements of soil moisture, water vapor, energy fluxes, and related environmental variables, and the network covers a large variety of ecosystem types including forests, grasslands, croplands, shrublands, wetlands, and savannas (Boden et al. 2013). Before estimating mutual information between AmeriFlux SM and EF observations, we averaged datasets of observed SM , low-pass-filtered SM [soil wetness index (SWI)], and EF into 2-day compositing values. While the decision to utilize a 2-day compositing window is somewhat arbitrary, sensitivity analysis (not shown) reveals that results are relatively robust to the use of any compositing window size less than or equal to 5 days.
To remove the seasonal cycle and obtain anomaly time series, climatological averages for , , SWI, and EF were calculated for each 2-day period of the annual cycle (within the entire multiyear observation period) and subtracted from each raw 2-day time series. However, since the removal of the seasonal cycle could potentially impact the estimated coupling, all analyses are conducted with and without removal of the seasonal cycle. Temporal periods (December–February) with potential frozen soil conditions and unstable surface flux observations are masked from the analysis. Additional processing details are given below.
a. Soil moisture measurements
In situ soil moisture observations at a half-hour time step were collected from the AmeriFlux Site and Data Exploration System (see http://ameriflux.ornl.gov/). The level 2 (L2) standardized soil moisture data without gap filling were selected for this study. Soil moisture measurements at AmeriFlux sites were generally taken at two layers, and the measurement depths of top-layer and bottom-layer soil moisture vary between sites (Table 1). The sensitivity of results to these depth variations will be examined further in section 3a. Soil moisture observations were masked if they were uncorrected or did not fall within the (physically feasible) range of 0–0.6 m3 m−3 (Xia et al. 2015).
The top-layer measurements were used to represent θS variations (as obtained, for example, from microwave remote sensing). In contrast, θV was calculated as the depth-weighted average of both and the linearly interpolated value using and the bottom-layer soil moisture measurements. As a result, captures vertically averaged soil moisture conditions between the top layer and the bottom layer. Note that this vertical integration does not consider local vegetation rooting depths and should not be assumed to represent a true root-zone soil moisture product. Instead, simply represents a vertically integrated soil moisture product that (should) provide a relatively better approximation of true root-zone conditions than .
To reduce temporal sampling errors, only 36 AmeriFlux sites (four sites outside the contiguous United States were excluded) with over 5 years of two-layer soil moisture and surface flux observations were considered. These sites are located in a variety of climate zones within the contiguous United States. Figure 1 shows a map of these AmeriFlux sites using different symbols to represent landscape-scale land cover surrounding each plot-scale site.
b. Surface fluxes and vegetation cover observations
The EF was calculated as the ratio between surface latent heat flux and the sum of latent heat and sensible heat fluxes:
Since evaporation is very strongly tied to net radiation (Rn; Koster et al. 2009), this normalization of LE [i.e., evapotransporation (ET)] removes the impact of non-soil-moisture influences on ET, such as net radiation, wind speed, and soil heat flux G. The L2 AmeriFlux LE and SH observations are based on high-frequency (typically >10 Hz) eddy covariance measurements that are processed into 30-min averages by individual AmeriFlux investigators. Using these L2 products, we averaged LE and SH over a daytime period between 0700 and 1800 local standard time (LST) and calculated daily EF using these averaged LE and SH quantities. To ensure reliable surface flux observations, we utilized temporally coincident 30-min precipitation accumulation observations to exclude rainy conditions and masked days lacking any observations between 1130 and 1330 LST. In addition, we excluded days with energy budget closure rates lower than 0.6 (unitless) or higher than 1/0.6 (1.67) (unitless). Data removal ratio (removed days divided by actual sample size before removal) varies from 0% to 24% for different sites. The impact of this masking on analysis results will be discussed in section 3b.
To describe the long-term vegetation cover conditions around each site, we derived 2001–14 multiyear leaf area index (LAI) time series from 4-day composites of the 1-km Moderate Resolution Imaging Spectroradiometer (MODIS) MCD15A3 product (version 5). The LAI dataset was quality checked by the LAI–Fraction of Photosynthetically Active Radiation (FPAR) general quality assessment field in the MCD15A3 product and only high-quality retrievals, that is, those categorized as “good quality,” “detectors apparently fine for up to 50% of channels 1,2,” “significant clouds NOT present,” or “main method used, best result possible” were utilized. In addition to the MCD15A3 product, an extra LAI dataset was obtained from 8-day composites of the 1-km Global Land Surface Satellite (GLASS) LAI product (http://glcf.umd.edu/data/lai/) that was generated from time series of AVHRR/MODIS reflectance data using general regression neural networks. The 2001–14 multiyear LAI mean for each site derived from both of these LAI datasets (i.e., MCD15A3 and GLASS) is listed in Table 1.
c. SWI from low-pass filtering
Low-pass filtering was applied to surface observation time series θS in order to create SWI to function as a proxy for vertically integrated soil moisture. Following Wagner et al. (1999), this filtering was based on the assumption that the water flux between two soil layers is proportional to the difference in soil moisture content between these two layers. The resulting solution takes the form of a semiempirical exponential filter that converts superficial acquired between t − M and t into a proxy SWI estimate at time t:
where T is a response time scale (parameter) and M is the measurement length over which past information is integrated. With increasing elapsed time, the historical time series exert exponentially decaying weight on current SWI conditions; therefore, a cutoff threshold M is used for computation efficiency. In this study, the threshold was taken to be 3T unless otherwise stated (Albergel et al. 2008; Qiu et al. 2014). The response time scale is controlled by regional hydroclimatic conditions, land-cover types, and the local soil characteristics, which impact the vertical hydraulic connectivity between various soil layers (Albergel et al. 2008; Ford et al. 2014a). Specific strategies for parameterizing T are described later in section 2e.
d. Information measures
Mutual information (Cover and Thomas 1991) is a nonparametric measure of correlation (here defined strictly as the lack of independence) between two random variables. It is a more rigorous measure compared to the commonly used metrics such as Spearman’s rank correlation coefficient and Pearson product-moment correlation coefficient, the latter being an approximation of MI under certain conditions (Nearing and Gupta 2015).
MI and the related Shannon-type entropy (Shannon 1948) are calculated as follows. Entropy about a random variable ζ can be interpreted as a measure of uncertainty according to its distribution pζ and is estimated as the expectation [E(⋅)] of information from the pζ sample:
MI between ζ and another variable ψ can be thought of as the expected amount of information about variable ζ contained in a realization of ψ and is measured by the expected Kullback–Leibler divergence D (Kullback and Leibler 1951) between the conditional and marginal distributions over ζ:
In this context, the generic random variables ζ and ψ represent EF and soil moisture, respectively. The observation space of the target random variable EF was discretized using a Gaussian-optimized histogram bin width w given by Scott (2004):
where is the standard deviation of the given target random variable and k is the sample size of soil moisture and EF pairs. Following Nearing et al. (2016), a bin width of 0.01 m3 m−3 (1% volumetric water content) was applied for soil moisture. Integrations required for MI calculation in Eq. (4) are then approximated as summations over the empirical probability distribution function bins.
The rationale behind the application of this approach for SM–EF coupling strength is that the MI content between various soil moisture products and current/near-future EF can be used to quantify the value of each product for estimating/forecasting variations in surface energy fluxes and thus the general value of various soil moisture representations for land–atmosphere coupling applications. Specifically, we applied this approach to calculate the MI content between current soil moisture conditions (as reflected by , , and SWI) and current/near-future energy flux (in particular, EF) at each AmeriFlux site. All the estimated site-specific MI is normalized by the entropy of the corresponding AmeriFlux-based EF measurements to remove the effect of intersite variation on the magnitude of difference. The normalized MI (NMI) using , , and SWI are denoted as , , and NMI(SWI), respectively. By definition, MI between two variables represents the amount of entropy (uncertainty) in either of the two variables that can be reduced by knowing the other. Therefore, the MI normalized by the entropy of the AmeriFlux-based EF measurements represents the fraction of uncertainty in EF that is resolvable given knowledge of the soil moisture state (Nearing et al. 2013). Unlike Pearson’s correlation coefficient, MI is insensitive to the impact of nonlinear variable transformations. Therefore, it is well suited to describe the strength of the (potentially nonlinear) relationship between SM and EF.
To quantify the standard error of NMI differences between various soil moisture products, we applied a nonparametric bootstrapping approach with 1000 replicates and calculated pooled average differences across all sites as an estimation of error bars, assuming spatially independent sampling error.
e. Parameterization of T
Calculating SWI requires an estimate of the time-scale parameter. The general lack of information on microscale soil characteristics available at AmeriFlux sites necessitates a more generalized method to parameterize T, which involves the following three steps. The T was first optimized (i.e., Topt) at each site so that the corresponding SWI time series had the highest possible MI content with EF time series at zero time lag (τ = 0). Then, the resulting (optimized) values of Topt were used to establish a global relationship based on readily available observations from each site, including long-term (2001–14) annual mean LAI from MCD15A3 product, long-term (2001–12) annual mean LAI from GLASS product, annual mean precipitation (mm), annual mean net solar radiation (W m−2), annual mean air temperature (°C), annual mean relative humidity (%), and the estimated De Martonne aridity index (De Martonne 1926). Stepwise-regression results indicate that long-term climatological air temperature is the only effective regressor of Topt among this entire list of potential candidate predictors. Consequently, a single global linear relationship was established between the long-term mean air temperature at each site and Topt. Estimations of T at all selected AmeriFlux sites based on this single fixed relationship were denoted as Treg. It should be stressed that this regression relationship is treated here as a purely empirical result and that obtaining a clear theoretical motivation for its particular form is outside the scope of this study. In addition, to prevent unrealistic large estimates from this regression equation, we imposed an upper bound of 60 (days) on T.
Given that the success of SWI is likely to be dependent on the parameterization of T, we conducted the analysis using three different approaches for acquiring T: 1) T was individually optimized (i.e., Topt) at each site, which represents an upper limit on the SWI performance; 2) T was globally regressed (i.e., Treg) using the single linear relationship with annual mean air temperature as mentioned above, which seems to be a more realistic parameterization scheme; and 3) T was a constant value (i.e., Tcon) for each site, which was taken to be the mean of Topt from all selected AmeriFlux sites. The sensitivity of NMI(SWI) results to each of these parameterization approaches is discussed in section 3b.
a. Comparison between NMI() and NMI()
Before examining the difference between and , we first show in Fig. 2 the overall magnitude of and to illustrate the degree to which surface soil moisture and vertically integrated soil moisture contain EF information. In particular, the results at zero time lag (τ = 0) reflect the skill of soil moisture product in diagnosing instantaneous land–atmosphere feedback, while positive lags (τ > 0) reflect the skill of prediction of future EF using current soil moisture conditions. It is seen that for anomaly time series of SM–EF pairs, the seasonal variation of and are relatively limited, although maximum values of NMI(θV) are notably higher in summer (June–August). For spring (March–May), summer (June–August), and autumn (September–November), the median of and are about 0.3, which indicates that approximately 30% of uncertainty in EF can be reduced given the knowledge of surface or vertically integrated soil moisture state.
Based on this, Fig. 3 plots the seasonal variation of NMI differences between and at different temporal lags and with 2σ sampling errors calculated from bootstrapping. Results in Fig. 3 demonstrate that, in general, vertically integrated soil moisture does not contain more information for diagnosing or predicting EF than surface soil moisture. This tendency is not sensitive to seasonal variations and is present regardless of whether the analysis is conducted on anomaly soil moisture and EF time series (Fig. 3a) or on the original time series containing seasonality (Fig. 3b).
Figure 4 plots the NMI difference between and across various land-use types and at different temporal lags using both the anomaly time series (Fig. 4a) and the original time series (Fig. 4b) of the entire year (winter excluded). The error bars in Fig. 4 also represent 2σ sampling error calculated using a 1000-replicate bootstrapping ensemble. Figure 4a illustrates that is not significantly different from for any land-cover type except grassland. For results using the original time series (Fig. 4b), is nearly indistinguishable from for all land-use types (i.e., the ±2σ confidence interval for plotted differences typically contains zero). Furthermore, despite significant differences in values across different land-use types, the multiple pairwise comparison of the land-cover-based difference between and show no significant difference based on results obtained from a nonparametric version of classical ANOVA test (i.e., Kruskal–Wallis test). Based on Figs. 3 and 4, NMI differences are qualitatively similar for both anomaly and the original time series. Therefore, unless otherwise stated, all future NMI results will be based on the original time series (with no seasonal cycle removal).
In terms of geographical distribution, symbol colors in Fig. 1 capture the z scores for minus differences (scaled by their corresponding 1σ sampling errors) for each of the selected AmeriFlux sites at lag τ = 0. In addition to confirming that is generally highly comparable to in terms of its mutual information with EF for most sites, Fig. 1 demonstrates that no clear spatial, climatic, or land-cover patterns emerge in the sampled differences between and .
It is noteworthy that in Qiu et al. (2014), SCAN observations taken from fixed soil depths are employed to estimate the utility of multiple soil moisture variables in forecasting near-future vegetation conditions. However, in this study, the soil moisture measurements of AmeriFlux were taken from different depths (Table 1). To examine the potential impact of these variations, Fig. 5 plots the relationship between minus and the vertical depths for measurements at the AmeriFlux sites listed in Table 1. The lack of any clear trend in Fig. 5 suggests that measurement depth variations between AmeriFlux sites do not exert a systematic bias on differences between and .
b. The added value of exponential filtering
In addition to the direct use of , we examined the impact of smoothing via the exponential filter in Eq. (2) to produce SWI. Using the three different approaches for estimating the SWI T parameter (i.e., Treg, Topt, and Tcon; see section 2e), NMI(SWI) results for all land-cover types using observations from spring through autumn are plotted in Fig. 6.
Under the optimal parameterization conditions, SWI better diagnoses and forecasts EF than either or . That is, NMI(SWI) is greater than both and in Fig. 6 for the Topt case. However, in practice, such optimal parameterization will likely require either extensive site-specific calibration or a high level of site-specific information and will therefore be impossible to apply over spatially distributed domains. In contrast, for the cases of applying more practical globally parameterized and constant T approaches (Treg and Tcon), SWI is a comparable diagnostic EF monitor than (at lag τ = 0), and SWI does not improve upon either. In addition, when mapped in space (Fig. 7), sites with positive z scores for NMI(SWI) minus differences at lag τ = 0 do not demonstrate any clear geographical pattern.
c. Sensitivity to temporal scales and evaluation metrics
1) Sensitivity to temporal scales
To examine whether our previous conclusions are robust at multiple time scales, we also conducted the above analysis at a monthly time scale in comparison to the original 2-day time scale. Specifically, we selected AmeriFlux sites with temporal coverage of more than 10 years to conduct the analysis at monthly time scale. Overall, as the criterion of 10-yr coverage reduces the site number and reduces the sampling power, error bars are generally wider than those shown in Figs. 3 and 4 for 2-day results. However, outside of these differences, modifying the analysis to be based on a monthly time scale does not lead to any qualitative changes in results.
2) Sensitivity to evaluation metrics
As discussed above, our use of NMI as an evaluation metric is based on the expectation that the relationship between soil moisture and EF will often be nonlinear and on the robustness of NMI in such cases. However, it is also instructive to consider the impact of other evaluation metrics. Therefore, in addition to NMI-based results in Fig. 6, we include results in Fig. 8 based on applying the Pearson correlation coefficient R to describe the strength of the linear relationship between soil moisture and current/near-future EF. The R between EF and three soil moisture time series are denoted as , , and R(SWI), respectively. To ensure consistency with NMI-based analysis in Fig. 6, sampled values of R(SWI) minus are also based on three different T parameterization schemes mentioned in section 2e, except that Topt is now determined by maximizing the Pearson R (instead of MI) between SWI and EF at lag τ = 0.
Results based on the Pearson R are qualitatively similar to earlier NMI-based results. In particular, they reflect that surface measurements θS are equally, or even slightly more, valuable than deeper measurements θV for diagnosing and/or predicting EF. As with earlier NMI-based results, SWI performance is significantly better than if the T parameter is individually optimized (i.e., using Topt) for each site. In contrast, for the constant T parameterization (i.e., using Tcon), the transformation of into SWI actually degrades the strength of the linear relationship between soil moisture and the (current) EF time series for short time lags. In addition, relative to NMI-based results in Fig. 6, R-based results in Fig. 8 suggest relatively better SWI results associated with the derivation of T from a global regression equation (i.e., using Treg).
In this study, we calculated the normalized mutual information (NMI) between time series of three soil moisture products (with various vertical supports) and current/near-future EF time series to identify the impact of vertical support on information content of soil moisture for latent heat flux estimation. Specifically, we examined the surface energy flux information content contained in 1) vertically integrated soil moisture observations , 2) superficial observations , and 3) a simple low-pass transformation of (i.e., SWI).
Overall, results demonstrate that (surface only) and (vertically integrated) contain comparable levels of MI relative to current and near-future EF. For the specific target of surface flux (EF) prediction, major differences in and are not observed over various land-use types, different seasons, and different hydroclimatic conditions (Figs. 1–4). This indicates that, contrary to expectations, the information content of soil moisture (for EF prediction) does not increase significantly as the soil moisture observations are integrated over deeper vertical depths. In addition, the degree to which the SWI time series can match in EF monitoring/forecasting is strongly dependent on T parameterization schemes. Using site-calibrated T, SWI is a significantly better monitor and forecaster of EF than for all nonnegative lags less than 10 days. However, when utilizing a globally parameterized T or/and constant T, SWI is comparable to θV and θS as a diagnostic EF monitor (Figs. 6, 7). Similar results are obtained when analyses are conducted at the monthly time scale and using evaluation metrics of Pearson R (Fig. 8).
These results are somewhat at odds with earlier results in Qiu et al. (2014) demonstrating that, using a parsimonious globally regressed T parameter, the application of exponential filtering (to convert into SWI) significantly increased the MI between soil moisture and near-future normalized difference vegetation index (NDVI). This inconsistency demonstrates that the particular application (i.e., EF monitoring vs NDVI forecasting) must be defined when describing the information content of soil moisture time series. Additionally, T parameterizations with different levels of complexity should be adopted when targeting particular application. In this case, different target applications and different T parameterization approaches lead to different conclusions regarding the value of exponential filtering.
In summary, we see no clear evidence that the information content for surface energy flux predictions is enhanced by the vertical integration of superficial soil moisture observations. The results of this study suggest that it is unclear how soil moisture variations impact the surface energy balance and, therefore, how the land surface wetness measurements can (and should) be used to better constrain surface energy balances. Nonetheless, at the very least, we should reconsider the common assumption that the vertical support of satellite-based surface soil moisture retrievals must be increased before they can meaningfully be available to surface energy balance studies. However, it should be noted that our results are limited to just EF monitoring and prediction and the specific datasets that were examined here. To better understand how the soil moisture variations could be used to constrain surface energy balances, similar analyses using various sources of datasets (e.g., from both land surface models and remote sensing) should be addressed in future research.
This work was supported by National Natural Science Foundation of China (Grant 41501450), Natural Science Foundation of Guangdong Province, China (Grant 2016A030310154), and the Fundamental Research Funds for the Central Universities (16lgpy06). We thank the anonymous reviewers for their helpful comments.