Proposal of method for calculating the volume of drying cracks in the soil

Clays and clay minerals are an important and permanent component of the natural environment. Their crystal-chemical properties react very sensitively to the hydrothermal environment and its physico-chemical properties. Their presence in the soil environment causes volumetric changes. Volume changes are three-dimensional process. They are accompanied by the fissures formation in horizontal plane and by the movement of the soil surface in vertical plane. Dynamic two-domain soil structure is created while fissure formation. It consists of the fracture domain and the soil matrix domain. In the natural space, desiccation cracks are located in the unsaturated zone of the soil. They significantly influence the rainfall-runoff process. Direct volume measurement of the fracture network is an extremely difficult task. Calculation procedure of fracture network retention properties is designed in the paper. It is based on the assumption that the cracks volume is dependent on the clay particles content in the soil and on the soil moisture volume. From a material point of view, the soil profile is divided according to the volume changes potential. In the calculation, the isotropy of the soil profile is assessed based on the geometric factor. Shrinkage characteristics, cracks geometric dimensions on the soil surface and volumetric moisture along the vertical of the soil profile are included in the calculation. It is advantageous to use pedo-transfer functions in the calculation.


Introduction
The water regime of heavy soils is significantly different from the water regime of rigid soils.The difference lies in the ability of heavy soils to swell and shrink (Novák, 1978).These processes are caused by soil volume changes, which are manifested by the formation of drying cracks and by the vertical movement of soil surface.The quantification of the mentioned processes and determination of the characteristics have a common basis in the measurement of heavy soils volume changes as a function of its moisture and clay mineral content (Amorim et al., 2007;Chen et al., 2022;Crescimanno et al., 1995).Volume changes research depending on humidity does not present special difficulties.The situation is different in connection with the volume changes quantification depending on the type and content of clay minerals.The ability of individual minerals to assimilate water into their structure and therefore actually increase their volume is different (Johnston, 2018;Velde, 1995).Clay minerals of the illite and montmorillonite groups have the greatest ability in this direction (Ross, 1978).Clay minerals do not occur in pure form but form mixed structures in natural conditions (Hovorka, 1994;Šucha, 2001).These are made up of different types of clay and other minerals.Term "clay" used in the literature is identified with the clay minerals.The term "clay" is used to denote material whose particles are smaller than 0.002 mm (Tall, 2002;Igaz et al., 2020) and also for material with a dominant representation of clay minerals (Velde, 1995).The need to know what types of clay minerals are present in the locality of interest is significant.The identification of types of clay minerals in soils is very difficult.The way out is through the quantification of the clay content without claims to distinguish their individual types.Dried cracks in heavy soils with a high content of clay particles are the main factor determining their hydraulic properties (Čubanová et al., 2022).The soil profile of heavy soils forms a two-component system consisting of cracks (i.e.macropores) and the soil matrix (which includes micropores) at the time of occurrence of desiccation cracks.Each of these components is characterized by diametrically different conditions for the transfer and retention of water (Vitková et al., 2022;Toková et al., 2023).When a large amount of water enters the soil in a short time and the cracks have not yet had time to close completely, rainwater moves in the cracks analogously to water in open gutters.This movement has an order of magnitude higher speed than the movement of water through micropores (Richards, 1931;Ross, 1990;Šimůnek et al., 2017;van Genuchten et al., 2000).The fissure system serves as a system of preferred paths for horizontal (horizontal conductivity of fissures) but also for vertical flow of water (vertical conductivity of fissures) (Chertkov, 2001;Chertkov, 2002).The geometry of cracks changes over time and the hydraulic parameters of the preferred paths change also.A free level can exist in cracks only if they are open.The presence of water in the fissures causes a rapid soil moisture rise and subsequent closing of fissures.It follows that the movement of water in cracks is limited by the speed with which the soil reacts to changes in humidity.When investigating the two-domain soil environment, it is necessary, in contrast to the soil environment of rigid soils, to additionally know and quantify phenomena that are conditioned by volume changes.These include spatial and temporal definition of fracture network geometry, investigation of vertical movements of the soil surface, determination of shrinkage characteristics, shrinkage potentials, flow into cracks, flow inside cracks and flow from cracks into the soil matrix.These facts point to the originality of these processes and their characteristics, implying that it is difficult to transfer experimental studies between different sites.The aim of the presented work is to propose a calculation method for determining the geometric characteristics of cracks, based on which it will be possible to calculate their volume in the investigated locality.Cracks volume calculation is based on the calculation and analysis of volume changes.The calculation is based on the assumption that the volume of cracks is dependent on the magnitude of volume changes with respect to the saturated state.The magnitude of volume changes depends on the content of clay particles in the soil and on the volume moisture of the soil (Beven and Germann, 1982).The calculation method is designed so that the range of input data is as small as possible and that it can be obtained in the simplest possible way (Čurlík and Bedrna, 1971).

