The Root-Soil Water Relationship Is Spatially Anisotropic in Shrub-Encroached Grassland in North China: Evidence from GPR Investigation

: To analyze the root-soil water relationship at the stand level, we integrated ground-penetrating radar (GPR), which characterized the distribution of lateral coarse roots (>2 mm in diameter) of shrubs ( Caragana microphylla Lam.), with soil core sampling, which mapped soil water content (SWC) distribution. GPR surveys and soil sampling were carried out in two plots (Plot 1 in 2017 and Plot 2 in 2018) with the same size (30 × 30 m 2 ) in the sandy soil of the semi-arid shrubland in northern China. First, the survey area was divided into ﬁve depth intervals, i.e., 0–20, 20–40, 40–60, 60–80, and 80–100 cm. Each depth interval was then divided into three zones in the horizontal direction, including root-rich canopy-covered area, root-rich canopy-free area, and root-poor area, to indicate different surface distances to the canopy. The generalized additive models (GAMs) were used to analyze the correlation between root distribution density and SWC after the spatial autocorrelation of each variable was eliminated. Results showed that the root-soil water relationship varies between the vertical and horizontal directions. Vertically, more roots are distributed in soil with high SWC and fewer roots in soil with low SWC. Namely, root distribution density is positively correlated with SWC in the vertical direction. Horizontally, the root-soil water relationship is, however, more complex. In the canopy-free area of Plot 1, the root-soil water relationship was signiﬁcant ( p < 0.05) and negatively correlated in the middle two depth intervals (20–40 cm and 40–60 cm). In the same two depth intervals in the canopy-free area of Plot 2, the root-soil water relationship was also signiﬁcant ( p < 0.01) but non-monotonic correlated, that is, with the root distribution density increasing, the mean SWC decreased ﬁrst and then increased. Moreover, we discussed possible mechanisms, e.g., root water uptake, 3D root distribution, preferential ﬂow along roots, and different growing stages, which might lead to the spatially anisotropic relationship between root distribution and SWC at the stand level. This study demonstrates the advantages of GPR in ecohydrology studies at the ﬁeld scale that is challenging for traditional methods. Results reported here complement existing knowledge about the root-soil water relationship in semi-arid environments and shed new insights on modeling the complex ecohydrological processes in the root zone.


Introduction
The root-soil water relationship has been increasingly recognized as a vital component of the Earth's surface ecosystem, determining plants' growth and succession in arid and semi-arid ecosystems [1][2][3][4]. It also influences the water cycle [5], critical zone processes [6], According to the Food and Agriculture Organization (FAO) and the United States Department of Agriculture (USDA) soil classification systems, this region's soils are classified as Calcic Kastanozems and Calcic Orthic Aridisol, respectively [46][47][48]. These soils are derived from loess and are mainly distributed in the eastern and central parts of the Mongolian Plateau. The soils of this region are excessively drained with high permeability and have relatively homogeneous physiochemical properties and low organic content. These properties are conducive to GPR-based root detection [35,40].
C. microphylla is the dominant shrub species in the study site. It is strongly tolerant of cold, dry, and barren conditions [42,49] and controls the sequestration of organic carbon, nitrogen accumulation, and the hydrological cycle in the study region [49].

Field Experiment Design
Field experiments were conducted on 12 July 2017, in Plot 1 and 16 August 2018, in Plot 2. There was no precipitation in the three days before or during the experimentations. In each field experiment, a relatively flat area was selected as the experimental site, and a square plot with sides of 30 m was established. Figure 1b shows satellite imagery of the two plots; Figure 1c,d are photographs of Plot 1 and Plot 2, respectively.

Field Experiment Design
Field experiments were conducted on 12 July 2017, in Plot 1 and 16 August 2018, in Plot 2. There was no precipitation in the three days before or during the experimentations. In each field experiment, a relatively flat area was selected as the experimental site, and a square plot with sides of 30 m was established. Figure 1b shows satellite imagery of the two plots; Figure 1c,d are photographs of Plot 1 and Plot 2, respectively. First, the above-ground parameters of each C. microphylla, including its spatial coordinates, the lengths of the long and short axes of its canopy, and its height, were measured in each plot. The parameters of the canopy's long diameter in the study are shown in Table 1. The C. microphylla shoots were removed, and the ground surface was cleaned. Orthogonal grid survey lines at 25-cm intervals were set up for the GPR survey ( Figure 2a). Each square plot included 240 GPR survey lines (120 along each axis). Soil auger sampling locations were laid out at 3-m intervals (Figure 2b), yielding a total of 100 evenly distributed SWC samples in each plot. in each plot. The parameters of the canopy's long diameter in the study are shown in Table  1. The C. microphylla shoots were removed, and the ground surface was cleaned. Orthogonal grid survey lines at 25-cm intervals were set up for the GPR survey ( Figure 2a). Each square plot included 240 GPR survey lines (120 along each axis). Soil auger sampling locations were laid out at 3-m intervals (Figure 2b), yielding a total of 100 evenly distributed SWC samples in each plot.  The distribution soil auger sampling locations, a total of 100 with the intervals of 3-m.

