THREE-DIMENSIONAL INTERPOLATION OF SOIL DATA: FERTILITY AND PEDOMORPHOLOGICAL FEATURES IN SOUTHERN BRAZIL

SUMMARY The graphical representation of spatial soil properties in a digital environment is complex because it requires a conversion of data collected in a discrete form onto a continuous surface. The objective of this study was to apply three-dimension techniques of interpolation and visualization on soil texture and fertility properties and establish relationships with pedogenetic factors and processes in a slope area. The GRASS Geographic Information System was used to generate three-dimensional models and ParaView software to visualize soil volumes. Samples of the A, AB, BA, and B horizons were collected in a regular 122-point grid in an area of 13 ha, in Pinhais, PR, in southern Brazil. Geoprocessing and graphic computing techniques were effective in identifying and delimiting soil volumes of distinct ranges of fertility properties confined within the soil matrix. Both three-dimensional interpolation and the visualization tool facilitated interpretation in a continuous space (volumes) of the cause-effect relationships between soil texture and fertility properties and pedological factors and processes, such as higher clay contents following the drainage lines of the area. The flattest part with more weathered soils (Oxisols) had the highest pH values and lower Al 3+ concentrations. These techniques of data interpolation and visualization have great potential for use in diverse areas of soil science, such as identification of soil volumes occurring side-by-side but that exhibit different physical, chemical, and mineralogical conditions for plant root growth, and monitoring of plumes of organic and inorganic pollutants in soils and sediments, among other applications. The methodological details for interpolation and a three-dimensional view of soil data are presented here. apresentou os maiores valores de pH e as concentrações mais baixas de Al 3+ . Essas técnicas de interpolação e visualização de dados têm grande potencial para ser utilizada em diversas áreas da ciência do solo, como a identificação de volumes de solos que ocorrem lado a lado, mas que oferecem diferentes condições físicas, químicas e mineralógicas para o crescimento das raízes das plantas; o monitoramento de plumas de poluentes orgânicos e inorgânicos em solos e sedimentos; entre outras aplicações. Os detalhes metodológicos para a interpolação e visualização em três dimensões de dados de solo foram apresentados neste trabalho.


SUMMARY
The graphical representation of spatial soil properties in a digital environment is complex because it requires a conversion of data collected in a discrete form onto a continuous surface. The objective of this study was to apply three-dimension techniques of interpolation and visualization on soil texture and fertility properties and establish relationships with pedogenetic factors and processes in a slope area. The GRASS Geographic Information System was used to generate threedimensional models and ParaView software to visualize soil volumes. Samples of the A, AB, BA, and B horizons were collected in a regular 122-point grid in an area of 13 ha, in Pinhais, PR, in southern Brazil. Geoprocessing and graphic computing techniques were effective in identifying and delimiting soil volumes of distinct ranges of fertility properties confined within the soil matrix. Both three-dimensional interpolation and the visualization tool facilitated interpretation in a continuous space (volumes) of the cause-effect relationships between soil texture and fertility properties and pedological factors and processes, such as higher clay contents following the drainage lines of the area. The flattest part with more weathered soils (Oxisols) had the highest pH values and lower Al 3+ concentrations. These techniques of data interpolation and visualization have great potential for use in diverse areas of soil science, such as identification of soil volumes occurring sideby-side but that exhibit different physical, chemical, and mineralogical conditions for plant root growth, and monitoring of plumes of organic and inorganic pollutants in soils and sediments, among other applications. The methodological details for interpolation and a three-dimensional view of soil data are presented here.

