# Machine learning of clod evolution under rain for numerical simulation of microtopographic variations by clod layout

Edwige Vannier, Richard Dusséaux

► To cite this version:

Edwige Vannier, Richard Dusséaux. Machine learning of clod evolution under rain for numerical simulation of microtopographic variations by clod layout. Biosystems Engineering, 2025, 253, pp.104130. 10.1016/j.biosystemseng.2025.104130. hal-05022989

HAL Id: hal-05022989

https://hal.science/hal-05022989v1

Submitted on 7 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Research Paper

# Machine learning of clod evolution under rain for numerical simulation of microtopographic variations by clod layout

E. Vannier* ID, R. Dusséaux ID

LATMOS/IPSL, Université de Versailles Saint-Quentin-en-Yvelines, Paris-Saclay, 11 boulevard d'Alembert, 78280, Guyancourt, France

<figure><img src="image_2.png" /></figure>

## ARTICLE INFO

### Keywords:

- Modelling

- Soil clods

- Digital elevation model

- Estimation

- Evolution under rainfall

## ABSTRACT

Soil surface roughness (SSR) is shaped by tillage operations and evolves with weathering. It is related to geomorphologic processes and to soil fertility. Thus it is an input of various models and the object of many studies. Numerical generation of rough surfaces is an alternative to measurement, which can be cumbersome.

This paper proposes a way to numerically generate soil surfaces resembling seedbeds with roughness that have evolved due to rain. As SSR is related to clod-size distribution, the principle is to set modelled clods on a planar surface, as a first approximation of the surface. An experiment was designed to get controlled surface roughness by setting pre-sieved clods on a nearly horizontal surface of loose soil, and then subjecting the surfaces to rainfall events performed by a rainfall simulator. Digital elevation models (DEMs) of each state of the surfaces were recorded by laser scanner to monitor the evolution of clods under rain. Clods were segmented and matched to form a data base of individual clods at each state of the surface. The evolution of clods under rain was modelled by Machine Learning. A set of DEMs was used for learning, the other for test, and several metrics were applied. Obtained results show the robustness of the model. This simple surface representation captures the main properties of the surface. Small scale surface is useful for various applications, such as rough surface scattering, and more generally, modelling where soil surface is an input.

## 1. Introduction

Conservation of water and soil resources requires an interdisciplinarity engagement to study and understand processes and their spatiotemporal dynamic in order to inform choices in agricultural practices. Agricultural soils consist of several levels of roughness, due to superposition of micro variations at grain scale, random small scale structure (clods and depressions), almost deterministic surface of greater scale (furrows), and some low undulation of even greater scale representing the landscape (also quasi-deterministic) (Takken, 2000; Vinci et al., 2020). Several definitions of surface roughness co-occur and it designs the variability of a topographic surface at the scale of the geomorphologic features of interest (Grohmann et al., 2011; Smith, 2014). At centimetre to metre scale, it is the widely called random roughness. At the scale of metres to tens metres it is the so called topographic roughness (Doane et al., 2024). Hereafter, Soil Surface Roughness (SSR) represents the undulation of the surface at small scale, due to the presence of small elevations and depressions. It results from tillage operations and changes over time due to weathering. Therefore, SSR can be related to clod-size distribution (Ajdadi et al., 2016; Arvidsson &

Bölenius, 2006; Bogrekci & Godwin, 2007; Garbout et al., 2013; Ishii et al., 2023; Jensen et al., 2016; Sandri et al., 1998; Schapel et al., 2019; Vannier et al., 2018a; 2018b; Wang et al., 2011).

The role of SSR in the geomorphologic processes at small scale is commonly acknowledged, for example Garbout et al. (2013); Thomsen et al. (2015); Martinez-Agirre et al. (2016); Gilliot et al. (2017); Schapel et al. (2019), Dusséaux and Vannier (2022), Wang et al. (2023). Millimetre-scale is significant for agricultural production because it determines soil erosion, deposition and water infiltration, along with seed germination and plant growth (Zhao et al., 2024). Centimetre scale is suitable for the radar signature of soils in L, C and X bands, for surveying soils by remote sensing (Dusséaux et al., 2012; Franceschetti & Riccio, 2007, chap 1 and 4; Ogilvy, 1991; Pan et al., 2016; Verhoest et al., 2008; Zheng et al., 2019).

The interpretation of radar measurements is a difficult task due to geometrical complexity of agricultural soils and moisture influence. Therefore, electromagnetic simulations of radar signal are an essential step to understand the influence of roughness and moisture on the backscattered signal. Rough surface scattering is the purpose of many studies including remote sensing application. The scattered field depends on the incident wave, the medium properties and the scatterer.

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: edwige.vannier@latmos.ipsl.fr (E. Vannier).</footnote>
Nomenclature

DEM Digital Elevation Model

SSR Soil Surface Roughness

x,y,z Coordinates in Cartesian system (mm)

1D One dimensional

2.5 D Two dimensional and a half for z(x,y) data

P Cumulated precipitation (mm)

O Barycentre of a clod

R = (Oxy) Barycentric coordinate system

Re = (OXY) Rotated coordinate system

a Semi-major axis length

b Semi-minor axis length

c Semi-height

θ Orientation angle

I<sub>Ox</sub> Second order horizontal moment

I<sub>Oy</sub> Second order vertical moment

I<sub>xy</sub> Total inertia

m<sub>1</sub>, m<sub>2</sub> Semi-axes length

u, v Angles for parametric equation

HES<sub>p</sub><sup>(k)</sup>(x,y) Half-ellipsoid surface

k State of a surface

p Surface number

HES<sub>lp</sub><sup>(k)</sup>(x,y) Set of learning surfaces

HES<sub>lp</sub><sup>(k)</sup>(x,y) Set of test surfaces

C<sub>p</sub><sup>(k)</sup>(x,y) Bi-directional autocorrelation function

(x<sub>i</sub>,y<sub>j</sub>) Sampling points

N<sub>p</sub> Number of samples

ε<sub>p</sub><sup>(k)</sup>(x<sub>i</sub>,y<sub>j</sub>) Error function

MSE Mean Square Error

MSE(C<sub>p</sub><sup>(k)</sup>) Mean Square Error of error function on the autocorrelation function C<sub>p</sub><sup>(k)</sup>(x,y)

V<sub>n</sub><sup>(k)</sup> Normalised volume of clods

V<sub>i</sub><sup>(k)</sup> Volume of clod i

MSE(V<sub>n</sub>) Mean Square Error of normalised volume

α Amplitude parameter

β Decrease parameter