GPR Scanning and Processing
A field-portable GPR system (RIS MF Hi-Mod; Ingegneria Dei Sistemi Inc., Pisa, Italy) with two pairs of antennae, shielded in the same antenna box, was used to scan the shrub root system along the survey lines. It could simultaneously collect data of two frequencies (900 and 400 MHz). A survey wheel attached to the antenna box was used to record the survey distance. There were 512 samples per trace at a time step of 0.0586 (or 0.1172) ns, for a total recorded length of 30 ns for the 900 MHz antenna and 60 ns for the 400 MHz antenna. The traces were triggered every 1.6 cm for the 900 MHz antenna and every 3.2 cm for the 400 MHz antenna along each survey line.
The GPR data processing software Reflex-Win 7.2 (Sandmeier Scientific Software, Germany) was used to process the field-collected GPR data. The primary processing steps were the first break correction, background noise removal, amplitude gaining, and travel time to depth conversion. First break correction was used to correct the vertical and horizontal scales on each radargram and ensure that all reflection two-way travel times were correct. Based on the significant difference of energy before and after the first arrival time point, the maximum ratio of the summed energy from two moving time windows was applied on each trace to detect the first break time [30,50]. The background noise on the

GPR Scanning and Processing
A field-portable GPR system (RIS MF Hi-Mod; Ingegneria Dei Sistemi Inc., Pisa, Italy) with two pairs of antennae, shielded in the same antenna box, was used to scan the shrub root system along the survey lines. It could simultaneously collect data of two frequencies (900 and 400 MHz). A survey wheel attached to the antenna box was used to record the survey distance. There were 512 samples per trace at a time step of 0.0586 (or 0.1172) ns, for a total recorded length of 30 ns for the 900 MHz antenna and 60 ns for the 400 MHz antenna. The traces were triggered every 1.6 cm for the 900 MHz antenna and every 3.2 cm for the 400 MHz antenna along each survey line.
The GPR data processing software Reflex-Win 7.2 (Sandmeier Scientific Software, Karlsruhe, Germany) was used to process the field-collected GPR data. The primary processing steps were the first break correction, background noise removal, amplitude gaining, and travel time to depth conversion. First break correction was used to correct the vertical and horizontal scales on each radargram and ensure that all reflection two-way travel times were correct. Based on the significant difference of energy before and after the first arrival time point, the maximum ratio of the summed energy from two moving time windows was applied on each trace to detect the first break time [30,50]. The background noise on the GPR radargram included horizontal bands caused by antenna-ground interactions or nearby radio interference, multiple reflections from some below-ground reflectors, and high-frequency spike events. Some of this noise was high-frequency and some was lowfrequency. The bandpass wave filtering method was used to eliminate the high-frequency and low-frequency noise in each profile. To compensate for the energy lost to medium attenuation, scattering, and dissipation, a combination of linear and exponential gain was applied on radargram to intensify the visibility of hyperbolas reflected by root. Finally, signal velocity was determined according to the hyperbola shape and then used to convert the travel time to depth, and the root depth could be obtained. All the data processing methods used were a standard GPR processing routine of REFLEX-WIN 7.2 software and can be found in the user manual. Therefore, the detailed processing procedures are not repeated here.
After radar data processing, all coarse roots were automatically identified from the post-processed GPR radargrams. The reflected signal from coarse roots has a characteristic hyperbolic shape ( Figure A1), and we used the randomized Hough transform (RHT) algorithm to automatically identify such hyperbolic reflections in the GPR images [40]. The peak of a hyperbola was considered the location of the lateral coarse root ( Figure A1), and each located root was counted as a root point. We summed the number of identified root points within each plot. These processes were performed with MATLAB (MathWorks, Natick, MA, USA). Detailed descriptions of automatic identification and location of coarse roots are provided in Li et al. [40] and Liu et al. [36]. Field validation of root detection accuracy by dual-frequency GPR in the study region was performed in our previous study [30].

Root Points Cloud Reconstruction from GPR Data
We combined the root points identified by the two frequencies (900 MHz and 400 MHz) in the 3D space within each plot. The two antenna frequencies had different target resolutions and detection depths, which influenced their performance in detecting roots. In some cases, one root might be errantly identified in different locations in the radargrams for the two frequencies. To prevent this, we carefully combined the root point data. First, we divided the root points into two groups: those identified only at the 900 MHz frequency and those identified only at the 400 MHz frequency. For the i th root point (R i ) obtained under the 900 MHz frequency, a spherical search volume, centered at R i and with a radius of 5 cm, was created in the 3D space ( Figure A2). If a root point (R j ) obtained under the 400 MHz frequency was found in the search volume, as determined by the distance between R i and R j (denoted by L in Figure A2) being less than 5 cm, it was considered the same root point (R i ). If no root point appeared in the search volume, R i was considered a root point only detected by the 900 MHz antenna. The same analysis was conducted for root points obtained under 400 MHz frequency. Finally, all root points were combined into a complete 3D root distribution map.