INTRODUCTION
Physical and chemical phenomena act on parent material leading to progressive transformations that determine the morphological, physical, chemical, and mineralogical characteristics of the soil formed. Topography is also an important factor in the amount of water retained and the transport of solids or solution materials, producing different soil types (Curi & Franzmeier, 1984). Soils from the same region and parent material may have distinct properties when located at different positions in the landscape and under diverse conditions of moisture, drainage, and relief (Schaefer et al., 2002), also providing them with distinct environmental and agricultural aspects (Reyniers et al., 2006;Chi et al., 2009).
Knowledge of soil spatial distribution, as well as the factors affecting such distribution, is fundamental for natural resource surveys and for the establishment and development of soil management strategies. Some studies have been conducted in soils derived from the same region and parent material as the current study (argillite of the Guabirotuba Formation, Curitiba, Parana, Brazil) aiming at determining the variability of the chemical and mineralogical characteristics in the landscape.
Oliveira Júnior et al. (2014) evaluated different soil sampling units on a slope (122 samples in georeferenced grid with a regular 30-m spacing) by using classical statistics and geoestatistics and concluded that the highest coefficients of variation occurred for Al 3+ , P, and K + (117, 99, and 78 % respectively). The occurrence of drainage lines on slopes led to increases in organic carbon. In the mineralogical data, the total content of crystalline-Fe oxides was the property that showed the least spatial variability in the clay fraction.
The authors cited and several other environmental studies used 2D tools to assess spatial data distribution, following the classical form of geographical representation.
A digital, three-dimensional, and interactive environment has been used in an incipient form in studies involving soil survey data (Grunwald & Barak, 2001;Delarue et al., 2009). Unlike other fields, such as medicine and those that study atmospheric phenomena, which make use of more sophisticated sensors able to collect data continuously, in soil science, interpolation of the data is necessary in order to produce continuous volumes (Neteler, 2001).
Since soil is a three-dimensional body (pedon), its spatial variation extrapolates the limits of the profile, either in horizontal or vertical directions. Analysis of the set of chemical properties can be made based on visualization techniques that allow such variations to be understood under a three-dimensional perspective.
As examples of application of 3D evaluation in environmental studies, Grunwald & Barak (2001) applied Virtual Reality Modeling Language (VRML) to create virtual 3D soil landscape and study relationships between soil horizons and terrain features. Delarue et al. (2009) also used a virtual environment to reconstruct the soil horizons in order to study spatial distribution in the landscape and, although they did not generate volumes, the threedimensional representation helped to interpret the evolution of pedogenetic soil horizons. Khan et al. (2012) estimated soil permeability and mapped the microstructures within the soil aggregates into threedimensional images in order to quantify the connectivity of the pore network, and Hu et al. (2013) developed a process to combine multiple horizontal layers and recreate a profile to view the paths of water flow in three dimensions in calcareous soil.
Three-dimensional data interpolations have also been used for environmental purposes. Ouyang et al. (2002) showed that Pb concentration decreased with depth in sediments of the Cedar River and Ortega River in Florida, USA, in a study to investigate the characteristics and 3D spatial distribution of heavy metals (Pb, Cu, Zn, and Cd).
The objective of this study was to apply threedimension techniques of interpolation and visualization to soil texture and fertility properties and establish their relationships with pedogenetic factors and processes in a slope area.