Smooth planar scatterers will concentrate the scattering energy at the specular angle, whereas roughness can influence the directivity of the scattered signal spreading to other angles (Choi et al., 2018; Ogilvy, 1991). The simulation of the backscattered signal can be performed by applying an exact or analytical electromagnetic model to rough surfaces presented as input. Within their domains of applicability, the analytical models allow a fast analysis of the electromagnetic signature by means of analytical formulae (Ogilvy, 1991). Exact methods associated with numerical techniques and Monte Carlo simulations demand a large number of rough surfaces (Bourlier et al., 2013). Numerical generation of rough surfaces is a widely used alternative to surface measurement, which is cumbersome (Ogilvy, 1991; Marchand & Brown, 1999; Franceschetti & Riccio, 2007, chap 1 and 4; Lievens et al., 2009; Dusséaux et al., 2012; Pan et al., 2016; Choi et al., 2018; Bykhovsky & Lyandres, 2018; Zheng et al., 2019; Dusséaux & Vannier, 2022). Small scale numerical surfaces are useful for various applications, such as rough surface scattering, and more generally, modelling where soil surface is an input. And so far, little attention has been paid to surface evolution due to weathering.

In remote sensing, the soil surface is generally modelled by its autocorrelation function, and a Gaussian distribution of elevations. It can be also modelled by geostatistical parameters (Grohmann et al., 2011; Trevisani & Rocca, 2015; Trevisani et al., 2023). Another approach is to model clods or depressions (Borselli & Torri, 2010; Cierniewski et al., 2002; Kamphorst et al., 2005; Vannier et al., 2014; Zribi et al., 2010). This approach is pertaining to Boolean algorithms, in which simple patterns are set randomly onto a background (Koltermann & Gorelick, 1996). The half circle (for 1D profile) or half sphere (for 2.5D surface) is commonly used. In Vannier et al. (2014), the clod base is modelled by an ellipse and the height by a semi-arch of cosine. In Barbosa and Gerke (2022), clods were modelled by cylinders to study clod friability. In Forooghi et al. (2017), the roughness geometry is generated by distributing roughness elements of random size and prescribed shape, with disc base and various profiles of heights, on a flat surface to reproduce sand roughness.

The present study follows this approach of surface numerical generation by setting objects of desired geometry onto a plane. An experiment was conducted to get different surfaces by setting pre-sieved clods on a nearly horizontal surface of loose soil and subjecting them to controlled rainfall. The objective is to reproduce numerically the

evolution of clods by Machine Learning. Machine Learning gathers several kinds of problems, mainly, prediction, detection, segmentation, partition and prediction problems are usually addressed in terms of regression, classification and clustering (Géron, 2022; James et al., 2021). If the prediction is to be a quantitative value, it is a problem of regression, whereas if the prediction is to be a label, it is a problem of classification for given labels, or of clustering for unknown labels. Machine Learning consists in designing an adapted model by statistical learning on a data set. The approach is to divide the data into a learning set used to design the model and estimate its parameters, and a test set used to validate the model. In an attempt to generate numerically cloddy surfaces that have evolve under rainfall, our goal is to segment the clods on the recorded surfaces, represent each clod by its geometrical model, and learn the evolution of the parameters of the geometrical models of the various clods. Given the correlations between parameters, the linear regression has proved to be suitable. The proposed global model of clod evolution has a deterministic part and a random part.

## 2. Data of soil surfaces

### 2.1. Surface preparation

The database used in this paper is made of two laboratory soil surfaces composed in trays with loose soil of silty loam structure (11 % clay, 60 % silt, 29 % sand) topped by pre-sieved clods collected in the field. Clods were first air-dried and then set randomly upon a nearly horizontal surface of soil, in order to look like a seedbed except for the furrows. The first tray and the second tray have standard deviations of elevations of 10.5 mm and 9 mm respectively. Clods stand out well from the soil support and are isolated on tray 1 and more gathered on tray 2. The clods of tray 2 are also generally smaller than those of tray 1 and are 1.5 times more numerous (see Fig. 1). The soil surfaces were eroded with a laboratory rainfall simulator similar to the one described in Foster et al. (1979) also described in Vannier et al. (2018a, 2018b), Vannier et al., 2014. Controlled rainfalls were performed with intensity of 33 mm h<sup>-1</sup> to reach swelling and then 42 mm h<sup>-1</sup> to see flattening and crusting. A millimetre DEM was recorded at each state of the surface with the laser-scanner described in Darboux and Huang (2003) and Vannier et al. (2018a, 2018b), Vannier et al., 2014. This laser-scanner had a grid spacing of 0.5 mm in x and y, and of 0.1 mm in z.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Photographs of two soil surfaces at different states: a) tray 1, and b) tray 2.</figcaption></figure>

Accuracies in x and y are 0.5 mm, but, the vertical accuracy in case of high roughness is only 1 mm. Thus, the elevations were re-scaled at a 1 mm sampling interval. And, as we adopted a geometrical approach to model the clods, it is important to have the same sampling interval in all directions, so a resampling was done also in x and y. Hence, the sampling is 1 mm in x, y, and z. Initial state is the state of original soil surfaces. State 2 is the state after a rainfall of 120 min at 33 mm h⁻¹, state 3 to 5 are the states after rainfalls at 42 mm h⁻¹ of respectively of 38 min, 51 min and 90 min. So the cumulated precipitation *P* for the 5 states is 0, 66 mm, 93 mm, 129 mm and 192 mm.

## 2.2. Surface segmentation

Wavelet-based clod segmentation (Vannier et al., 2019) was used to delineate individual clods on the DEMs. The result of segmentation can be seen in Fig. 2. A clod to clod matching was performed among states of the surfaces in order to form a data base containing the individual clods at each state of the surface. Each clod volume was calculated by summing the elevations within the clod contour, and a survey of clod volumes was done to select only stable individual clods by control of

volume variation. If some clods merged or split under rainfall, they were discarded from the data base. Also, the clods located too close to the edges of the trays were discarded because of unstable segmentation. As a result, 172 out of the 185 clods from the first tray were recorded in the database as the initial state, while for the second tray, 241 out of 288 clods were retained.

## 2.3. Clod modelling

Clods were modelled by adjusted half ellipsoids as in Chimi, 2012. The principle is developed hereafter.

### 2.3.1. Equivalent ellipse modelling clod contour

The equivalent ellipse is the best fitting ellipse of a clod contour. It has three parameters: semi-major axis length *a*, semi-minor axis length *b* and orientation angle *θ*. They are estimated by minimization of the moments of inertia as in Feral et al. (2000). Let *O* be the barycentre of the considered clod, and *R* = (*Oxy*) its barycentric reference, with unitary vectors *e*<sub>x</sub> and *e*<sub>y</sub>. The orientation angle *θ* is the angle between the major axis of the ellipse, supported by (*Oxy*), and the horizontal axis

<figure><img src="image_4.png" /><figcaption>Fig. 2. Segmentation of clods on DEMs for two soil surfaces at different states: a) tray 1, and b) tray 2.</figcaption></figure>
supported by (Ox). The central second order moments in the barycentric reference are given by:

