Why We Should Include Soil Structural Dynamics of Agricultural Soils in Hydrological Models

Abstract

Surface soil structure is sensitive to natural and anthropogenic impacts that alter soil hydraulic properties (SHP). These alterations have distinct consequences on the water cycle. In this review, we summarized published findings on the quantitative effects of different agricultural management practices on SHP and the subsequent response of the water balance components. Generally, immediately after tillage, soils show a high abundance of large pores, which are temporally unstable and collapse due to environmental factors like rainfall. Nevertheless, most hydrological modeling studies consider SHP as temporally constant when predicting the flow of water and solutes in the atmosphere-plant-soil system. There have been some developments in mathematical approaches to capture the temporal dynamics of soil pore space. We applied one such pore evolution model to two datasets to evaluate its suitability to predict soil pore space dynamics after disturbance. Lack of knowledge on how dispersion of pore size distribution behaves after tillage may have led to over-estimation of some values predicted by the model. Nevertheless, we found that the model predicted the evolution of soil pore space reasonably well (r² > 0.80 in most cases). The limiting factor to efficiently calibrate and apply such modeling tools is not in the theoretical part but rather the lack of adequate soil structural and hydrologic data.

1. Introduction

The need for sustainable intensification of agriculture to feed our growing population against the backdrop of climate change has led to worldwide attention on agricultural management practices (AMP) that are adaptive towards these challenges. AMP can denote multiple measures such as tillage, fertilization, plantation, fallows, etc. Such practices are often not separate but are part of the broad spectrum of sustainable land-management practices [1]. Although AMP are applied as instantaneous interventions, their effects on soil structure are active over various time scales (Figure 1, [2]). On long-term scales (years to decades), land-use changes such as deforestation or conversion from pasture to cropland result in long-term adaptation of soil structure [3,4,5,6]. On the other hand, tillage is the most analyzed AMP and, along with crop management (e.g., cover crops, crop rotation), affects the soil structure on a seasonal scale [6,7,8]. As a consequence, how much water infiltrates, evaporates, or leaks depends not only on the climate and vegetation but also on the geometry of the soil pore space, pore size distribution (PoSD), and soil structure.

While the textural or intra-aggregate pores in agricultural soils are not expected to vary over time, the structural or inter-aggregate pores, which drain at saturated and near-saturated conditions, are most affected by AMPs and natural stresses like rainfall, freeze-thaw (FT) and wetting-drying (WD) cycles [9,10]. The retention and transport of water and nutrients in the soil which are governed by soil hydraulic properties (SHP) are also affected as a result. Nevertheless, the temporal variations of soil pore space and SHP are often not addressed in hydrological modeling studies. In fact, an invariant pore system is assumed when modeling the transport of water and solutes in the atmosphere-plant-soil system [11]. Ignoring the alterations in soil pore space and SHP may lead to incorrect predictions of water and solute fluxes and therefore to invalid statements, misrepresentations and finally to poor decision making, for instance by policymakers and land-use planners.

In recent years, there have been advancements in using mathematical models to evaluate the impacts of AMP, climate and land cover on water balance components and erosion risks [12,13] as well as the temporal evolution of soil PoSD [14]. The few studies that have included these dynamics in mechanistic hydrological models have shown that it results in better and more reliable modeling outcomes [6,8,15,16,17,18]. Therefore, efforts have been made in the direction of investigating and modeling the temporal evolution of soil PoSD. In this context, Or et al. [19] proposed to use the Fokker–Planck equation (FPE) to capture the changes in soil pore space geometry following tillage with time and pore radius. Nonetheless, Vereecken et al. [20] identified that ‘approaches to simulating temporal changes of soil structure … are relatively rare and at an early stage of development’. Quantifying the temporal dynamics of soil structure is essential to gain a holistic perspective on the physical support services of soil and to understand the response of crop systems to changing soil water contents.

As a basis to understand the changes in soil hydrology due to AMP or climate change, a better understanding of dynamics in soil pore space is required [19,21,22,23]. This knowledge may help us develop tools for sustainable food production and integrated use of natural resources [6,16,24,25]. Based on these assertions, the objective of this paper is to review the temporal dynamics of SHP due to AMP and environmental impacts and the response of soil water simulations to the incorporation of time variant SHP. We then applied the model of Or et al. [19] to two water retention data sets to evaluate its suitability to predict the evolution of soil PoSD following tillage as well as when there is a change in tillage regime. Through this review article, we intend to contribute to the development of tools to model temporal changes in soil structure and incorporate them in hydrological models.

2. Impacts of Anthropogenic and Environmental Influences on SHP

2.1. Overview of SHP

Generally, SHP refer to the soil water retention function (θ(h)) and the hydraulic conductivity function (K(h)) [26,27]. θ(h) is the relationship between volumetric water content (θ) and soil water head (h) and is theoretically based on the soil PoSD via the law of capillarity. θ(h) reflects the ability of the soil to hold and retain water whereas K(h) governs the movement of water and nutrients in the soil. SHP influence the plant available water capacity (PAWC), plant rooting depth and pattern as well as the dynamics of hydrological-biogeochemical processes. Both θ(h) and K(h) are influenced by soil structure, texture, organic matter, chemical and biological properties. Accurate modeling of flow and transport processes in the soil requires good understanding of pore space and SHP dynamics and their hydrological impacts.

