INTERPOLATION AND 3 D VISUALIZATION OF SOIL MOISTURE

Adaptation to climate change demands the optimal and sustainable water management in agriculture, with an inevitable focus on soil moisture conditions. In the current study we developed an ArcGIS 10.4. platform-based application (software) to model spatial and temporal changes in soil moisture in a soy field. Six SENTEK Drill & Drop soil moisture sensors were deployed in an experimental field of 4.3 hectares by the contribution of Elcom Ltd. Soil moisture measurement at each location were taken at six depths (5, 15, 25, 35, 45 and 55 cm) in 60-minute intervals. The model is capable to spatially interpolate monitored soil moisture using the technique. The time sequence change of soil moistures can be tracked by a Time Slider for both the 2D and 3D visualization. Soil moisture temporal changes can be visualized in either daily or hourly time intervals, and can be shown as a motion figure. Horizon average, maximum and minimum values of soil moisture data can be identified with the builtin tool of ArcGIS. Soil moisture spatial distribution can be obtained and plotted at any cross sections, whereas an alarm function has also been developed for tension values of 250, 1,000 and 1,500 kPa.


Introduction
Preparation for climate change demands optimal and sustainable water supplies for cost-efficient and sustainable agricultural productivity (Makó et al. 2010), which requires the thorough understanding and sound knowledge on vadose zone soil moisture conditions.Soil moisture is a highly variable environmental parameter, both spatially and temporally (Loew -Schlenz 2011), and knowledge on its spatial heterogeneity may provide useful information for precision farming and costefficient crop production (Paul -Speckmann 2004).
Soil moisture spatial distribution varies both vertically and horizontally in a small scale largely due to topography (Anderson -Kneale 1980;Zhu -Lin 2011), soil texture, soil organic matter content and vegetation (Novák 2005;Novák et al. 2013).Principal options to obtain sufficient knowledge on spatial distribution of soil moisture include high resolution ground monitoring (Zhu -Lin 2011;Bárdossy -Lehmann 1998) or satellite remote sensing applications (Jia et al. 2013;Penna et al. 2009;Mohanty -Skaggs 2001).A relatively few papers are available on in situ measurement with sensors, whereas a large number of studies tested remote sensing data on large-scale soil moisture products.One drawback of the remote technique is that it only obtains data from the very shallow subsurface in limited depth (max.: 15-20 cm) (Brocca et al. 2009;Heathman et al. 2012) and is more applicable for larger area then on plot-scale or small watersheds (Hegedüs et al. 2015).Other factors that limit radar technology measurement comprise bulk density, surface roughness and high vegetation density.Radar soil moisture measurements may also be biased due to frequency shifts and loss in power of a radar signal (Paul -Speckmann 2004).Precision can be further reduced by background signal disturbances (spreading, attenuation, reflection and scattering) and flight density (Zhou et al. 2016).
Nevertheless, the use of ground soil moisture measurements, at least in high spatial density, is costly and labor-intensive.If the correlation is sufficiently robust, a single soil moisture sensor could provide the necessary input to predict soil moisture across the watershed (Wang -Qu 2009;Venkatesh et al. 2011;Cosh et al. 2004), assuming that the environmental variables (i.e.land use and soil characteristics) remain unchanged (Lacava et al. 2010).Such predicted datasets within the watershed may provide sufficiently accurate soil moisture (Fu et al. 2003;Loew -Schlenz 2011;Bárdossy -Lehmann 1998).When multiple sensors are needed for the studied area, to overcome the challenge of cost-efficiency, measured data should be interpolated among the measurement points through correlation functions with various GIS applications (Gravalos et al. 2013).A large number of interpolations techniques, methods and functions have been used to estimate soil moisture among the measurement points (e.g.: Yao et al. 2013).Yao et al. (2013) tested four different interpolation methods, namely the ordinary kriging, inverse distance weighting, linear regression and regression kriging.According to their results the latter performed the best.Perry and Niemann (2008) used empirical orthogonal functions to estimate soil moisture among measurement points.
Nonetheless, most interpolation models focused on 2-dimensional distribution of soil moisture, and only a limited number of publications are available on 3D-modelling, most commonly using Hydrus-3D, which, in most cases performed well (Honar et al. 2011).Gravalos et al. (2013) developed a soil moisture interpolation model under laboratory circumstances.In their soil tank they measured soil moisture in a depth of 15 cm.Although their model performed well at laboratory-scale, data interpretation at field scale is challenging.
Our objective was to develop a relatively simple and cost-efficient 3-dimensional soil moisture visualization software which is suitable to present the moisture dynamics of multiple soil layers and horizons in real time.We aimed to indicate the temporal changes of soil moisture in each soil master horizons and project and interpolate data at plot scale (few hectares).We also planned to create a plotting function that is capable to show and visualize temporal changes of soil moisture in any arbitrarily selected point is space (including interpolated points).The result of the research is the property of Elcom Ltd.

