Integrating Invasive and Non-invasive Monitoring Sensors to Detect Field-Scale Soil Hydrological Behavior

In recent decades, while great emphasis has been given to the monitoring of point-scale soil moisture patterns and field-scale integrated soil moisture, the measurement of matric potential has attracted little attention. Information on the soil matric potential is available in point-scale measurements but is still missing at field-scale. This state variable is necessary to understand hydrological fluxes and to determine the soil water retention function (WRF) for field-scale applications. In this study, we combine data from cosmic-ray neutron probes (CRNP, non-invasive proximal soil moisture sensors) and SoilNet wireless sensor networks (invasive ground-based soil moisture and matric potential sensors) installed in two sub-catchments with contrasting land-use (agroforestry vs. near-natural forest) to derive a field-scale WRF. We investigate the hypothesis that both sensor types provide effective measurements that are representative for the entire sub-catchment, as well as the drawbacks of integrating the different measurement scales of the sensor types (i.e., spatial-mean of distributed point-scale data vs. an integrated field-scale measurement). We found discrepancies in the data of the two sensor types related to the effects of the time-varying vertical measurement footprint of the CRNP, which induces a scale mismatch between CRNP-based soil moisture (referring mostly to near-surface depths) and the spatially averaged soil matric potential data measured at soil depths of 0.15 and 0.30 m. To remove the offsets, we opted to use the soil moisture index (SMI) based on the estimation of field capacity and wilting point, retrieved from the knowledge of the field-scale WRF. We found that the bimodality of SMI calculated with SoilNet-based soil moisture induced by Mediterranean rainfall seasonal behavior is not well-captured by CRNP-based soil moisture, except in a particularly dry year like 2017. The contrasts in SMI values between the two test sites were associated with differences in the spatial variability of soil moisture patterns explained by soil texture or terrain characteristics. We argue that field-scale WRFs are useful for the analysis of hydrological processes at the sub-catchment (field) scale and the application of distributed models.