Field area, soil sampling, and analysis
The experimental plot (13 ha The experimental station where the study was developed has been part of the environmental protection area of the Irai River since 1996. This is a protected area because it is within the bounds of a source of water supply for some cities of the Curitiba Metropolitan Region. Thus, activities involving soil management and farming in general have become quite restricted, the use of fertilizers and pesticides not being allowed. This study was conducted on a hillslope near a site with the presence of cattle before 1996, with the predominance of vegetation for grazing for over 30 years. There is no record of liming and animals in the area after 1996. The hillslope orientation is SW-NE, with altitude ranging from 905 to 940 m. The climate is designated as Cfb (Köppen classification), characterized as humid temperate with well-defined seasons and monthly rainfall ranging from 73 mm in August to 172 mm in January.
The two main soil types were Inceptisols, on the upper and most dissected portion of the landscape (Figure 1c), and Oxisols, on the footslope.
The geological lithostratigraphy included smectiterich argillite followed by sandy arkose deposits and small caliche horizons, typical of the local Guabirotuba geological formation (Salamuni & Stellfeld, 2001). These sedimentary materials are derived from crystalline rocks which weathered under wet conditions and settled during semi-arid periods, over humid-dry cycles.
In accordance with relief and field observations, Oliveira Júnior et al. (2010) identified four distinct geomorphological surfaces in the study area: S01, S02, S03, and S04 at the 940, 935, 930, and 905 m altitudes, respectively. Surfaces S01 to S03 were delimited by small water dividers, which remained as evidence of previous erosion in the form of small crests, with convex divergent relief. In contrast, S04 represented a flat area which received the colluvial sediments ( Figure 1c).
Soil was sampled following a regular 122-point grid, with 30 m spacing between points (Figure 1b). An undisturbed core-sample (cylinder) was taken from the 0 to 110 cm depth (10 cm diameter) at each sampling point using a PVC tube enclosed in a mechanical auger apparatus driven by a tractor power take-off.
After being removed, the PVC tube was cut lengthwise, and the sample was divided into two monoliths, where the A, AB, BA, and B horizons where identified and morphologically described. Samples of the four horizons were separated, dried at 40 o C for 24 h, passed through a 2 mm mesh, and stored for further analysis. A total of 488 soil samples were used in the study (122 sampling points × 4 horizons).

Volume modeling of soil and visualization
Soil fertility data were organized on a spreadsheet according to identification of the sampling point, the UTM geographic coordinates (x, y), and sampling depth (z). The z parameter of each horizon had to be represented as a single value (discrete position), so it was considered as the depth of the upper limit of the horizon. Data was exported into the ASCII (American standard code for information interchange) format.
To generate three-dimensional soil models, we used the GRASS GIS 6.3.0 software (Geographic Resources Analysis Support System) (GRASS Development Team, 2012). Before data import, a georeferenced environment was created in order to limit processing to only within the sampling site or bounding box (Figure 1b,c). This environment was set in a UTM coordinate system referred to as SAD 69 Datum, and the central meridian of 51 o W.G.
After data import, a mask (raster map) was created to cover and limit processing to only the sampling site. Soil fertility data were interpolated and represented in three-dimensional space through the regularized spline with tension (RST) method (Mitasova & Hofierka, 1993;Mitasova & Mitas, 1993). The spatial resolution was 2 m and, by default, the system adopts a tension of 40 and smoothness of 0.1.
After interpolation, the resulting representation of soil volume was exported into the VTK (visualization toolkit) format for viewing with ParaView 3.6.0 software (Henderson, 2001). This procedure required the use of a mask for the 3D volume, which ensured that only data within the perimeter of the area was exported. This was done by transforming the mask raster map from a 2D to 3D mask.
Only three display filters were used in Paraview: contour, slice, and threshold. Contour filters data through isovalues determined by the user; slice is a filter that allows cuts into the volume, enabling visualization inside the volume as cutting profiles; threshold is a filter that extracts scalar values defined by the upper and lower limits specified by the user.
The information contained in volumes was obtained through the surface extraction procedure. Each surface was interactively determined according to the values of interest for each property, e.g., pH above 5.4, with is considered the lower limit for precipitation of exchangeable Al 3+ (Pionke & Corey, 1967).

Soil texture
Heavy clay soils (>600 g kg -1 of clay) occurred in the lower and flat portion of the landscape, between the S03 and S04 levels, where Oxisols predominate, especially associated with drainage lines located between erosion planes (Figure 2). Isolated soil volumes with clay content >700 g kg -1 were also found in this portion (red volumes, Figure 3). Sand content, for its part, was higher mainly on ridges between S01 and S03, in accord with lower clay contents.
The higher clay content (>600 g kg -1 ) in the lower portion of the area, where Oxisols predominate (Figure 3), may be associated with a preferential deposition of fine colluvial material and with a more advanced weathering stage of the subsurface horizons.
The increase in clay in drainage lines (Figure 2) may be attributed to two processes, which are not mutually exclusive: i) dissection of the area, in which soil removal by erosion may have exposed deeper horizons closer to the parent material (rich in clay), mainly in the valley bottoms of the drainage lines, as can be seen by the differential accumulation of clay in these positions in the longitudinal cross section ( Figure 2b); ii) lateral migration of clay, as seen by the increase in clay content towards the drainage lines (transversal cross section, Figure 2b). Supporting this evidence, Oliveira Júnior et al. (2014), in a previous study, found lower water content, divergent convex relief forms, and higher index values of soil loss on the ridges of the terrain, defined as S01, S02, and S03 (Figure 1c).
Soil acidity: pH, exchangeable Al 3+ , and nonexchangeable H All soil samples had an acidic character, with pH ranging from 4.8 to 6.0 (Figure 4a), with 5.2 being the modal pH value (Figure 4b). There was no coincidence in the spatial distribution of pH and Al 3+ ( Figure 5). The highest pH values were observed up to the 20-cm depth, with the occurrence of some isolated volumes with pH > 5.4 in the soil subsurface ( Figure 5). The highest Al 3+ concentrations were observed from the 20 to 80 cm depth, with some isolated soil volumes > 4 cmol c kg -1 (red, Figure 5).
The main Al 3+ source in these soils was possibly smectite, a major constituent of the argillite parent material (Guabirotuba Formation) (Salamuni & Stellfeld, 2001). Melo et al. (2009) studied the mineralogical composition of the clay fraction of Inceptisols derived from this same argillite and found montmorillonite and beidellite (dioctahedral smectites with Al in the octahedral sheet) in the C2, C3, and C4 horizons, up to 4.5 m deep. Through diagenesis and neogenesis processes during weathering, smectite is transformed preferentially into 1:1 phyllosilicate (kaolinite) in subsurface horizons, with intense release of Al 3+ (Juo & KampratH, 1979;Thomas & Hargrove, 1984).
In Oxisols of the footslope, the occurrence of soil volumes with pH < 5.4 or Al 3+ > 2.8 cmol c kg -1 were not as expressive as in the upper land ( Figure 5). The low Al 3+ concentrations in the footslope are attributed to more intense weathering, characterized by silica removal and the formation of gibbsite (Oliveira Júnior et al., 2014), and higher pH values. The low Al 3+ concentration in highly weathered soils (Fox, 1982) has positive practical implications in terms of less lime requirements for crop production.
The lower Al 3+ concentrations in the top 0-20 cm layer are attributed to three simultaneous factors: (i) the horizon is more weathered, with Al 3+ being transformed into gibbsite (Motta & Melo, 2009); (ii) higher contents of organic matter, which complexes Al 3+ (Simas et al., 2007); and (iii) higher pH (>5.4), inducing Al precipitation (Pionke & Corey, 1967). In the top layer, there are extensive soil volumes with pH >5.4 (Figures 4 and 5) and OC contents > 25 g kg -1 (Figure 6) that do not coincide with surfaces that  Figure 3. Volume representation of the spatial distribution of clay content from 600 to 700 g kg -1 (light green), and >700 g kg -1 (red). Inceptisols occur in the upper portion of the landscape, and Oxisols in the lower.   Figure 6. Volume representation of the spatial distribution of organic carbon >25 g kg -1 (blue), Al 3+ between 2.8 and 3.9 cmol c kg -1 (yellow) and >4.0 cmol c kg -1 (red) (a); and of Al 3+ >4.0 cmol c kg -1 (red) and the horizontal distribution of the clay fraction CEC at 80-cm depth (b). Inceptisols occur in the upper portion of the landscape, and Oxisols in the lower. represent the highest Al 3+ concentrations. According to Conyers (1990) and Rodeja et al. (2004), the extraction of Al 3+ specifically bound to functional groups of the organic matter (inner sphere complex) requires molar solutions with a stronger extraction capacity than KCl, like CuCl 2 .
There was a predominance of soil volumes with Al 3+ > 4.0 cmol c kg -1 in drainage lines of the upper land (Inceptisols) (Figure 6b). The subsurface layer (80 cm depth) in these drainage lines accumulated Si from lateral and vertical flows, which led to favorable conditions for maintaining the smectite derived from the parent material (argillite), and thus the formation of less developed soils. The higher cation exchange capacity (CEC) of the clay fraction in drainage lines (Figure 6b) suggests the occurrence of 2:1 expansive minerals. Therefore, smectite would buffer the higher Al 3+ concentrations in these parts of the landscape (Turner & Brydon, 1967;Pratt et al., 1969;Thomas & Hargrove, 1984;Motta & Melo, 2009).
There was a close association between OC and H contents (Figure 7a). In the upper parts (better drained areas and lower OC content), the volume of lower H concentrations extended up to the surface soil. The functional groups of organic matter and pH-dependent variable charges of the clay fraction minerals (kaolinite and Fe-and Al-oxides) are main sources of covalent H -dissociation of carboxylic and phenol groups from humic substances (Dick et al., 2002;Novotny et al., 2006), and of silanol (-SiOH), aluminol (-AlOH) and the iron group (-FeOH) from clay fraction minerals (Tarì et al., 1999).
Soil cation exchange capacity (CEC) at pH 7.0 Higher soil CEC was observed near the soil surface and landscape portions close to the drainage lines (Figure 7b), which coincides with portions of higher OC contents (Figure 7c). Pei et al. (2010) and Hancock et al. (2010) found higher organic matter accumulation in areas with a higher wetness index, a factor that delays the decomposition process. Furthermore, Moore et al. (1993), using the wetness index and slope gradient, found high correlations between these topographic indices and organic matter content, sand content, and A horizon depth.
At greater depth, with little influence from organic matter, an expressive effect of the quality of the clay fraction was observed -higher soil CEC at pH 7.0 (figure not shown) and higher clay fraction activity (clay CEC) in drainage lines (plane at 80 cm depth, Figure 6b).