SHP are defined by parameters that are obtained from direct measurements of θ(h) and K(h) or through indirect methods. Direct measurements of SHP enable capturing site-specific information of soil structure with high spatial and temporal resolutions. Here, water retention is generally measured in the laboratory using undisturbed soil samples from pressure heads close to saturation till h = −1000 cm which is a methodologically based limit. After that, usually, the next few measurement points are below this value till h = −15,200 cm (wilting point) using disturbed samples. Retention data are interpolated and the goodness of SHP parameterization is a function of the number of available measurement points and the required fitting of points coming from different techniques of measurements. With respect to the conductivity measurements, unsaturated hydraulic conductivity (K_u) is generally measured in the field using disk or tension infiltrometers (e.g., [6,8]) while in the lab, evaporation method is employed (e.g., [10]). Saturated hydraulic conductivity (K_s) is often determined in the field using single or double ring infiltrometers or using hood infiltrometer [28]. Recently, Weninger et al. [29] evaluated the combination of different measurement methods to parameterize SHP to predict soil water dynamics. The authors concluded that a combination of field and lab measurements (e.g., hood infiltrometer in combination with the evaporation method) may provide a basis for representative predictions of soil water dynamics.

There are also developments in visual soil evaluation methods that help to assess the structural state of a soil (for review and future directions in this technique, see [30,31]). Proceeding, indirect methods estimate SHP from readily and easily available data that are more convenient and less expensive. For instance, pedotransfer functions (PTF) are empirical relationships between models of soil parameters and easily available proxy variables on soil properties. Though they have become an indispensable tool for modeling soil processes, their ability to capture spatio-temporal dynamics of soil properties have been under debate [32].

The estimated parameters, either through direct or indirect methods, are then used as inputs which define SHP in numerical modeling techniques. Model functions for θ(h) and K(h) have been subject to extensive research and have resulted in a large number of combined water-retention-hydraulic-conductivity models (for a recent example, see [33]) for a performance evaluation of the most used θ(h) models, see [34]). A widely used example is based on the assumption of a lognormal PoSD (f(r)) [35]:

$$f\left(r\right)=\frac{{\varphi }_0}{r\sigma \sqrt{2\pi }}{e}^{-\frac{{\left[\ln \left(\frac{r}{r_m}\right)\right]}^2}{2{\sigma }^2}}\mathrm{where}{{\displaystyle \int }}_0^{\infty }{f}_0\left(r\right)dr={\varphi }_0,0<r<\infty$$

(1)

where ϕ₀ [-] is the total initial porosity, r_m [L] is the initial median pore radius or geometric mean and σ [-] is the standard deviation of the log-transformed pore radius. Integrating Equation (1) yields an expression for θ(h) which may be combined with corresponding models to derive an expression for the relative hydraulic conductivity. For example, description of hydrophysical properties of the soil using θ(h) and K(h) by van Genuchten–Mualem method [33,36,37] is frequently used. θ(h) obtained from Equation (1) shows great flexibility in representing θ in the wet and dry ranges for all soil types.

2.2. Dynamics of SHP

For agricultural soils, tillage is the most widely used management practice that has the greatest effect on soil structure [7,38]. Tillage is defined as mechanical manipulation of soil for crop production and it alters the soil moisture regime as well as infiltration and evapotranspiration processes [39]. Based on the depth and intensity of soil disturbance, tillage systems can be classified into conventional tillage (CT), reduced tillage (RT), and no tillage (NT). Generally, CT represents the most intensive tillage treatment and includes manual or mechanized plowing using chisels, moldboards, and disk cultivators to a depth of about 30 cm [7]. Under RT, at least 30% of the crop residues are left in the field and it involves treatments such as mulch, ridge and contour tillage [39,40]. NT involves little or no soil disturbance and physical impact is usually limited to seedbed preparation. In addition, cover crops and crop rotation practices are undertaken to maintain or enhance soil fertility and improve soil quality. They improve environmental benefits (i.e., in soil and water conservation) and provide ecosystem services that address food security and overall agricultural sustainability [41]. The effects of AMPs on SHP have been quantified by several studies [7,42,43,44,45]. Strudley et al. [7] presented a state of the science review on the tillage effects on SHP in space and time. The authors found that most tillage practices have considerable effects on SHP immediately after tillage which then diminish with time. However, it was seen that the effects of tillage (e.g., comparison of CT vs. NT) did not produce consistent results which was attributed to factors such as soil type, climate, location as well as the challenge to separate the effects of temporal and spatial variations.

Table 1. Effects on soil hydraulic properties (SHP) due to agricultural management practices and atmospheric impacts published after 2008.