Site description
The study area (Koplaló Plot) is located on the southern foothills of the SW part of the Mecsek Hills (Jakab-hegy) in SW Hungary (Fig. 1.).The plot is owned by B-Aranykorona Ltd., and was fallowed from mid-July to mid-September, when rape was planted in the plot.The study site covers a land area of 4.3 hectares and is characterized by Ramann-type brown forest soils (WRB: Stagnic Cambisol).
The elevation ranges between 171.85 and 182.76 meters a.s.l. with a maximum slope of 5.41% on the southern, south-western steep, a former foothill colluvial deposits (Fig. 2.).The long-term  average annual precipitation total is 700 mm with most rainfall measured between May to mid-July (295 mm), when the precipitation regime is characterized by high-intensity convective rainfall events.Based on our field measurements, using a TOPCON HiPer Pro (Topcon Positioning Systems, Inc., Tokyo, Japan) RTK GPS device, a digital elevation model (DEM) of a meter horizontal resolution was developed in ArcGIS 10.4.software environment.

Soil moisture monitoring and soil sampling
In July, 2016 a soil profile of 150 cm depth was excavated to identify master horizons and for soil sampling purposes about 30 meters north of the northern borderline of the studied plot.Soil samples were taken in three repetitions form each master horizon (Ap, A, B, BC and C) with Veér-type brass rings with a volume of 100 cm 3 and a diameter of 48 mm at depths of 15, 25, 50 and 115 cm to determine bulk densities and soil moisture with the gravimetric method according to Flint -Flint (2002).Four Decagon 10HS sensors (Decagon Devices Inc., Pullman, WA, United States) and four MPS-2 in one profile at the norther edge of the experimental plot (Fig. 3.) were inserted into depths of 15, 25, 50 and 115 cm.Sensors were deployed on July 18, 2017 and data have been collected since then continuously in 30-minute intervals with a Decagon EM50 datalogger to obtain the water retention curves for the aforementioned soil depths Soils of the studied plot were sampled with an AMS hand sampling tool at 35 locations to a depth of maximum 290 cm and at one point, at the northern border of the plot a soil profile was excavated to a depth of 120 cm to identify the major soil horizon and to insert the Decagon 5TM and MPS-2 sensors.The soil sampling size were first positioned according to a standardized grid method with equal distances (2×8 drilling holes in North-South direction), second, at the margins (2×4 holes) in order to get an interpolation surface and then 11 randomly located intermediate points were sampled.The depth and thickness of each master soil horizons (Ap, A, B, BC and C) were determined at all sampling sites.Soil samples were taken at each sampling site from all master horizons, therefore altogether 175 soil samples were collected.
Six SENTEK Drill & Drop (Campbell scientific, Edmonton, Alberta, Canada) soil moisture sensors were deployed in an experimental field of 4.3 hectares.Soil moisture measurement at each location were taken at six depths (5, 15, 25, 35, 45 and 55 cm) in 60-minute intervals.

Laboratory analyses
To determine soil textural classes, the liquid limit according to Arany were determined according to Buzás (1993).Soil texture was also determined with a dynamic light scattering method using a Malvern MasterSizer 3000 (Malvern Inc.Malvern, England, UK) particle size analyzer.Samples were pretreated with 10% HCl for CaCO 3 and with 10% H 2 O 2 for organic matter removal.
Walkley-Black titration for organic carbon with 0.1n K 2 Cr 2 O 7 as an oxidizing agent was used, while 96% H 2 SO 4 was used to extract organic matter (quantification).Color intensity was determined with a Biochrom Libra Spectophotometer (Biochrom Ltd.Cambridge, England) at the wavelength of 585 nm.

Development of the DEM for each master horizon
After generating the surface DEM, we delineated the surfaces of each soil horizons (Ap, A, B, BC, C).Master horizon depths were imported line-by-line into an ArcGIS geodatabase table.Then string-like stitches were strung on soil drilling point feature classes (by the Make Route Event Layer geoprocess function) from the different depth data of the table.From the Z values of the points, raster surfaces were made with the Spline geoprocess function.Since the raster surfaces cannot be extruded into three-dimensional layers (multipatch), TIN surfaces were generated with the Raster To Tin geoprocess function from the raster maps generated by Spline.
Between the adjacent TIN layers, each soil horizons were generated with the ArcGIS Extrude Between geoprocess function.The extruded three-dimensional layers (multipatch) from the TIN generated irregularities (iridescence) at their edge.This phenomenon occurs when the edges of the soil layers are not perfectly smooth.To resolve this error, we used the Buffer geoprocess function to make the TIN layers 5 percent wider.Using the functionality of Extrude Between geoprocess function, we gave the original basemap (Area of Interest -AOI) as boundary.Thus, each soil layer extent from top view deepens according to the base map beneath the different surfaces (Fig. 7.).For the three-dimensional display, however, we had to vertically exaggerate each layer downward to make the layers visually interpretable.