SWC Measurement
One soil sample was extracted from each sampling location in each plot using a soil auger (5 cm inner diameter and 1 m length; Figure 2b). Each soil sample was divided into five parts after extraction according to the depth ranges (0-20, 20-40, 40-60, 60-80, and 80-100 cm), and each section was quickly packed into an aluminum box and weighed in the field. The soil samples were then brought back to the laboratory and dried in the oven at a temperature of 105 • C until a constant weight was reached. The SWC was calculated based on the mass loss: where θ v is the SWC, W f is the weight of the fresh soil sample, W d is the weight of the dried soil sample, and ρ is the bulk density of soil. For the convenience of subsequent analysis, the Ordinary Kriging method was used to estimate the planar distribution of SWC within each depth interval separately. In order to ensure the interpolation accuracy, we checked the original SWC data of each depth interval in the following steps: firstly, we used a histogram and Quantile-quantile plot to verify whether the SWC data of each depth interval satisfied the normal distribution. Secondly, the semi-variogram was used to test whether there was spatial autocorrelation in the SWC data of each depth interval. In our case, the SWC data of each interval followed a normal distribution, and the spatial autocorrelation existed. Finally, by comparing different semi-variance models, the Gaussian function was selected for interpolation. The results were divided into 600 × 600 grids. The side length of the grid cell is 5 cm. Spatial interpolation was performed in ArcGIS 10.2.2 [51].

Exploring Root-Soil Water Relationship on Different Axes
As shown in Figure 3b,d, each plot volume was divided into five depth intervals (0-20, 20-40, 40-60, 60-80, and 80-100 cm), corresponding to soil sampling. The number of root points and the mean SWC of each depth interval were calculated. The root distribution density and the mean SWC of each canopy gradient zone in the same vertical interval were calculated to reveal the root distribution patterns and the spatial relationship between root and SWC. Specifically, the root distribution density of a canopy gradient zone was represented by the quotient of the total number of root points in each grid cell and the total area of the grid cells with the same gradient value. This was calculated as: where is the root distribution density of the m th canopy gradient zone, is the number of grid cells in the m th canopy gradient zone, and is the number of root points of the i th grid in the m th canopy gradient zone. The mean SWC of a canopy gradient zone was represented by the mean value of the SWC in the grid cells with the same gradient value: where is the mean SWC in the m th canopy gradient zone, is the number of the grid cells in the m th canopy gradient zone, and is the SWC of the i th grid of the m th canopy gradient zone.
The entire study area was divided into three regions based on their canopy gradient zones, from near to far: root-rich and canopy-covered, root-rich and canopy-free, and rootpoor and canopy-free. The 10% of gradient zones with the fewest root points were cate- The exploration of the horizontal relationship was complicated by the range of ages of the shrubs since root scope generally increases with age, and the ratio of the above-ground crown scope to the underground root scope of natural shrubs is relatively stable [52,53]. To better explore the relationship between roots and soil water in different physiological portions of the root structures, each interval obtained by the depth stratification method was divided into canopy gradient zones. These zones were defined by their distance from the taproot collar of the nearest shrub, based on the gradient stratification method [24,54]. The distance of each canopy gradient zone was weighted based on shrub size. To obtain the horizontal canopy gradient zones, we discretized each plot into a 600 × 600 grid, such that the grid spacing was 5 cm × 5 cm. Each grid cell was assigned a value (k) for the canopy gradient zone in which it was located. A higher value of k indicated a greater distance from the shrub's root collar. The value of k in the i th grid cell was calculated as: where d i,j represents the distance from the center of the ith grid cell to the jth shrub, and C j represents the long axis dimension of the canopy of the jth shrub. The 0.1 constant in the denominator indicates that we took 10% of the long axis of the shrub canopy as unit distance and extended outward from there. The "round" function removed any decimal component, thereby rounding down to the nearest integer. With this formula, regions with the same k value in each plot will cover the same canopy gradient zone. The root distribution density and the mean SWC of each canopy gradient zone in the same vertical interval were calculated to reveal the root distribution patterns and the spatial relationship between root and SWC. Specifically, the root distribution density of a canopy gradient zone was represented by the quotient of the total number of root points in each grid cell and the total area of the grid cells with the same gradient value. This was calculated as: where ρ m is the root distribution density of the mth canopy gradient zone, n m is the number of grid cells in the mth canopy gradient zone, and N m i is the number of root points of the ith grid in the mth canopy gradient zone. The mean SWC of a canopy gradient zone was represented by the mean value of the SWC in the grid cells with the same gradient value: where θ m is the mean SWC in the mth canopy gradient zone, n m is the number of the grid cells in the mth canopy gradient zone, and θ m i is the SWC of the ith grid of the mth canopy gradient zone. The entire study area was divided into three regions based on their canopy gradient zones, from near to far: root-rich and canopy-covered, root-rich and canopy-free, and root-poor and canopy-free. The 10% of gradient zones with the fewest root points were categorized as root-poor areas, while the rest was considered root-rich.
The root-rich area was used to explore the root-soil water relationship, which was subsequently divided into two parts: the canopy-covered area, which comprised the first five gradient zones, and the canopy-free area, which took the rest of the gradient zones. The root-soil water relationships were examined and discussed separately for these categories because the distribution of canopy is associated with wind protection, shelter, and precipitation pooling [55,56] that further impacts soil water patterns.