Study	Study Period	Location	Soil Type	Properties Measured	Land-Use/Cover	Effects on SHP
[10]	2006	Oldenburg, Germany	Podzol and Stagnosol	ρb, KsCH, KuEM, ICDRI, θ(h)	Forest, grassland and cropland with CT	Low ρb immediately after CT followed by gradual increase through a season as well as low Ks, Ku. Both forest and grassland show the opposite trend. IC, Ku and Ks was highest in forest.
[46]	2008	Prague, Czech Republic	Haplic Luvisol	ρb, KsCH, θ(h)	CT (after harvest of winter barley) and grassland	Lower θ(h), Ku under CT compared to grassland. Higher Ks values in CT indicating greater fraction of larger pores and pore connectivity.
[6]	2008–2009	Canterbury, New Zealand	Haplic Cambisol	KsHI, KuHI,TI, θ(h)	CT and pasture	Ks and Ku values under CT were 4 times higher than under pasture and flow active macropores were 80% larger. Pasture showed a more balanced PoSD.
[8].	2008–2010	Raasdorf, Lower Austria	Typic Vermudoll	KsTI, KuTI, θ(h)	CT, RT, NT under winter wheat	Ks, Ku increased in the order CT < RT < NT. Higher connectivity and smaller tortuosity of macropores under NT.
[47]	2005	Piedmont, NW Italy	Typic Endoaquept	ρb, KuTI	CT under rice and maize with different water management systems	Progressive soil compaction and increase in ρb, destruction of macro- and mesopores due to submerging water and rainfall.
[48]	2007–2009	Liudaogou watershed, China	Calcaric Arenosol and Calcaric Regosol	ρb, KsCH	4 different mixed shrub landscapes	Over-winter increase of Ks and decrease of Ks in the drying season due to FT and WD cycles. Reverse trend was observed in ρb.
[49]	2011	Saskatchewan, Canada	Haploboroll	θ(h) and infiltration rate	Prairie	Lower infiltration rates because of increased FT cycles, importance of antecedent soil moisture conditions for potential ice formation in soil pores.
[50]	2009–2012	Raasdorf, Lower Austria	Chernozem	θ(h)	Post harvest agricultural soils under NT	WD cycles enhanced macroporosity, decreased pore heterogeneity. Prediction of temporal changes in soil PoSD from environmental variables.
[51]	2007–2010	Hněvčeves, Čáslav & Humpolec, Czech Republic	Haplic Luvisol, Greyic Phaeozem, Haplic Cambisol	Aggregate stability, ρb, Ks, KuTI, θ(h)	CT under winter barley and winter and spring wheat	Variable trends in SHP for all the years depending on rainfall, plant growth and tillage. Different Ks results were observed when different measurement methods were used
[17]	2012	Lamothe, France	Stagnic Luvisol	ρb, KsTI, KuTI, θ(h)	CT under maize	Increase in ρb and decrease in Ks over the growing season of the crop from seeding to harvest. Unstable porosity after tillage
[52]	-	Zhonggou catchment, China	Calcaric Regosol	KsHI, KuHI, θ(h)	Forest (black locus), grassland and cropland (under CT with rape)	Creation of temporary macropores after tillage after which there is a decrease in PAWC. Conversion from CT to grassland or forest increases IC and Ks. Forest demonstrated increases macropore connectivity and water storage capacity.
[53]	2010	Ceará, Brazil	Typic Hapludult	Aggregate stability, θ(h)	Natural vegetation forest converted to irrigated pasture for 10 years	Lesser total porosity and macroporosity in irrigated pasture compared to natural vegetation. Reduction in θ(h) under pasture.
[54]	2010–2013	Attica, Greece	Typic Xerofluvent	KsCH, θ(h)	Rototillage and NT with bare soil and NT with weed cover	Rototillage caused a decrease in the number of mesopores and an increase in the micropores. Smaller Ks values for all three management practices after rain events and an increase in Ks during dry periods.
[55]	2013	Southern Chile	Histi-Silandic Andosol	θ(h)	Land-use change from native forest to naturalized grassland	Soils under grassland showed smaller water storage capacity, air capacity and permeability compared to forest. Land-use changes affected spatial and temporal variability of soil physical functions.
[42]	2011–2012	Zaragoza, Spain	Hypercalcic Calcisol	ρb, θ(h)	CT, RT and NT on fallow land	Tillage operations caused a decrease in ρb. Rainfall was a major factor that affected the θ(h) curve due to disintegration and deformation of soil aggregates.
[43]	2010–2013	Ponta Grossa, Brazil	Typic Haplorthox	θ(h)	CT and NT with crop rotation	Better pore connectivity in NT. Larger values of air filled porosity for NT in relation to CT. Tri-modal and bi-modal PoSDs for CT and NT respectively.
[56]	-	NW Iran	Calciexerept	ρb, θ(h), infiltration rate	Raindrop impact on CT soils	Strong increase in aggregate breakdown, ρb and compaction with high tendency to form crusts due to raindrop impact. Decrease in infiltration rate and enhanced runoff
[44]	2015	Nebraska, USA	Typic Argiudoll	KsCH, θ(h), total porosity, infiltration	Chisel plow, disk, moldboard plow and NT	Ponded infiltration was higher under moldboard plow tillage plow than for NT, disk Long-term till had no significant impact on soil porosity, Ks, θ(h) and PAWC.
[57]	2014–2015	Lower Saxony, Germany	Luvisol	KsFH, θ(h), ρb, KuTI	Moldboard plow (CT), chisel plow (RT1) and disc harrow (RT2) under winter wheat	Higher Ks and K(h) in RT1 and RT2. Higher θ(h) in CT immediately after tillage both in untrafficked areas and tramlines.
[58]	2013–2014	Aquidauana, Brazil	Ultisol	ρb, porosity, aggregate stability, penetration resistance	CT with bare soil and soybean, NT with soybean and pasture	Highest macroporosity in pasture. Higher aggregate stability in pasture and NT due to organic matter. High ρb and penetration resistance in NT. In the initial stages, soil tillage influenced infiltration the most while in the final stages, vegetation cover predominantly influenced infiltraition.

Note 1: ρ_b: bulk density; θ(h): water retention function; K_s: saturated hydraulic conductivity; K_u: unsaturated hydraulic conductivity; IC: infiltration capacity; NT: no-tillage; RT: reduced tillage; CT: conventional tillage; FT: freezing-thawing; WD: wetting-drying; PoSD: pore size distribution; PAWC: plant available water capacity; SHP: soil hydraulic properties. Note 2: Properties measured: the superscripts denote the following measurement method—CH: constant head; FH: falling head; EM: evaporation method, DRI: double ring infiltrometer; HI: hood infiltrometer; TI: tension infiltrometer.