INTRODUCTION
Soil moisture (θ ) is an important state variable in environmental systems and can be inferred from ground-based devices, proximal sensors, and remote-sensing platforms enabling observations to be performed across different spatial and temporal scales. Satellite systems are non-invasive measurement methods, with a spatial resolution ranging from about 10 1 m (e.g., Sentinel-1) to 10 3 m (e.g., SMAP, SMOS), which roughly determine near-surface soil moisture when canopy disturbance is not significant, but are unable to provide information in the entire soil rooting zone (Brocca et al., 2017). Direct measurement based on the thermogravimetric technique (Topp and Ferré, 2002), can operationally provide reliable data in sparse locations and is commonly used during sporadic field campaigns. This technique enables the soil water content distribution within the root zone to be measured occasionally and is destructive, timeconsuming, expensive, and therefore unfeasible for large-scale applications (Entin et al., 2000;Romano, 2014). Unattended and automated in-situ monitoring networks for monitoring soil moisture are designed to overcome most of these drawbacks and comprise invasive ground-based instrumentation or noninvasive proximal sensors. The former include point-scale sensors installed in multiple positions and soil depths, thus providing localized information about soil moisture dynamics in a field. The latter consist of stationary "passive" cosmicray neutron probes (CRNPs) that monitor areal soil moisture over a footprint of hundreds of meters in diameter and of several decimeters of soil depth (Ochsner et al., 2013;Vereecken et al., 2015;Babaeian et al., 2019). The issue of coupling rootzone sensor networks and CRNP-based observations has been investigated only recently (e.g., Franz et al., 2012;Peterson et al., 2016;Nguyen et al., 2019). Other proximal sensors are the thermal or spectral cameras carried by unmanned aerial systems (UASs), although this relatively recent sensing technique can provide only sporadic measurements of soil moisture patterns.
At the field scale, soil moisture strongly depends on the local dynamics of the hydrological processes and therefore varies considerably in space and time. Soil moisture is greatly affected by complex and interwoven time-variant (e.g., climate, vegetation) and (mostly) time-invariant (e.g., topography, soil characteristics) factors. Especially in croplands, a variety of human activities and management practices induce changes in soil physical and hydraulic properties that have quite direct repercussions on the spatial-temporal evolutions of soil moisture (Hébrard et al., 2006;Price et al., 2010;Jonard et al., 2013). Grayson et al. (1997) identified the controls exerted by local (such as soil properties, vegetation, etc.) and non-local (such as terrain attributes) characteristics on soil moisture patterns: the former exerts a major influence under drierthan-normal soil conditions, when vertical water fluxes, such as evapotranspiration and infiltration, become the dominating hydrological processes; the latter, instead, is linked more to lateral fluxes and has a greater influence when wetter-than-normal conditions establish in the soil (Orth and Destouni, 2018).
Over the last two decades, monitoring of soil moisture has entered a stage of unprecedented growth, mainly because it plays an important role in controlling the exchange of water and energy between land and atmosphere and partly because it is increasingly employed to validate distributed hydrological models of different complexity over different spatial scales (Orth and Seneviratne, 2015;Nasta et al., 2019). Nevertheless, another fundamental state variable is the soil matric potential (ψ) that enables the total potential gradient to be determined, and hence the water fluxes in the soil domain. Low-cost sensors provide point-scale measurements of soil matric potential that can be coupled to sensors that indirectly measure soil moisture. Data pairs of ψ and θ provide the point-scale water retention function θ (ψ) (WRF) that is a necessary input soil characteristic to solve a Richards-based distributed hydrological model, but also represents valuable information to parameterize a bucket-type hydrological model. Knowledge of the fieldscale WRF allows the estimation of the soil moisture values at the conditions of "field capacity" and "permanent wilting." These two points are commonly employed for the purpose of irrigation scheduling.
In an ideal situation, large scale distributed modeling of hydrological processes should rely on measurements of both field-scale soil moisture and matric potential values. While fieldscale soil moisture can be provided by non-invasive proximal sensors (e.g., CRNPs), to our knowledge there is still a lack of sensing techniques enabling field-scale soil matric potentials to be monitored.
To offer a step forward, a major goal of this study is to explore the feasibility of assessing field-scale soil hydrological behavior of two experimental fields by coupling the measurements provided by non-invasive and invasive sensors. The areaaverage soil moisture monitored by CRNP is integrated with the spatial-mean of the point-scale soil moisture and soil matric potential values measured by a network of multiple low-cost sensors (SoilNet) installed in two sub-catchments of the Alento observatory with different physiographic characteristics. The results section is organized into four sub-sections. The first part concerns the correspondence between the vertical measuring soil volume of the CRNP (i.e., the CRNP's support depth) and the relative positions of the SoilNet sensors below the soil surface (section Assessing the CRNP Vertical Footprint at MFC2 and GOR1). The second part of the results is devoted to identifying the field-scale soil water retention functions at the two experimental sites (section Field-Scale Soil Water Retention Characteristics). The paper then proceeds by analyzing the temporal variability of CRNP-based and SoilNet-based soil moisture data by questioning their ability to respond to the typical rainfall seasonality of the Mediterranean climate (section Temporal Variability of the Soil Moisture Index). Finally, we quantify the spatial variability of soil moisture data explained by more easily retrievable information, such as soil texture and terrain features (section Spatial Variability of Soil Moisture). To remove the offset between CRNP-based and SoilNet-based soil moisture data and the impact of physiographic characteristics on the soil moisture temporal variability, we employed the soil moisture index (SMI) that is a transformation of soil moisture data by referring to the soil moisture values at the conditions of "field capacity" and "permanent wilting" retrieved from the field-scale WRF. The obtained results are firstly discussed to shed some more light on the estimation of a field-scale WRF (section Drawbacks in Setting Up the Field-Scale Water Retention Function) while critically evaluating the drawback associated with the identification of the field capacity and permanent wilting soil moisture values (section Shortcomings Related to the Calculation of the Soil Moisture Index). Then, we explore whether, and to what extent, noninvasive datasets obtained from the stationary CRNPs can be assumed as representative of areal soil moisture values whose spatial variations can be explained, in turn, by easily retrievable controlling factors, such as the soil textural classes and terrain attributes (section Explaining the Spatial Variability of Soil Moisture Data).

Site Descriptions and Monitoring Devices
The Alento River Catchment (ARC) is located in the region of Campania in southern Italy (Nasta et al., 2020a) (Figure 1) and forms part of the TERENO (TERrestrial ENvironmental Observatories) long-term ecosystem infrastructure network . For the purposes of the present study, we selected two sub-catchments, namely MFC2 and GOR1, located in the Upper Alento River Catchment (UARC), which is a hilly upper part of ARC, with a drainage area of ∼102 km 2 and delimited downstream by the Piano della Rocca earthen dam (Nasta et al., 2017). The climate and environmental characteristics of these two experimental sub-catchments are described in Romano et al. (2018). MFC2 is located near the village of Monteforte Cilento on the south-facing hillslope of UARC, and has a drainage area of ∼8.0 hectares (Figure 1a). This site is representative of the cropland zone of UARC with a coexistence of relatively sparse horticultural crops, olive, walnut, and cherry trees. MFC2 exhibits a typical V-shaped topography to form an ephemeral creek in the valley-bottom; hence it has a classic signature from a hydrological viewpoint. GOR1 is the other experimental sub-catchment located near the village of Gorga, on the north-facing hillslope of UARC (Figure 1b). It has a drainage area of ∼22.8 hectares and is representative of the woodland zone of UARC, characterized by chestnut and oak trees with brushwood made up of ferns and brambles growing during summer. This forested site is on average steeper than MFC2 and has a bedrock mainly consisting of turbidite sandstones, with medium permeability and mantled by a regolith zone of sand-silt mixtures.
In 2016 and at each of the above-mentioned experimental sites, a SoilNet wireless sensor network (Forschungszentrum Jülich, Germany) was installed comprising twenty end-devices connected to sensors positioned at the soil depths of 0.15 and 0.30 m. At each soil depth, the apparent soil dielectric permittivity (which is used to estimate θ through an empirical calibration relation), soil temperature, and soil electrical conductivity are measured by the GS3 capacitance sensors (METER Group, Inc., Pullman, WA, USA), whereas the soil matric pressure potential, ψ, is determined by the MPS-6 sensor (METER Group, Inc., Pullman, WA, USA). Note that the MPS-6 sensor can measure ψ values only in the range from −90 hPa (i.e., a matric suction head of about 0.92 × 10 0 m of H 2 O) to −10 6 hPa (i.e., a matric suction head of about 1.02 × 10 4 m of H 2 O) and is, therefore, unable to provide measurements when the soils are near saturation. The SoilNet data are transmitted wirelessly to a local gateway and then via Global System for Mobile Communications (GSM) modem to a central data server in near-real time (Bogena et al., 2010). The SoilNet enddevices were installed around a stationary cosmic-ray neutron probe (CRS2000/B by Hydroinnova LLC, Albuquerque, USA) by covering the experimental sub-catchments as much as possible for future application and validation of hydrological models (Figures 1a,b). The CRNP is a particle detector that measures the neutron intensity in the well-mixed neutron pool above the land surface, which is mainly determined by the amount of hydrogen atoms in the soil. The resulting inverse relationship between soil moisture and neutron intensity is described by the equation proposed by Desilets et al. (2010), where an empirical parameter requires site-specific calibration. The soil volume probed by the CRNP (i.e., in terms of radial footprint and penetration depth) is still a matter of debate, as new theories, additional data, and further interpretations are constantly being added (Köhli et al., 2015). According to Schrön et al. (2017), the CRNP provides indirect measurements of soil moisture over a circular footprint with an effective radius ranging from ∼150 to 210 m (i.e., from about 7 to 14 hectares) depending on various factors, e.g., soil moisture, atmospheric pressure, air humidity, vegetation biomass, etc. The CRNP is hyper-sensitive to soil moisture within the immediate vicinity and the sensitivity to soil moisture decreases non-linearly with radial distance (Schrön et al., 2017). The CRNP is most sensitive to soil moisture in the upper soil horizon, and this sensitivity decreases exponentially to a penetration depth of about 0.3-0.8 m depending on the soil moisture content. Therefore, we used the weighting procedure proposed from Schrön et al. (2017) to calculate appropriate mean values of our point-scale soil moisture measurements with the GS3 capacitance sensors installed at soil depths of 0.15 m and 0.30 m for comparisons with our CRNP-based soil moisture measurements.
One automatic weather station is located near MFC2, at 400 m a.s.l., and another near GOR1, at 711 m a.s.l. (Figure 1). Both weather stations are equipped with the same types of sensors for measuring the following variables at hourly time-steps: rainfall (R), air temperature, air relative humidity, wind speed, and direction, and net solar radiation using four-component net radiation sensors (NR01 net radiometer, Hukseflux Thermal Sensors, The Netherlands). According to World Meteorological Organization (WMO) standards, wind speed and air temperature are measured at a height of 3.0 m, whereas solar radiation is measured with the sensors positioned at a height of 2.0 m above the soil surface. Reference evapotranspiration (ET 0 ) was calculated by using the Penman-Monteith equation according to the protocol proposed by Allen et al. (1998).
At each position of the SoilNet end-devices, disturbed soil samples were collected at the two sensing depths for the laboratory determination of the soil particle-size distribution (PSD). Undisturbed soil cores (steel cylinder of 0.072 m inner diameter and 0.070 m height) were also collected at the soil depth of 0.15 m (vertical sampling at a soil depth of 0.115-0.185 m) of each measuring position to determine oven-dry soil bulk density (ρ b ). Table 1 lists the descriptive statistics of United States Department of Agriculture (USDA) soil particle classes (sand, silt, and clay contents expressed in percentages) at the two soil depths (0.15 and 0.30 m) as well as the oven-dry soil bulk density and soil porosity at a soil depth of 0.15 m for the 20 SoilNet units in both MFC2 and GOR1 (see also Tables A1, A2 in the Appendix reporting measurements of soil physical properties at each location). The individual triplet of sand-siltclay percent is inserted in the USDA soil textural triangle of Figure 2, showing that the MFC2 cropland site is dominated by the clay and silty-clay-loam textural classes, whereas the predominant soil textural class of the GOR1 woodland site is loamy. At MFC2 and GOR1, the mean and standard deviation values of the textural particles are virtually the same for both sensing soil depths. Therefore, the shallow soil layer of both sites can be considered quite uniform in terms of soil properties. However, the mean sand content at MFC2 is almost half that at GOR1, thus offering a first glimpse of the diversity between the two experimental sites. Mean oven-dry soil bulk density (ρ b ) values measured at the soil depth of 0.15 m are also virtually the same at both MFC2 and GOR1, but the spatial variability of ρ b is slightly higher in the forested site. Mean soil porosity is computed from the knowledge of ρ b , assuming the soil particle density always equal to 2.65 g cm −3 . Table 2 reports the descriptive statistics of the following terrain attributes: elevation (ζ ), slope tangent, surface curvature (κ), aspect (in degrees), and upslope contributing area (A c ). The elevation of the forest site is almost twice that of the agricultural site. The high steepness at GOR1 induces a wide range of elevation (almost 100 m) of the SoilNet end-device positions, whilst MFC2 has gentle slopes and hence a very narrow elevation range (about 20 m). The two sites have similar curvature, but different slope aspect as described above.

Calibration Procedure for Invasive and Non-invasive Soil Moisture Sensors
A calibration step is necessary for ensuring good accuracy and precision when estimating point-scale and field-scale soil moisture values from the soil dielectric properties measured by capacitance sensors Gasch et al., 2017;Domínguez-Niño et al., 2019) and from neutron counts measured by CRNP (Franz et al., 2013;Baroni et al., 2018), respectively. Calibration of the CRNP requires simultaneous measurements of soil moisture using the thermogravimetric method and neutron intensity measured by the CRNP to estimate the area-wide soil moisture, θ CRNP . We followed the calibration protocol suggested by Heidbüchel et al. (2016) and conducted three field campaigns to cover both wet and dry climate conditions to account for local climate seasonality. Soil sampling was carried out using a stainless steel core sampler with a plastic liner inside, with a length of 0.30 m and an inner diameter of 0.05 m. The soil samples were collected in 18 positions around the CRNP (six locations along radial distances of 1.0, 10.0, and 110.0 m). Each plastic liner was then cut into six pieces (each of 0.05 m length) to measure the oven-dry soil bulk density and soilwater content with the thermogravimetric method. Therefore, for each experimental sub-catchment, a total of 108 undisturbed soil cores were collected in each field campaign. The annex to this paper provides the Excel file employed for the calculations of the calibration procedure (Heidbüchel et al., 2016).
For calibrating the GS3 capacitance sensors, we carried out five field campaigns to collect undisturbed soil cores with the same core sampler used for calibrating the CRNP. The soil samples were collected over the twenty positions of the SoilNet end-devices. After soil sampling, we cut each core into six subcores (each of 0.05 m in length) to measure the soil moisture value by the thermogravimetric method and also the particlesize distribution. For calibration purposes, only the two sub-cores relating to the soil depths of 0.15 and 0.30 m were considered. To convert the apparent soil dielectric permittivity, ε a , into the GS3based soil moisture, θ GS3 , we used the relation proposed by the METER company (Ferrarezi et al., 2020): which was validated with the soil moisture values determined in our laboratory with the thermogravimetric method (not shown in this paper). Two MPS-6 sensors were tested and calibrated in the laboratory using a pressure plate apparatus in a way similar to that described by Malazian et al. (2011). Finally, the empirical equation was used to compensate for the temperature effect on soil matric potential measured by the MPS-6 sensor (see Equation (3) in Walthert and Schleppi, 2018).
For both sub-catchments, the monitoring program started in spring 2016. Daily (indicated by the subscript d) values of rainfall, R d , and reference evapotranspiration, ET 0,d , together with the soil moisture content time series are shown in Figure 3 for MFC2 and in Figure 4 for GOR1.
The presence of sporadic sensor malfunctioning of the GS3 sensors appears in these graphs due to the abrupt voltage drops of some batteries. Unfortunately, the CRNP at GOR1 (from June 2018 to January 2019) and at MFC2 (from April 2019 to July 2019) experienced a fault because of the modem interruption. The CRNP measurements, θ CRNP , were affected  by episodic spikes during the rainy season probably because of water ponding, low air temperature on the soil surface, or water being retained by leaf interception. Therefore, such values were removed from the data analysis. Table 3 reports the descriptive statistics of SoilNet-based (θ GS3 ) and CRNP-based (θ CRNP ) soil moisture, rainfall (R a ), cumulative reference potential evapotranspiration (ET 0,a ) depths on an annual (indicated by subscript a) basis over the years 2017, 2018, and 2019. It is worth noting that the year 2017 was characterized by a notable meteorological drought, with R a being about <41.5% in 2018 and <33.2% in 2019 for the cropland area and <31.4% in 2018 and <8.5% in 2019 for the forested area. During these 3 years, the mean annual rainfall in GOR1 was constantly and noticeably less than in MFC2, whereas the mean annual potential evapotranspiration ET 0,a in GOR1 was only slightly less than in MFC2.