$$ 
\left\{ 
\begin{array}{l} 
I_{0x} = \frac{1}{N} \sum_{i=1}^{N} y_i^2 \\ 
I_{0y} = \frac{1}{N} \sum_{i=1}^{N} x_i^2 \\ 
I_{xy} = \frac{1}{N} \sum_{i=1}^{N} x_i y_i 
\end{array} 
\right. 
\qquad (1) 
$$

where the summation extends over all the pixels located within the clod contour.

Now, the inertia principal moments $I_{0x}$ and $I_{0y}$ can be expressed in the principal axis (OX) and (OY), which are supported by unitary vectors $\vec{e}_X$ and $\vec{e}_Y$:

$$ \vec{e}_X = \cos \theta \vec{e}_x + \sin \theta \vec{e}_y \quad \text{and} \quad \vec{e}_Y = -\sin \theta \vec{e}_x + \cos \theta \vec{e}_y \qquad (2) $$

$$ I_{0X} = f(\theta) = I_{0x} \cos^2 \theta + I_{0y} \sin^2 \theta - I_{xy} \sin 2\theta \qquad (3a) $$

$$ I_{0Y} = g(\theta) = I_{0x} \sin^2 \theta + I_{0y} \cos^2 \theta + I_{xy} \sin 2\theta \qquad (3b) $$

$I_{0X}$ and $I_{0Y}$ minimize the f and g functions for:

$$ \tan 2\theta = \frac{2I_{xy}}{I_{0y} - I_{0x}} \qquad (4) $$

Then the orientation angle $\theta$ and the length of the semi-axes (a, b) are calculated as follows:

$$ \theta = \frac{1}{2} \tan^{-1} \left( \frac{2I_{xy}}{I_{0y} - I_{0x}} \right) \pm k \frac{\pi}{2} (k \in \mathbb{N}) \qquad (5) $$

where only one orientation angle is kept so that $\theta$ varies within the open interval $]-\pi/2, \pi/2]$, and:

$$ m_1 = 2\sqrt{I_{0x} \cos^2 \theta + I_{0y} \sin^2 \theta - I_{xy} \sin 2\theta} \qquad (6) $$

$$ m_2 = 2\sqrt{I_{0x} \sin^2 \theta + I_{0y} \cos^2 \theta + I_{xy} \sin 2\theta} \qquad (7) $$

$$ \left\{ 
\begin{array}{l} 
a = \max(m_1, m_2) \\ 
b = \min(m_1, m_2) 
\end{array} 
\right. 
\qquad (8) $$

### 2.3.2. Half-ellipsoid modelling entire clod

Parametric equation of half ellipsoid in barycentric reference is given by:

$$ x = a \cos(\nu)\cos(u) \qquad (9a) $$

$$ y = b \cos(\nu)\sin(u) \qquad (9b) $$

$$ z = c \sin(\nu) \qquad (9c) $$

with $u \in ]-\pi, \pi[$ and $\nu \in [0, \pi/2]$, a and b, semi-axes of the equivalent ellipse and c, height of the half-ellipsoid. In our case, the height c is estimated by matching the volume of half-ellipsoid to the volume of clod.

## 2.4. Modelled surface

The modelled surface is obtained by representing each individual clod as its adjusted half-ellipsoid and positioning the half-ellipsoid at the clod's location on a planar surface. It will be called half-ellipsoid surface in the following.

A half-ellipsoid surface is denoted $HES_p^{(k)}(x,y)$. Subscript p numbers

surfaces, $p=0$ is reference clod surface estimated from DEM (see Fig. 3), and $p>1$ designates numerical evolution of initial state for simulation p. At last, $k$ is the state of surface. A subscript l or t is added to designate the learning data or the test data. The clods of the first tray were used for learning and that of second tray for test. This leads to two sets: $\{HES_{lp}^{(k)}(x,y)\}$ and $\{HES_{tp}^{(k)}(x,y)\}$.

## 2.5. Learning of clod evolution

This model of clod evolution has two components, a deterministic part and a random part. The global algorithm of surface evolution can be seen in Fig. 4.

For clods present since initial state, the basic trend of parameters of half-ellipsoid a, b, c and $\theta$ was estimated by linear regression from one state towards the next, and $\theta$ was bounded in the $]-\pi/2, \pi/2]$ interval (see Table 1). An assumption that the centres of gravity of clods do not change was made for convenience. Then, probabilities for clods to disappear or to appear at state k were found size-dependent and inferred from the histograms of semi axis a.

In case of appearance, a value of semi-axis a was picked randomly, and then the correlated parameters b and c were predicted from a by linear regression with noise, and the last parameter $\theta$ was set randomly within proper range. A position of the half-ellipsoid was then set randomly on the plane where it was possible without overlay. Indeed, until state 4 of the surface, about 2-3 % of the clods may create a new small clod by splash detachment. The random positioning is a simplification.

A clod would disappear at state k either if semi axis a is in within some small span, with a probability that is greater the smaller a is (see Fig. 5), or if height c is getting null. Therefore, there are two causes of disappearance, hence the probabilities from a needed to be scaled down, and the basic trend of c relearned.

For learning, the evolution of clods was validated by comparing the number of clods and the maximum of autocorrelation function at each state of the surface. Parameter c has been tuned to get a balance. Indeed, deterministic linear regressions for evolution of c were shaped by the bigger clods, so the intercepts were too much for small clods that disappeared too fast. Table 1 shows the relearned regressions of parameter c.

# 3. Metrics of goodness

## 3.1. Autocorrelation function

The bi-dimensional autocorrelation function $C_p^{(k)}(x,y)$ of the surface $HES_p^{(k)}(x,y)$ measures the similarity between the function $HES_p^{(k)}(x,y)$ and its lagged copy along the directions x and y as a function of the two lags. It is a characteristic of surface widely used in remote sensing modelling (Dusséaux et al., 2012). Let us notice that the variogram, which provides geostatistics, was also created to describe the spatial correlation between observations, and is broadly used in Earth Science (Chilès & Delfiner, 1999, p. 695). In discrete form, the autocorrelation function is computed at each sampling points $(x_i, y_j)$ and scaled by the total number of samples $N_p^2$ (see Dusséaux & Vannier, 2022). It is thus the biased estimate. Therefore, the value at the origin $C_p^{(k)}(0,0)$ is the maximum of the autocorrelation function and is the sample variance.

For measuring the proximity between a surface and its reference based on autocorrelation function, the error function $\epsilon_p^{(k)}(x_i, y_j)$ is computed at each sample point $(x_i, y_j)$,

$$ \epsilon_p^{(k)}(x_i, y_j) = \frac{C_p^{(k)}(x_i, y_j) - C_0^{(k)}(x_i, y_j)}{C_0^{(k)}(0, 0)} \qquad (10) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3. Half ellipsoid surface for two soil surfaces at initial state: a) tray 1, and b) tray 2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Algorithm of surface evolution.</figcaption></figure>