lists selected studies on the impact of AMPs on SHP published after 2008 from which some studies are discussed in detail below.

2.2.1. Effects of Tillage

Kargas et al. [54] investigated the temporal dynamics of θ(h) and K_s for sites under annual rototillage (RoT) and NT on surface soil kept bare using glyphosate as well as NT under weed cover in Attica, Greece. Undisturbed soil samples were taken in the months following tillage and water retention curves were determined in the laboratory using a sand-kaolin box in the pressure head range h = 0 to −200 cm (h = 0, −40, −80, −120, −160, −200 cm). The authors observed reduced K_s values and decreasing θ(h) under RoT in relation to NT due to destruction of the pore sequences and their continuity that existed before application of tillage. Additionally, raindrop impact rapidly disintegrated soil aggregates into single particles and reduced soil pore space which has also been observed by [2,8,42]. This caused a further reduction in the K_s values for all treatments. However, K_s recovered for all treatments during spring and summer which is advocated to changes in the soil pore space conditions created by the intensity and frequency of WD cycles [2]. Aggregate disintegration due to rainfall may occur through slaking and/or mechanical breakdown related to the kinetic energy of a specific event [59].

Schwen et al. [8] investigated the temporal variability of SHP under CT, RT and NT on SHP for two years in a wheat field. The infiltration measurements were performed using a disc infiltrometer for h values of −10, −4, −1 and 0 cm. Both K_s and saturated water content (θ_s) strongly decreased after tillage in both CT and RT due to raindrop-induced pore sealing and settling. Furthermore, the near-saturated hydraulic conductivity varied in the order of CT > RT > NT and a lower volume of hydraulically active pores were observed in winter due to FT cycles and raindrop impact. FT cycles lead to expansion of water volumes as it freezes in the soil and exerts considerable pressure on the soil pore walls leading to reduced K_s and infiltration rates. Furthermore, Xu and Mermond [15] studied the effects of CT, NT and sub-soiling tillage (conducted using a sub-soiler till 40 cm depth to structurally disrupt compacted subsurface layer and rototilled till a depth of 15 cm to break down large clods). In wet years, subsoiling promoted infiltration when compared to NT and CT and reduced the risk of water logging in the fields. Meurer et al. [60] reviewed the effects of CT, intermediate tillage (chisel, disk, harrow, mulch plough) and NT on soil organic carbon. The authors found significant differences in bulk density (ρ_b) between the three tillage treatments. For soil organic carbon, they found that intermediate tillage and NT contributed to higher organic carbon content in the top 30 cm of the soil when compared to CT. Organic matter content has implications for soil macropore characteristics, water transmission and K(h) and may impact soil structural resilience towards rainfall, WD and FT cycles [61,62].

2.2.2. Effects of Cover Crops

The effects of cover crops on SHP were studied by Blanco-Canqui et al. [41]. Cover crops promoted aggregate stability by means of roots generating organic binding agents. In the long-term and in combination with NT, cover crops promoted K_s by increasing macroporosity and pore connectivity and increased water infiltration and storage. A review of Ball et al. [63] on the role of crop-rotations indicated that the amount of organic matter incorporated into the soil was more important to preserve soil structure than the farming technique employed. Crop-rotation led to development of macropores due to root growth and faunal activity which influenced the water and nutrient availability in the soil. The study by Bhattacharyya et al. [64] showed that a crop-rotation of soybean–pea in combination with NT resulted in higher unsaturated hydraulic conductivity values at suctions of −2 and −4 kPa. A 3.5% decrease in ρ_b was observed for cover crop plots when compared to no-cover plots by Haruna and Nkongolo [65]. Furthermore, management decisions such as the number of crops and crop type in a rotation also have significant effects on soil organic matter and biota. For instance, McDaniel et al. [66] found that increasing the number of crops from two to three increases the total carbon by almost 3.6%. Similarly, an increase of 33% in organic carbon was observed by Tiemann et al. [67] when the number of crops were increased from one (monoculture) to five. Organic matter content accelerates the formation of macroaggregates which promotes soil water retention and transmission [68]. Blanco-Canqui and Ruis [69] highlighted the need for accompanying AMPs, such as cover crops, diversified cropping systems and manure, together with NT in order to achieve a beneficial soil physical environment for crop growth.

2.2.3. Effects of Land-Use Change

The evolution of PoSD and SHP after a shift in tillage regime from CT to NT or to pasture and grassland has also been evaluated. For instance, Schwärzel et al. [6] characterized changes in soil structure and SHP due to a crop-pasture rotation cycle in Canterbury, New Zealand in the period from 2008–2009. Undisturbed soil cores were used to determine water retention curves by means of dewatering using ceramic plates connected to a hanging water column (h = −10, −31, −100, −316, −1000, 15,800 cm) and fitted to the Kosugi lognormal distribution model. The authors observed a reduction in mean pore sizes and an increase in the spread of pores (σ in Equation (1)) as well as a loss of large pores. This indicates that macropore rich structure created by tillage is rather unstable and subject to rapid collapse [70]. Hu et al. [71], who investigated the contribution of different pore classes to flow under four different land-uses, found that K(h) showed a significant change only at h = −15 cm for all four land-uses. Furthermore, a study by Yu et al. [52] investigated the effects of different kinds of land-use (cropland and grassland) on soil structure in the Zhongghou catchment in the Loess Plateau region of China in the period 2012–2013. Undisturbed soil cores were used to determine the water retention curves by means of evaporation method and pressure cell (h = −10, −20, −31, −63, −100, −200, −316, −15,000 cm) and fitted to the model of van Genuchten [37]. The authors observed that conversion from CT to grassland or forest tended to promote infiltration capacity, K_s due to creation of meso- and macropores that also increases the soil’s water retention capacity.