Ca 2+ and Mg 2+ contents
The highest Ca 2+ and Mg 2+ concentrations, both in the surface and depth layers (Figure 8), occurred in portions of higher organic matter accumulation and higher soil CEC at pH 7.0 values in drainage lines (Figure 7b,c). Possibly, these exchangeable bases were leached from higher and convex portions towards lower portions of higher adsorption capacity. In spite of the same trend in the surface, the Ca 2+ and Mg 2+ concentrations were lower at the 100 cm depth (Figure 8).
Due to the higher adsorption strength of Ca 2+ compared to Mg 2+ , the Ca:Mg ratio in the surface soil, at the higher frequency of each one (4.0 cmol c kg -1 /2.1 cmol c kg -1 ), was higher than that observed at 100 cm (1.3 cmol c kg -1 /1.0 cmol c kg -1 ) (Figure 8), suggesting a greater leaching and accumulation of Mg 2+ in the subsurface.
The higher Mg 2+ concentrations at the end of the main drainage route at the 100 cm depth, which divide the region of Oxisols from the regions of Inceptisols (Figure 8g) could also be attributed to preferential leaching of this cation in relation to Ca 2+ . The volumes with higher Ca 2+ +Mg 2+ concentrations and higher base saturation (V) values did not coincide with the higher Al 3+ concentrations (Figure 9a) or the higher Al 3+ saturation values (Figure 9b), respectively. This trend opposes the common generalization that highly weathered soils of the tropics are very acid (Lopes & Cox, 1977). However, the effect of the form of relief on these relationships should also be considered (Kämpf & Curi, 2012), where the higher contents of Ca 2+ + Mg 2+ tend to follow the drainage lines of the area (Figure 8).
For the purpose of agricultural use, the flat lowland area with greater occurrence of more weathered soils (basal view-Oxisols) would require less lime. The ringshaped volume with higher Al 3+ concentrations (yellow) was similar to the volume with higher Al 3+ saturation (m) (red) within the Ca 2+ +Mg 2+ concentrations and base saturation surfaces, respectively (Figure 9a,b).
Three-dimensional distribution allowed visualization of volumes with better soil chemical conditions, as seen in the basal view of the distribution of base saturation (V) (Figure 9b). Near the most central line of the area (limit between Oxisols and Inceptisols), there are better growing conditions due to the greater projection of the volume of higher base saturation values (23 to 68 %) towards the deeper horizons. Nevertheless, even in this region, plants growing side-by-side can find different fertility conditions due to the variation in the volumes of base saturation values.
The use of geoprocessing and graphic computing enabled visualization of the distribution of fertility properties in three-dimensional space, allowing observation of the extension of internal volumes with differentiated values/concentrations confined within the soil matrix. These techniques have great potential for more extensive use in diverse areas of soil science,   Figure 8. Horizontal distribution of Ca 2+ concentration in the surface soil (a) and its frequency distribution (b); and at the 100-cm depth (c) and its frequency distribution (d); and distribution of Mg 2+ concentration in the surface soil (e) and its frequency distribution (f); and at the 100-cm depth (g) and its frequency distribution (h). Inceptisols occur in the upper portion of the landscape, and Oxisols in the lower position.
for example, identification of soil volumes occurring side-by-side but which offer different conditions for plant root growth; and monitoring of plumes of organic and inorganic pollutants in soils and sediments.

CONCLUSIONS
1. The 3D tools enabled visualization and interpretation in a continuous space (volumes) of the cause-effect relations between chemical properties and pedogenetic factors and processes.
2. Higher clay content in drainage lines and in flat portions of the terrain, and higher sand content on ridge positions are attributed to internal clay migration and to erosion processes.
3. The preferential flow of water towards the drainage lines favored the accumulation of organic matter in the surface layer and of more clay activity in the subsurface layer (confined volumes of soils with higher clay CEC).  . Volume representation of the spatial distribution of the highest Ca 2+ +Mg 2+ concentrations (2.1 to 10.8 cmol c kg -1 ) (orange) and of the highest Al 3+ (>2.8 cmol c kg -1 ) concentrations (yellow) (a); and of base saturation (V) from 23% to the highest value of 68% (green), and of Al 3+ saturation (m) from 60% to the highest value of 83% (red) (b). Inceptisols occur in the upper portion of the landscape, and Oxisols in the lower.