**Table 1**

Evolution of initial clod parameters as a function of previous state: A matrix of semi-major axes, B matrix of semi-minor axes, C matrix of semi-heights, and Θ matrix of orientation angles.

<table><thead><tr><th>Matrices of semi-major axes</th><th>Matrices of semi-minor axes</th><th>Matrices of semi-heights</th><th>Matrices of orientation angles</th></tr></thead><tbody><tr><td>A(2) = A(1)×<br>1.007 + 2.107</td><td>B(2) = B(1)×<br>1.071 + 0.596</td><td>C(2) = C(1)×<br>0.985-1.5</td><td>Θ(2) = Θ(1)×<br>1.207 + 0.004</td></tr><tr><td>A(3) = A(2)×<br>0.996-0.321</td><td>B(3) = B(2)×<br>0.969 + 0.132</td><td>C(3) = C(2)×<br>0.97-1.3</td><td>Θ(3) = Θ(2)×<br>0.945 + 0.025</td></tr><tr><td>A(4) = A(3)×<br>0.980 + 0.061</td><td>B(4) = B(3)×<br>0.994-0.153</td><td>C(4) = C(3)×<br>1.04-2</td><td>Θ(4) = Θ(3)×<br>1.207 + 0.004</td></tr><tr><td>A(5) = A(4)×<br>1.003-0.798</td><td>B(5) = B(4)×<br>0.956 + 0.176</td><td>C(5) = C(4)×<br>1.01-2.5</td><td>Θ(5) = Θ(4)×<br>0.916 + 0.048</td></tr></tbody></table>

and, the mean squared error (MSE) of the error function is designated as $MSE(C_p^{(k)})$.

## 3.2. Volume of clods

Another characteristic more specific to DEMs of soil surface is the sum of elevations which represent the volume delimited by the surface. With a clod made surface, the volume represents well the size of clods. The total volume occupied by clods for a surface at state k normalised by the volume of clods at initial state is considered:

<figure><img src="image_7.png" /><figcaption>Fig. 5. Histograms of clod disappearing.</figcaption></figure>

$$V_n^{(k)} = \frac{\sum_i V_i^{(k)}}{\sum_i V_i^{(0)}} \qquad (11)$$

with $V_i^{(k)}$ the volume enclosed in clod i at state k of the surface. $V_n^{(k)}$ corresponds to the normalised volume variation $VV(k)$ introduced in Vannier et al. (2018a). It is a measure of the volume variation due to evolution of clod shape under rain. The MSE of the normalised volume is designated as $MSE(V_n)$.

Therefore, volume decrease since state 2 can be modelled by exponential function of cumulated precipitation P (Zobeck and Onstad (1987); Vannier et al. (2018a; 2018b), Vannier et al., 2014):

$$V_n^{(k)} = \alpha e^{-\beta P^{(k)}} \quad k > 1 \qquad (12)$$

## 4. Results and discussion

Several studies have established empirical equations to relate change of surface roughness over time as a function of other parameters, such as slope, cumulated precipitation, rainfall energy (Bauer et al., 2015; De Oro & Buschiazzo, 2011; Vinci et al., 2020; Zhang & Xie, 2019; Zheng et al., 2019). In De Oro and Buschiazzo (2011), the decrease of random roughness due to rain can be related to previous roughness value or to cumulated precipitation and rain energy. The approach is the same in present study, with the special feature of characterising the roughness by structural elements that are clods. This is suitable because SSR is related to clod-size distribution. Also, the model was validated on a separate data set that was not used to estimate its parameters.
## 4.1. Learning phase

Two simulations of surface evolution were performed to validate the model of clod evolution under rain learned with tray 1. So, for these simulations, the initial state is the surface.

$HES_{l0}^{(1)}(x,y)$ composed of half-ellipsoids estimated from DEM of tray 1 at initial state, set at their place on a plane, and further states are simulated numerically to reproduce clod evolution under rain, according to the laws of evolution learned with tray 1. The difference between simulations lies in the random aspect of disappearance or appearance of clods. Table 2 gathers the metrics used to characterise a soil surface, that is to say, the number of clods, the volume occupied by clods and the maximum of autocorrelation function, which is equal to the variance of elevations. The last column is the measure of difference between autocorrelation functions of the reference surface ($HES_{l0}^{(k)}(x, y)$) and the simulation ($HES_{l1}^{(k)}(x,y)$ or $HES_{l2}^{(k)}(x, y)$). For both simulations, the number of clods is comparable to the reference at each state, and the volume of clods and the maximum of autocorrelation function are very comparable. This shows the goodness of the learning phase, as expressed also by the very low MSEs, which are $5.2 \times 10^{-5}$ in average and less than $7.2 \times 10^{-5}$. Therefore, the simulations of evolution reflect well the actual surface evolution as measured by DEMs, clod segmentation, and clod modelling by half ellipsoids. This evolution shows that the clods globally increase in volume during the first rain and then decrease during the following rains. Also, some clods disappear all along the rains.

## 4.2. Test phase

The same process of two simulations was performed with tray 2 for testing the robustness of the model of evolution when applied to a new surface. In that case, initial state is the surface.

$HES_{l0}^{(1)}(x,y)$ composed of half-ellipsoids estimated from DEM of tray 2 at initial state, set at their place on a plane, and further states are simulated numerically to reproduce clod evolution under rain, according to the laws of evolution learned with tray 1. Table 3 gathers the same metrics as Table 2 in order to follow the actual and the simulated evolution of tray 2. The number of clods are comparable at each states. Concerning the volume occupied by clods and the maximum of autocorrelation function, the values are comparable but are globally higher in the simulations than in the reference. This can be due to the little difference in clod size and layout between tray1 and tray 2. Clods globally decrease in volume since the first rain for actual evolution of tray 2, because they are smaller. The simulation reproduces the volume increase learned with tray 1. However, regarding the MSEs, the average is $9.8 \times 10^{-5}$ and the maximum $1.4 \times 10^{-4}$, that is to say that the mean