Data Analysis for Assessing Temporal Variability of the Soil Moisture Index
The soil moisture index (SMI j ) at time j provides a rough quantification of soil wetness and was computed in this study according to Hunt et al. (2009): where θ j is soil moisture content at time j, θ FC is soil moisture at the condition of "field capacity" in the soil profile, and θ WP is soil moisture when on average a plant wilts permanently, being unable to recover its turgor. Note that the difference θ FC -θ WP is commonly defined as the Plant Available Soil Water Holding Capacity (Romano and Santini, 2002). SMI takes on negative and positive values indicating the presence of soil conditions from relatively dry to driest (−5 ≤ SMI < 0) or from relatively wet to wettest (0 ≤ SMI ≤ +5). The values SMI = −5 correspond to the soil moisture content at "permanent wilting, " whereas SMI = +5 is indicative of the soil moisture content at "field capacity" (Hunt et al., 2009). The permanent wilting point (θ WP ) is commonly computed from the knowledge of the soil WRF, θ (ψ), as θ at the matric pressure potential ψ = −15,300 hPa [i.e., 15 bars; Romano and Santini (2002)], and this criterion is employed in the present study. Instead, the condition of "field capacity" in a soil profile  does not have a clear and shared definition and is still subjected to slightly different meanings. Actually, the soil moisture content at "field capacity, " θ FC , is a process-dependent parameter that can be defined as the average soil moisture measured in the entire soil profile when the water flux at its lower boundary becomes virtually nil during a drainage process starting with the initial condition of full saturation and with evaporation prevented at the upper boundary of soil surface (Romano and Nasta and Romano, 2016). This parameter is considered a critical threshold of soil water-holding capacity and commonly employed in bucketing-type models to control hydrological processes, such as overland flow and drainage below the root zone of a uniform soil profile (Romano et al., 2011). The value of θ FC should be obtained using a specifically designed insitu drainage experiment, but can be conveniently retrieved also through numerical simulations, especially in the case of layered soil profiles that are the rule rather than an exception in the real world (Nasta and Romano, 2016). Under simplified assumptions and in the case of a uniform soil profile, the value of θ FC can instead be roughly estimated from the knowledge of the soil WRF using a static criterion that assumes, on average, that θ FC is equal to the soil-water content at the matric pressure potential of −330 hPa. Allowing for the effect of soil texture on the "field capacity" value, Romano and Santini (2002) suggested setting a soil matric pressure potential of ∼-100 hPa for coarse soils and −500 hPa for fine soils. The matric pressure value of −330 hPa should be mostly used in the cases of medium-textured soils (Romano et al., 2011). Given the simplified picture of soil condition offered by the SMI, in this study we estimated the θ FC value through the analytical equation proposed by Assouline and Or (2014) that computes the soil matric potential at "field capacity" (ψ FC ) as follows:

Santini
and hence determines the soil moisture content at field capacity as follows: where α vG (hPa −1 ) and n vG (-) are the two shape parameters featuring in van Genuchten (1980) analytical soil-water retention relationship θ (ψ): with the condition m vG = 1-1/n vG . We use the acronym "vG" as a subscript to refer to the above parameters. The θ s (m 3 m −3 ) and θ r (m 3 m −3 ) parameters are the saturated and residual soil water contents, respectively. While θ s has a clear physical meaning and is measured with laboratory or field tests, θ r is often set at zero or assumed as an additional unknown parameter to be estimated by the fitting procedure together with the other two unknown parameters α vG and n vG . By substituting Equation (3) into Equation (4), the knowledge of θ s , θ r , and n vG , allows to determine θ FC through the following closed-form expression: To test whether the SMI distributions are unimodal or bimodal, in this study we refer to the empirical method that evaluates a bimodality coefficient, termed BC (Pfister et al., 2013), which ranges between 0 and 1 and is computed as follows: where n is the sample size, m 3 is the skewness and m 4 is kurtosis of the distribution. To avoid sample bias, the skewness and kurtosis values are calculated as follows: where x is the arithmetic mean of the set of data. A bimodal distribution is characterized by high skewness, low kurtosis, or both. Specifically, if BC is >0.555 (BC > 0.555), then the hypothesis that the data follow a bimodal distribution cannot be rejected (Kang and Noh, 2019).