Spatial Correlation between Root Distribution Density and SWC
To account for the potential non-monotonic relationships between variables, we used the generalized additive models (GAMs) with the spatial term to analyze the complicated root-soil water relationship [57]. The GAMs have the advantage of finding the "segmentation point" automatically and getting different or the same root-soil water relationship on the left and right side of the "segmentation point", then making a smooth connection to form a continuous line to represent the root-soil water relationship. Therefore, this model can capture the unknown root-soil water relationship more accurately.
The GPR scanning space was divided horizontally and vertically in 2D to better explore the spatially anisotropic root-soil water relationship. Vertically, the correlation between the number of the root points and the mean SWC was analyzed, and the spatial location was added to the GAMs as a variable: where θ i is mean SWC of the ith interval, α represents the intercept of the model, f (I i ) is the spatial smooth function to account for the spatial pattern of θ that eliminates the spatial autocorrelation, and I i is the relative spatial position of the intervals used to divide the GPR scanning space. f (N i ) is the smooth function of the number of the root points, N i is the number of the root points of the ith interval, and ε i is the error term of the model. The spatial term can reveal the correlation between SWC and spatial location, that is, the Remote Sens. 2021, 13, 1137 9 of 24 inherent spatial autocorrelation of SWC. Compared with the GAMs models without spatial term, the GAMs with the spatial term can eliminate the effect of the spatial autocorrelation and highlight the real "signal" of the correlation between the number of the root points and SWC at the same time.
Horizontally, the correlations between root distribution density and mean SWC in the canopy-covered and canopy-free areas of five intervals were analyzed in the similar way as did for the vertical direction: where θ i is mean SWC of the ith canopy gradient zone. f (Z i ) is the spatial smooth function to account for the spatial pattern of the θ and eliminates the spatial autocorrelation, and Z i is the relative spatial position of the canopy gradient zones. f (ρ i ) is the smooth function of the root density and ρ i is the root density of the ith canopy gradient zone. The thin plate regression splines were selected as the smooth functions because they have the advantages of reducing the computational requirements of the smoothing splines and avoiding the problems of the knot placement of the regression splines [58]. After the GAMs were established, the P-values were calculated to judge whether there was a significant correlation between root distribution density and SWC. When the P-value of the root independent variable (N i in Equation (5) and ρ i in Equation (6)) is less than 0.05, it is considered a significant correlation between the root distribution density and SWC, and the specific relationships were analyzed by drawing the smooth function fitting diagrams. Finally, we plotted the semi-variograms of the residuals for each model to check whether the spatial autocorrelation was eliminated. All statistical analyses were carried out in the R package "mgcv" [59]. The distribution of the number of root points in different depth intervals within the two plots is shown in the left panel of Figure 4. Both plots showed a trend of increasing from the surface to the mid-depth and then decreasing with additional depth. The number of root points in the surface interval (0-20 cm) was the smallest, with only 398 in Plot 1 and 355 in Plot 2; but the number peaked in the third interval (40-60 cm), with more than 3800 in Plot 1 and 4000 in Plot 2.

Vertical Distribution and Correlation of Roots and Soil Water
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 25

Distribution Pattern of Roots
After processing, 10,973 and 12,113 coarse root points were obtained from the GPR scans of Plot 1 and Plot 2, respectively. Figure 3a,c show the 3D distribution of lateral coarse roots in Plot 1 and Plot 2, based on the combined GPR root detection.
The distribution of the number of root points in different depth intervals within the two plots is shown in the left panel of Figure 4. Both plots showed a trend of increasing from the surface to the mid-depth and then decreasing with additional depth. The number of root points in the surface interval (0-20 cm) was the smallest, with only 398 in Plot 1 and 355 in Plot 2; but the number peaked in the third interval (40-60 cm), with more than 3800 in Plot 1 and 4000 in Plot 2.

Distribution Pattern of Soil Water
We collected a total of 500 soil samples (100 sampling points with five depths for each sampling point) inside each plot. The right panel of Figure 4 is the distribution of the mean SWC in the different depth intervals of the two plots. The plots had similar trends: SWC increased from the surface to the moderate depth before decreasing in the deep soil. However, there were differences between them. For Plot 1, the range of the mean SWC was

Distribution Pattern of Soil Water
We collected a total of 500 soil samples (100 sampling points with five depths for each sampling point) inside each plot. The right panel of Figure 4 is the distribution of the mean SWC in the different depth intervals of the two plots. The plots had similar trends: SWC increased from the surface to the moderate depth before decreasing in the deep soil. However, there were differences between them. For Plot 1, the range of the mean SWC was from 4.57 to 6.87%, with the minimum value from 0 to 20 cm and the maximum value between 40 and 60 cm. In contrast, the SWC in Plot 2 ranged from 4.56 to 7.65%, and the maximum value appeared closer to the surface (20-40 cm).

Correlation between Roots and Soil Water in the Vertical Direction
In the vertical direction, the number of root points was positively correlated with mean SWC in Plot 1 and Plot 2 ( Figure 5). Namely, more roots are distributed in soil with high SWC and fewer roots in soil with low SWC. However, the above correlations were only significant in Plot 2 ( Table 2). The percentage of deviance is explained by the GAMs ranges from 98.3% to 93.4%. The spatial term eliminated the spatial autocorrelation of the residuals ( Figure A4).