<table><caption>**Table 2**
Metrics of surface to compare two simulations with reference: learning case.</caption><thead><tr><th>State</th><th>Surface</th><th>N<sub>clods</sub></th><th>V<sub>n</sub><sup>(k)</sup></th><th>C<sub>p</sub><sup>(k)</sup>(0,0)</th><th>MSE(C<sub>p</sub><sup>(k)</sup>)</th></tr></thead><tbody><tr><td>n°1</td><td>HES<sub>l0</sub><sup>(1)</sup>(x,y)</td><td>172</td><td>1</td><td>85.1</td><td>–</td></tr><tr><td>n°2</td><td>HES<sub>l0</sub><sup>(2)</sup>(x,y)</td><td>111</td><td>1.119</td><td>85.4</td><td>–</td></tr><tr><td>n°3</td><td>HES<sub>l0</sub><sup>(3)</sup>(x,y)</td><td>106</td><td>0.975</td><td>71.8</td><td>–</td></tr><tr><td>n°4</td><td>HES<sub>l0</sub><sup>(4)</sup>(x,y)</td><td>93</td><td>0.878</td><td>65.8</td><td>–</td></tr><tr><td>n°5</td><td>HES<sub>l0</sub><sup>(5)</sup>(x,y)</td><td>81</td><td>0.715</td><td>54.0</td><td>–</td></tr><tr><td>n°2</td><td>HES<sub>l1</sub><sup>(2)</sup>(x,y)</td><td>126</td><td>1.112</td><td>85.9</td><td>3.9 &times; 10<sup>-5</sup></td></tr><tr><td>n°3</td><td>HES<sub>l1</sub><sup>(3)</sup>(x,y)</td><td>119</td><td>0.976</td><td>72.2</td><td>3.7 &times; 10<sup>-5</sup></td></tr><tr><td>n°4</td><td>HES<sub>l1</sub><sup>(4)</sup>(x,y)</td><td>99</td><td>0.890</td><td>66.7</td><td>4.1 &times; 10<sup>-5</sup></td></tr><tr><td>n°5</td><td>HES<sub>l1</sub><sup>(5)</sup>(x,y)</td><td>85</td><td>0.724</td><td>54.2</td><td>5.6 &times; 10<sup>-5</sup></td></tr><tr><td>n°2</td><td>HES<sub>l2</sub><sup>(2)</sup>(x,y)</td><td>119</td><td>1.084</td><td>85.1</td><td>5.8 &times; 10<sup>-5</sup></td></tr><tr><td>n°3</td><td>HES<sub>l2</sub><sup>(3)</sup>(x,y)</td><td>109</td><td>0.949</td><td>71.3</td><td>5.5 &times; 10<sup>-5</sup></td></tr><tr><td>n°4</td><td>HES<sub>l2</sub><sup>(4)</sup>(x,y)</td><td>95</td><td>0.870</td><td>65.9</td><td>5.6 &times; 10<sup>-5</sup></td></tr><tr><td>n°5</td><td>HES<sub>l2</sub><sup>(5)</sup>(x,y)</td><td>82</td><td>0.715</td><td>53.8</td><td>7.2 &times; 10<sup>-5</sup></td></tr></tbody></table>

<table><caption>**Table 3**
Metrics of surface to compare two simulations with reference: test case.</caption><thead><tr><th>State</th><th>Surface</th><th>N<sub>clods</sub></th><th>V<sub>n</sub><sup>(k)</sup></th><th>C<sub>p</sub><sup>(k)</sup>(0,0)</th><th>MSE(C<sub>p</sub><sup>(k)</sup>)</th></tr></thead><tbody><tr><td>n°1</td><td>HES<sub>l0</sub><sup>(1)</sup>(x,y)</td><td>241</td><td>1</td><td>70.8</td><td>–</td></tr><tr><td>n°2</td><td>HES<sub>l0</sub><sup>(2)</sup>(x,y)</td><td>150</td><td>0.981</td><td>64.8</td><td>–</td></tr><tr><td>n°3</td><td>HES<sub>l0</sub><sup>(3)</sup>(x,y)</td><td>131</td><td>0.880</td><td>56.2</td><td>–</td></tr><tr><td>n°4</td><td>HES<sub>l0</sub><sup>(4)</sup>(x,y)</td><td>114</td><td>0.769</td><td>49.0</td><td>–</td></tr><tr><td>n°5</td><td>HES<sub>l0</sub><sup>(5)</sup>(x,y)</td><td>84</td><td>0.634</td><td>39.5</td><td>–</td></tr><tr><td>n°2</td><td>HES<sub>l1</sub><sup>(2)</sup>(x,y)</td><td>155</td><td>1.073</td><td>70.6</td><td>9.0 &times; 10<sup>-5</sup></td></tr><tr><td>n°3</td><td>HES<sub>l1</sub><sup>(3)</sup>(x,y)</td><td>144</td><td>0.933</td><td>58.9</td><td>9.0 &times; 10<sup>-5</sup></td></tr><tr><td>n°4</td><td>HES<sub>l1</sub><sup>(4)</sup>(x,y)</td><td>115</td><td>0.850</td><td>54.0</td><td>9.8 &times; 10<sup>-5</sup></td></tr><tr><td>n°5</td><td>HES<sub>l1</sub><sup>(5)</sup>(x,y)</td><td>84</td><td>0.663</td><td>42.3</td><td>1.4 &times; 10<sup>-4</sup></td></tr><tr><td>n°2</td><td>HES<sub>l2</sub><sup>(2)</sup>(x,y)</td><td>156</td><td>1.070</td><td>71.2</td><td>8.0 &times; 10<sup>-5</sup></td></tr><tr><td>n°3</td><td>HES<sub>l2</sub><sup>(3)</sup>(x,y)</td><td>140</td><td>0.923</td><td>59.3</td><td>8.2 &times; 10<sup>-5</sup></td></tr><tr><td>n°4</td><td>HES<sub>l2</sub><sup>(4)</sup>(x,y)</td><td>105</td><td>0.837</td><td>54.1</td><td>9.0 &times; 10<sup>-5</sup></td></tr><tr><td>n°5</td><td>HES<sub>l2</sub><sup>(5)</sup>(x,y)</td><td>88</td><td>0.677</td><td>43.1</td><td>1.1 &times; 10<sup>-4</sup></td></tr></tbody></table>

difference between MSEs of learning and test phases is $4.6 \times 10^{-5}$ which is very good.

