Multifractal analysis of soil penetration resistance under sugarcane cultivation

: Soil resistance to penetration (PR) is an indirect measure of the state of soil compaction. Thus, the objective of this study was to characterize PR in vertical profiles in an area cultivated with sugarcane using multifractal models for different relief units. The experiment was carried out in an Oxisol with a clay texture, with 6.85 ha in the municipality of Coelho Neto (Maranhão state, Brazil), where 60 sampling points were demarcated. The area was divided into four relief units (Type A > 74 m, Type B from 71 to 74 m, Type C from 68 to 71 m and Type D from 65 to 68 m). The PR was measured at the 60 sampling points using an impact penetrometer, and the PR determined in the 0-0.60 m depth layer every 0.01 m. The multifractal analysis was performed considering the scale property of each profile and typified the singularity and Rènyi spectra estimated using the current method. Multifractal analysis allowed the identification of patterns at different scales and with high heterogeneity. The multifractal behavior was represented by the singularity spectrum (α), versus f(α), and the generalized dimension (Dq). The multifractal analysis allowed the differentiation between the profiles of the relief units (Types A, B, C and D), resulting in an important tool for studies of soil resistance to penetration.


Introduction
Fractal and multifractal models can be used to elucidate the complexity of processes that determine soil characteristics (Halsey et al., 1986;Caniego et al., 2006). According to Evertsz & Mandelbrot (1992), a monofractal system is characterized by a simple scale, while a multifractal system involves a continuous spectrum of scale exponents which are related to complex processes at different levels of intensity. According to Banerjee et al. (2011), a multifractal analysis can reveal more information about data heterogeneity than other spatial techniques.
In this sense, research on soil compaction should consider mathematical models that take soil heterogeneity into account, even at small scales. The degree of compaction of soil can be characterized by measuring its resistance to penetration based on its dynamic nature and its interaction with other factors such as soil composition, structure, bulk density, porosity and water content (Lipiec & Hatano, 2003;Siqueira et al., 2015;Lima et al., 2017;Tavares et al., 2017). Folorunso et al. (1994) examined the soil penetration resistance (PR) in parallel transects and found that fractal analysis allowed a qualitative discrimination of soils, despite small variations across scales. Roisin (2007) found that soil management affected PR multifractality. Siqueira et al. (2013) used multifractal and geostatistic techniques to study vertical profiles of PR and concluded that PR is variable on a vertical and horizontal scale. Studying the effects of soil bulk density and water potential on PR multifractality, Paz-Ferreiro et al. (2013) found that increasing bulk density caused a decrease in PR multifractality. Wilson et al. (2016) used multifractal analysis techniques to study PR due to the increase of water deficit and observed that PR values become increasingly heterogeneous with decreasing soil water content.
Thus, the aims of this study were a) to analyze the vertical PR profiles under sugarcane cultivation by means of multifractal techniques, b) to test whether the multifractality of vertical PR profiles was influenced by the relief, and c) to characterize the spatial variability of PR and multifractal parameters.

Material and Methods
The study area was conducted at the Itajubara Mill in the municipality of Coelho Neto, state of Maranhão, Brazil (4° 19' 5'' S, 43° 0' 22'' W, at an altitude of 71 m). The experimental area covered 6.85 ha and had been cultivated with sugarcane (Saccharum officinarum), variety RB92579, since 2014. The soil of the area was classified as an Oxisol.
In total, 60 sampling points were randomly selected to determine soil penetration resistance (PR) (MPa), and they were classified into four relief units ( Figure 1): type A, > 74 m; type B, from 71 to 74 m; type C, from 68 to 71 m; and type D, from 65 to 68 m.
Soil PR was measured using an impact penetrometer with a 0.0125-m diameter steel cone, at a 30° angle. PR measurements were performed on April 22, 2016 with soil at field capacity (θ = 0.14 m 3 m -3 ), which was determined in the laboratory by placing the samples in soil retaining rings in a Richards pressure plate device. The PR values were determined in the 0-0.60 m depth layer at 0.01 m intervals, and 60 vertical PR profiles were generated.
Descriptive statistics were obtained using R software version 3.3.1 (R Development Core Team, 2017). The multifractal analysis was performed using the method of moments to obtain multiple PR scales, considering a one-dimensional profile of length L, this involves successive partitioning into k stages (k = 1, 2, 3 …), and successive partitioning was produced at each scale δ, a number of segments N(δ) = 2 k are obtained with a characteristic size length, δ = L × 2 k , covering the whole extent of the support L (Evertsz & Mandelbrot, 1992;Caniego et al., 2006;Vidal-Vázquez et al., 2013).
The expression of the data as a function of mass p i (δ), at time resolution δ was estimated according to Eq. 1 and the size of the segments was calculated using Eq. 2: where p i (δ) -probability mass distribution; N i (δ) -value of the measure in a given segment I; N t -sum of the measure in the whole transect; N(δ) -number of segments with size δ. The statistical moment q is defined as -∞ < q < ∞; and, χ(q,δ) -scale of partition function, in the moments (q) and generated segments of order (δ).
To estimate the exponent of mass function, τ q , with the mass scaling function of order q, using Eq. 3, and the distribution probability was given as p i (δ) = δ αi , where α i is the singularity or Hölder exponent characterizing density in the i th box (Halsey et al., 1986). The Hölder exponent was calculated as α i = log μ i (δ)/log δ, which can be used to evaluate the degree of concentration of the measure μ.
where τ(q) -mass exponent function; μ i -region of concentration of a measure in the Hölder exponent; and, δ -generated segment length.
The relation between the exponent τ q and f(α) can be obtained using the Legendre transformation. However, the Legendre transformation is unsatisfactory in some situations, mainly due to errors in the estimation of f(α). Thus, in the present study, τ q and f(α) were estimated using the direct method of Chhabra & Jensen (1989) using Eqs. 4 and 5.
where ∝ -is proportional to.
The generalized dimension was defined according to Eq. 6, where D 1 is indeterminate because the denominator value is zero. Therefore, for the particular case in which q = 1 was used, according to Eq. 7, considering a stationary stochastic process defined by c = (c i : i = 0,1,2, ...), with constant mean μ and finite variance α 2 .
In a monofractal, D q is a constant function of q, and in multifractals, the relationship between D q and q is not constant. The Hurst exponent (H) was determined using the expression D 2 = 2H -1.
Geostatistical analyses were performed considering the intrinsic hypothesis and after adjustment of the semivariograms, obtaining the following parameters: C 0 (nugget effect), C 0 +C 1 (sill) and a (range, m). The software Surfer 11 (Golden Software, 2014) was used to construct the isoline maps of soil PR and multifractal parameters.