Initial assumptions and input data
The proposed calculation procedure is based on a hypothesis based on the following assumptions: -Cracks are created due to volume changes.
-Volume changes are caused by changes of water content in soils with a high content of clay particles.-Volume changes are a three-dimensional process.-In the space, the desiccated cracks are located in the water-unsaturated soil zone between the soil surface and the groundwater level.
The requirements for the structure and extent of the database are based on the above assumptions.From this point of view, the database can be divided into: a) Information on the hydrophysical characteristics of the soil profile.Lubrication characteristics, soil profile isotropy in terms of volume changes, texture.b) Information on volume moisture along the vertical of the soil profile.c) Information on the geometric characteristics of the fracture network on the soil surface.

Calculation procedure
The volume of cracks in the soil profile most often refers to the volume   =    , where S is the horizontal crosssection of the soil and zl is the height of the investigated soil monolith.In natural conditions, it is the thickness of the unsaturated soil layer.The formation and development of a fracture network is a complex process (Brown, 1977;Pekárová et al., 2022;Chertkov and Ravina, 2004).It is a dynamic system whose geometric dimensions, and hence volume, change over time (Chertkov, 2000).The volume of cracks is calculated layer by layer in the proposed calculation procedure.For the investigated soil monolith, the equation of the volume balance of cracks Vc in a unsaturated environment at time ti can be expressed in the form (Fig. 1): Vc,iis the volume of cracks in the i-th calculation layer, while if  , ≤  ⟹  , = 0 where ε is the chosen value of the crack width, which can already be neglected.It is the width ε at which gravity and capillary forces are in balance. , is the volume of cracks in the investigated soil profile in the i-th time interval.
All terms of equation ( 1) are related to the selected time ti.If the volumes of cracks in the soil profile are related to the time interval ∆ =   −  −1 then the time course of the volume of cracks Vc,t during k-time intervals is: The calculation process of the volume and depth of cracks along the vertical of the soil monolith can be divided into the following steps.1) Measurement of the width and length of the fracture network on the soil surface.2) Division of the studied monolith into material layers.
3) Determination of shrinkage and isotropic characteristics and humidity ratios in individual material layers.4) Breakdown of individual material layers into computational layers.5) Calculation of the volume of cracks in individual calculation layers.6) Calculation of the volume of cracks in the entire monolith.

Measurement of the width and length of the fracture network on the soil surface
At the beginning of the calculation period, the length of cracks Lc is measured (Novák, 2000).For simplicity, it is assumed that the length of cracks does not change over time and only their width  , changes.If the length changes, it is necessary to adjust the frequency of measurements according to the variability of the cracks.The width of cracks  ,1 on the soil surface can be measured during the entire calculation period.The measurement of the width of the cracks is used for the verification of the calculation and for the estimation of the calculation error of the cracks volume.
In the calculation, the thickness of the calculated layers zc >10 cm is not assumed, and for that reason the width of the cracks in the individual calculated layers is assumed to be  , = .Tortuosity is neglected along the vertical of cracks (Fig. 1).

Material layers
It is necessary to determine the number of homogeneous soil layers in the studied monolith.The homogeneity of the soil layer is assessed on the basis of the potential for volumetric changes, which is characterized by the coefficient of linear extensibility (COLE) (Yan et al., 2022).The latter is defined by the equation: where Vsis the volume of soil in a saturated state, Vdryis the volume of dried soil.
If direct COLE measurements are not available, this parameter can be estimated using the clay content in the soil, i.e. from the content of particles with fraction   < 0.002 .