Table 4 shows the MSEs of volumes and the parameters $\alpha$ and $\beta$ of the volume decrease exponential model (Eq. (12)) for learning and test reference surfaces and respective simulations. As expected, the MSEs are a little better for simulation of learning surface evolution than for test surface evolution because the small individual errors on clods cumulate for the volume occupied by clods. Yet, the difference of average MSEs between learning and test cases is $3.9 \times 10^{-3}$ which remains good. The values of $\alpha$ are similar for the simulations despite the original difference of 1.389 and 1.228 for the reference surfaces, because the volume occupied by clods depends mainly on the greater clods that do not disappear. This is a signature of the evolution model, which reproduces the volume increase of tray 1 at state 2. The parameter $\beta$ in the exponential model shows the decrease of volume as a function of cumulated precipitation. We can notice that it has the same value of $3.5 \times 10^{-3}$ for the reference surfaces $HES_{l0}^{(k=2.5)}(x,y)$ (tray1) and $HES_{l0}^{(k=2.5)}(x,y)$ (tray2). Then, the simulations make the clod volume decrease a little faster in test case (initial state of tray 2) than in learning case (initial state of tray 1). This can be explained, for learning case, by the fact that parameter c (height of the half ellipsoid) was tuned not too decrease too fast so that the number of clods is close to target. In test case, this phenomenon is compensated by the fact that smaller clods decrease faster. For all cases, and both parameters $\alpha$ and $\beta$, the relative difference between simulation and reference is less than about 10%, which is satisfactory. Moreover, let us notice that the order of magnitude of parameters $\alpha$ and $\beta$ are comparable to that published by Zhang & Xie, 2019 ($\alpha \in [0.97; 1]$ and $\beta \in [0.1; 0.8] \times 10^{-3}$). In the cited paper, the soil surfaces have larger standard deviation of elevations ranging from 9.9 to 25.3 mm (against 9.2 and 8.4 mm for half-ellipsoid surfaces in present study), and the rainfall intensities are also lower than in present study (0.64 1.29 and 1.75 mm h$^{-1}$ against 33 mm h$^{-1}$ and 42 mm h$^{-1}$). Therefore, for comparable amounts of cumulated precipitation, the erosion is slower in their study, traduced by somewhat smaller values of $\alpha$ and $\beta$.

<table><caption>**Table 4**
MSE and model parameters of normalised volume for reference and two simulations in learning and test cases.</caption><thead><tr><th>Evolved Surfaces</th><th>MSE(V<sub>n</sub>)</th><th>α</th><th>β</th></tr></thead><tbody><tr><td>HES<sub>l0</sub><sup>(k=2.5)</sup>(x,y)</td><td>–</td><td>1.389</td><td>3.5 &times; 10<sup>-3</sup></td></tr><tr><td>HES<sub>l1</sub><sup>(k=2.5)</sup>(x,y)</td><td>6.7 &times; 10<sup>-5</sup></td><td>1.366</td><td>3.3 &times; 10<sup>-3</sup></td></tr><tr><td>HES<sub>l2</sub><sup>(k=2.5)</sup>(x,y)</td><td>5.1 &times; 10<sup>-4</sup></td><td>1.319</td><td>3.2 &times; 10<sup>-3</sup></td></tr><tr><td>HES<sub>l0</sub><sup>(k=2.5)</sup>(x,y)</td><td>–</td><td>1.228</td><td>3.5 &times; 10<sup>-3</sup></td></tr><tr><td>HES<sub>l1</sub><sup>(k=2.5)</sup>(x,y)</td><td>4.6 &times; 10<sup>-3</sup></td><td>1.356</td><td>3.7 &times; 10<sup>-3</sup></td></tr><tr><td>HES<sub>l2</sub><sup>(k=2.5)</sup>(x,y)</td><td>4.0 &times; 10<sup>-3</sup></td><td>1.328</td><td>3.6 &times; 10<sup>-3</sup></td></tr></tbody></table>
than in this paper. In Zobeck and Onstad (1987), in average $\beta = 2.6$
$10^{-3}$, and this parameter may also be higher than this average.

This shows, that the numerical model of clod evolution under rain is
suitable to reproduce clod evolution of surfaces close to that of the
learning case, including with a small variation. If more precision were
required with a slightly different surface, it would be possible to make
relearning to enhance the range of surfaces.

4.3. Representing soil surface by clod layout

Agricultural soil surfaces show several levels of roughness, due to
superposition of random small scale structure (clods and depressions),
and almost quasi-deterministic components of greater scale corre-
sponding to furrows and landscape. On a small sample of a square metre
range, the higher level of roughness is negligible. For convenience, the
deterministic component was not taken into account in this study. The
laboratory soil made surfaces were composed of clods set upon an
approximate plane. Future work will have to deal with the clod support.
In remote sensing studies, numerical simulations of soil surfaces can be
made with a filter reproducing the autocorrelation function model. In
Dusséaux and Vannier (2022), a five parameter model of bi-dimensional
autocorrelation function was introduced.

In order to make comparisons, the metric of autocorrelation function
was retained. Fig. 6 shows autocorrelation functions along (Oy) axis for
the laboratory surface, the surface obtained by settings segmented clods
onto a horizontal plane, the half-ellipsoid surface, and the 5-parameter
model of the autocorrelation estimated on the laboratory surface, at
initial state and final state. We can see that the clod surface and even
more the half-ellipsoid surface overestimate the maximum and departs
from the main lobe of autocorrelation function of laboratory surface
(from 0 to 50 mm). However, these surfaces follow the laboratory sur-
face for the rest of autocorrelation function (oscillations from 50 to 500
mm). The 5-parameter model of autocorrelation function is designed to
fit well the main lobe and neglect further oscillations. The MSEs between
autocorrelation functions of clod surface, half ellipsoid surface, 5-
parameter model and that of the laboratory surface (reference), are
respectively of 1.05 10⁻⁴, 1.52 10⁻⁴, 3.31 10⁻⁴ at initial state, and 1.80
10⁻⁴, 2.02 10⁻⁴, 5.29 10⁻⁴ at final state. This means that in our

simplified case of laboratory made surface, representing the surface by
clod layout is relevant and modelling clods by half-ellipsoids is a good
model.

The purpose of this study was centred on clod evolution model, but it
is also possible to generate an initial state of a soil surface, by adapting
the algorithm presented in Vannier et al., 2014. To generate a cloddy
surface, by setting half ellipsoids onto a plane, the statistical laws of the
six parameters have first to be learned. For example, with tray 1, co-
ordinates of the centres of gravity as well as orientation angles θ follow
uniform laws (p-values >70 %). Semi-axe a follows a Weibull distribu-
tion with parameters 7.25 and 1.21 plus a minimum value of 4.15 mm
(p-value = 51 %). Semi-axes b and c are correlated to semi-axe a (R² >
85 %) and can be derived according to linear regressions with noise.
Estimation of these parameters from surfaces recorded in the field
should be carried on to refine these values. This would have first to
tackle the challenging problem of clod segmentation. Then, the principle
of generation is to set the number of clods and the size of the surface,
next to draw the parameters of each clod and last, to place the
half-ellipsoids on the plane without overlapping. This point is specific to
the experimental design of this study. In Kamphorst et al. (2005), a
Boolean model was proposed to superimpose half-spheres. So far, rep-
resenting a soil surface by clod layout is a good approximation. With the
properties of clods at initial state of surface, and with the model of
evolution, the whole process of numerical surface generation is raised. It
will have to be enhanced by modelling the support of clods and studying
the effect of clod spacing. Moreover, the model parameters were set from
an experiment in the laboratory. As the soil surfaces were composed of
sieved clods extracted from the field, the evolution of clods is indicative.
A soil surface in the field would show also depressions that will have to
be studied.

5. Conclusion