Summing up, from the results of the above studies and the outcomes of the studies listed in Table 1, it is observed that immediately after tillage, soils generally show an abundance of large pores. However, these pores are unstable and the PoSD changes with time due to biological activity, rainfall, WD and FT cycles as well as shrink–swell phenomena [10,72]. Raindrop impact on agricultural soils generally leads to aggregate breakdown, compaction, soil crusting and reduced infiltration capacity [56]. As a result, the variance of the PoSD (σ² in the case of [35]) increases indicating a heterogenization of the pore system. Frequently, a reduction in the transmission pores or macropores (50 < r < 500 µm) and an increase in the storage pores or micro and mesopores (0.5 < r < 50 µm) is observed in the months following tillage or when there is a shift in tillage regime from CT to NT or pasture. K_s decreases in periods of rainfall due to sealing of pores and settlement [73,74] and increases in spring due to biological activity and WD cycles [50,74]. The presence of cover crops and crop rotation cycles lead to development of macropores, pore connectivity and increase of organic carbon which enhances soil water infiltration and storage.

2.3. Inconsistencies in the Outcomes on Dynamics of SHP

From the previous section, it can be concluded that SHP are subject to temporal variations post-tillage through a season or after a shift in the tillage regime. However, studies that have reported alterations in SHP are inconsistent in their outcomes and they vary across experimental designs, locations and soil types [7] (see Table 1). For instance, Blanco-Canqui et al. [44] did not observe any significant impact on K_s and θ after long-term CT while Schwen et al. [16] observed a strong decrease in K_s and θ under CT. These contradictory findings may be explained by differences in land-cover (continuous corn vs. winter-wheat), use of different equipment for CT (chisel plow, tandem disk vs. moldboard plow) and soil composition.

In some cases, NT and RT exhibited an increase in near-zero infiltration rates and hydraulic conductivity which is associated with macropore connectivity while in others water infiltration rates may be greater for soils under CT [7,44]. Most studies have also reported an increase in ρ_b under NT due to soil densification in the absence of soil turnover by tillage operations [69,75]. The contradictory results may be attributed to temporal variability of soil infiltration rates; Water infiltration is higher immediately following tillage because tillage creates a macropore rich structure after which there is generally a collapse of large pores. This is associated with a decrease in infiltration in the weeks following the event.

Despite progress in methods to quantify SHP, there exist inconsistencies in the results from different management practices on SHP. These inconsistencies may be due to the use of different methods for the characterization of SHP. Different measurement techniques may result in different initial and boundary conditions being applied so that the behavior of SHP may actually depend on the technique used rather than on the intrinsic nature of soils [76]. Some authors used laboratory methods for the characterization of K(h) (e.g., multi-step outflow method) while others used field methods (e.g., tension infiltrometer). Field techniques are often inconclusive owing to the variability of soil structure, especially in the upper layers of the soil. Moreover, some field methods (such as disc or hood infiltrometer) are able to account for the contribution of macropores to the infiltration rate and soil hydraulic conductivity at saturation [28] whereas other methods (such as double ring infiltrometer) cannot differentiate the contribution of macro-, meso- and micropores to the water flow. Laboratory measurements are generally seen as more precise and straightforward than those in the field [76,77] because they can be conducted in a controlled environment whereas field measurements are subject to variations and conditions of the field. In this context, Basile et al. [78] reported that by using water content measured at zero head and air entry value, both measured in the field, it was possible to derive field hydraulic functions from laboratory functions with adequate approximation. However, reliable measurements of near-saturated hydraulic conductivity in the laboratory depend on the soil core size where the soil sample has to be large enough to represent the soil horizon at the pedon scale [76]. The relationship between the sample volume and variability in SHP measurement is an example of the concept of representative elementary volume (REV). If a REV does not have approximately 20 elementary units of soil structure, issues such as non-detection of large-range spatial structures of SHP and observation of inaccurate values of SHP may emerge [79].

Different stages of root growth also affect the outcome of SHP measurements. Murphy et al. [80] found a decrease in K_s in wheat from sowing to tillering phases of plant growth. However, almost 17% of the roots produced from emergence decayed when the plants reached the phase of stem elongation while 37% of the roots decayed by the ripening stage which led to an increase in K_s. Bodner et al. [81] studied the root effects on macroscopic pore parameters and found that rooting densities greater than 0.5% of pore space stabilized soil structure against pore loss. The authors also observed that coarse root systems increased macroporosity by 30%. Ajayi et al. [82] observed strong contribution of root development to improved pore networks and pore continuity for conversion of CT to grassland over a 19-year period. It is often difficult, however, to separate land-use induced changes from changes caused by other factors such as climate.

Strudley et al. [7] highlighted that many studies usually investigate contrasted “snapshots” of different treatments; i.e., an instantaneous capture of soil conditions immediately after AMPs without considering the evolution of the soil pore space at later stages. Furthermore, most studies dealing with the impact of management practices focus only on the uppermost layers of the soil profile [9,19,83]. This may be appropriate for soil dynamics on a rather short temporal scale such as changes induced by weather or tillage. However, mid- and long-term changes of management practices and long-term climate change effects have a significant influence in deeper soil layers [5,11,84].