𝐶𝑂𝐿𝐸 = 𝑓(𝑓𝑟 𝑐 )
[-] (4) For the conditions of the East Slovak Lowland (ESL), this dependence was expressed in the form: Moreover, if frc is expressed as a proportional number then theoretically   ∈ 〈0.16; 0.7〉 or in the case of a percentage expression   ∈ 〈16%, 70%〉.It follows from the above in East Slovak Lowland (ESL) conditions COLE ∈ 0.00; 0.16.To divide the soil profile into material layers according to COLE, it is possible to use the division according to (Parker et al., 1977) extended by one category for extremely high shrinkageswelling potential (COLE>0.12)(Table 1).

Determination of shrinkage and isotropic characteristics and humidity ratios in individual material layers
Soil volume changes are a three-dimensional process, manifested in the horizontal plane by the formation of a fissure network and in the vertical direction by the movement of the soil surface in the course of moisture changes (Bronswijk, 1988;Bronswijk, 1990).
The mutual relations between these elements can be mathematically formulated on the basis of a simple reasoning.The saturated cube-shaped soil element has dimension zs and volume Vs (Fig. 2).After isotropic shrinkage by Δz, the soil element has dimension z and volume V.  Then the volume of the saturated soil element is   =   3 , the edge of the cube after shrinkage is z = zs -∆z, and the volume of the cube after shrinkage is V = z 3 .Based on the above, equations ( 6) and ( 7) can be formulated.
where Vris the change in total volume during drying compared to the saturated state.
The exponent in equation ( 7) has a value of 3 for isotropic shrinkage.In general, the specified exponent can take values from the interval ⟨1, ∞).The exponent is denoted rs and is called the geometric factor.By combining equations ( 6) and ( 7) and introducing rs, it is possible to formulate equation ( 8) and then equation ( 9) for calculating the vertical movement of the soil surface.
The total volume change Vr is formed by the volume Vv which is projected into the vertical movement and the volume Vh which is projected into the horizontal plane.Then equations ( 10) and ( 11) apply.
By combining equations ( 10) and ( 11), relation ( 12) is obtained for calculating the volume that is reflected in the formation of soil cracks.
The geometric factor   can be expressed from equation ( 8) in the form of equation ( 13).It is possible to measure the values of V, Vs, ∆z, zs and thereby calculate the value of the geometric factor rs according to equation ( 13).12).After taking an undisturbed soil sample from the material layer, calculating rs and measuring the volume in the saturated state Vs, the variable V is unknown in equation ( 12).Equation ( 14) applies to V.
In equation ( 14), it is necessary to know the variable Vr for different volumetric moistures that may occur in the given material layer.In other words, each volumetric humidity corresponds to a Vr value.For this reason, it is necessary to introduce the variable  = 100     ⁄ and to know its course according to equation ( 15).