Data Analysis for Assessing Spatial Variability of Soil Moisture Data
The partial least-squares regression (PLSR) model is used to reveal the factor explaining the spatial variance of daily soil moisture data. As discussed in (Romano and Chirico, 2004) and Nasta et al. (2018b), most soil moisture variations are cross-correlated with soil physical properties, such as the percentages of sand, silt, and clay contents (see Table 1 and Tables A1, A2), as well as terrain attributes, such as elevation above mean sea level, slope aspect, slope gradient, tangential curvature, and upslope contributing area (see Table 2). The PLSR model generates new predictor variables (components) as combinations of the original predictors and finds combinations of the predictors that have a large covariance with the soil moisture data. The predictive ability of PLSR is expressed in terms of the coefficient of determination and root mean squared error (RMSE) which, respectively, evaluate the scatter of the data points around the fitted regression curve and the bias between observed and modeled soil moisture values. PLSR also provides the percent of variance explained by the components, distinguishing between soil textural characteristics and topographic attributes.

Assessing the CRNP Vertical Footprint
In this study, we investigated the vertical footprint of the CRNP in the two experimental areas (MFC2 and GOR1). Among the variety of weighting procedures provided in the literature, we opted to use the equal (unweighted) spatialmean and a reliable physically-based weighted spatial-mean of point-scale soil moisture measurements (Schrön et al., 2017). The first step is to compare the CRNP-based soil moisture with the unweighted soil moisture spatial-mean, as depicted in Figures 3, 4. We also computed the physicallybased weighted spatial-mean of the SoilNet-based soil moisture data as recommended by Schrön et al. (2017) by using the MATLAB script provided by these authors. The comparisons between the CRNP-based soil moisture and all the unweighted or weighted spatial-mean SoilNet-based soil moisture are depicted in the scatter plots of Figure 5 for both the MFC2 and GOR1 sites. Although RMSE takes on very similar values for the four cases considered, the Schrön et al. (2017) weighted procedure seems to provide a slightly better correlation than the unweighted one especially at the higher soil moisture contents, as evident from the fact that the scatter cloud tends to be closer to the identity line as θ increases. This situation is more pronounced for MFC2 than GOR1. It is also worth noting that the scatter cloud increasingly diverges from the identity line when soil moisture values decrease. Overall, the impact of employing a weighting procedure for our case studies seems to provide only a scant improvement. This situation can be explained if one considers that, allowing for the limited number of available devices, the positioning of the sensor nodes was primarily designed to detect the spatial patterns of the soil hydrological variables as well as possible and not for comparison purposes between the two sensing systems (i.e., a relatively low number of SoilNet enddevices is located in the direct vicinity of the CRNP). That said, in the following we assumed that both vertical and radial weights also apply to obtain the averaged soil matric potential values. Because the soil moisture values observed by the two sensing systems show different means and ranges of variation, the comparison among the time series is more robust when these data are standardized using the Z-scores, which can take on positive or negative values (for example, positive values indicate that the data are above the arithmetic mean). As evident from a perusal of Figure 6, both unweighted and weighted time series of θ GS3 overlap almost perfectly over the CRNPbased soil moisture time series. Nevertheless, the CRNP data show significantly higher temporal variability, thus indicating a stronger response to climate forcing. This feature reflects the higher sensitivity of the CRNP to near-surface soil moisture due to the decreasing sensitivity with soil depth (Schrön et al., 2017).

