2 Materials and Methods

2.1 Study selection

A literature search was conducted using keywords “soil erosion”, “experiment” and “productivity” in Web of Science, Science Direct and others. This resulted in a sample of 69 studies published from 1969 to March 2019. Different research methods may result in different yield reduction rates due to erosion (Bakker et al., 2004). Most studies explored erosion effects on productivity by observing crop yields in field plots with different erosion levels (light, moderate, severe) classified according to various criteria; in this case, erosion effects on productivity cannot be accurately quantified (Duan et al., 2016; Lin et al., 2019). We used only those studies that clearly specified erosion depth. Then, crop yield variability was analyzed for six equidistant soil erosion depths (i.e., 5, 10, 15, 20, 25, and 30 cm). More details of the publication selection process can be found in Figure 1. The geographical distribution of the studies selected for meta-analysis included different areas across the globe as shown in Figure 2. Studies were published between 1995 and 2015, and regarded as eligible if they included the changing patterns of yield along erosion gradient, or yield states of eroded and non-eroded soils. Each study reported data of one or more relationships between yield and erosion depth were selected for the following meta-analysis. Subsequently, only 13 out of the 41 studies were selected after filtering titles, abstracts, and results. The selected 13 studies explicitly quantified crop yield response along erosion gradients.

2.2 Data processing

Crop yield data from different erosion depths at each site were extracted. In this process, WebPlotDigitizer software (Burda et al., 2017) was used to extract data from graphs, or data were directly retrieved from tables and main text. In addition, soil types, grain types, measure and geographical locations of each study were recorded. Two soil types (i.e., clay loam and sandy clay loams) involved in the dataset were standardized with China’s soil classification; and undetermined soil types were set to “test”. Three grain species (i.e., maize, soybeans, and wheat) were involved. Other important information that was expected to influence results was also collected, including monitoring year, and management practices (i.e., fertilization, fertilization add manure, irrigated site). Pearson’s correlation coefficients (r) for the relationships of erosion depth-crop reduction were either directly taken from the published studies, or calculated using soil erosion depth and crop yield reduction if reported for multiple plots.

2.3 Data analysis

Pearson’s correlation coefficient (r) for each case was normalized using Fisher’s z transformation for an effect size, sample size was the product of the number of repetitions and the number of erosion layers (Peng et al., 2019). Subsequently, the log-response ratio (lnRR) of crop yield in non-eroded and eroded treatments was used as the effect size in our meta-analysis, and was calculated as follows:
\begin{equation} lnRR=\ln{\ \left(T_{A}\right)}-ln\ (T_{B})\nonumber \\ \end{equation}
where TA was the average crop yield under erosion conditions, and TB was the average crop yield without erosion (erosion depth equal to 0 cm) in the same environment. All parameters in the meta-analytical models were estimated using maximum likelihood method (Zuur et al., 2009). The possibility of publication bias and temporal changes in effect size were examined using the fail-safe number. The procedure of Fisher’s z and lnRR were conducted using OpenMEE software (Wallace et al., 2017). Full analyses were performed using the Metafor package (Viechtbauer, 2010) in R software (version 3.5.2).
Meta-analysis assumes that individual studies are statistically independent; thus, obtaining crop yield for multiple erosion levels or several observations (e.g., cases in different site and year) from one publication could violate the assumption of independence and create a hierarchical dependence structure among the effect size estimates (Stevens and Taylor, 2008). Hence, data were analyzed with multilevel linear mixed-effect models using the ‘rma.mv’ function in Metafor (Viechtbauer, 2010). Additionally, models were fitted with nested random effect terms as follows: (ID |reference) (see Code meta-analysis), which performs noise processing on similar or non-independent cases, increasing the reliability of the analysis.