This paper proposes a way to generate numerically soil surfaces close
to seedbeds with roughness evolved with rain. As soil surface roughness
is related to clod-size distribution, the principle is to set modelled clods
on a planar surface, as a first approximation of the surface. DEM
recording and Machine Learning allowed to learn clod evolution under

<figure><img src="image_8.png" /><figcaption>Fig. 6. Autocorrelation functions of tray1 with respect to (Oy): a) initial state, and b) final state.</figcaption></figure>
rain and to simulate a seedbed-like surface evolution with very good accuracy. Numerical simulation of surfaces is an alternative to surface measuring that enables to get quickly a large number of surfaces of desired properties. The small scale is useful for various applications, such as rough surface scattering, and more generally, modelling where soil surface or other like surface is an input. This simple surface representation will have to be enhanced by adding depressions, studying the objects spacing, and adding a proper support.

An interesting point is that setting modelled clods on a flat surface is representative of the autocorrelation function of the surface. That is to say that it captures the main properties of the surface. Also, as soil surface roughness is shaped by tillage operations and evolves due to weathering, it was important to learn how to generate an evolved roughness. It may benefit also for modelling signals of natural origin that have an exponential like autocorrelation function.

## CRediT authorship contribution statement

E. Vannier: Writing – original draft, Visualization, Validation, Supervision, Software, Resources, Methodology, Investigation, Funding acquisition, Formal analysis, Conceptualization. R. Dusséaux: Writing – review & editing, Validation, Resources, Project administration, Methodology, Funding acquisition, Formal analysis, Conceptualization.

## Funding