Field-Scale Soil Water Retention Characteristics
The simultaneous availability of soil moisture, θ GS3,j , and matric pressure head, ψ j , values, at the same time, j, and soil depth, z, from the wireless sensor networks, allows us to identify the SoilNet-based soil water retention functions (WRFs) at the two sites. The WRFs of Figures 7A,B, shown for the two sensing depths of 0.15 and 0.30 m, were determined by coupling the daily unweighted values of the spatial-mean of θ GS3,j and log 10transformed ψ j . Because the MPS-6 sensor is unable to measure matric pressures >-90 hPa, data pairs |ψ j |-θ GS3,j are not available in the near-saturated region of the WRF. Consequently, we set the saturated soil water content, θ s , as 90% of the spatialmean soil porosity obtained from the direct measurements of ρ b (see Table 1). Figures 7A,B also show van Genuchten's water retention relations whose unknown parameters α vG and n vG were obtained by non-linear regression, while setting θ r at zero. All parameters, including the prescribed (θ s ), optimized (α vG , n vG ), and those derived from the knowledge of the WRFs (θ FC , θ WP ) are reported for each site in Table 4 with the corresponding RMSE values. Overall, the two SoilNet-based WRFs belonging to the MFC2 and GOR1 sites (solid black line for GOR1 and dashed black line for MFC2) are characterized by different shapes, with the solid curve of GOR1 showing a shape typical for medium-or coarser-textured soils, with a quite rapid decrease in soil moisture contents in the wet region of the diagram. The dashed curve for MFC2 desaturates with a smooth decay which certainly reflects the typical behavior of finer-textured soils.
According to our method and as depicted in Figure 7C, an "effective" (field-scale) soil WRF over the CRNP footprint can be suitably obtained by coupling, for the same time j, the areal-based soil moisture (θ CRNP,j ) contents with the weighted spatial-mean matric pressure potential (ψ j ) values. The weights applied to ψ j -values were obtained through the weighting procedure recommended by Schrön et al. (2017).
In all the soil water retention curves of Figure 7, the black circles identify the soil moisture contents at "field capacity, " as computed by Equation (6), and the black squares identify the soil moisture contents at "permanent wilting, " corresponding to |ψ| = 1.50 × 10 4 hPa. The picture provided by the θ CRNP (|ψ|) retention curves of Figure 7C allows us to frame in a new perspective the potential of using the areal soil moisture contents offered by a CRNP, and also gives us some preliminary indications on possible drawbacks related to a scale  mismatch between the areal θ and the point-scale ψ data. For both MFC2 and GOR1 sites, the observed scatter of the retention points (see RMSE values in Table 4) is likely due to the different spatial scales associated with the CRNP-based soil moisture contents and the point-measured matric potential pressure values. CRNP-based soil moisture (θ CRNP,j ) and spatial-mean SoilNet-based soil matric potential (absolute value |ψ j |) at MFC2 (red circles) and GOR1 (green circles). Dashed black (MFC2) and solid black (GOR1) lines represent the soil-water retention curves fitted using van Genuchten's parametric relation. Black circles and black squares denote the θ FC and θ WP values at the conditions of "field capacity" and "permanent wilting," respectively.

Temporal Variability of the Soil Moisture Index
The soil moisture index (SMI) is a transformation of soil moisture data and was employed in this study to remove (i) the offset between CRNP-based and SoilNet-based soil moisture observed in Figures 3-5, and (ii) the impact of time-invariant environmental controls, such as the soil textural characteristics and topographical features, on the soil moisture temporal variability (Mittelbach and Seneviratne, 2012). The SMI is computed by using the soil moisture contents at "field capacity" (θ FC ) and "permanent wilting" (θ WP ) as determined in section Field-Scale Soil Water Retention Characteristics and reported in Table 4 (see also the plots in Figure 7) following the procedure described in section Data Analysis for Assessing Temporal Variability of the Soil Moisture Index. Based on the results shown by Romano et al. (2018) regarding the soil moisture distributions, we further processed our datasets by analyzing the histograms of daily SMI values (SMI d ) that are depicted in Figure 8. The individual graphs of this figure enable us to suitably compare the hydrological responses of the cropland (red bars) and woodland (green bars) areas. The SMI d values mostly fluctuate around the zero value with the occurrence of wetter-than-normal conditions induced by rainfall surplus and drier-than-normal conditions determined by rainfall deficits. Note that the period from March to September 2017 was characterized by an exceptionally long dry spell. For almost all months over the investigated 4-year period, the rainfall at MFC2 exceeded that at GOR1, which is reflected in the amount of water present in the soil.
The graphs of Figure 8 show that the histograms of SMI d computed from the θ GS3 data, at the two soil depths of 0.15 and 0.30 m, are characterized by an evident bimodal shape and this visual evidence is also quantified by the bimodal coefficient BC, which always takes on values >0.555. In contrast, the histograms of SMI d obtained from the θ CRNP data follow a more Gaussian-like unimodal distribution, as indicated by BC values lower than 0.555. A perusal of these plots provides a useful indication of the soil hydrological behavior of these two experimental sites. Inspection of Figures 3, 4 highlights that the soil moisture time series measured by the GS3 sensors show the presence of persistent wet and dry conditions that are interrupted by sharp wetting pulses and smooth soil desaturation. The two observed long-term plateaus induce the bimodal shape of the corresponding distribution of SoilNetbased SMI d values that reflect the typical rainfall seasonality of the Mediterranean climate (Vilasa et al., 2017). By contrast, CRNP-based soil moisture data show a constant decrease during desaturation in the growing season and sharp wetting-drying pulses during the rainy season. This temporal variability is typical of the near-surface soil moisture dynamics and is reflected in the unimodal frequency distribution of CRNPbased SMI d values in Figure 8. Nevertheless, a further step is to verify whether these patterns hold within a year-byyear analysis given the remarkable variability of annual rainfall sums reported in Table 3. Figures 9, 10 report the frequency distributions of SMI d values calculated by using SoilNet-based soil moisture at z = 0.15 m and CRNP-based soil moisture for each year at MFC2 and GOR1, respectively. Interestingly, CRNP-based soil moisture is able to capture bimodality only in 2017 (characterized by a long drought) in MFC2 (BC = 0.57). Extremely wet conditions recorded in 2018 induce a unimodal normal distribution in CRNP-based SMI d values in MFC2. Unfortunately, a similar situation could not be explored in GOR1 due to a large amount of missing data at this site.

