Mapping Soil Clay Content and Hydraulic Properties over an Agricultural Semiarid Plain Using Remote Sensing and Interpolation Techniques

Clay content is an important parameter governing hydrodynamics property of soils and consequently crucial to environmental management and agricultural development. The present study aims to use the textural middle infra-red index (MID index) product of Landsat-8 Operational Land Images to map clay content over the Haouz plain (Central Morocco). The clay content was mapped at 30 m grid spatial resolution based on the relationship between the MID index and a large set of soil samples. Over the areas covered by green vegetation, the clay content was predicted using the ordinary cokriging technique. Then, this information was used to derive soil hydraulic proper - ties such as the field capacity ( θ fc ), the wilting point ( θ wp ) and the saturated hydraulic conductivity ( K sat ) by using different pedotransfer functions. The validation of the maps was performed by using independent soil samples and measurements. The results showed that the clay content is significantly correlated to MID index. The ordinary cokriging improved mapping of clay content over the Haouz plain (R 2 = 0.70, RMSE = 3.5%). The obtained maps of θ fc , θ wp and K sat revealed a good correlation between the simulated values and the measured values.


INTRODUCTION
Spatial variability of soil characteristics is one of the main requirements of agricultural development and environmental management [Blanco-Canqui and Lal, 2008;Mulder et al., 2011].Furthermore, several models that address hydrological processes, climate change and land degradation require soil input parameters [Bastiaanssen et al., 2005;Anderson, 2008;Lal and Stewart, 2013].Soil texture is one of the main characteristics that affect many physical properties [Greve et al., 2012;Wang et al., 2015].More specifically, the clay content allows predicting soil hydraulic parameters [Shabou et al., 2015].
Mapping the spatial distribution of clay content over large areas using ground measurements is very expensive, time consuming and requires extensive field work [Vidhya Lakshmi et al., 2015;Wang et al., 2015].During the last decades, the remote sensing methods offer the potential of direct or derived soil properties mapping [Kaihua et al., 2013;Shabou et al., 2015].Several studies highlighted that the reflectance of bare soils was successfully used to estimate the soil texture and moisture content [Kaihua et al.,

ECOLOGICAL ENGINEERING & ENVIRONMENTAL TECHNOLOGY
2013].Clay contents were mapped using hyper spectral data [Ouerghemmi et al. 2011, Gomez et al. 2012), or using the visible and near-infrared spectra (400-2500 nm) [Zeynal et al., 2019].The latter is based on the significant correlation between the clay content and the soil reflectance [Kaihua et al. 2013] since the presence of clay in the soil enhances the absorption in the domain of 2.2 μm [Shabou et al., 2015].Because hyperspectral data are not easily available, multispectral data are useful to derive soil surface characteristics from satellites products combining radiometric bands in the visible (0.45-0.69 μm), nearinfrared (0.76-0.9 μm) and shortwave infrared (1.55-2.35μm) spectra [Camacho-Tamayo et al. 2014; Vohland et al. 2014].Consequently, Landsat Thematic Mapper (TM), Landsat-7 Enhanced Thematic Mapper Plus (ETM+), Landsat-8 Operational Land Imager (OLI) and new satellites as Sentinel-2 from ESA are adequate satellite data sources.The two shortwave infrared bands are particularly interesting to determine clay content using the textural indice middle infrared index (MID-infrared) [Shabou et al., 2015].This index requires dry bare soil, which makes this method usable in semi-arid conditions where vegetation is generally scarce.When the soil surface is partly covered with dry or green vegetation, the MID index is not applicable.Therefore, interpolation geostatistical methods are needed [Webster and Oliver, 2001;Nielsen and Wendroth, 2003].Ordinary Cokriging is a powerful technique of interpolation to estimate and predict values at unsampled locations using the sampled locations [Reza et al., 2010].Recently, soil scientists have used this approach to map soil properties from small to large scales using pedotransfer functions (PTFs) [Kaihua et al., 2013 Various PTFs have been proposed based on clay and silts contents [Giarola et al., 2002;Oliveira et al., 2002), sand and clay contents [Reichert et al., 2009;Saxton and Rawls, 2006] or clay content as the main variable [Boulet et al., 2009;Merlin el al., 2016].The PTFs are generally established by multivariate regression based on selected laboratory soil samples or field measurements [Reichert et al;2009;Boulet et al., 2009].
The main objective of the present work was to map clay content on the topsoil of the Haouz plain in the Central Morocco, using remote sensing techniques, and then convert this information into soil hydraulic properties based on PTFs.The study used Landsat-8 Operational Land Imager (OLI) images and a training soil dataset of 200 samples collected over the area.The approach consists in four steps: (1) studying relationships between the MID index and the measured soil clay content over the plain, (2) applying Ordinary cokriging (OK) as interpolation technique to obtain complete coverage of topsoil clay even in the vegetation zones, (3) validating the resulting map based on independent soil samples sets, (4) using PTFs based only on clay content information to derive the hydraulic conductivity, the field capacity and the wilting point.

Study area
The Haouz plain is located at the central Morocco, between 7˚2' and 9˚1' W and 31˚5' and 32˚ N. It extends over an area of about 6000 km 2 .The plain is limited by the Tensfit Wadi in the north and High-Atlas Mountains in the south (Fig. 1).This region has a semi-arid climate, characterized by low and irregular rainfall (250 mm/year) with rainy season concentrated between November-May, and dry one from May to September [Chehbouni et al., 2012;Khabba et al., 2013].The Haouz plain encompasses irrigated and rainfed agricultural zones.The major agricultural crops are winter wheat, olive, citrus and apricot.In the rainfed zone, the main crop is wheat [Khabba et al., 2013].The Haouz plain is filled with alluvial plioquaternary deposits [Boukhari, 2015], mainly constituted by clay, silt, sand, and conglomerate (Fig. 2).
The soil texture map of irrigated areas in Haouz plain (Fig. 3) is obtained by compiling existing soil maps from 1970 to 1990 [Benabdelouahab, 2005].Nine soil classes of the twelve USDA texture classes are observed in the Haouz plain.The FAO textural classes allow simplifying these classes in three classes: fine, moderate/ medium and coarse.The fine class (clay and silty clay) covers the eastern irrigated part (Tassaout) and from the center to the southwest of irrigated zone as shown in (Fig. 3).The moderate texture (loamy and silt loam) covers a large area of these irrigated perimeters in northeast, central and western.Whereas the coarse texture (sandy and loamy sand) appears clearly on the wadis beds of (Ghdat, Zat, Ghmat, Ourika and Rheraya) as well as near the western part of the irrigated area between Ourika and Rheraya wadis.

Soil sampling and analyses
A total of 200 samples were taken from various pedological formations at the topsoil (0-15 cm); 140 in the agricultural zones and 60 over the bare soil (Fig. 1).The samples were packed in plastic bags and properly marked for identification and analyzes.20 g of the mixed soil was sampled for analyzing grain size distribution.These soil samples were air-dried and sieved in two fractions: (0.05-2 mm) to calculate the percentage of sand content.While the smaller fraction passed through (0.05 mm) sieve was recovered and collected in vial then analyzed using pipette and/or granulometry laser methods to measure coarse loam (20-50 μm), fine loam (2-20 μm) and clay content (< 2 μm).The soil samples were classified according to the classification system of the USDA, using the Talwin 42 software.

Selection of the appropriate Landsat-8 images
The images download from http://earthexplorer.usgs.gov/were processed geometrically and atmospherically for 30 m spatial resolution by the MUSCATE software using CNES computing center (add reference for this software).The products are ortho-rectified surface reflectance after atmospheric correction, along with a mask of clouds and their shadows, as well as a mask of water and snow.The atmospheric correction and the cloud detection were processed using the MACCS processor, developed at CESBIO [Hagolle et al., 2010].The top-of-atmosphere reflectance was converted to surface reflectance then used to calculate normalized difference vegetation index (NDVI) and mid infrared index (MID).
The use of remote sensing for soil studies faces 2 issues: (i) the presence of green or dry vegetation that covers the soil surface [Wang et al., 2015], pixels with green vegetation fraction cover more than 20% and crop residue can prevent direct visualization of bare soil and hamper the soil properties estimation [Bartholomeus et al., 2007]; (ii) the high soil moisture that decreases soil reflectance and increases soil absorption in middle infrared [Muller et Décamps, 2000].Therefore, to guarantee a maximum of bare soil we selected four images corresponding to dry periods, the dates selected are (22 June 2013, 08 July 2013, 24 July 2013 and 09 August 2013) (Fig. 4).
It was necessary to extract only bare soil zones and apply the MID index to retrieve soil clay fraction.We applied a method proposed by Shabou et al. (2015) in similar conditions in Tunisia.The green vegetation areas for each four selected images were masked based on the NDVI values by applying a threshold of 0.14.To distinguish bare soil from soil with crop residue and agricultural practices, we tested the evolution of the Red band, NDVI and MID-infrared profiles over a field of cereal crop (  The dates at the beginning of the summer season (22 June 2013, 08 July 2013, 24 July 2013 and 09 August 2013) are considered the best period to map a large fraction of bare soil.In the example of the figure 3, red band recorded a stable value of about 0.2.We cannot conclude if the field was ploughed, but the straw coverage is very low with a slight decrease of MID infrareds index.Straw is completely collected by farmers and finally grazed by sheeps.
Afterwards, the pixels corresponding to the 200 soil samples were extracted from the selected Landsat images (22 June 2013, 08 July 2013, 24 July 2013 and 09 August 2013).

Over bare soil using the MID index
The MID-infrared index is calculated as Shabou et al. (1) () = (1) MID-infrared is a normalized difference index between TM5 and TM7 shortwave infrared bands.The index is applied over the bare soil of the selected images.After masking the vegetation zones, only 100 samples (55 in irrigated zone and 45 in non-irrigated zones) were selected, 50 were used calibration and 50 for validation (Fig. 5).
A simple linear regression model was conducted between measured topsoil clay content and the MID index to produce clay content maps at 30 m spatial resolution over the bare soil zones of the four selected images.

Interpolation over the covered soil using ordinary cokriging
The ordinary cokriging (OCK) was applied to predict the clay content information at 100 m spatial resolution by forming the best linear unbiased estimator between variables with minimum prediction variance [Wu et al., 2009] ( where: γ(u) -represents the variogram for a separation distance (lag) u between two locations Z i (u αi ) and Z i (u αi +u), n i (u) -the number of pairs separated by u.
The geostatistical computations and spatial analysis were performed in a GIS environment, using the geostatistical wizard.
The predicted topsoil clay content was compared to the on 50 field samples dedicated to the validation process.The determination coefficient (R 2 ) (Eq. 9) and the Root Mean Square Error (RMSE) (Eq.10) were used to evaluate the performance of the prediction models [Martin et al., 2020].The first one measure the effectiveness of a variable to predict another variable, and the second one measures the discrepancy of predicts values around observed ones [Dulakshi, 2022] MID − Infrared = TM6−TM7 TM6+TM (1) () = (3) (1) θ wp = 0.0076 *   − 0.0035 θ wp = 0.037 *   0.51 (14) MID − Infrared = TM6−TM7 TM6+TM (1) () = (3)  (10) where: Z(x i ) and Z * (x i ) -the observed and the predicted value of clay content, respectively and n is the number of observations.

Estimation and validation of the soil hydraulic properties
The clay content information was used to derive the soil moisture at field capacity (θ fc ), the wilting point (θ wp ) and the saturated conductivity hydraulic (K sat ) at 100 m spatial resolution in the Haouz plain, by applying different pedotransfer functions (PTFs) that are based only on clay content.
The calculation of θ fc and θ wp makes it possible to estimate the total available water (TAW) which is defined as the capacity of a soil to retain water available to plants which plays an important role for irrigation management.This parameter is estimated as the difference between (θ fc ) and (θ wp ) as mentioned in [Kirkham, 1972]: (1)   (17) To validate the prediction of θ fc , θ wp and K sat , and to check the accuracy of the applied formulas, two different datasets of test fields were used (4) (Fig. 6).The first one corresponds to a soil sampling on the experiment sites of Chichaoua, Aitimour, Saada, Agafay, Agdal, R3 and Sidi Rahal (Fig. 1).

Soil classification using clay content of the soil samples
The sampled textures in the non-irrigated zones consist of high sand fractions (35-85%), moderate clay fraction (25-60%) and low silt fractions (15-40%); the dominant soil type is sandy loam (SaLo).In the irrigated zones the sampled textures consist of high silt fractions (40-60%), lower clay fractions  and a wide range of sand fractions (0-80%); the dominant soil types

Clay content maps over bare soil using the MID index
The MID index is extracted on a 3×3 windows (90×90m) corresponding to the 50 selected samples for each of the four selected images (22 June 2013, 08 July 2013, 24 July 2013 and 09 August 2013).There is a positive correlation between the MID index and the clay content, with significant coefficients of determination (R 2 ) of 0.41, 0.63, 0.56 and 0.72 (Fig. 8).The best correlation (R 2 = 0.72) is obtained for the image of 09 August 2013, according to the following relationship: clay content = 380.8MID-infrared + 5.90.The other images could be affected by the presence of few crops residue, especially in June 2013 that coincided with the harvest period of wheat.
The obtained relationships between the MIDinfrared index and the topsoil clay content are then used to produce the clay content maps at 30 m spatial resolution over bare soil (Fig. 9).The clay content varies from 5 to 40%.

Clay content prediction over the covered soil
The values of clay content obtained for the bare soil at 30 m of spatial resolution were extended to the covered zones using the ordinary cokriging to elaborate the clay content maps at 100 m resolution.
To calculate the variogram, we apply the relationship of 09 August 2013 as primary variable and those achieved for 24 July 2013, 08 July 2013 and 22 June 2013 as secondary variables.The obtained variogram and cross variograms of clay content are well fitted by spherical model (Fig. 10), which corresponds to the best geostatistical model recording the highest R 2 values of about 0.79 (Table 1).In this fitting procedure, the range is 69588 m and the relative nugget effect (RNE) of the variogram is 28.9%, while of the cross variograms are 41.4%, 33.4% and 33.1%, respectively (Table 1).
The interpolated data were validated using the 50 validation soil samples by comparing with the observed topsoil clay content.The obtained results recorded a strong correlation R 2 (0.70) with an RMSE equal to 3.5% (Fig. 11), attesting   The predicted topsoil clay content map at 100 m spatial resolution (Fig. 12) shows six classes of clay content according to the USDA texture, varying from 5 to 50%.The entire Haouz plain is characterized by moderate percentage of clay content varying from 20 to 30% (Fig. 12, 13).
Low clay content (5-15%) soils follow the wadis of Ghdat, Zat, Ghmat, Ourika and Rheraya.This could be explained by the nature of the stream sediment materials characterized by coarse textures [Cochet et al., 1998].The zones with a moderate clay content (20-30%) correspond to the central and northeast of the Haouz plain.The highest clay content (>30%) appears in the eastern part of plain (Tassaout), in the piedmont of the High Atlas, and the western part of the plain along the Chichaoua wadi.This might be inherited from clay and marls material largely present around the plain.

Mapping of derived soil hydraulic properties
The maps of the soil moisture at field capacity (θ fc ) and wilting point (θ wp ) at 100 m spatial resolution over the Haouz plain (Fig. 14, 15) were derived from the map of clay content by using the different pedotransfer functions (PTFs), (Eqs.11 to 14).θ fc varies from 0.16 to 0.35 (mm 3 /mm 3 ) based on the PTFs of Noilhan and Mahfouf (1996) (Fig. 14a) and from 0.12 to 0.51 (mm 3 /mm 3 ) using the one derived from the Beerkan test [Benabdelouahab, 2005] (Fig. 14b).θ wp varies from 0.03 to 0.38 (mm 3 /mm 3 ) by using the Eq.13 (Fig. 15a), and from 0.08 to 0.27 (mm 3 /mm 3 ) using Eq.14 (Fig. 15b).The variation of θ fc and θ wp over the plain is related to the clay content.The low values of θ fc and θ wp are extended along the streambeds, and the high values in the eastern part of Haouz plain (Tassaout).The map of the saturated hydraulic conductivity, which is calculated with (Eq.15) shows a large variation, from 0.14 to 1300 (mm/day) (Fig. 16).Since K sat depends essentially on pore size related to the soil texture [Haghnazari et al., 2015;Eck et al., 2016], therefore the high values of K sat are recorded for coarse textured soils unlike clay and silty clay soils where the K sat values are low.
Figure 17 presents the validation results of θ fc , θ wp and K sat .The results showed that the prediction of θ fc , θ wp over Haouz plain using the (Eq.12) and (Eq.14) is considered reasonable (R 2 = 0.75 and 0.8, RMSE = 0.035 and 0.043 mm 3 /mm 3 respectively).The results also showed that those PTFs are more suitable for our study area compared to the equations derived from the Beerkan test (Eq.11 and Eq.13).Furthermore, the (Eq.15) is suitable for our study area to predict K sat (R 2 = 0.9 and RMSE = 88.8 mm/day).
Diallo and Mariko (2013) stated that a ratio between θ fc and θ wp for a clayey soil equals 1.6.This value is higher compared to our case study with values from 1.39 to 1.44.This underestimation is explained by the fact that the Haouz plain is characterized by a large distribution of moderate clay content class and a small extent of clayey soils.
Using the values θ fc (Eq.12) and θ wp (Eq.14), we calculated the plant available water (TAW) based on Eq. 17.The resulting map (Fig. 18) shows values between 121 and 135 (mm/ mm) with a high frequency ranged from 132 to 134 (mm).This range is appropriate for soils

CONCLUSIONS
The present study used remote sensing data, Landsat-8 Operational Land Imagery (OLI), for mapping topsoil clay content over the agricultural Haouz plain.The obtained map was used to derive soil hydraulic properties such as soil moisture at field capacity (θ fc ), wilting point (θ wp ) and saturated hydraulic conductivity (K sat ).
The results have shown that topsoil clay content over bare soil zones is positively correlated with MID infrared index (R 2 = 0.72).This relationship allowed predicting clay content over the bare soil zones.For zones covered by green vegetation and crop residues, where the MID index is not applicable, the ordinary cokriging was used to estimate clay content distribution.The final map of clay content at 100 m resolution is characterized by the high frequency of moderate clay content.It was validated with independent soil samples (R 2 = 0.70 and RMSE = 3.5%).This map was converted into soil hydraulic properties based on pedotransfer functions.The obtained maps can be used as input data to hydrological and ecological modeling for agricultural and environmental management.

Fig. 1 .Fig. 2 .Fig. 3 .
Fig. 1.The Haouz plain, its land use and location of the soil samples Fig 3).Based on the Red band time series, we could identify the period of agricultural practices such as plowing period by the decrease in Red band reflectance values [Shabou et al., 2015].It was possible then to select the period when there is no crop residue by choosing a period that records stable evolution of MID index [Shabou et al., 2015].

Fig. 6 .Fig. 7 .
Fig. 6.Location of validation points the hydrodynamic properties of the soil

Fig. 9 .
Fig. 9. Clay content maps at 30 m spatial resolution inferred from MID index for (a) 22 June 2013, (b) 08 July 2013, (c) 24 July 2013 and (d) 09 August 2013.The grey pixels correspond to the clay content values, the white one corresponds to the masked zones

Fig. 11 .Fig. 12 .
Fig. 11.Comparison between estimated clay content by using the ordinary cokriging with the observed topsoil clay content at 50 independent sampling points