This work has been supported by the Programme National de Télédétection Spatiale (PNTS, http://programmes.insu.cnsr.fr/pnts/), grant no PNTS-2020.20).

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## References

Ajdadi, F. R., Abbaspour-Gilandeh, Y., Mollazade, K., & Hasanzadeh, R. P. R. (2016). Application of machine vision for classification of soil aggregate size. *Soil and Tillage Research, **162***, 8–17.

Arvidsson, J., & Bölenius, E. (2006). Effects of soil water content during primary tillage – laser measurements of soil surface changes. *Soil and Tillage Research, **90***, 222–229.

Barbosa, L. A. P., & Gerke, H. H. (2022). Structural heterogeneity of soil clods: Correlating Weibull parameters to fracture surface topography. *Geoderma, **428***, Article 116161.

Bauer, T., Strauss, P., Grims, M., Kamptner, E., Mansberger, R., & Spiegel, H. (2015). Longterm agricultural management effects on surface roughness and consolidation of soils. *Soil and Tillage Research, **151***, 28–38.

Bogrekci, I., & Godwin, R. J. (2007). Development of an image-processing technique for soil tilth sensing. *Biosystems Engineering, **97***, 323–331.

Borselli, L., & Torri, D. (2010). Soil roughness, slope and surface storage relationship for impervious areas. *Journal of Hydrology, **393***, 389–400.

Bourlier, C., Pinel, N., & Kubické, G. (2013). *Method of moments for 2D scattering problems: Basic concepts and applications*. Hoboken, NJ, USA: Wiley.

Bykhovsky, D., & Lyandres, V. (2018). On the numerical generation of positive-axis-defined distributions with an exponential autocorrelation function. *Digital Signal Processing, **77***, 43–47.

Chilès, J. P., & Delfiner, P. (1999). *Geostatistics (modeling spatial uncertainty)*. New York: John Wiley & Sons INC.

Chimi, O. (2012). *Probabilistic characterization and synthesis of agricultural surfaces by structuring objects from high resolution images. (Caractérisation probabiliste et synthèse de surfaces agricoles par objets structurants à partir d'images haute résolution.)*. PhD thesis (pp. 102–109). France: Université de Versailles Saint-Quentin-en-Yvelines.

Choi, W., Shi, F., Lowe, M. J. S., Skelton, E. A., Craster, R. V., & Daniels, W. L. (2018). Rough surface reconstruction of real surfaces for numerical simulations of ultrasonic wave scattering. *NDT & E International, **98***, 27–36.

Cierniewski, J., Verbrugghe, M., & Marlewski, A. (2002). Effects of farming works on soil surface bidirectional reflectance measurements and modelling. *International Journal of Remote Sensing, **23(6)**, 1075–1094.

Darboux, F., & Huang, C. (2003). An instantaneous-profile laser scanner to measure soil surface micro topography. *Science Society of America Journal, **67***, 92–99.

De Oro, L. A., & Buschiazzo, D. E. (2011). Degradation of the soil surface roughness by rainfall in two loess soils. *Geoderma, **164***, 46–53.

Doane, T. H., Gearon, J. H., Martin, H. K., Yanites, B. J., & Edmonds, D. A. (2024). Topographic roughness as an emergent property of geomorphic processes and events. *AGU Advances, **5***, https://doi.org/10.1029/2024AV001264

Dusséaux, R., & Vannier, E. (2022). Soil surface roughness modelling with the bidirectional autocorrelation function. *Biosystems Engineering, **220***, 87–102.

Dusséaux, R., Vannier, E., Taconet, O., & Granet, G. (2012). Study of backscatter signature for seedbed surface evolution under rainfall - influence of radar precision. *Progress in Electromagnetics Research, **125***, 415–437.

Feral, L., Mesnard, F., Sauvageot, H., Castanet, L., & Lemorton, J. (2000). Rain cells shape and orientation distribution in south-west of France. *Physics and Chemistry of the Earth - Part B: Hydrology, Oceans and Atmosphere, **25(10)**, 1073–1078.

Foroogh, P., Stroh, A., Magagnato, F., Jakirlić, S., & Frohnapfel, B. (2017). Towards a universal roughness correlation. *Journal of Fluids Engineering, **139(12)**, Article 121201.

Foster, G. R., Eppert, F. P., & Meyer, L. D. (1979). A programmable rainfall simulator for field plots. In *Proceedings of the rainfall simulator workshop. Tucson, Arizona. March 7-9, 1979, pages 45-59, Sidney, Montana. Agricultural reviews and manuals, ARMW-10*. Oakland, CA: United States Department of Agriculture - Science and Education Administration.

Franceschetti, G., & Riccio, D. (2007). *Scattering, natural surfaces and fractals*. Elsevier. https://doi.org/10.1016/B978-0-12-265655-2.X5000-4

Garbout, A., Munkholm, L. J., & Hansen, S. B. (2013). Temporal dynamics for soil aggregates determined using X-ray CT scanning. *Geoderma, **204-205***, 15–22.

Géron, A. (2022). *Hands-on machine learning with scikit-learn, keras, and TensorFlow* (3rd ed.). O'Reilly Media, Inc. ISBN: 9781098125974.

Gilliot, J. M., Vaudour, E., & Michelin, J. (2017). Soil surface roughness measurement: A new fully automatic photogrammetric approach applied to agricultural bare fields. *Computers and Electronics in Agriculture, **134***, 63–78.

Grohmann, C. H., Smith, M. J., & Riccomini, C. (2011). Multiscale analysis of topographic surface roughness in the Midland Valley, Scotland. *IEEE Transactions on Geoscience and Remote Sensing, **49***, 1200–1213.

Ishii, S., Suto, I., Tabe, H., Aono, K., Zhao, M., Ueshige, Y., Matsushita, K., Iritani, T., Hanamoto, T., Nakao, M., & Nagato, K. (2023). Construction of a soil clod recognition bench-scale experiment for discrete element method modelling of tilling phenomena. *Journal of Terramechanics, **109***, 63–71.

James, G., Witten, D., Hastie, T., & Tibshirani, R. (2021). An introduction to statistical learning (with applications in R). *Springer texts in statistics* (Second.). New York: Springer.

Jensen, T., Green, O., Munkholm, L. J., & Karstoft, H. (2016). Fourier and granulometry methods on 3D images of soil surfaces for evaluating soil aggregate size distribution. *Applied Engineering in Agriculture, **32***, 609–615.

Kamphorst, E. C., Chadkeuf, J., Jetten, V., & Guérif, J. (2005). Generating 3D soil surfaces from 2D height measurements to determine depression storage. *Catena, **62***, 189–205.

Koltermann, C. E., & Gorelick, S. (1996). Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating, and descriptive approaches. *Water Resources Research, **32(9)**, 2617–2658.

Lievens, H., Vernieuwe, H., Alvarez-Mozos, J., De Baets, B., & Verhoest, N. E. C. (2009). Error in radar-derived soil moisture due to roughness parameterization: An analysis based on synthetic surface profiles. *Sensors, **9***, 1067–1093.

Marchand, R. T., & Brown, G. S. (1999). On the use of finite surfaces in the numerical prediction of rough surface scattering. *IEEE Transactions on Antennas and Propagation, **47(4)**, 600–604.

Martinez-Agirre, A., Alvarez-Mozos, J., & Gimenez, R. (2016). Evaluation of surface roughness parameters in agricultural soils with different tillage conditions using a laser profile meter. *Soil and Tillage Research, **19***, 30–30.

Ogilvy, J. A. (1991). *Theory of wave scattering from random rough surfaces*. Bristol: Adam Hilger. https://doi.org/10.1088/0034-4885/50/12/001

Pan, C., Dusséaux, R., & Emad, N. (2016). A complex mix-shifted parallel QR algorithm for the C-method. *Progress in Electromagnetics Research, **B68***, 159–171.

Sandri, R., Anken, T., Hilfiker, T., Sartori, L., & Bollhalder, H. (1998). Comparison of methods for determining cloddiness in seedbed preparation. *Soil and Tillage Research, **45***, 75–90.

Schapel, A., Marschner, P., & Churchman, J. (2019). Influence of clay clod size and number for organic carbon distribution in sandy soil with clay addition. *Geoderma, **335***, 123–132.

Smith, M. W. (2014). Roughness in the Earth sciences. *Earth-Science Reviews, **136***, 202–225.

Takken, I. (2000). *Effects of roughness on overland flow and erosion* (p. 218). Belgium: Katholieke Universiteit Leuven. Thesis.

Thomsen, L. M., Baartman, J. E., Barneveld, R. J., Starkloff, T., & Stolte, J. (2015). Soil surface roughness: Comparing old and new measuring methods and application in a soil erosion model. *Soils, **1***, 399–410.

Trevisani, S., & Rocca, M. (2015). Mad: Robust image texture analysis for applications in high resolution geomorphology. *Computers & Geosciences, **81***, 78–92.

Trevisani, S., Teza, G., & Guth, P. (2023). A simplified geostatistical approach for characterizing key aspects of short-range roughness. *Catena, **223***, Article 106927.

Vannier, E., Dusséaux, R., Taconet, O., & Darboux, F. (2019). Wavelet-based clod segmentation on digital elevation models of a soil surface with or without furrows. *Geoderma, **356***, Article 113933.

Vannier, E., Taconet, O., Dusséaux, R., & Chimi-Chiadjeu, O. (2014). Statistical characterisation of bare soil surface micro relief. In *M. Marghany (Ed.), Advance geoscience remote sensing* (pp. 207–228). Rijeka, Croatia: Edition Intech.

Vannier, E., Taconet, O., Dusséaux, R., & Darboux, F. (2018a). A study of clod evolution in simulated rain on the basis of digital elevation models. *Catena, **160***, 212–221.
Vannier, E., Taconet, O., Dusséaux, R., & Darboux, F. (2018b). Using digital elevation models and image processing to follow clod evolution under rainfall. *Journal of toxicology and ecology, 2*, 113.

Verhoest, N. E. C., Lievens, H., Wagner, W., Alvarez-Mozos, J., Susan Moran, M., & Mattia, F. (2008). On the soil roughness parameterization problem in soil moisture retrieval of bare surfaces from synthetic aperture radar. *Sensors, 8*, 4213-4248.

Vinci, A., Todisco, F., Vergni, L., & Torri, D. (2020). A comparative evaluation of random roughness indices by rainfall simulator and photogrammetry. *Catena, 188*, Article 104468.

Wang, W., Kravchenko, A. N., Smucker, A. J. M., & Rivers, M. L. (2011). Comparison of image segmentation methods in simulated 2D and 3D micro tomographic images of soil aggregates. *Geoderma, 162*, 231-241.

Wang, M., Liu, H., Rezanezhad, F., Zak, D., & Lennartz, B. (2023). The influence of micro topography on soil carbon accumulation and nutrient release from a rewetted coastal peatland. *Geoderma, 438*, Article 116637.

Zhang, G.-H., & Xie, Z.-F. (2019). Soil surface roughness decay under different topographic conditions. *Soil and Tillage Research, 187*, 92-101.

Zhao, Z., Wei, H., Liu, S., & Xue, Z. (2024). Estimation of agricultural soil surface roughness based on ultrasonic echo signal characteristics. *Soil & Tillage Research, 239*, Article 106038.

Zheng, X., Li, L., Chen, S., Jiang, T., Li, X., & Zhao, K. (2019). Temporal evolution characteristics and prediction methods of spatial correlation function shape of rough soil surfaces. *Soil and Tillage Research, 195*, Article 104417.

Zobeck, T. M., & Onstad, C. A. (1987). Tillage and rainfall effects on random roughness: A review. *Soil and Tillage Research, 9*, 1-20.

Zribi, M., Le Morvan, A., Dechambre, M., & Baghdadi, N. (2010). Numerical backscattering analysis for rough surfaces including a cloddy structure. *IEEE Transactions on Geoscience and Remote Sensing, 48*(5), 2367-2374.