*Article* **Simulating Soil–Disc Plough Interaction Using Discrete Element Method–Multi-Body Dynamic Coupling**

**Mustafa Ucgul 1,2**


**Abstract:** Due to their (a) lower draught force requirements and (b) ability to work at deeper operation depths and faster operation speeds, disc ploughs have gained interest in Australia. A modified version of the disc plough that involves removing every second disc and fitting larger and often more concave discs has become popular. However, the development of the one-way modified disc plough is in its infancy, and a detailed analysis is required, particularly on soil movement. Historically, the soil movement analysis of the soil–tool interactions is conducted using empirical methods. However, the experimental tests are resource and labour intensive. When the soil and tool interaction can be accurately modelled, more efficient tools can be designed without performing expensive field tests, which may only be undertaken at certain times of the year. This study modelled the interaction between soil and a one-way modified disc plough using the discrete element method (DEM). As the disc plough is a passive-driven tool, the rotational speed of the disc plough was modelled using DEM-MBD (multi-body dynamic) coupling. The results of the study show that DEM-MBD coupling can predict the rotational speed of the disc plough with a maximum relative error of 6.9%, and a good correlation was obtained between the DEM-predicted and actual soil movement (R2 = 0.68).

**Keywords:** discrete element method (DEM); multi-body dynamics (MBD); DEM-MBD coupling; topsoil burial

#### **1. Introduction**

A disc plough is a primary tillage tool used when a mouldboard plough is unsuitable, such as in stony, stiff, and dry soils [1]. Although a mouldboard plough provides complete soil inversion, disc ploughs are cheaper and can be used in broader soil conditions. Additionally, the draught force requirement of the disc ploughs is lower than that of the mouldboard plough, which also reduces fuel consumption [2]. In addition, disc ploughs can be used at faster operating speeds. In Australia, disc ploughs are used to bury (a) non-wetting topsoil, (b) weed seeds, or (c) surface-applied amendments (organic matter, nutrients, and lime) from the surface to 250–350 mm depth. A modified version of the disc plough that involves removing every second disc and fitting larger and often more concave discs has also become popular. Using larger discs creates a larger break-out pressure and allows the disc plough to be used deeper, whereas removing every second disc provides more space for soil to turn over and a greater degree of inversion [3]. However, to improve the effectiveness of the burial performance of the modified disc ploughing, understanding the soil movement process is essential. Soil movement due to tillage operations are generally investigated using empirical methods. Empirical methods create relationships from specific test results and therefore have no extrapolative capabilities outside this context [4,5]. Although the finite element method (FEM) and computational fluid dynamics (CFD) can be used to model soil flow around a tool, these methods are limited, as they cannot model mixing between soil layers or crack propagation [6]. Modelling the soil–tool interaction is

**Citation:** Ucgul, M. Simulating Soil–Disc Plough Interaction Using Discrete Element Method–Multi-Body Dynamic Coupling. *Agriculture* **2023**, *13*, 305. https://doi.org/10.3390/ agriculture13020305

Academic Editor: Francisco J. Castillo Ruiz

Received: 15 December 2022 Revised: 13 January 2023 Accepted: 23 January 2023 Published: 27 January 2023

**Copyright:** © 2023 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

a highly complex process due to the non-linear behaviour of the soil. The discrete element method (DEM) can be used to simulate the soil–disc plough interaction. Recent studies have shown that DEM can simulate the soil–implement interaction (considering tillage forces and soil movement) [7–12].

To model the performance of the soil–disc plough interaction, the determination of the disc's rotational speed is essential. A disc plough is a passive-driven tool whose rotation is driven by the soil reaction forces. Therefore, using DEM simulation is not sufficient to simulate soil–disc plough interaction. However, DEM simulation can be coupled with multi-body dynamic (MBD) software to predict the rotational speed of the disc plough. Nevertheless, using DEM-MBD coupling in real-scale simulations is computationally expensive due to the longer simulation time.

This study aims to model a disc plough's topsoil burial performance using DEM-MBD coupling, which was not investigated in prior studies. The DEM parameters of the field soil were determined using the angle of repose test. Then, the rotational speed of the disc plough was predicted using DEM-MBD coupling. After that, the topsoil burial performance of the disc plough was simulated using DEM and then validated with field tests. The colour difference between the deep- and topsoil was used along with a digital image processing algorithm to validate the DEM simulation results. The proposed simulation approach will help designers and researchers to evaluate the performance of any passive-driven soil-engaging tool without performing cost- and time-intensive field tests.