Horizontal Distribution and Correlation of Roots and Soil Water
The values of k ranged from 1 to 54 in Plot 1 and from 1 to 37 in Plot 2, so the e mental plots were divided into canopy gradient zones named Z1, Z2, ..., Z54 and Z ..., Z37, respectively ( Figure 6). The specific range of the gradient zone for three areas rich canopy-covered area, root-rich canopy-free area, and root-poor area) in di depth intervals is shown in Table 3. Table 3. The range of the gradient zone for three areas of different depth intervals in Plot 1 a 2. Zone 1 represents the root-rich canopy-covered area, Zone 2 represents the root-rich cano area, and Zone 3 represents the root-poor area.

Horizontal Distribution and Correlation of Roots and Soil Water
The values of k ranged from 1 to 54 in Plot 1 and from 1 to 37 in Plot 2, so the experimental plots were divided into canopy gradient zones named Z1, Z2, ..., Z54 and Z1, Z2, ..., Z37, respectively ( Figure 6). The specific range of the gradient zone for three areas (root-rich canopy-covered area, root-rich canopy-free area, and root-poor area) in different depth intervals is shown in Table 3. Combined with the horizontal canopy gradient zone, it is clear that areas close to the plant's taproot collar had denser root point detection, which is in good accordance with the basic principles of root growth. The pattern is more evident in each soil depth interval (as shown in Figure A3). This phenomenon is also shown in Figure 6c,d, i.e., that 90% of the root points were distributed in the first 24 gradient zones and 20 gradient zones for Plot 1 and Plot 2. However, root aggregation varied across intervals, decreasing from shallow to deep soil. The statistical results in Figure 6c,d illustrate this pattern: in Plot 1, the first 30% of the gradient zones (Z1 to Z15) contained 88%, 77%, 70%, 68%, and 63% of the total number of root points in the first through fifth soil depth intervals, respectively. In Plot 2, the first 30% of gradient zones (Z1 to Z11) accounted for 83%, 59%, 56%, 54%, and 52% ( Figure A3). Furthermore, the C. microphylla root system's horizontal extension range was far greater than the coverage range of the above-ground canopy, suggesting that the roots tend to explore the largest soil area possible to exploit soil resources and provide anchorage and stability. Figures 7 and 8 show the variation in root distribution density with canopy gradient zones in the root-rich area (Zone 1 and Zone 2) of Plots 1 and 2, respectively. Although  Table 3. The range of the gradient zone for three areas of different depth intervals in Plot 1 and Plot 2. Zone 1 represents the root-rich canopy-covered area, Zone 2 represents the root-rich canopy-free area, and Zone 3 represents the root-poor area.  Figure 6a,b depict the distribution of all of the detected root points in Plot 1 and Plot 2. Combined with the horizontal canopy gradient zone, it is clear that areas close to the plant's taproot collar had denser root point detection, which is in good accordance with the basic principles of root growth. The pattern is more evident in each soil depth interval (as shown in Figure A3). This phenomenon is also shown in Figure 6c,d, i.e., that 90% of the root points were distributed in the first 24 gradient zones and 20 gradient zones for Plot 1 and Plot 2. However, root aggregation varied across intervals, decreasing from shallow to deep soil. The statistical results in Figure 6c,d illustrate this pattern: in Plot 1, the first 30% of the gradient zones (Z1 to Z15) contained 88%, 77%, 70%, 68%, and 63% of the total number of root points in the first through fifth soil depth intervals, respectively. In Plot 2, the first 30% of gradient zones (Z1 to Z11) accounted for 83%, 59%, 56%, 54%, and 52% ( Figure A3). Furthermore, the C. microphylla root system's horizontal extension range was far greater than the coverage range of the above-ground canopy, suggesting that the roots tend to explore the largest soil area possible to exploit soil resources and provide anchorage and stability. Figures 7 and 8 show the variation in root distribution density with canopy gradient zones in the root-rich area (Zone 1 and Zone 2) of Plots 1 and 2, respectively. Although the values of root distribution density in the two plots were different, the trend in root distribution density's variation with the gradient of canopy amplitude was almost identical for any given soil depth interval. For different soil depth intervals, the pattern in root distribution density was different. In the near-surface interval (0-20 cm), the roots grew concentrated in the root-rich canopy-covered area (Zone 1), whereas there were almost no roots in the root-rich canopy-free area (Zone 2). In intervals 2 and 3 (20-40 and 40-60 cm), the root distribution density decreased significantly in Zone 1, whereas it declined gradually in Zone 2, indicating a more gradual dispersal of roots with distance from the root collar. In interval 4 (60-80 cm), the root distribution density slowly decreased from the first gradient zone to the last gradient zone. Unlike the other depth intervals, the distribution of root distribution density in interval 5 (80-100 cm) had no noticeable trend across canopy gradient zones.

Distribution Pattern of Soil Water
The variation of SWC with canopy gradient zones at different intervals in Zone 1 and Zone 2 of Plot 1 and Plot 2 are also shown in Figures 7 and 8. For Plot 1, the mean SWC gradually increased with the canopy gradient zone, except for the surface interval. For Plot 2, the mean SWC slightly decreased from the Z1 to Z5 gradient zones (Zone 1) and then gradually increased from Z6 to Zmax (Zone 2), except for the surface interval. In the surface interval of both plots, SWC continued to decrease as the canopy gradient zone increased.