𝑑 = 𝑓(𝜃 𝑣 )
[%] ( where dis the shrinkage characteristic that expresses the total change as a function of the volume humidity for the given texture, expressed in [%] of the saturated volume Vs.
It is advantageous to express the measured empirical dependence analytically.The solution in heavy soils often leads to a second degree polynomial of the form: In that case, theoretically  ∈ 〈0, 〉, where c is the parameter from equation ( 16).Soil moisture monitoring does not currently represent a significant problem.There are a number of methods for determining it, often with remote data transmission (Kuban et al., 2022).

Computing layers
The thickness zc,i and thus the number of "n" calculation layers depends on the number of material layers and on the density of volumetric moisture measurements along the vertical soil profile (Fig. 1).

Calculation of the volume of cracks in individual calculation layers
At the beginning of the calculation, the measured length of the cracks on the surface of the soil monolith and the volumetric soil moisture θv are entered.Volumetric humidity is the input data with the greatest variability and impact on the dynamics of the crack network for a given texture.The influence of θv is reflected in the volume of cracks through equations ( 15), ( 16) and ( 12).In next calculation, it is convenient to express d as a relative number n in the form: If the F is the volume of saturated soil sample in [%] than it applies that F = 100%.Subsequently, the volume of shrinken soil sample B can be expressed as B = Fd or B = Fn F. If the volume changes in horizontal plane expressed in [%] of the F volume are designated as C than the equation ( 12) is: For simplicity, C is expressed as a relative number m in the form: If the volume of the calculation layer of the investigated soil monolith is Vl with the area  =     where   ,   ,   are the dimensions of the calculation layer (Fig. 1) then: where Vcis the volume of cracks in the calculation layer, Scis the area of cracks in the horizontal section of the calculation layer then: where wb-is the average width of cracks in the calculation layer, where Lcmis the specific length of cracks, Scmis the specific area of cracks, Vcmis the specific volume of cracks.

Calculation of the volume of cracks in the entire monolith
After calculating the volumes of cracks in all calculation layers Fig. 1 in the i-th time interval and substituting them into equation ( 1), the value of the total volume of cracks in the examined soil monolith is obtained.The time course of Vc,t in k-time intervals is expressed in equation (2).By analogy to the single calculation layer (26), it is also advantageous to calculate the total volume of cracks per unit area.The dynamics of fissures is influenced by the increased infiltration surface, which is formed by the side walls of the fissures.Surface runoff occurs only after the fissures are filled with water.However, the water in the cracks quickly infiltrates through the side walls and the cracks close.If the number of calculation layers is n, then the thickness of the monolith zl is (Fig. 1): If n is the last calculation layer where   = 0 then the crack depth ℎ  =   −  , and the crack infiltration area can be expressed by equation (28).
Also in this case, as in the case of one layer (25), it is advantageous to convert Sinf into a unit of area.
In precipitation events with a total rainfall P, part of the precipitation Pc falls directly into the crack through the area Sc,1, which is in the surface calculation layer.
After converting to unit area: where Pcmis the precipitation that fell on the area unit directly into the crack.

Discussion
The calculation was designed on the basis of our own knowledge obtained during research at ESL and knowledge published in domestic and foreign literature.
The presented calculation method is intended for expert estimation of the volume of the fracture network and other related parameters.These are the infiltration area through the side walls of the cracks (Novák and Šoltész, 1994), the width of the cracks and its change along the vertical of the soil profile, and the depth of the cracks.
The model makes it possible to estimate the amount of precipitation that will fall directly into the fracture network.Knowledge of the retention volume of the fracture network will allow to assess their influence on runoff conditions during sudden, intense rainfall.For the model, it is necessary to sensitively determine the shrinkage characteristics and volume moisture along the vertical of the soil profile.The benefit is that the isotropy of volume changes and the distribution of the soil profile according to the potential of volume changes are taken into account.It is important that the results of calculations can be verified on the width of cracks that can be monitored.
When designing the calculation method, there were some simplifications that may be a source of inaccuracies.The tortuosity of cracks along the vertical soil profile was neglected (Chertkov and Ravina, 1999).A constant crack width is assumed in the calculation layer.The calculation process can begin only when a measurable fracture network is developed in the investigated location.

Conclusion
The presence of clay minerals in soils causes their volume changes.The manifestation of volume changes during drying is the formation of a crack network.This significantly changes the hydrodynamics of the soil profile.The retention capacity of the fracture network affects rainfall-runoff relationships.Measuring the volume of cracks in the natural environment is very difficult.The presented calculation method for calculating the volume of cracks was developed for that reason.The calculation method is based on the analysis of volume changes, shrinkage characteristics and humidity regime of the investigated area.Simplifying assumptions were introduced in the proposed calculation to simplify the input data.Nevertheless, it is assumed that the modelled results will be in sufficient agreement with the verified ones.
value of this factor, the isotropy of volume changes of each material layer in the soil monolith is assessed.Theoretically,   can have the following values:   = 1 shrinkage cracks are not formed; the entire volume change is reflected in vertical movement; 1 <   < 3 vertical movement predominates over formation of cracks;   = 3, isotropic shrinkage;   > 3, cracks formation predominates over vertical movement;   → ∞, the entire volume change is reflected in the formation of cracks; there is no vertical movement.The quantification of the shrinkage characteristics for the calculation of the crack volume is based on equation (