The described inconsistencies hamper the generalization of results on the impacts of AMP on SHP. Often, experimental protocols (e.g., experiment designs, measurement times, sample sizes, etc.) are not rigorous enough to reveal the differences between different AMP [85]. Particularly at high and low values of h, there is a lack of reliable data owing to the difficulty in measurement in these ranges. For a better generalization of results from different measurements, there is a need for accurate and harmonized description of methods used to quantify alterations in SHP. Science-based quantitative data and standardization of research methodologies are critical towards improving our interpretation of results from different experiments on management practices. Derpsch et al. [85] suggest certain steps towards standardizing research protocols, particularly when designing NT experiments. An effort to standardize most aspects of NT including weed cover, soil moisture at seeding, equipment details and crop-rotation sequences among other factors is made. Besides, it is important that all AMP are described in a detailed manner and be inclusive of relevant antecedent conditions as well as all actions taken from beginning till the end of the experiment.

3. Time Variable SHP in Soil Water Simulations

The alterations in SHP have distinct effects on evapotranspiration, groundwater recharge, water storage capacity, surface runoff and groundwater flow in a watershed [86]. The sensitivity of soils and the response of the water cycle to soil structure alterations are site-specific and driven by environmental interactions, i.e., between climate, topography and geological substrate that control soil formation [40] and intensity and direction of water fluxes. This section addresses the response of using temporally variant SHP in hydrological and land-use models.

Simulations that consider time variant SHP in hydrological modeling studies are challenging and require an adequate set of soil hydraulic parameters. Though limited in number, the studies that have considered these variations show promising results to improve the accuracy of our numerical simulations for soil water balance. Xu and Mermoud [15] used a simulation model with time-dependent hydraulic conductivity parameters to predict the effect of CT, RT and subsoiling tillage on the components of the soil water balance during the summer maize growing season. They found that using temporally variant hydraulic conductivity values improved the model’s predictive capability in estimating percolation (losses of more than 6% in wet years by subsoiling), water storage (13–14% lesser storage under subsoiling), transpiration and evaporation (lower actual evaporation in subsoiled plots).

Schwen et al. [8] studied the impacts of CT, RT and NT on SHP (derived from inverse simulation of frequent tension infiltrometer measurements) for two seasons and implemented the changes interpolated to a daily base in water balance simulations. The application of time-variable SHP significantly improved the agreement of near-surface θ and water storage in the profile to a depth of 30 cm for all tillage treatments in both seasons. The root mean square error (RMSE) values for storage was reduced by 93% while the average relative error was halved for both θ and water storage. Similarly, Alletto et al. [17] quantified the temporal dynamics of SHP in a conventional maize monoculture and included the results of the modified SHP due to CT in modeling the water dynamics of the area. The efficiency of the whole simulation period increased by 40% when using varying θ(h) and the results had a low RMSE value when compared to the simulation with constant soil parameters. The results of this study indicate that the assumption of time-invariable, i.e., constant SHP leads to a poor description of water dynamics and to an over- or under-estimation of water drainage.

An effort to include the modified SHP for rainfall-runoff simulations was undertaken by Wahren et al. [5]. The authors studied the effects of afforestation on runoff and soil water storage in foothill catchments in NE Germany for arable land, 6- and 50-year old afforestation sites and ancient natural forest. It was seen that changes in SHP were already detectable in the 6-year site. Two cases were compared in the rainfall-runoff simulations: one with constant SHP and the other with altered SHP after land-use change. It was seen that inclusion of altered SHP in modeling led to an increase in the simulated soil water storage by more than 50% and resulted in a more realistic description of the land-use effects on the soil water balance. The authors recommend considering changes in both vegetation and soil properties in rainfall-runoff models.

Finally, Bormann et al. [87] did a comparison of three hydrological models (SWAT, TOPLATS, WASIM) to assess the importance of considering changes in soil properties as part of land-use scenario calculation. PTF were used to obtain the parameters of the hydrological model. The water balance was simulated based on land-use scenarios using both changed and unchanged SHP. All three models were considerably sensitive to the changes in soil properties. It was also observed that if the effects of land-use and of an additional change in SHP are set into relation, both the changes demonstrate effects in the same order of magnitude.

Despite evidence that soil structure is subject to temporal variations and inclusion of these dynamics lead to better modeling results, SHP are set constant over time in most simulation studies. Most hydrological models that predict water fluxes under changing environmental conditions do not consider that alterations in land-use will also induce changes such as SHP [87]. Such studies usually assume only a change in the vegetation cover and overlook the dynamics of the soil structure in the process. However, there are some steps that have been taken in the direction of capturing the complex dynamics of tilled soils which are discussed in the following section.

4. Modeling Techniques to Capture Evolution of PoSD Post-Tillage

In general, few mathematical expressions were found to be appropriate for a continuous expression of the temporal dynamics of PoSD following management activities. The approach presented by Xu and Mermoud [15] used an empirical time-dependent decay function for K_s for soil water simulations. The described temporal variations of K_s were introduced in the van Genuchten-Mualem equation [36,37] to obtain time-dependent hydraulic conductivity values of partly saturated soils. As a result, the capability of the model for predicting the effects of different tillage components on soil water balance improved, particularly for subsoiling tillage.

A stochastic approach to simulate rainfall-driven soil water dynamics in the unsaturated zone using a non-linear, space domain random walk of water particles was proposed by Zehe and Jackisch [88]. Soil water particles are constant in mass and volume and their movement can be expressed by a divergent form of the Richards equation, an equivalent of the Fokker–Planck equation (FPE). Ultimately, the model accounted for the varying mobility of soil particles in different soil pore sizes. The approach was a feasible alternative to the Richards equation for simulated wetting and drying phases while also accounting for preferential and non-equilibrium flow. However, the model by Zehe and Jackisch [88] did not accurately represent long-term drainage phases in fast draining soils.