Correlation between Roots and Soil Water in the Horizontal Direction
In the canopy-covered area, there was no significant correlation between root distribution density and mean SWC in each depth interval of Plot 1 and Plot 2. However, in the canopy-free area, there were significant correlations between the root density and mean SWC in the middle two depth intervals (20-40 cm, 40-60 cm) of the two plots (Table 4), with more than 67% of the root points distributed in this depth range of Plot 1 and more than 69% in Plot 2 ( Figure 4). There was no significance in other intervals (Table 4). Interestingly, there were differences observed between the two plots: Table 4. p-values of the model components and percentage of deviance explained by the generalized additive models (GAMs) for each depth interval of the canopy-free area. The meaning of each model component is described in Equation (6).

Plots & Zones
Depth ( ually in Zone 2, indicating a more gradual dispersal of roots with distance from the root collar. In interval 4 (60-80 cm), the root distribution density slowly decreased from the first gradient zone to the last gradient zone. Unlike the other depth intervals, the distribution of root distribution density in interval 5 (80-100 cm) had no noticeable trend across canopy gradient zones.

Distribution Pattern of Soil Water
The variation of SWC with canopy gradient zones at different intervals in Zone 1 and Zone 2 of Plot 1 and Plot 2 are also shown in Figures 7 and 8. For Plot 1, the mean SWC gradually increased with the canopy gradient zone, except for the surface interval. For Plot 2, the mean SWC slightly decreased from the Z1 to Z5 gradient zones (Zone 1) and For Plot 1, the root distribution density was negatively correlated with mean SWC in the middle two depth intervals (20-40 cm and 40-60 cm), i.e., that as root distribution density increased, the mean SWC decreased (Figure 9). That is to say, with the increase of the gradient zones, the root distribution density decreased, and the mean SWC increased (Figures 7 and 8). than 69% in Plot 2 ( Figure 4). There was no significance in other intervals (Table 4). Interestingly, there were differences observed between the two plots: For Plot 1, the root distribution density was negatively correlated with mean SWC in the middle two depth intervals (20-40 cm and 40-60 cm), i.e., that as root distribution density increased, the mean SWC decreased (Figure 9). That is to say, with the increase of the gradient zones, the root distribution density decreased, and the mean SWC increased (Figures 7 and 8). Figure 9. Estimated smooth functions (the solid lines) and 95% confidence intervals of the root density variables (the dashed lines) for the generalized additive models (GAMs) of the area with significant correlations between root distribution density and mean SWC. The above four images all represent the canopy-free areas. Figure 9. Estimated smooth functions (the solid lines) and 95% confidence intervals of the root density variables (the dashed lines) for the generalized additive models (GAMs) of the area with significant correlations between root distribution density and mean SWC. The above four images all represent the canopy-free areas.
For Plot 2, the relationship between root distribution density and mean SWC showed a non-monotonic trend in the middle two depth intervals (20-40 cm and 40-60 cm), i.e., that as root distribution density increased, the mean SWC decreased first and then increased (Figure 9). That is to say, with the increase of the gradient zones, the root distribution density decreased, and the mean SWC decreased first and then increased (Figures 7 and 8). The percentage of deviance explained by the GAMs ranged from 94.8% to 99.9%. The spatial term eliminated the spatial autocorrelation of the residuals ( Figure A5).

Root-Soil Water Relationship at the Study Site
Previous studies revealed that the root-soil water relationship is a comprehensive expression of a variety of factors, including root uptake [8][9][10][11][12], root interception [21,22,60], preferential flow redistribution [13,14,19], and surface evaporation [22]. This study confirms that the relationship between C. microphylla roots and soil water is spatially heterogeneous, both vertically and horizontally. In the vertical direction, the distribution of roots has a good agreement with mean SWC, and roots were more distributed in the soil interval with higher SWC (Figure 4). This pattern is consistent with the conclusions of previous studies [21,22,60], which suggested that the vertical distribution of plant root biomass was significantly correlated with soil water in arid and semi-arid regions. Soil water in arid and semi-arid areas mainly comes from rainfall [61], and the abundance of root and soil water distribution often show a positive correlation [60]. As a result, in the vertical direction, the shrub keeps the soil water around the root system sufficient by root interception [21,22,60], which improves the adaptation of shrubs to the arid climate.
In the horizontal direction, the root-water relationship in the vicinity of individual shrubs is rather complex, likely regulated by the tradeoff between root uptake and preferential flow-induced SWC redistribution [9,12,22,60]. On the one hand, the stem flow concentrates to lateral roots to form preferential infiltration [15,[17][18][19], resulting in a higher SWC in the first few gradients. The root system could increase soil porosity and thus the frequency and magnitude of preferential flow [62], which transports water to far and deep soil [13,14,19]. Therefore, the soil water gradually decreases with the root distribution density decrease. On the other hand, root water uptake could reduce the SWC in the area of high root distribution density [63,64]. Therefore, such a correlation was negative when root uptake predominated but positive when preferential flow predominated, and not significant when the two major effects were equivalent.
For Plot 1, in the middle two depth intervals (20-40 cm and 40-60 cm) of Zone 2, there were significant negative correlations between root distribution density and mean SWC (Table 4), where the mean SWC showed a trend of decrease with the increase of the root distribution density (Figure 9). According to the conclusions of previous studies, affected by the "fertile islands" of shrubs, the redistribution of preferential flow under the canopy is mainly related to the root distribution density, and the SWC is often higher where the root distribution density is higher [65]. However, in the canopy-free area, the redistribution of preferential flow is weaker, the root water uptake becomes stronger [9,10], and the correlation between root distribution density and SWC was negative [24]. There is no canopy influence in this region, and rainfall replenishes soil water relatively evenly in all intervals [61]. Therefore, in the areas with lower root distribution density, the root water uptake is weaker, and the SWC is relatively higher. On the contrary, in areas with higher root distribution density, the roots uptake is more substantial, resulting in a lower SWC. Therefore, the mean SWC and root distribution density showed opposite trends in this area.
For Plot 2, the root distribution density and mean SWC showed a non-monotonic change in the middle two depth intervals (20-40 cm and 40-60 cm) of Zone 2 (Table 4 and Figure 9). This may be due to the fact that in the area closer to the stem, the replenishment of soil water by preferential flow is stronger than the root water uptake [14,19], whereas, in the area away from the stem, the opposite is true. The difference in the root-water relationship between the two plots may be due to the different soil water conditions. The precipitation during July 2017 and August 2018 in the study area was 81.5 mm and 132.8 mm. Moreover, during 11-13 August 2018, there was 12.5 mm rainfall that occurred just three days before the experiment (from China Meteorological Data Center). Thus, the rainfall enriched the soil water of Plot 2, and the water was blocked by the root system during infiltration, making the highest SWC in the root-rich area in the 20-40 cm depth interval in Plot 2. With the infiltration of water, the soil water in the root-rich area in the 40-60 cm interval was replenished accordingly. Therefore, for Plot 2, the supplementation of soil water by the preferential flow along the roots in the middle two intervals (20-40 cm and 40-60 cm) may have been greater than the uptake of soil water by the roots [14,19], and the root-soil water relationship exhibited a different situation to that in Plot 1. In addition, in the other depth intervals of Plot 1 and Plot 2, the correlations between the root distribution density and mean SWC are not significant, which may be a result of the tradeoff between root water uptake and the supplement of soil water by the preferential flow [9,12,22,60]. The spatial heterogeneity of the root-soil water relationship reflects the survival strategy of shrubs when adapting to the arid and semi-arid environments.

Potentials and Limitations of the GPR Method in the Ecohydrological Study at Field Scales
Previous studies mostly use whole-root excavation or stratified excavation to collect root distribution, which is destructive to the soil environment and limits the study area and the number of root samples [9][10][11][12]. In this study, the spatial distribution information of roots was obtained through GPR non-destructive detection. Adopting the GPR method to study the root-soil water relationship has two advantages: GPR can be used to complete a large area of root detection in a short time [28,31,34,36,39] and can detect the root system in the whole experimental space, which eliminates the errors caused by the bias in sampling [24]. Investigating the root-soil relationship over a large area is often challenging via traditional sampling methods. Moreover, using dual-frequency pairing GPR to obtain the spatial distribution of the underground root system is an emerging method to overcome the low detection accuracy of roots detected by a single frequency [30]. The spatial distribution of roots ( Figure 3) obtained using dual-frequency pairing GPR helps us know more about the developing pattern of the C. microphylla root system on the scale of the shrub population, which facilitates the investigation of the root-soil water relationship.
Further, the canopy gradient zoning is a presentable method for converting the absolute distance between the root point and the taproot collar on the horizontal direction into the canopy gradient zones with ecological significance. This method has been previously used to study spatial patterns and morphological characteristics of plants [24,54]. Using the canopy gradient zones to represent the absolute distance to the shrub's taproot collar has two advantages. Firstly, for shrubs of different canopy sizes, roots within the same gradient zone have a similar physiological significance. As shown in Figure A3, the gradient zones of the same background color represent the similar physiological significance of the root system. Therefore, the influence caused by the difference in shrub size can be avoided in examining root distribution in the horizontal direction. Secondly, this method is helpful to analyze the actual extension range of the crown of underground roots and is of great significance in analyzing the below-ground soil volume occupied and the utilization of soil resources by this species, as well as the further analysis of its survival strategies in arid environments [66]. Specifically, we found that the gradient zones reaching 100% cumulative root points of Plot 1 and Plot 2 were Z40 and Z30, respectively (Figure 6c,d). It can be inferred that the lateral extension length of the root system of the C. microphylla was around 6-8 times the diameter of the above-ground crown in the study area, suggesting that the C. microphylla can absorb limited water over a large area in arid habitats, and maximize the use of soil water, so as to obtain more survival opportunities. However, such data are difficult to obtain using traditional root sampling methods, which led to such important information often being overlooked in previous studies. As the GPR technique has been more frequently used at larger scales for ecohydrological applications, the canopy gradient zoning method used here will help the spatial analysis of massive GPR data.
In this study, we mainly obtained the information of the lateral coarse roots of shrubs without considering the influence of taproots and fine roots on the distribution of SWC. As an essential organ for plants to absorb soil water and nutrients [67,68], fine roots strongly impact the spatial pattern of SWC. However, fine roots grow on the basis of coarse roots, and their spatial distribution is determined by the spatial distribution of coarse roots to a large extent [23,69]. Therefore, we believe that it is reasonable to use the coarse root as a proxy for the entire root system to analyze its relationship to SWC in terms of spatial correlation only. Further, affected by the upper root system, the taproots are difficult to detect by GPR in the deep intervals [70], which is embodied in the root-rich canopy-covered area of the fifth interval, the root distribution density is lower than expected (Figures 7e and 8e), which is a concrete manifestation of the missing detection of the taproot. Last but not least, plant roots are highly dynamic, and different soil environmental factors largely drive the spatial distribution of plant roots in vertical and horizontal directions. In addition to SWC, soil nutrients are also important factors influencing root growth and placement [71][72][73]. Moreover, roots also respond to soil biota and neighboring plants [74,75]. In this study, the measured root spatial distribution was the result of the joint action of various factors.
However, we only analyzed the relationship between SWC and root distribution without accounting for the influence of other factors, thus limiting further interpretation of the root distribution pattern of C. microphylla. Results of this study suggest that current ecohydrological models may not be able to reconstruct the complex relationship between root distribution and SWC distribution. Therefore, we advocate conducting more field studies to disentangle such a complex relationship to boost future models. For example, modeling SWC via root distribution pattern (or vice versa) needs to take into account factors like direction (vertical or horizontal), inside or outside the canopy area, and root density.

Conclusions
A dual-frequency GPR system (operating at 900 and 400 MHz) was used to detect the distribution of the coarse roots of C. microphylla at the stand level in the semi-arid shrubland in China. The spatial variations of root distribution density of the shrub and SWC along canopy gradient zones in different soil depth intervals and their relationships were analyzed. Vertically, the mean SWC and the number of root points in different soil depth intervals of the two plots showed the same trends of increasing at first and then decreasing from shallow to deep soil. Horizontally, the correlations between root distribution density and the mean SWC varied between plots. For the canopy-free area of Plot 1, the root-soil water relationship was significantly negatively correlated in the middle two depth intervals (20-40 cm and 40-60 cm). For the canopy-free area of Plot 2, the root-soil water relationship was non-monotonic in the middle two depth intervals (20-40 cm and 40-60 cm), that is, with the root distribution density increased, the mean SWC decreased first and then increased. The complexity of the root-soil water relationship was likely a result of a combination of factors, such as root water uptake, replenishment of SWC by preferential infiltration along roots, and evapotranspiration during the growing season. More field experiments and observational data are needed to enhance the in-depth understanding of the root-soil water relationship in arid and semi-arid environments. This study improves the application of GPR in the investigation of the root-soil relationship at the stand level, which is often challenging for traditional methods. Results of this study help describe and model the ecohydrological processes in shrub-encroached grassland.

Data Availability Statement:
The GPR raw data and the soil moisture data obtained by soil core sampling will be archived in ZENODO (www.zenodo.org, accessed on 11 March 2021). The source codes of GPR data processing, spatial interpolation, and data visualization are available from the authors upon request.

Appendix A
cessing, spatial interpolation, and data visualization are available from the authors upon request.
Acknowledgments: This study was supported by the National Natural Science Foundation of China (Grant No. 41571404) and the Fundamental Research Funds for the Central Universities at Sichuan University (Grant No. YJ202093. The authors are thankful for the assistance in field data collection from Junxiong Zhou, Qing Li, Liqin Gan, Luyun Zhang, Qixin Liu, Qi Dong, Chishan Zhang, and Wenqing Wang.     Figure A4. The semi-variogram of residuals from the generalized additive models (GAMs) for Plot 2 in the vertical direction. The blue scatter indicates the value of the semi-variance of the residuals, the red "X" is the mean value of the semi-variance of the same interval distance, and the red line is the quadratic fit curve. Figure A5. The semi-variogram of residuals from the generalized additive models (GAMs) for the area with significant correlations between root distribution density and mean SWC. The blue scatter is the value of the semi-variance of the Figure A4. The semi-variogram of residuals from the generalized additive models (GAMs) for Plot 2 in the vertical direction. The blue scatter indicates the value of the semi-variance of the residuals, the red "X" is the mean value of the semi-variance of the same interval distance, and the red line is the quadratic fit curve.
Remote Sens. 2021, 13, x FOR PEER REVIEW 22 of 25 Figure A4. The semi-variogram of residuals from the generalized additive models (GAMs) for Plot 2 in the vertical direction. The blue scatter indicates the value of the semi-variance of the residuals, the red "X" is the mean value of the semi-variance of the same interval distance, and the red line is the quadratic fit curve. Figure A5. The semi-variogram of residuals from the generalized additive models (GAMs) for the area with significant correlations between root distribution density and mean SWC. The blue scatter is the value of the semi-variance of the Figure A5. The semi-variogram of residuals from the generalized additive models (GAMs) for the area with significant correlations between root distribution density and mean SWC. The blue scatter is the value of the semi-variance of the residuals, and the red "X" is the mean value of the semi-variance of the same interval distance, and the red line is the quadratic fit curve.