#### **2. Materials and Methods**

The Material and Methods section is two-fold: simulating the soil–disc plough interaction using DEM and DEM-MBD coupling and performing field tests to validate the results.

#### *2.1. DEM and DEM-MBD Coupling Simulations*

This study conducted the DEM simulations using commercial DEM software, EDEM 2019TM (developed by DEM Solutions Ltd. In Edinburgh (United Kingdom). A DELL Precision T7610 Dual Xeon E5-2680 v3 @ 2.5 GHz and 128 GB RAM—workstation computer was used to run the EDEM 2019TM software. A hysteretic spring contact model (HSCM) integrated with a linear cohesion model suggested by [13] was used as the mathematical model for contact calculations.

Compressible materials such as soil can be modelled using the HSCM. In the HSCM, the particles behave as a linear elastic material until their stress level reaches a predefined stress (yield) (Figure 1). Once the total stress on the contact area exceeds the yield stress, the particles are deformed plastically [14]. On the other hand, the linear cohesion model adds a cohesive force through normal contact direction (for modelling soil cohesion).

The HSCM uses the following equations to calculate the total normal (Fn) and tangential (Ft) forces [14].

$$\mathbf{F\_n = F\_n^s + F\_n^d} \tag{1}$$

$$\mathbf{F\_n = F\_n^s + F\_n^d} \tag{2}$$

$$\mathbf{F\_n^s} = -\begin{cases} \mathbf{K\_1} \cdot \mathbf{U\_{abn}} & \text{loading} \\ \mathbf{K\_2} \cdot (\mathbf{U\_{abn}} - \mathbf{U\_0}) & \text{unloading/reloading} \\ 0 & \text{unloading} \end{cases} \tag{3}$$

where Fs n, Fd n, F<sup>s</sup> <sup>t</sup> , and F<sup>d</sup> <sup>t</sup> are the normal contact force, normal damping force, tangential contact force, and tangential damping force, respectively. Uabn is the normal component of the relative displacement, and U0 is the residual overlap. K1 and K2 are the loading and unloading stiffnesses, respectively, and [15,16] suggest the following equation for calculating K1 and K2.

$$\mathbf{K}\_1 = \mathbf{5} \times \mathbf{r}\_{\text{eq}} \times \min(\mathbf{Y}\_{\text{a}}, \mathbf{Y}\_{\text{b}}) \tag{4}$$

$$\mathbf{K}\_2 = \frac{\mathbf{K}\_1}{\mathbf{e}^2} \tag{5}$$

where Y is the yield strength and e is the coefficient of restitution; req is the equivalent radius, and [14] calculates req as

$$\frac{1}{\mathbf{r\_{eq}}} = \frac{1}{\mathbf{r\_a}} + \frac{1}{\mathbf{r\_b}} \tag{6}$$

For each time step, the U0 is updated as

$$\mathbf{U}\_0 = \begin{cases} \mathbf{U}\_{\text{absn}} \cdot \left( 1 - \frac{\mathbf{K}\_1}{\mathbf{K}\_2} \right) & \text{loading} \\ \mathbf{U}\_0 & \text{unloading/reloading} \\ \mathbf{U}\_{\text{absn}} & \text{unloading} \end{cases} \tag{7}$$

HSCM calculates the F<sup>s</sup> <sup>t</sup> , Fd n, and F<sup>d</sup> <sup>t</sup> using the following formulas [14]:

$$\mathbf{F}\_{\mathbf{t}}^{\mathbf{s}} = \mathbf{n}\_{\mathbf{k}} \cdot \mathbf{K}\_{\mathbf{l}} \cdot \mathbf{U}\_{\mathbf{a}\mathbf{b}\mathbf{t}} \tag{8}$$

$$\mathbf{F\_n^d} = -\mathbf{n\_c} \times \sqrt{\frac{4 \times \mathbf{m\_{eq}} \times \mathbf{K\_1}}{1 + \left(\frac{\pi}{\text{Im } \mathbf{e}}\right)^2} \cdot \dot{\mathbf{U}}\_{\text{abn}}} \tag{9}$$

$$\mathbf{F\_t^d} = -\sqrt{\frac{4 \times \mathbf{m\_{eq}} \times \mathbf{n\_k} \times \mathbf{K\_1}}{1 + \left(\frac{\pi}{\text{Im } \mathbf{e}}\right)^2} \cdot \dot{\mathbf{U\_{abt}}}} \tag{10}$$

where Uabt is the tangential component of the relative displacement, nk is the stiffness factor, defined as the ratio of tangential stiffness to normal loading stiffness, U˚ abn and U˚ abt are the normal and tangential components of the relative velocity, respectively, nc is the damping factor that controls the amount of velocity-dependent damping, and meq is the equivalent mass and calculated as suggested by [14]:

$$\frac{1}{\mathbf{m}\_{\rm eq}} = \frac{1}{\mathbf{m}\_{\rm a}} + \frac{1}{\mathbf{m}\_{\rm b}} \tag{11}$$

The following should also be considered when calculating Ft:

$$\mathbf{F\_{t}} = -\min\left(\mathbf{n\_{k}} \cdot \mathbf{K\_{l}} \cdot \mathbf{U\_{abt}} + \mathbf{F\_{t}^{d}}, \boldsymbol{\mu} \cdot \mathbf{F\_{n}^{s}}\right) \tag{12}$$

The magnitude of the moments caused by total tangential force (M) and the rolling resistance (Mr) were computed as per [17]:

$$\mathbf{M} = \mathbf{r}\_{\text{con}} \times \mathbf{F}\_{\text{t}} \tag{13}$$

$$\mathbf{M\_{f}} = -\mathbf{r\_{con}} \times \boldsymbol{\mu\_{r}} \cdot \mathbf{F\_{n}^{\boldsymbol{\theta}}} \cdot \boldsymbol{\lambda\_{\boldsymbol{\theta}}} \tag{14}$$

where rcon is the perpendicular distance of the contact point from the centre of mass, μ<sup>r</sup> is the coefficient of rolling friction, and λ<sup>θ</sup> is the unit vector of angular velocity at the contact point. The position of the particle is calculated by integrating Equations (15) and (16).

$$
\ddot{\mathbf{U}} = \frac{(\mathbf{F\_n} + \mathbf{F\_t})}{\mathbf{m}} \tag{15}
$$

$$\mathbf{R} = \frac{(\mathbf{M} + \mathbf{M}\_{\mathbf{r}})}{\mathbf{I}} \tag{16}$$

Cohesion is calculated and added to Equation (1) as [14]:

$$\mathbf{F}\_{\mathbf{c}} = \boldsymbol{\xi} \times \mathbf{A}\_{\mathbf{c}} \tag{17}$$

where ξ is the cohesion energy density and Ac is the contact area. Hence, Equation (1) becomes

$$\mathbf{F\_n = F\_n^s + F\_n^d + F\_c} \tag{18}$$

Using actual particulate sizes and shapes in 3D DEM simulations increases computation time (due to the increased number of contacts) and computation costs. Therefore, larger particle sizes are used in 3D DEM simulations. In addition, spherical particles are generally preferred due to their computational simplicities. However, when larger particle sizes are used, the parameters of these larger particles should be determined via a calibration process, so that these larger particles act like the simulated granular material.

In this study, randomly generated 7.5 mm nominal radii spherical particles (between nominal radii x (0.5 to 2)) were used. This particle size distribution was based on the actual particle size distribution (Table 1). A 0.3 mm sieve size was considered as the base value (7.5 mm nominal radii), particle sizes less than 0.3 mm sieve size were considered 0.5 × 7.5 mm nominal radii, and particle sizes greater than 0.3 mm sieve size were considered 2 × 7.5 mm nominal radii. In other words, very small and large particle sizes were grouped in ×0.5 and ×2 nominal radii.

**Table 1.** Particle size distribution of the test soil.


Note that considering the number of particles created in the DEM simulations and available computational software licenses, 7.5 mm nominal radii was the smallest particle size that could be selected.

The DEM parameters were determined from the literature or some physical tests (Table 2). The friction coefficient of soil–soil and soil–steel (using a piece of polished steel placed in the upper portion of the shear box that was sheared over the soil) were determined using direct shear box tests. The coefficient of rolling friction between soil and steel was measured by performing inclined plane tests using a steel sphere (9.5 mm radius). Firstly, a container was filled with test soil, and the steel ball was placed. After that, the tilt angle of the container gradually increased. The tilt angle was measured when the ball began to roll over the soil. The coefficient of rolling friction was then calculated from the measured tilt angle.

However, the coefficient of rolling friction and cohesive energy density of soil particles were calibrated using an angle of repose test. This calibration process is based on matching simulation results to actual measured results using a soil sample taken from the field [18]. The soil sample from the field was collected using a sealed bag (to keep the soil's moisture), and the angle of repose test was conducted on the same day in the soil laboratory. Note that the soil samples collected from the field were used for the abovementioned soil tests.

A soil sample was placed in a pipe (100 mm diameter and 300 mm long) to measure the angle of repose. After that, the pipe was lifted upward at a constant speed of 500 mm s−<sup>1</sup> using a Hounsfield Tensiometer. Soil flowed into a cylindrical tray (200 mm diameter and 22.5 mm long) until the soil overflowed and formed a pile. When at rest, images of the angle of repose were captured and processed as per [18]. The flow chart of the process for the measurement of the angle of repose is given in Figure 2. The method to measure the angle of repose was based on an automatic search of the best-fitted line formula along the boundary of the material pile using linear regression analysis. After the images of the soil pile were taken, digital image processing techniques were used to segment the conical shape of the pile from the original image. For this process, a global batch clustering algorithm developed by [19] was employed: "Morphological operations were then applied to segment the boundary pixel points of the pile. In the final step, the angle of repose was determined by searching the maximum correlation coefficients using linear regression, using a series of gradually increasing sizes of windows along the entire boundary of the material pile" [18].

**Figure 2.** Image processing scheme for automatic measurement of the angle of repose.


**Table 2.** EDEM parameters used for simulation.

After measuring the angle of repose, by varying the coefficient of rolling friction of soil–soil and cohesive energy density, the measured angle of repose of 31.48◦ (±0.5◦) was achieved by using a trial-and-error process in the DEM simulation (Figure 3).

**Figure 3.** (**a**) Angle of repose test and (**b**) angle of repose simulation.

Simulating the rotational speed of the disc plough cannot be carried out without using multi-body dynamic coupling, as the disc is rotated due to the force from soil particles (as the disc plough is a passive-driven tool). To model disc rotation, MSC ADAMS software was coupled with EDEM. However, using EDEM-ADAMS coupling for real-scale simulation is not viable due to the available computational power. Therefore, this study determined the rotational speed of the disc plough with one disc using a shorter soil bin (10,000 mm long × 3250 mm wide × 500 mm deep).

To simulate the actual field operation, an open furrow was produced at 300 mm depth (Figure 4). Into this open furrow, a furrow was ploughed using the disc plough at 300 mm depth.

**Figure 4.** EDEM-ADAMS co-simulation to calculate the rotational speeds of the disc.

The CAD model of the disc plough assembly was created using SolidWorksTM software as an assembly (which consists of a hub and a disc). The designed assembly was then imported to MSC ADAMSTM software. In MSC ADAMSTM software, a cylindrical joint was added between the hub's centre and the disc plough's centre. A translational connection was also added between the ground and the disc hub, which was also restrained from all other motions except the direction of travel. After that, a translational motion was assigned to the disc hub through the direction of travel. Finally, the EDEMTM and MSC ADAMS TM were coupled, and the simulations were carried out. The rotational speeds of the soil–disc plough interactions were predicted at 300 mm depth and 2.5, 5, 7.5, and 10 km h−<sup>1</sup> forward speeds for a single disc.

In order to simulate the real-scale simulation of the soil–disc plough interaction, a virtual field soil with dimensions of 20,000 mm long × 3250 mm wide × 500 mm deep was generated to predict topsoil burial. The DEM soil bulk density was set to 1481 kg/m3 (the bulk density used in the field test). To achieve the desired bulk density, particles (8,819,576 in total) were compressed using an upper physical plane until the desired bulk density was achieved. After that, the disc plough assembly was imported, the depth, speed, and rotational speed (predicted in the first set of simulations) were set, and simulations were carried out (Figure 5). The simulation was run at 300 mm operation depth and 5 and 7.5 km h−<sup>1</sup> forward speeds, the same operational depth and speed used in the field experiment.

**Figure 5.** DEM simulation of soil–disc plough interaction to predict topsoil burial.

In order to compare the DEM-predicted soil movement with the field experiment, the colour of the DEM particles in a 2000 mm wide × 1000 mm long × 40 mm deep trench was changed from brown to blue (Figure 5). After the simulations were carried out, the brown-coloured sand particles were removed from the simulation, and the coordinates of the blue-coloured particles were exported to Microsoft Excel. Then, the percentage of the volumes of blue-coloured particles at each 50 mm increment was calculated.

A 5-furrow disc plough was designed and used to reduce the simulation size (in opposition to the 13-furrow disc plough used in the field tests). As the first disc creates an open furrow, the soil burial of this furrow is more excessive than that of the other discs. Therefore, the topsoil burial results of the first furrow were not considered in the analysis.

#### *2.2. Field Tests*

In order to validate the DEM simulation results of the modified one-way disc plough– soil interaction, a field experiment was performed at Malinong (35.5106◦ S, 139.5157◦ E), South Australia, in May 2019. The soil was sandy loam with a bulk density, and moisture contents were 1481 kg m−<sup>3</sup> and 1.6% (dry basis), respectively. The moisture content of the soil was measured by the oven drying method (ASABE Standards, 2008), and the bulk density of the soil was measured using a core sample (50 mm diameter and 77 mm long). A thirteen-furrow modified disc plough was used in the experiment (Figure 6).

**Figure 6.** One-way disc plough used in the field tests.

Field experiments were undertaken at 300 mm depth and 5 and 7.5 km h−<sup>1</sup> forward speeds. Due to the long-term no-till practice, the top 40 mm of the topsoil had a darker colour than deep soil. This colour difference was used in the tests to validate the DEM simulation results of topsoil burial. After ploughing, three vertical excavations were conducted (using a shovel) at 750 mm increments across the direction of travel. The digital photos of each slice were taken. A metal frame was designed to ensure the camera was always the same distance from the cut when taking the pictures to assist in the accuracy of the later image processing. Note that the topsoil burial profile of a plough (mouldboard or disc) is continuous and uniform, multiple slices were investigated in the field test, and the quantification method is not subjective; therefore, no repeated tests were performed.

The images were then digitally analysed to determine the pixel locations of the darkcoloured topsoil. For this process, the photos were cropped to remove undesired features such as excessive background and shadows to help the software detect the target features. A clustering algorithm developed by [25] was then run over the image to highlight the location of the dark-coloured topsoil and produce a binary plot of the pixels.

In the algorithm, after the region of interest (ROI) is selected, the original RGB (Red, Green, Blue) components are transformed into the perceptual colour components (i.e., a\* and b\* from Lab colour space and Hue component from CIE standard). After that, the PCA solution from the image data set or subset is used in a hierarchical manner along with the termination measure under an ellipsoidal property. Then, the initial value of centroids of the number of colour clusters is found in the first part of the algorithm (Figure 7). In the second part, the initial parameters are refined further by a modified global clustering algorithm. Then, the refined centroids are used to label the different colour clusters (Figure 7).

**Figure 7.** Digital image processing (DIP) method for acquiring the featured coloured soil cluster image.

The pixel coordinates were then processed in MS Excel, where all the pixel coordinates from the three pictures were combined for each test condition, and the percentages of the topsoil burial at each 50 mm interval were calculated.

#### **3. Results**

Measured and DEM-MBD coupled predicted rotational speeds are given in Figure 8. As shown in Figure 8, increasing the forward speed increases the rotational speed. Increasing the forward speed from 3 km h−<sup>1</sup> to 9 km h−<sup>1</sup> increases the rotational speed from 13.07 rpm to 42.73 rpm. The results show that the measured and DEM-MBD-predicted rotational speed values are in good agreement, with a maximum relative error of 6.9%. This finding indicates that the proposed modelling approach can be effectively used to predict the rotational speed of the disc plough.

**Figure 8.** Measured (field test) and DEM-MBD predicted rotational speed vs. forward speed.

The predicted rotational speed results in Figure 8 were used in full-scale DEM simulations to evaluate the topsoil burial performance of the one-way modified disc plough.

Examples of images captured from the field experiments for each test are given in Figure 9. The combined cross-sectional locations of the buried topsoil observed in the excavated profiles (sum of three images) are represented in Figure 10. The results of the DEM simulations are also shown in Figure 10. The field test and DEM simulated results using a deep working modified one-way plough showed similar visual trends. At slower speeds, more distinct and vertical burial profiles were obtained, while burial profiles became horizontal at higher speeds [26].

300 mm depth and 7.5 km h<sup>−</sup>1 speed

**Figure 9.** An example of digital pictures taken from each excavation.

**Figure 10.** Comparison of the experiment and DEM results.

The percentage of particles in each layer computed for both field soil and DEM simulations is shown in Figure 11. Results show that increasing forward speed decreases the topsoil buried to 200–300 mm layers. Decreases of 16.1% and 27.9% were found in field tests and DEM simulations, respectively. It was also found that increasing speed creates a less uniform topsoil burial and creates a bulge in concentration between 100–200 mm depths. The concentration of the bulge between 100–200 mm increases with speed. It was found that the bulge concentration between 100–200 mm increases by 21.3% and 15.3% in the field test and DEM simulation results, respectively. It was also determined that only a small portion of the topsoil was buried in 250–300 mm layers at both speeds (12.4% and 7.4% at 5 and 7.5 km h<sup>−</sup>1, respectively).

(**b**)

**Figure 11.** Percentage of particles in each layer for an experiment vs. DEM simulation for 300 mm depth and (**a**) 5 km h−<sup>1</sup> and (**b**) 7.5 km h−<sup>1</sup> forward speeds.

As shown in Figure 11, the percentage of the topsoil buried at 5 and 7.5 km h−<sup>1</sup> forward speeds in the field was comparable to topsoil buried at 5 and 7.5 km h−<sup>1</sup> forward speeds in the DEM. The good correlation between the test and DEM results (R2 = 0.68) also validates this finding (Figure 12).

**Figure 12.** Correlation between experiment vs, DEM simulation. The cross dots represent the test vs. DEM predicted percentage of burial (as (x,y)), and, the black line represents the perfect prediction (y = x or R2 = 1).

The results summarised above indicate that there is a need to improve the design of the disc plough for improved topsoil burial. The results also showed that DEM-MBD coupling could be used to simulate the soil–disc plough interaction, which also helps to optimise the operational and design characteristics of the one-way plough.

The differences between the DEM and field results can be attributed to the experimental error during the excavation process when taking the images (each time a slice was dug, there was always a chance that some soil particles could fall into the sliced section) and the larger DEM particles used in the simulations. Using a more powerful computer and more software licences might help to improve the accuracy of the results. The simulation results may also be able to be improved by further refining the DEM-calibrated parameters to better suit the larger than actual particles.

#### **4. Conclusions**

This study simulated the interaction between soil and a modified one-way disc plough using DEM to investigate the topsoil burial. The rotational speed of the disc plough (as a passively driven tool) was precited using DEM-MBD co-simulation. The results were validated using field experiments. The results showed that the rotational speed of the disc plough was predicted with a maximum relative error of 6.9%. In contrast, the coefficient of determination between the test and DEM-predicted topsoil burial was R<sup>2</sup> = 0.68. This shows that the developed DEM and DEM-MBD models can be used to model the soil–disc plough interaction.

It was also found from the test and DEM results that the used, modified one-way disc plough does not (1) effectively bury the topsoil into a deeper layer or (2) provide uniform topsoil burial throughout the surface. Therefore, future work is required to improve the design and geometry of the disc plough to achieve a more effective topsoil burial using DEM.

**Author Contributions:** Conceptualisation, M.U.; methodology, M.U.; validation, M.U.; formal analysis, M.U.; writing—original draft preparation, M.U.; writing—review and editing, M.U. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Australia's Grains Research and Development Corporation (GRDC) and the Department of Agriculture and Food Western Australia (DAWFA) (Project Number: DAW 00244).

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon request from the first author, M.U.

**Acknowledgments:** The research was funded by Australia's Grains Research and Development Corporation (GRDC) and the Department of Agriculture and Food Western Australia (DAWFA) DAW 00244. The technical and field work assistance of Dr Chris Saunders is gratefully acknowledged.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