Post-tillage changes in soil PoSD tend to evolve in rather predictable pathways such as (1) loss of large pores or reduction in the mean pore size, (2) reduction in overall porosity, and (3) change in the spread of pore radii [19]. In this context, Or et al. [19] proposed to use the FPE to capture the changes in PoSD following tillage:

$$\frac{\partial f}{\partial t}=\frac{\partial }{\partial r}\left(D\left(r,t\right)\frac{\partial f}{\partial r}\right)- \frac{\partial }{\partial r}\left(V\left(r,t\right)f\right)-M\left(t\right)f$$

(2)

where f is the PoSD or frequency [L⁻¹] of pores as a function of time t [T] and pore radius r [L], D the dispersion coefficient [L²T⁻¹], V the drift coefficient [LT^-1] and M the degradation coefficient [T⁻¹]. D and V quantify the changes with time of the variance of the PoSD and mean pore radius, respectively, while M is a first-order degradation factor representing instantaneous pore loss, i.e., the fraction of pores that are lost due to instantaneous collapse. An analytical solution for the FPE was obtained by Leij et al. [9,83] when subjecting the FPE to an initial condition such as Equation (1) and a zero-probability flux in the upper and lower boundaries. The coefficients can be determined from moment analyses or from independent models. For more details on the model and its analytical solution, see [9,19,83]. The coefficients of the FPE encompass the physical processes that occur during the evolution of PoSD as a response to external stresses.

We applied this pore space evolution model and its analytical solution to results presented in two studies involving water retention data for different tillage treatments. Comprehensive measurements for the evolution of PoSD with time are limited. From the limited studies available, we selected two studies—one from USA and one from Spain—on the basis that they had data sets for temporal stages following tillage as well as for land-use change. We assessed the suitability of the model to predict the dynamics of soil PoSD due to seasonal changes following tillage as well as for evolution of PoSD when there is a change in tillage regimes. The Kosugi parameters according to Equation (1) and the r² and RMSE values for the studies are listed in

Table 2. Water retention parameter and model fits.

Study	Practice	Days after Tillage	θs	θr	σ	rm (µm)	r2	RMSE
[89]	CT (0–3 cm)	42	0.56	0.1500	1.98	30.47	-	-
		186	0.60	0.0500	1.86	11.31	0.88	0.005
		250	0.61	0.0300	1.83	8.82	0.88	0.006
		335	0.57	0.0010	1.93	6.19	0.84	0.009
	CT to NT	-	0.50	0.0010	2.18	2.43	0.89	0.01
	CT (3–6 cm)	42	0.53	0.0010	2.34	1.73	-	-
		186	0.53	0.1250	1.87	7.79	0.73	0.009
		250	0.61	0.0120	1.95	5.78	0.84	0.010
		335	0.53	0.0010	2.14	3.53	0.89	0.005
	CT to NT	-	0.50	0.0010	2.18	2.43	0.75	0.01
[42]	CT	1	0.49	0.2000	0.98	89.36	-
		62	0.50	0.2100	1.31	37.03	0.89	0.001
		124	0.46	0.2100	1.25	15.52	0.84	0.009
		220	0.44	0.2700	1.03	17.71	0.83	0.011
	CT to NT	-	0.40	0.2200	1.37	11.82	0.90	0.004

θ_s: saturated water content; θ_r: residual water content; σ: standard deviation of the log-transformed pore radius; r_m: median pore radius; RMSE: root mean square error; NT: no tillage; CT: conventional tillage.

.

4.1. Model Application: Data Sets

4.1.1. Washington, USA

Fuentes et al. [89] analyzed the temporal dynamics of SHP – both seasonal and under CT and NT—from 2001 to 2002 at sites in Washington, USA. The CT during sampling was under a spring wheat (Triticum aestivum L.) and spring pea (Pisum sativum) rotation while NT was under winter wheat. Primary tillage operations were performed in October 2000 and 2001 and soil samples were taken on 16 May (46 days after tillage) and 30 November in 2001 and on 22 April, 26 June and 18 September in 2002. Water retention curves for h = 0 to −40 cm were measured using the hanging water column method while a pressure plate extractor was used for h = −100 to −1000 cm (h = −1, −2, −4, −8, −10.5, −12.5, −30, −40, −100, −200, −300, −400, −600, −1000 cm) and fitted to the van Genuchten [37] model. We predicted the seasonal evolution of CT following tillage only after the 2001 tillage as there was only one measurement after tillage in 2000. For the change in tillage regime from CT to NT, we used the dataset immediately following tillage for CT and the corresponding measurement for NT in 2001. The predictions were done for two depths 0–3 cm and 3–6 cm.

4.1.2. Zaragoza, Spain

Peña-Sancho et al. [42] studied the dynamics of SHP under CT and NT for a fallow period (in a cereal—fallow rotation) in Zaragoza, Spain from 2011–2012. Tillage took place in March 2012. The different stages of sampling were—pre-tillage (S1), after tillage (S2), after tillage + rain (S3), late fallow (S4) and just before tillage operations (S5). Soil volumetric water content was measured in the laboratory using Time Doman Reflectometry (h = −5, −15, −30, −100, −500, −1000, −5000 and −15,000 cm) and fitted to the van Genuchten [37] model. We used the retention data for CT and NT for 10–20 cm soil depth to predict the seasonal evolution of PoSD from S2 to S5 as well as for the change in tillage regime.