Interpolation of soil moisture values
Based on its suitability on soil moisture interpolation the ordinary kriging technique for interpolation was used.As long as deterministic functions interpolate with one exact value of the measured point, the geostatistical functions utilize the statistical properties of the measured point.Deterministic interpolation techniques create surface from the measured points, based on either the extent of the similarity (inverse distance weighted) or the degree of the smoothing (spline).
As a first step soil moisture can be interpolated with IDW technology at each time step (hour and day) in each depth (5, 15, 25, 35, 45, 55 cm) based on soil moisture values measured by six SENTEK Drill&Drop sensors (Fig. 4.).
The result raster map has a quite bad resolution.However, we know the liquid limit value according to Arany of the given 6 different depths of the 6 sensor (Fig. 5.).From these values a raster map can also be created still with IDW type interpolation.
Soil samples were taken at 35 locations in the field and the liquid limit value according to Arany at each genetic soil horizon of the 35 points were determined under laboratory conditions (Fig. 6.).Samples were collected at the measurement depths of the SENTEK Fig. 4. Soil moisture interpolation process sensors (5, 15, 25, 35, 45, 55 cm).We also read the liquid limit IDW raster value at each soil drilling point and added it to their feature class as an attribute.
Soil moisture data were collected with a manual soil moisture meter (Spectrum TDR-300) in the sample plot from the soil horizons with different liquid limit values seven times between the summers of 2016 and 2017.The data were inserted into a table (Fig. 4.).In each depth of soil drillings, the soil moisture IDW interpolated values were also saved.Based on our table we changed the soil moisture values.
For the development of the interpolated raster, as an example, if the IDW calculated liquid limit value of the H3 borehole at a depth of 5 cm (Fig. 6.) was 46 and the moisture value was 21.4% then we looked at the percentage of moisture at liquid limit value 43 in our table based on the predefined functions.This new moisture value (23.9%) was also saved as an attribute for each depth of the given drilling.We measured moisture values for 41 points (at six different depths).These 41 points were interpolated again by ordinary kriging technique, which requires a minimum number of 30 data points (Or-Wraith 2000).
An alarm function has also been developed for the model, with four categories, with threshold values at tension values of 1,500 (permanent wilting point), 1,000 and 250 kPa.The moisture content -water potential relations were generated with water retention curves based on measured data and fitted with RETC 6.02 using the van Genuchten -Mualem relationship (van Genuchten 1980;Maulem 1976).

Transformation of raster maps into vector maps
Raster maps are easily used for visualization, whereas date attribute cannot be as easily assigned as to vector files.Since we aimed to visualize temporal changes of soil moisture, small vector files were needed (80-90 kB).Firstly, we created the contour lines of the time interval raster maps with the help of Contour geoprocess function then we bounded around polygons from contour lines with the Feature To Polygon geoprocess function.We projected the center of polygons onto the original raster and read out the corresponding soil moisture value.In other words, we generalized soil moisture spatial data by converting raster-based maps to vector maps that are 20 times smaller in size than the corresponding raster files.The completed polygon map allows the realtime motion of time sequence of long-term soil moisture.The soil moisture polygons of different depth were then merged for the entire measurement period to obtain real- time vector files that are movable both in time and space.

2D view of soil moisture values in ArcMap
We also created the top view time sequence for any single master horizons in ArcMap in order to be able to follow each point of each depth of soil moisture values.Soil moisture vector maps calculated for different depths were placed spatially shifted by 300 meters in to the east (EOV Easting 300) while 600 and 1200 meters to the south (EOV Northing -600, -1200).This way moisture pattern of each master horizons was visualized at the same time, however only the map in the upper left corner kept its actual position (Fig. 10.).

Soil data
According to our field soil survey and the subsequent laboratory analyses silt fraction dominated the soils of the studied plot.Genetic soil type was classified as Ramann-type brown forest soils (WRB: Stagnic Cambisol), with clay illuviation in the B-horizon and with occasional stagnic properties.Ap horizon typically ranged between 7 to 12 cm indicating only shallow tillage operations.Depth of the parent material (Pleistocene loess) from the ground surface varied between 50 and 290 cm, depending on hill position, identified as Colluvic Regosol (Table 1).The measured Arany values well correspond with the former measurements (KA = 34-47) carried out by the MINERAG Ltd. and prepared for landowner B-Aranykorona Ltd.For the 35 soil samples taken with borehole drillings, the liquid limits according to Arany ranged between 37.65 and 50.9, with a mean value of 43.8, indicating loamy and clay loam texture for most of the soil samples.
This west margin of the pediment is clearly an accumulation ground where soil sediment is continuously deposited from the uphill part  Soil moisture values over the monitoring period ranged between 6.07 and 43.72% whereas water potential reached its lowest value (-2163 kPa) on October, 3, 2016.The field and laboratory measured soil physical properties formed an integral part of the soil moisture interpolation model, and all measured data have been inserted and employed in the model.