Results and Discussion
The vertical PR profiles of different relief units are shown in Figure 2 Tavares et al. (2017) suggested that PR is influenced by soil clay content and bulk density, which is associated with the volumetric content of water. Therefore, the greater PR homogeneity over depth in the type D relief unit is likely associated with the dynamics of water in soil since this relief unit had the lowest elevation (between 65 and 68 m). The mean and standard deviation of PR in the vertical profiles of the 0-0.60 m layer are shown in Figure 2E. The increase in standard deviation with depth is because soil management under sugarcane cultivation typically results in homogenization of the surface layer, as reported by Tavares et al. (2017). The plot shows that the lowest PR (4) , log 1 value was 0.594 MPa, and the highest value was 4.277 MPa. Several studies reported that PR values above 2 MPa limit the development of sugarcane crop (Otto et al., 2011;Sá et al., 2016), and other studies suggested different threshold values to limit the development of most crops, e.g. 2.5 MPa (Gonçalves et al., 2014), 3.0 MPa (Souza et al., 2015), or ranging from 2.8 to 3.2 MPa (Vepraskas & Miner, 1986). The lowest mean PR occurred in profile 18 (1.79 MPa) and the maximum value was observed in profile 13 (2.67 MPa) in types C and B relief units, respectively. Sá et al. (2016) studied soil PR limiting the development of sugarcane and reported values of 3.8 MPa as limiting for crop growth; in the present study, however, the PR values in the representative profiles of each relief unit (Table 1; Figure 2) were below 3.8 MPa. Lipiec & Hatano (2003) suggested that PR values above 4 MPa prevented root development in sugarcane. However, Siqueira et al. (2013) found that under favorable moisture conditions in the field, root development was not compromised, even at high PR values.
The partition functions were constructed for successive segments of different sizes, i.e., 2 k , with k = 0 to k = 5, and moments in the interval q = + 5 to q = -5. The partition functions of the profiles 22 (type C) and 60 (type D) and the graphs that presented the best and worst coefficient of determination (R 2 ) which was R 2 = 0.99993 and R 2 = 0.99615, respectively are shown in Figure 3 (A and B). According to Siqueira et al. (2018), R 2 above 0.9 indicates strong scales for the moments of order q [-5, 5]. Thus, the data followed a power law characterizing multifractal scales, as described by Banerjee et al. (2011). Vidal-Vázquez et al. (2013) reported that the presence of multiple scales indicates a multifractal, whereas single scales characterize a monofractal.
The multifractal dimensions for PR in the dimensions (D -5 -D 5 , D -5 , D 5 , D 0 , D 1 and D 2 ) is presented in Table 2. The in the type B relief unit were more heterogeneous compared to the profiles of the relief unit types A, C, and D ( Table 2). The multifractal description of PR data on the type B relief unit is important information, as this relief unit represents 3.27 ha (47.73%) of the study area, thus this information can be used for the delimitation of specific management areas using multifractal parameters.
The values of α and f(α) for the singularity spectra of the four patterns were calculated for the moment q, with R 2 values above 0.90. For the four relief units (types A, B, C and D), q + was constant (-5) and qranged from -2 to -5 among the profiles ( Table 2).
The singularity spectra for the PR profiles (types A, B, C and D) are shown in Figure 3D. The singularity spectra had a parabola shape with variation in the skewness of the branches. According to , Paz-Ferreiro et al. (2018) and Siqueira et al. (2018), the skewness of spectra indicates the dominance of outliers in the data. High coefficients of variation were observed in the present study (Table 1). The singularity spectra were broad and stretched to the right for all analyzed profiles ( Figure 3D), albeit with small variations in type D. The width and amplitude (α max -α min ) of the singularity spectra of PR profiles varied between the different relief units ( Table 2).
Spatial distribution maps and models of the semivariogram for the mean and maximum PR at 0-0.60 m depth, the parameters PR D-5-D5 , PR α0-α5 and the Hurst exponent (H) produced by ordinary kriging are shown in Figure 4. The best fit to the data set was observed in the exponential model ( Figure  4), apart from PR maximum that was adjusted to the Gaussian model, which may be explained by its range value (a = 32 m). This was below PR mean (a = 38.7 m), indicating that PR maximum values were more heterogeneous in the study area, whereas PR mean values varied on a smaller scale in the study area. The parameters PR D -5-D5 , PR α0-α5 , and H showed range values of 41, 31.7, and 41 m, respectively.
The spatial maps for the mean and maximum PR suggest an inverse relationship between the parameters PR D-5-D5 and PR α0-α5 . In this sense, linear correlations were calculated for all parameters involved in spatial analysis in order to confirm the spatial pattern of the contour lines. Positive correlations of PR mean and PR maximun (|r| = 0.565), PR D-5-D5 and PR α0-α5 (|r| = 0.486), and of PR maximun and PR α0-α5 (|r| = 0.333) were observed. Negative correlations were observed in PR D-5-D5 and H (|r| = -0.498), PR α0-α5 and H (|r| = -0.183), and in PR mean and H (|r| = -0.102). The negative correlations of H and PR mean observed in the present study corroborate the results of Siqueira et al. (2013) who examined the multifractality of vertical profiles of PR and found a negative correlation of H and PR mean (|r| = -0.316).  (D -5 -D 5 , D -5 , D 5 , D 0 , D 1 and D 2 ) and the singularity spectra (q + , q -, α 0 , α max and α min ) for profile of penetration resistance (PR) *± Standard error of mean The maps of spatial variability for mean and maximum PR at depths of 0-0.60 m ( Figures 4A and B) showed a pattern of distribution of contour lines similar to that of the topographic map (Figure 1), which may explain higher averages and higher variability in the type B relief unit (Table 1), compared with the other relief units. As mentioned above, the type B relief unit showed the highest multifractality of the data ( Figure 3C; Table 2) and the highest values of Dq (D -5 -D 5 ). This suggests that management zones can be delimited based on relief units. The spatial variability map of PR D-5-D5 ( Figure 4C) thus allows the distinction of two management zones, one in the upper part of the area comprising types A and B, and the other in the lower part including types C and D.
The map of PR α0-α5 ( Figure 4D) facilitates the delimitation of two specific management zones, following the same pattern observed in the PR D-5-D5 map ( Figure 4C), with two specific management zones, which mainly follow the contour in the PR maximum map ( Figure 4B). This is in line with the results of Siqueira et al. (2013) who assessed multifractality of vertical PR profiles and observed the same pattern of distribution. The relation of the spatial variability maps of PR maximum and PR α0-α5 ( Figures 4B and D) confirms the values of the observed positive linear correlation of these parameters (|r| = 0.333). Using the PR maximum modeled by means of geostatistical tools and multifractal analysis allows inferences on the state of soil compaction. In this sense, it is verified that most of the study area showed PR maximum values of approximately 3.6 MPa and showed no problems due to soil compaction, which is in line with the studies of Lipiec & Hatano (2003) The spatial pattern of H in the study area was above 0.91 ( Figure 4E). The values of H describe the persistence or anti-persistence of a series following a law of scale. Thus, the spatial variability map of H indicates the multifractality and non-randomness of the vertical PR profiles, suggesting that independent of the relief unit the PR measured in the 0-0.60 m depth layer is autocorrelated.

Conclusions
1. The vertical profiles of PR in the 0-0.6 m depth layer at 0.01 m intervals was multifractal as indicated by the singularity spectrum (α) versus f(α) and by the generalized dimension (Dq).
2. The vertical PR profiles in the type B relief unit showed greater heterogeneity in the Dq (D -5 -D 5 ) dimension compared to the other relief units (types A, C and D), demonstrating greater multifractality of the PR series.
3. The use of geostatistics and multifractal analysis allowed the identification of zones for homogeneous management of PR, considering variation in the relief. PR maximum and PR D-5-D5 allowed the delimitation of two management zones.