4.2. Results and Model Suitability

For Section 4.1.1, Figure 2a,b show the observed and predicted frequencies for both the seasonal evolution of PoSD following CT and the PoSD evolution from CT to NT respectively for the 0–3 cm depth while Figure 2c,d show the same for 3–6 cm depth. For Section 4.1.2, Figure 2e,f show the observed and predicted frequencies for both the seasonal evolution of PoSD following CT and the PoSD evolution from CT to NT respectively for the 10–20 cm depth. Data points show f(r) values that correspond to measured retention data. In general, the results indicate that the pore evolution model is able to capture to a certain extent the evolution of PoSD for both seasonal variations following tillage and for a change in the tillage regime. The coefficients of the FPE, which were calculated from moment analyses as well as from existing models, capture the physical processes following tillage. The r² and RMSE values indicate good fits of the observed and predicted values in the pressure head range between 0 to 330 cm (r ≥ 5 µm), where the evolution of larger or inter-aggregate pores mainly take place.

For the study from Fuentes et al. [89], the data from 0–3 cm depth show more temporal variations than for the 3–6 cm depth. This is attributed to the surface layer being more prone to environmental and natural disturbances. The smaller pores (r < 5 µm) become more prevalent through the season. The same is the case when there is a shift in the tillage regime from CT to NT. For the study from Peña-Sancho et al. [42], r² value are close to 90% in the 10–20 cm depth where the soil is not directly affected by rainfall events. We also investigated the pore evolution model for this study on the surface soil at 0–10 cm (results not displayed) and we found that lower r² and RMSE values indicating poor conformation of the model to the observed values. This implies that the model may have to be calibrated for high intensity rainfalls on the surface as the coefficients may not be able to capture the rapid collapse of the soil structure in such cases.

A closer look at the results reveals that the model exhibits better fits for the change in tillage regimes when compared to the seasonal evolution of soil PoSD. The change in the mean pore radius denoted by V and the degradation term M have been studied over the years and they adequately capture the evolution of the newly formed larger pores after tillage which are structurally unstable towards a more stable PoSD [19]. However, the dispersion coefficient or the instantaneous variance D poses a challenge. Or et al. [19] suggest that though some experimental evidence indicates that D decreases with a reduction in V, there is virtually no information on the rate of these changes and their dependency on the pore size. From Figure 2, it can be observed that there is an over-estimation of the predicted values in some cases which may come down to how we estimate dispersivity. Perhaps, a new term or a factor that takes into account the time-variable nature of dispersion may improve the modeling outcomes for seasonal evolution of PoSD. Nonetheless, the model in its current state still captures the temporal dynamics of PoSD reasonably well as is seen from the r² and RMSE values (Table 2).

Ultimately, this approach is promising as it incorporates the changes in soil pore parameters directly into the governing relations of the model. It can be further extended to predict temporal changes in saturated and unsaturated hydraulic conductivities, for instance, by using the Kozeny–Carman relationship [90,91]. Calibration and assessment of the model performance will only be valid if the model is systematically applied to data sets in different locations under different soil types and AMP. However, the lack of sufficient studies so far in this regard will make evaluation of the model challenging.

5. Summary and Outlook

Surface structural pore space continuously changes on spatial and temporal scales due to a wide range of natural and anthropogenic factors. These variations also alter SHP as well as the water budget components. Herein, we reviewed the impacts of AMP and changing atmospheric boundary conditions on PoSD and consequently SHP. Most of the studies indicate that tillage is followed by creation of large pores which collapse rapidly because of rainfall events, FT and WD cycles. Moreover, there is a heterogenization (increase in the variance of the PoSD curve) of soil PoSD for both seasonal evolution of PoSD following tillage as well as when there is a change in tillage regime. The studies, however, are not consistent in their outcomes due to different methods for characterization of SHP (field and laboratory methods), different experimental designs and measurements taken at different stages of plant growth. A need for harmonization of research methodologies for generalization of results of the impacts of different management practices on SHP and the water cycle was identified.

The changes in SHP have consequences for the soil water budget components. There have been very few efforts to include time variant SHP in hydrological models. The limited studies that have incorporated the temporal dynamics have shown that it increases the model’s accuracy and reliability. In this context, we applied a sophisticated modeling approach with an analytical solution [9,19,83] that describes changes in PoSD with time and pore radius to two case studies. The model was able to capture the dynamics of soil PoSD following CT as well as when there was a change from CT to NT well. Parameterizing the model using long-term SHP measurements under different AMP and locations will help to better evaluate the suitability of the model and to calibrate its coefficients as well as improve model performance, especially for the seasonal evolution of PoSD. To our knowledge, the approach by Or et al. [19] has not been applied to a complex time series of measured hydraulic parameters in a systematic manner. The combination of the soil pore space evolution model with predicting changes in the soil hydraulic conductivity function has also not yet been tested against measured data. Future research on these topics may pave the way to include time variant SHP in hydrological modeling studies.

There arise several questions with respect to inclusion of temporal variations of SHP in measurement strategies and hydrological modeling. How often or frequently should we take measurements to account for temporal variations in SHP? How long will it take for the soil structure to reach a quasi-static equilibrium after a land-use change or management practice? How would we separate the changes induced by a combination of agricultural management practices and how can we standardize the measurement process? How complex and feasible will a process-based model be that takes into account the interaction between different environmental and anthropogenic factors to account for temporal variation of SHP? To answer these questions using modeling tools, the limiting factor is not in the theoretical part but rather the lack of adequate soil structural and hydrologic data.