Spatial Variability of Soil Moisture
This section illustrates the spatial variability of SoilNet-based soil moisture explained by easily-available environmental variables such as soil textural classes and terrain features. Figure 11 illustrates the CRNP-based SMI d values against the spatialstandard deviation of θ GS3,j values measured at the soil depth of 0.15 m in both MFC2 and GOR1. The spatial variability of soil moisture data is related to SMI d values with a somewhat concave downward shape that is more pronounced for the forested site. This shape is similar to that reported by, for example, Teuling and Troch (2005), Famiglietti et al. (2008), and Rosenbaum et al. (2012), although different patterns of this relationship have been reported elsewhere (Mittelbach and Seneviratne, 2012). Because of the larger spatial variability of soil properties in the cropland site, one immediately notices the larger scatter of the points for MFC2 (Figures 11A,B) with respect to GOR1 (Figures 11C,D) in the central range of the catchment's spatially averaged soil moisture measured by the CRNP. This may explain the large scatter observed in the fieldscale WRF belonging to MFC2 (Figure 7C). Further insights are gained by evaluating the role exerted by easily-retrievable local and non-local factors, comprising the five topographical attributes (elevation, aspect, slope, tangential curvature, upslope contributing area) and the three soil textural classes (percentages of sand, silt, and clay content), to help describe the spatial organization of the observed soil moisture datasets. However, apart from a few general comments, we caution the reader that the impact of vegetation was neglected in this study due to a lack of direct robust measurements.
The variability explained by the soil terrain attributes (top panels, Figure 11A for MFC2 and Figure 11C for GOR1) and textural characteristics (bottom panels, Figure 11B for MFC2 and Figure 11D for GOR1) is quantified by using the PLSR method (described in section Data Analysis for Assessing Spatial Variability of Soil Moisture Data) and represented by the color bar on the right side of the plots in Figure 11. It is worth noting that, on average, the terrain attributes explain almost 50% of soil moisture spatial variability in the case of the cropland site (i.e., MFC2), under both wet and dry conditions. In the forested site (i.e., GOR1), the considered terrain attributes explain on average about 40% of soil moisture spatial variability. Nonetheless, this percentage of explained variability increases significantly under wet conditions (bluish colors in Figure 11C) and decreases considerably under dry conditions (reddish colors in Figure 11C). Figures 11B,D show that soil texture has a scant ability to explain the spatial organization of soil moisture data in both sites at the soil depth of 0.15 m.
The five terrain attributes (elevation, aspect, slope, tangential curvature, upslope contributing area) explain most near-surface soil moisture variability in both sites. Interesting aspects emerge if we iteratively repeat the PLSR exercise by removing one terrain attribute at a time: aspect-induced variability of soil moisture plays a major role in both experimental sites; the secondranked topographic attribute is elevation for MFC2, and upslope contributing area for GOR1.

Drawbacks in Setting Up the Field-Scale Water Retention Function
In this study, we proposed to construct an "effective" (fieldscale) soil WRF by coupling the areal-based soil moisture (θ CRNP,j ) contents with the weighted spatial-mean matric pressure potential (ψ j ) values with weights obtained through the weighting procedure recommended by Schrön et al. (2017). On the one hand, a persistent offset occurs between the spatialmean θ GS3 and θ CRNP (see Figures 3, 4) mostly to be attributed to calibration effects; on the other, we observed a relatively large scatter among the water retention data pairs of Figure 7 which may be attributed to several drawbacks. Enhancing the calibration procedure for the CRNPs is certainly to be considered a priority and more field campaigns would probably be required under extremely dry or wet conditions. However, the onset of a sort of scale mismatch represents another important issue which, however, can be tempered by installing some MPS-6 sensors at depths much closer to the soil surface and over more locations, especially near the CRNP (within a radial distance of 5 m). The biomass of dense understory vegetation during  the summer season may have influenced the neutron counts determined by the CRNP (Baatz et al., 2015;Jakobi et al., 2018), and the corresponding impact on soil moisture estimates needs more attention in the near future.
To some extent, we used the SMI values to remove the effect of sensor calibration on both the offset and impact of environmental conditions between the two experimental sites. An analysis of the time series of SMI d values highlights the different behavior between the CRNP-based and SoilNetbased soil moisture contents. The SMI from the sensor networks at soil depths of 0.15 and 0.30 m was able to reflect the impact of Mediterranean rainfall seasonality (Nasta et al., 2020b), although the possibility of having additional capacitance sensors installed at a shallower soil depth would have helped us to better compare our two sensing systems (Franz et al., 2012).