Main characteristics of the obtained software
The current version of the software has a 2D and 3D visualization module, that can be viewed in ArcGIS ArcMap (v.10.4) and ArcScene (v.10.4), respectively.Mandatory input data include name and depths of master soil horizons, number, location (in X, Y Hungarian Unified National Projection System) and depth of soil moisture sensors, liquid limits according to Arany, water The alarm function has four categories in the model denoted with red, yellow, orange and green colors.Threshold volumetric water contents are shown in Table 3.

Discussion
In the current study we have developed a measurement-based simple numeric soil moisture calculation and 2D and 3D visualization model.The model is applicable for many nature conservation and crop yield optimization projects, without complex calculations (e.g.: Richards equation).The software interpolates soil moisture data based on clay content of the soils and is capable to visualize data in two-and threedimensional GIS software environment.We are awere that applicability of described method is limited as it requires spatially detailed measurement of soil texture and soil type.Soil moisture interpolation models may also serve as a substitute for spaceborne radar analyses or rather as a verification tool to validate remotely sensed data.Although a large number of uncertainties are related to soil moisture spatial interpolation, this field of environmental research has significantly matured over the past decades.Despite these considerable developments, point ground soil moisture measurements only provide very limited information both on small-(a few m2) and large-scale (hectares or square kms) soil moisture conditions, especially on terrain of high relief, or uneven vegetation cover.
Besides developing adaptation techniques to global climate changes and altered ecosystem services (Palomo et al. 2017) point measurements and interpolated soil moisture data can be used for many practical applications (Bojovic et al. 2017).Among various uses, some of the most important include (a) agricultural water and irrigation management (b) weather forecasting (Klug -Oana 2015) (c) prediction on fog formation by integrating regional-scale soil moisture into numerical weather models, (d) water balance calculations and estimation of water and solute fluxes and storage in the vadose and root zones, (e) drought forecasting, (f) runoff estimation for flood forecasting, (g) rehabilitation and afforestation projects.
Although remotely sensed data offer unparalleled advantages compared to point ground measurement, it may pose scientific challenges in terms of coarse space and time resolution and shallow penetration depth.Similarly, remotely sensed and groundmeasured soil moisture data should be used in a combined manner in order to minimize ambiguity and improve accuracy related to soil moisture measurements.
The findings of our research can be utilized at plot, and small watershed and floodplain scale, in water balance studies.However, since the model presented is still in its initial state and has been developed on soils of silty and clayey-loam texture, calibration and validation runs are needed to test model applicability for other soil physical soil types.
Soil moisture interpolation accuracy depends on site properties and is also influenced by topography, soil and vegetation type as well as the selected interpolation method.Here, in the current study, we selected the ordinary kriging function for interpolation, whereas in the future we plan to test the model with other interpolation methods suggested in multiple former studies (e.g.Zhang et al. 2016) and the one that proves to be the best fit for our validation data.Our results may be inspiring in terms of understanding water dynamics of soils with complex structure.Also, many soil parameters need to be considered in the model, therefore future development is indispensable.Such parameter addition should include organic matter content, and applicability for sandy soils should also be considered in the future.In the current paper, we described homogeneous flow; whereas, in field cases preferential flow may also contribute to subsurface moisture conditions and solute pulses during individual rainstorms as indicated by former studies (e.g.Hendrickx -Flury 2001;Mullane et al. 2015).

Fig. 1 .
Fig. 1.Location of the study site

Fig. 5 .
Fig. 5. Volumetric Water Content IDW interpolation in 5 cm depth based on value measured by SENTEK sondes of the plot, which is typically characterized by shallow topsoils with a depth less than 60 cm to the top of the C-horizon.Intense rainfall events, due to gully erosional processes, has carved 50-100 cm deep erosion gully into the valley-bottom.Variations in topsoil depths were identified by drilling data.Horizon depth are shown and visualized in the interpolated 3D model (Fig.4.).

Fig. 7 .
Fig. 7. Depth of the master horizons (50× vertical exaggeration) and the location of the soil sampling boreholes

Fig. 10 .
Fig. 10.A cross-sectional view of soil moisture at a selected time

Table 1 .
Characteristics of soil layersHorizon name