Shortcomings Related to the Calculation of the Soil Moisture Index
In this study, we assumed that the soil moisture index (SMI) can give a useful picture of the soil hydrological behavior of zones with quite different physiographic characteristics, such as our MFC2 and GOR1 experimental sites. Yet the calculation of SMI values requires knowledge of soil moisture contents at field capacity and permanent wilting which, in turn, should be at least retrieved from the soil WRF. If an analytical θ (ψ) relation cannot be obtained from direct measurements, one can resort to a pedotransfer function (PTF) (van Looy et al., 2017) and then compute the values of θ FC and θ WP for example by following the static approaches proposed by Assouline and Or (2014) or Reynolds (2018). Alternatively, Guber et al. (2006) reported a list of PTFs to directly estimate the soil moisture contents at the prescribed matric suction pressure |ψ| of 330 hPa (field capacity) and 15,300 hPa (wilting point).
Another question concerns the lack of direct measurement of the soil moisture content at full saturation, θ s . We assumed θ s as 90% of the spatial-mean soil porosity, which was calculated from the spatial-mean soil bulk density as suggested by Pollacco et al. (2013). To better support non-invasive measurement methods, such as cosmic-ray neutron sensing, we will shortly carry out field campaigns to measure θ s at both experimental sub-catchments.
The observed area-wide θ CRNP values are on average lower than θ GS3 values (see Table 3). The sensitivity of the CRNP exponentially decreases with depth, with most information on θ CRNP being concentrated in the first 0.05 m of soil depth, which is not covered by the SoilNet sensors. The differences between SoilNet-based and CRNP-based soil moisture data can definitely affect the determination of the WRF, hence yielding different θ FC and θ WP values and consequently different SMIvalues (Table 4). In MFC2 the SoilNet-based WRFs lead to field capacity and wilting point values other than those obtained by the CRNP-based WRF (θ FC = 0.338 m 3 m −3 and θ WP = 0.180 m 3 m −3 ). Similarly, at GOR1 the SoilNet-based WRFs lead to field capacity and wilting point values different from those obtained by the CRNP-based WRF (θ FC = 0.356 m 3 m −3 and θ WP = 0.147 m 3 m −3 ). SMI was computed over a period of time (from 2016 to 2019) which may be considered as short for drawing sound conclusions. Therefore, this paper presented only some preliminary results of an ongoing long-term monitoring program underway in the Upper Alento River Catchment. As more data become available, a subsequent paper will explore the relationships between long-term time-series of SMI and climate-based standardized indices, such as the Standardized Precipitation Index (SPI) or Standardized Precipitation-Evapotranspiration Index (SPEI), obtaining robust correlations between meteorological and agricultural/hydrological droughts (Mozny et al., 2012;Martínez-Fernández et al., 2015;Barker et al., 2016). Capturing such relationships is of paramount importance in the decision-making process applied to the management of water resources.

Explaining the Spatial Variability of Soil Moisture Data
Assessing the spatial variability of soil moisture, even over relatively small areas, is challenging and commonly obtained using a large number of sensors sparsely deployed in the uppermost soil layer of a study area and along the soil profiles. This task is prohibitively expensive and laborious, and therefore motivated the present investigation, i.e., to explore whether, and to what extent, non-invasive soil moisture datasets obtained from stationary CRNPs can be assumed as representative of areal values whose spatial variations can be explained, in turn, by easily retrievable controlling factors, such as the soil textural classes and terrain attributes. The practical implication of this outcome concerns, for example, either physically-based or statisticallybased downscaling methods enabling high-resolution spatial maps of near-surface soil moisture (at meter grid-size) to be estimated from the cosmic-ray areal measurements (i.e., at hectometer grid-size) (Qu et al., 2015;Nasta et al., 2018a).
An increase in spatial variability of soil moisture as its spatial-mean increases was reported in several studies (Martínez-Fernández and Ceballos, 2003;Molina et al., 2014), but increasing spatial variability with decreasing spatial-mean soil moisture was observed by others (e.g., Famiglietti et al., 2008). Most of the investigations reported in the literature have detected a concave downward shape for the relationship that links the spatial soil moisture variability and the spatial-mean soil moisture (e.g., Rosenbaum et al., 2012;Fatichi et al., 2015) as we also found at GOR1 for the soil depth of 0.15 m.
By exploiting the results from the PLSR method, Figure 11 shows the spatial variance of soil moisture explained by terrain attributes and soil texture. In MFC2, almost 50% of the observed spatial variance of topsoil soil moisture (i.e., that measured at the soil depth of 0.15 m) is explained by non-local controls (terrain attributes); soil texture is instead able to explain only 5% of the total variance.
In the forested site (GOR1), the spatial variance of soil moisture data explained by easily-available soil and terrain factors increases when moving from dry to wet conditions. The maximum amount explained reaches 92% in the wet season. Beaudette et al. (2013) reported that terrain attributes were able to consistently account for 30-70% of the total variance of soil moisture at 0.10 m soil depth, and 10-40% at the soil depths of 0.30 and 0.50 m. Wilson et al. (2004) were able to explain a small component (<30%) of the soil moisture spatial variability in six experimental sites using terrain attributes.

CONCLUDING REMARKS
This study explores the feasibility of assessing the soil hydrological behavior of two experimental fields by integrating the area-average soil moisture monitored by CRNP and the spatial-mean of the point-scale soil moisture and soil matric potential values measured by a network of multiple low-cost sensors (SoilNet). An in-depth understanding of the hydrological behavior of a field helps improve the application of suitable hydrological models to run numerical simulations. Nonetheless, the measurement of field-scale near-surface soil moisture is sometimes not sufficient to calibrate and validate hydrological models which, in turn, require the knowledge of the soil water retention function. In the present study, we have proposed a technique for obtaining a field-scale soil water retention function as an alternative to laboratory-based methods that are time-consuming, expensive, and rather unfeasible for large-scale applications. However, some drawbacks need to be considered and better understood at least before a proximal sensor for monitoring soil matric potential comes out in the future. The scale mismatch affecting the scatter in the field-scale WRF can be reduced only at the cost of installing a large number of MPS-6 sensors, especially in the near-surface positions. Yet this would increase costs and might become unpractical. Improving calibration procedures can help reduce the offset between CRNPbased soil moisture and weighted spatial-mean SoilNet-based soil moisture as demonstrated by the strong agreement of their Zscores. The CRNP is not able to capture the bimodal distribution reflected by the SoilNet-based soil moisture data, except in a very dry year like 2017. Last but not least, the large scatter observed in the field-scale WRF in MFC2 can also be induced by the large spatial variability of soil moisture which, in turn, is largely explained by terrain attributes. If the CRNP-based soil moisture is assumed to be representative of areal values, the practical implication of this outcome opens room for using topographybased downscaling methods enabling high-resolution spatial maps of near-surface soil moisture.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
NR and PN made the conception and design of this study. All the authors participate in the analysis and interpretation of data.