Hybrid centrifugal spreading model to study the fertiliser spatial distribution and its assessment using the transverse coefficient of variation

A virtual twin-disc spreader is modelled combining theoretical motion equations and statistical information.Statistical distributions of input parameters are deduced from experimental measurements.A Monte Carlo process reproduces the random variability of fertiliser spread pattern deposition.Simulations demonstrate the influence of the application rate and the collection tray size on transverse CV. Studying centrifugal spreading by carrying out field or in-door experiments using fertiliser collection trays is tedious and labour intensive. This is particularly true when several implementation methods need to be compared, numerous replications are required or fertiliser sample characterisation is required. To circumvent cumbersome experiments, an alternative approach consists in performing in silico studies. In order to reach this objective, a hybrid centrifugal spreading model is designed by combining theoretical fertiliser motion equations with statistical information. The use of experimental measurements to characterise fertiliser properties, outlet velocity, angular mass flow distribution and spread pattern deposition, ensure a realistic calibration of the model. Based on this model, static spread patterns and transverse distributions are computed for a virtual twin-disc spreader. The number of fertiliser granules used to compute a spread pattern is deduced from the target application rate while the granule properties and their motion parameters are randomly selected from pre-established statistical distributions. This Monte Carlo process reproduces the random variability of fertiliser spread pattern depositions. Using this model, simulations demonstrate the mean and standard deviation of CV value decrease with the application rate. The CV mean value also decreases with the collection tray surface, while the standard deviation decreases with the collection tray length. Mathematical relationships are deduced from simulation results to express the mean and standard deviation of the CV as functions of the application rate and collection tray surface or length. The simulation model is also used to compare spreader test methods and study the influence of some fertiliser particles properties on the transverse distribution.


Introduction
In agriculture, the objective of mineral fertiliser supplies is to provide the right rate of nutrients to cultivated plants. Because of their low cost and high productivity, centrifugal spreaders are widely used for this application aiming to spread fertiliser at a target rate with an acceptable uniformity in the field. For 50 years, several works have demonstrated the negative effects of nonuniform spatial distributions concerning environmental impacts (Tissot et al., 2002) and yield or economical losses (Horrell et al., 1999;Jensen and Pesek, 1962;Miller et al., 2009;Richards and Hobson, 2013;Søgaard and Kierkegaard, 1994;Tissot et al., 1999). For the same decades, numerous works have been devoted to the measurement of fertiliser distributions, the assessment of distribution quality and the understanding of spread patterns. Throughout the world, transverse tray tests are traditionally performed to measure the spreading uniformity according to various standards such as: ISO Standard 5690/1 (1985), ASAE Standards S341.2 (1999); EN 13739-2 (2003); Spreadmark code of practice (New Zealand Fertiliser Quality Council (2015)) or ACCU Spread (Australian Fertiliser Services Association, 2001). The experimental transverse distribution is then used to compute the coefficient of variation CV after overlapping. This CV value is used to quantify the spreading quality, define the appropriate swath spacing according to the fertiliser and spreader setting, and thus certify the spreader bout width.
Some studies have addressed the comparison of transverse distribution measurement methods. Several works investigated the influence of the collection systems. Parish (1986)  manually-operated rotary spreader and two granular materials. The maximal effective swath width of this spreader was 4.3 m. Each test run consisted of three passes and three replications where carried out. Using the results obtained in this previous work, Parish and de Visser (1989) analysed the effect of the collection tray width on the CV value. In field, Parish et al. (1987) compared the crop response quality assessed by a horticulturist with fertiliser rates deduced from transverse distribution measurements. Three collection methods were compared using three replications for each test. All these studies demonstrated that major differences occurred in the measurement of transverse distributions depending on test methods. Therefore, the authors highlighted the importance of using the same test method for comparisons of spreader performances. Moreover regarding the low throwing distance of the spreader chosen for these studies and the low number of replications, these works illustrate the difficulties of carrying out such experiments.
To perform statistical comparisons of six international spreader tests, Jones et al. (2008) carried out a huge experimental work by using 18 transverse rows of 80 trays each. The experiments were carried out with urea, for three application rates and two replications so that 36 transverse distributions were obtained for each spreading situation. The bout width of the spreader was 15 m. Concerning the prediction of the certifiable working width, the authors concluded that the ACCU Spread test method (Australian Fertiliser Services Association, 2001) was superior to the other tested standards because it uses two rows of collector trays and multiple passes. Jones et al. (2008) concluded multiple rows of trays, multiple passes of the spreader and long trays can improve the accuracy of transverse tests.
Since the transverse distribution results from the combination of numerous parameters, it only provides a limited piece of information concerning the spread pattern. Thus, transverse tests are not efficient to study how mechanical parameters or fertiliser characteristics affect the 2D spread pattern deposition. This was illustrated by Piron and Miclet (2005) who showed that different 2D static spread patterns can yield to similar transverse patterns. Unfortunately, the measurement of the 2D static spread pattern is very tedious when a grid of collection trays is used, because of the wide size of spreader footprints and the high number of trays required to cover this area. Moreover, for indoor test, the high throwing distance of recent spreaders would require very expensive infrastructures. To circumvent these difficulties, Piron and Miclet (2005) developed a rotating test bench called CEMIB. With this method, the spreader is rotated during the spreading and a radial row of collection trays equipped with load cells records the cumulated mass of fertiliser according to the angular orientation of the spreader. The static spread pattern is then derived from Nomenclature a regression parameter A p particle frontal area, m 2 b regression parameter c regression parameter C d drag coefficient CV transverse coefficient of variation, % CV geom geometrical component of the CV, % CV k value of the CV obtained when the collection tray width is w k , % D continuous random variable, m d p fertiliser granule diameter, m d pi diameter of the ith fertiliser granule, m F D (d p ) cumulative frequency function of the granule diameter f D (d p ) probability density function of the granule diameter g acceleration due to gravity, m s À2 G D (d p ) cumulative mass distribution function of the granule diameter G M (h vane ) cumulative mass flow distribution with respect to the vane location g M (h vane ) mass flow distribution with respect to the vane location h vane height of the outer extremity of the vane, m K constant, m 3 K a aerodynamic coefficient, m À1 l tray length of the collection tray, m L w swath spacing, m m particle mass, kg m(d p ) mass of a granule of diameter d p , kg m i mass of the ith fertiliser granule, kg m tot total mass of fertiliser ejected by the two discs of the virtual spreader, kg n disc number of granules ejected by one disc of the virtual spreader (O, i, j, k) Cartesian frame centred on the disc centre, with j oriented in the travel direction q t target application rate, kg/ha q f in-field target rate, kg/ha r Pearson correlation coefficient the cumulated mass and the CV of the transverse distribution can be deduced (EN 13739-2, 2011). The measurement of the 2D spread pattern is of particular interest to improve the understanding of the spread pattern formation, the understanding of mechanical parameter effects and more generally to design new spreader. It is also useful to calibrate or validate spreading models. Recently, Cool et al. (2015) addressed the design of a simplified measurement technique to estimate the 2D static spread pattern in field, using a limited number of collection trays placed on a square or polar grid. The results of these two sampling techniques were compared with the results obtained with a transverse test. Tests were carried out with a spreader whose setting corresponded to 15 m bout width for ammonium nitrate fertiliser. Tests were performed for three fertilisers without replication. The authors observed large differences in the CV values deduced from the three measurement techniques and highlighted the importance of using the same measurement techniques to compare spread patterns. These experiments are tedious and do not make possible a sufficient number of replications to compare significant values. They illustrate the need of alternative approaches when the design or the assessment of new spreading quality measurement techniques is required.
The complexity and the labour-intensive nature of experimental measurements further increase when the study is not limited to the fertiliser mass distribution but aims to analyse the size or the nutrient formulation of the granules with respect to their spatial distribution. For example, very few studies investigated the effect of fertiliser particle size on spread distribution. Pettersen et al. (1991) studied the spatial distribution of fertiliser particle size using a twin disc spreader. Experiments were limited to the choice of one fertiliser, one spreader setting, one feeding flow rate and one measurement of the 2D stationary spread pattern. A set of 884 collected samples was analysed by image processing technique to draw the spatial distribution of the particle size. Thirty years later, Yule (2011) attempted to study the effect of fertiliser particle size on spread distribution. The transverse distribution of the percentage of particle size was drawn for two loads of superphosphate having different granule size distributions. Yule (2011) used these experimental results to simulate the transverse distribution of other materials with other particle size characteristics. Nevertheless, the author underlined the work is limited to representing only one particular spreading situation: one spreader with one fertiliser and one setting. Yule (2011) concluded that further work would be required but analysing each tray from field testing was too time-consuming and no laboratory measurement techniques were adapted for this kind of study at the present time.
The study of the spatial distribution of fertiliser granules according to their physical properties is of particular interest in the case of blended fertilisers. As these materials are produced by mixing mechanically single products, their components differ in the physical properties (size, shape and density). These differences can involve segregation of the fertiliser components during handling and spreading. This problem was already addressed by Hoffmeister et al. (1964). When fertiliser is applied with a centrifugal spreader, the differences in physical properties can affect the granule behaviour during the ballistic flight. Then, ballistic segregation can occur and yield heterogeneous spatial distribution of chemical elements. Several works have examined the ballistic segregation of blended fertilisers by carrying out field experiments (Miserque and Pirard, 2004;Tissot et al., 1999;Virk et al., 2013). All these studies required huge field tests to evaluate the mass and nutrient distribution.
In order to circumvent cumbersome experiments or reduce the number of experiments, some author attempted to develop new approaches based on modelling fertiliser granule motion. Recently, Antille et al. (2013) suggested the particle size range of new fertilisers could be designed to meet a target bout width and the authors proposed to model and simulate the ballistic flight to assess whether granule physical properties suited the spreading objective. Grafton et al. (2015) also suggested the use of a ballistic model to provide information reducing the risk of crop striping. In these recent works, the models were limited to the prediction of the landing distance of some individual granules projected by a spinning disc. Therefore, no spread pattern was computed so that no transverse distribution can be determined.
Numerous works have attempted to model the motion of fertiliser granules in the spreading process. Various mechanical models have been proposed to describe the motion of individual fertiliser granules on a spinner disc (Cunningham, 1963;Cunningham and Chao, 1967;Hofstee, 1995;Inns and Reece, 1962;Olieslagers, 1997;Patterson and Reece, 1962;Villette et al., 2005) and through the air (Antille et al., 2015;Mennel and Reece, 1963;Pitt et al., 1982). Concerning the motion on the disc, models using the discrete element method (DEM) have also been developed to take into account particle interactions (Casas et al., 2015;Coetzee and Lombard, 2011;Tijskens et al., 2005;Van Liedekerke et al., 2006). One drawback of DEM models is that they required input parameters which are difficult to obtain to characterise the physical behaviour of fertilisers. When results of simulations are compared with actual spread pattern depositions, they reach moderate success even when spreading distances are lower than 3 m (Coetzee and Lombard, 2011;Van Liedekerke et al., 2009). Consequently, at the present time, no model appears sufficiently advanced to correctly simulate actual spread pattern depositions. Moreover, to the best of our knowledge, no model reproduces the random variability observed from run to run in spreader test.
Despite some experimental studies, the comparison of spreader transverse tests using different collection trays or different test protocols is still difficult and the main conclusion is limited to the recommendation of using the same test to give sense to comparisons. For example, there is a lack of knowledge concerning the quantitative effect of the surface or shape of collection trays on the CV measurement. Similarly, the effect on the CV of increasing the number of runs or reducing the speed travel of some standard tests has not been studied. In addition, the effect of the application rate on the measurement of the CV value has never been studied. This lack of information results from the difficulty, not to say impossibility, of carrying out adapted experiments with enough replications. The same difficulty limited the production of knowledge on the spatial distribution of fertiliser particle size and on the ballistic segregation for blended fertilisers. An alternative solution lies in the use of models to simulate the physical phenomena and carry out in silico experimental studies. The main advantage of this approach is that it avoids practical and time limitations so that statistical parameters can be deduced from replications. Nevertheless, this implies that simulation models have to reproduce the stochastic nature of fertiliser dispersal processes.
The aim of this paper is to design such a model for simulating realistic fertiliser spread patterns and providing new solutions to carry out numerical experiments. This hybrid model combines a mechanistic approach based on the use of mechanical relationships and a stochastic approach based on the use of statistical distributions of input parameters. The simulation model is used to study the meaning of the CV value deduced from transverse tests according to the target application rate and the test method. The paper also presents an insight into the influence of particle size distributions and particle drag coefficients on transverse distributions.

Materials and methods
The particularity of the Hybrid Centrifugal Spreading Model HCSM lies in combining some theoretical motion models with experimental data obtained at various steps of the spreading process. This model assumes the spread pattern deposition is affected by the outlet velocity of the particles when they leave the spinning disc, the angular mass flow distribution around the disc and the fertiliser particle properties (specific density, size-distribution, drag coefficient). Other parameters such as wind, landform or guidance error are not taken into account in the present paper.

Model of granule motion on and off the spinning disc
Concerning the motion of granules on the spinning disc, the HCSM considers kinematic relationships between the disc configuration, the outlet angles and the outlet velocity components. In this section, the motion model is described for a clockwise spinning disc. Fig. 1 presents the main geometrical parameters used to describe the disc and the motion of the granules when they leave the vane.
Let (O, i, j, k) be a three dimensional right-handed Cartesian coordinate system having its origin O on the rotational axle of the disc and with j pointing in the travel direction.
In this coordinate system, the location (x out , y out , z out ) of the granule when it leaves the vane is: where r vane is the radius of the vane, h vane is the angular location of the vane, h vane is the height of the outer extremity of the vane. For a concave disc, the relationship between the vertical outlet angle X out of the particle when it leaves the vane, the horizontal outlet angle h out , the vertical angle of the vane X vane and the pitch angle of the vane a lv , is as follows: The horizontal component of the outlet velocity is also deduced from h out as follows: where x is the rotational speed of the disc.
Then, the outlet velocity is deduced: According to Fig. 1, for a clockwise rotating disc, the expression of the horizontal orientation h traj of the outlet velocity with respect to i is as follows: The components of the outlet velocity in the Cartesian coordinate system (O, i, j, k) are: In this model, owing to Eq. (2), the vertical outlet angle is not taken as the vertical angle of the vane unlike some other models suggested in the literature (Gomez-Gil et al., 2009;Olieslagers et al., 1996). This avoids coarse approximations in computing the initial conditions of the ballistic flight.
During the ballistic flight the model considers that the forces acting on the granule are only the gravity force and the drag force due to the motion of the granule through immobile air. This simple ballistic model had been used in numerous works such as Mennel and Reece (1963), Pitt et al. (1982), Griffis et al. (1983), Olieslagers et al. (1996), Grift and Hofstee (2002), Reumers et al. (2003), Aphale et al. (2003), Bradley and Farnish (2005). In the three dimensional Cartesian coordinate system, the motion in the air is described by the following differential equations: where x, y, z are the coordinates of the granule; vx, vy, vz are the velocity components of the granule, g is the acceleration due to gravity and K a is as follows: where m is the granule mass, C d is the drag coefficient, A p is the projected surface area of the granule, q air is the air density. In this study, A p is computed for spherical shapes.

Spreading process parameters
The HCSM uses experimental data measured at the beginning and at the end of the ballistic flight. This provides the initial parameters of the ballistic flight (i.e. outlet velocity) but also makes possible the estimation of the drag coefficient during the flight. The combination of these experimental measurements with mechanical models is of particular interest to take into account the actual behaviour of the fertiliser in the spreading process Reumers et al., 2003) and provide realistic simulations. In the spreading simulations, the motion of a high number of fertiliser granules is computed. To compute one simulation some variables are taken constant for all granules while some other variables assign random values for each granule. These last variables are associated with probability distributions. Thus, some probability distribution functions or the corresponding cumulative distribution functions need to be defined.
This section describes the input parameters used in the simulations and their measurement methods.

Experimental spreading device and spread pattern deposition
A custom-made spreader was used for the experimental measurements. This spreader consisted of a single clockwise rotating disc. This concave disc was equipped with two radial vanes (X vane was 13.5°; a lv was 0°; r vane was 0.395 m) and was spinning at 810 rpm. The height h vane of the outer extremity of the vanes was 0.8 m. The feeding mass flow was 0.97 kg/s. The stationary spread pattern obtained with this spreader and with ammonium nitrate was measured using a rotating test bench called CEMIB. This measurement device consisted in a rotating carrier and a motionless line of 80 collection trays. Each tray was equipped of a load cell and had a square collection area of 0.5 Â 0.5 m. The design and the advantages of this test bench are detailed in Piron and Miclet (2005) and Piron et al. (2010).
During the spreading, the spreader carrier rotated at a constant speed of 3.1°/s, so that the whole spread pattern passed above the collection tray row. During the rotation, the cumulated mass collected by each tray was recorded and the weight values were stored with the corresponding orientation angle of the carrier (measured by an angular optical encoder). The acquisition frequency was 10 Hz. During the whole carrier rotation, the total fertiliser mass ejected by the spreader was approximately 38 kg and the total mass collected by all the trays was approximately 0.45 kg.
Using this measurement device, the resulting raw data was a matrix in which each line corresponded to an orientation angle and each column corresponded to the cumulative fertiliser mass collected in each tray. The fertiliser mass collected at each angular location was then derived from cumulative measurements. Since the collection areas of all trays are the same, the spatial density of the fertiliser deposition was directly deduced from previous data in polar coordinates. Then, the spread pattern deposition was computed by the CEMIB algorithm, with respect to the disc centre in Cartesian coordinates, using a mathematical interpolation and a sampling interval of 0.25 Â 0.25 m. Fig. 2 presents the spread pattern deposition obtained with the experimental spreader and the ammonium nitrate fertiliser used in this study. This spread pattern is taken as reference data for the study.
The analysis of the spread pattern deposition shows that the mean radius of the spread pattern slightly increases with the rotation of the vane (i.e. from the beginning to the end of the spreading angular sector).

Horizontal outlet angle distribution
The horizontal outlet angles h out were measured using an imaging system based on the processing of motion-blurred images. In this acquisition technique, the exposure time is long relative to the velocity of fertiliser granules so that the granule displacements appear as streaks across the image. The horizontal outlet angles were derived from the distance between these streaks and the disc axle. The imaging system and the image processing are detailed in Villette et al. (2008).
Images were captured with a monochrome CCD camera (Sony XCD-SX910), equipped with a 6 mm lens. The camera was approximately placed at 0.7 m above the upper corner of the vane and approximately above the central part of the spreading angular sector. The optical axis of the camera was set parallel to the disc axle at a distance of approximately 0.5 m from this axle. Using a set of 300 images, the horizontal outlet angle was measured for trajectories selected near the principal point of the image (i.e. the point corresponding to the view axis in the image) to improve the measurement accuracy and avoid geometrical bias. Thus, 2280 trajectories lying in a 10°spreading angular sector were used to estimate the horizontal outlet angle. The angular location of the vane h vane corresponding to the middle of this sector was À20°.
The mean value lh out was 40.2°and the standard deviation rh out was 0.85°.
Considering the histogram of the measured value (Fig. 3), the probability density function was chosen as a normal distribution defined by the two parameters: lh out and rh out .
In order to take into account, the slight increase of the spread pattern radius with respect to the rotation of the vane, the mean value lh out was modelled by the following linear relationship: where lh out and h vane are expressed in degrees.

Vertical outlet angle distribution
For a given value of the horizontal outlet angle h out , the mean value of the corresponding vertical outlet angle X out was provided by Eq. (2). Nevertheless, this relationship did not model the dispersion of X out around its mean value. Consequently this dispersion was measured separately by analysing the vertical distribution of the mass flow. The experimental method consisted in recording granule impacts on a vertical screen placed in the vicinity of the spinning disc. The method is the same than the one described in Villette et al. (2013), except a simple flat screen was used instead of a cylindrical screen. Moreover a shutter system was added to control the exposure time of the screen to fertiliser impacts. The screen was covered with a paper of A4 size, a carbon film and a protective film, so that granules hitting the screen produced impact marks on the recording paper (Fig. 4). After the exposition of the recording paper to granule shocks, the paper was digitalised and a dedicated image processing was carried out to analyse the vertical distribution of the impacts. The details of the mathematical model and algorithms used to process impact records can be found in Villette et al. (2013). Considering the curve of the impact vertical distribution (Fig. 4), the probability density function was chosen as a normal distribution. Placing the recording screen at 1.08 m and 2.09 m from the axle of the spinning disc, 5 replications of impact recording were carried out at each distance. The mean values of the standard deviations of the impact heights were respectively: 13.3 mm and 25.9 mm. Then, the standard deviation rX out of the vertical outlet angle was estimated at 0.7°.

Angular mass flow distribution
The angular mass flow distribution was computed at the outer extremity of the vane, as a function of the angular location of the vane. The distribution was deduced from the spread pattern deposition and the horizontal outlet angle. The whole spreading angular sector was sampled each degree. For each angular location of the vane, the theoretical horizontal direction was computed using Eq. (5) and the relative fertiliser quantity was computed for each sampled angular sector from the outer extremity of the vane to a range of 30 m with a sampling interval of 0.25 m. Thus, the spread pattern deposition was computed as a function of the spreading distance (from the extremity of the vane to the landing point) and the angular location of the vane (Fig. 5). The angular mass flow dis-tribution g M (h vane ) at the extremity of the vane (Fig. 5) was obtained by summing the relative mass of fertiliser obtained for each vane location (whatever the spreading distance). The cumulative mass flow distribution with respect to the vane location G M (h vane ) was deduced from g M (h vane ) as follows: This distribution was used in the HCSM to compute the probability of the fertiliser mass ejected for each angular position of the vane for the clockwise spinning disc of the virtual spreader.

Fertiliser parameters
In this study, ammonium nitrate was used for actual experiments. The same fertiliser characteristics were used in numerical simulations computed with the HSCM. The shape of the granules was assumed to be spherical.

Specific density
The density of the fertiliser granules was deduced from weighing a material bulk volume and weighing anew the same bulk volume after completing with a liquid of known density. Then, the volume of the granules is deduced and the granule density is calculated. For the ammonium nitrate used in this study, the specific density q was 1563 kg m À3 .

Granule diameter distribution
The particle size analysis was performed with a sieving test according to the European standard EN 1235/A1 (2003). This provided the cumulative mass distribution function. Then, a twoparameter lognormal distribution was used to describe the distribution.
Thus, the normalised cumulative mass distribution (value range from 0 to 1) was fit with the following function: where d p is the granule diameter, erf() is the error function, l ln and r ln are the two fitting parameters (corresponding to the mean and standard deviation of the variable's natural logarithm). The derivative function of G D (d p ) is: Considering the probability density function f D (d p ) and the cumulative frequency function F D (d p ) of the random variable D, the probability of having the granule diameter D lower than d p is: Using the probability density function f D (d p ), the cumulative mass distribution G D (d p ) is also expressed as follows: where m(d p ) is the mass of granules of diameter d p .
Assuming the mass of the granule is proportional to d 3 p , Eq. (14) yields: where K ¼ R þ1 0 f D ðnÞ Â n 3 dn. This provides: Combining Eq. (13) and Eq. (16), the cumulative frequency function of the granule diameter is finally obtained as follows: In practice, the integration of Eq. (17) is computed numerically and the constant K is determined so that F D (d p ) is 1 when d p tends to infinity.
For the ammonium nitrate fertiliser used in this study, Fig. 6 presents the results of the sieve test with the corresponding cumu-lative mass distribution G D (d p ) and the cumulative frequency function F D (d p ).
The establishment of the cumulative frequency function F D (d p ) is of particular interest to select efficiently random diameter values corresponding to a random set of granules whose the mass distribution respects the fertiliser sieve test.

Drag coefficient in the air
To compute the ballistic flight of fertiliser particles, many researchers assumed the drag coefficient C d to be constant (Coetzee and Lombard, 2011;Grafton et al., 2015;Grift and Hofstee, 2002;Olieslagers et al., 1996;Pitt et al., 1982), while some other authors tried to improve the description of the ballistic flight by considering changes of the C d value during the motion (Antille et al., 2015).  In some works, the fertiliser particles were approximated as perfect spheres. Consequently, considering a turbulent flow regime, the C d value of the fertiliser granules was chosen at 0.44 (Coetzee and Lombard, 2011;Olieslagers et al., 1996;Walker et al., 1997). Taking into account the influence of shape and texture of fertiliser granules on their aerodynamic behaviour, some other authors chose higher C d values. For instance, Pitt et al. (1982) chose 0.46 for ammonium nitrate. Comparing modelled and measured fall time, Grift and Hofstee (2002) suggested multiplying the diameter of the equivalent sphere by a correction factor (named ''qfactor") ranged from 0 to 1.
In the present study, the drag coefficient was assumed constant during the ballistic flight. The value of C d was chosen by comparing the reference spread pattern (i.e. obtained with the CEMIB test bench) with simulated spread patterns computed for various C d values. Thus, for the ammonium nitrate used in this study, the value of the drag coefficient was estimated to 0.47. Moreover, the air density q air was assumed to be 1.21 kg m À3 in the spreading condition.

The virtual spreader
The virtual spreader considered for the simulations was a twin disc spreader for which the spacing between the two disc axles s disc was 1 m. Both discs had the same angular speed. The right disc rotated in the counter-clockwise direction while the left one rotated in the clockwise direction. Each disc of the spreader was fed by the same mass flow of fertiliser.
The setting of the virtual spreader consisted in modifying the angular location of the feeding point on each disc. With this setting mechanism, rotating the angular location of a set for the left disc involves the rotation of the left spread pattern in the same direction and with the same angle a set .

Static spread pattern simulation
Considering the virtual spreader, static spread patterns were computed using the HCSM and a Monte Carlo process. This approach consisted in computing the motion of a high number of fertiliser granules for which several characteristics were randomly drawn from pre-established statistical distributions. Simulations were implemented with Matlab (2005), and used the random number generator of this software. For normally distributed variables, the values were obtained using the randn function. In the case of other arbitrary distributions, the selection of random values was performed in two steps. First, random numbers were generated with a uniform distribution using the rand function on the range 0-1. Second, final random values were deduced from these random numbers by inversing the cumulative frequency function of the specified distribution.
For a given mass m tot of fertiliser, the computation of the spread pattern was decomposed in computing the left and the right spread patterns independently. For the left disc, the Monte Carlo simulation consisted of the following steps.
First, a set of virtual granules was generated by drawing a set of diameter values from the fertiliser diameter distribution using the cumulative frequency function F D . Then, the mass m i of each granule was computed as follows: where d pi is the diameter of the granule of mass m i . The total number n disc of granules ejected by the disc was adjusted so that: Second, the initial conditions of the ballistic flight were assigned to each granule independently from its diameter. For each granule, the values of the different variables were assigned as follows: (1) The angular location of the vane h vane corresponding to the granule ejection was randomly selected using the cumulative mass flow distribution G D .
(2) The corresponding coordinates of the ejection point (x out , y out , z out ) were deduced from h vane using Eq. (1). (3) The corresponding horizontal outlet angle h out was drawn from the normal distribution parametrized by lh out and rh out , where lh out was deduced from the vane location h vane using Eq. (9). (4) The corresponding vertical outlet angle X out was drawn from the normal distribution parametrized by lX out and rX out , where lX out was the vertical outlet angle deduced from h out using Eq. (2).
(5) The corresponding outlet velocity v out and its 3Dcomponents (vx out , vy out , vz out ) were deduced from h vane , h out and X out using successively Eqs. (3)-(6).
Third, the coordinates of the landing point of each granule were computed by solving Eq. (7) with the initial conditions of flight x out , y out , z out and vx out , vy out , vz out . Then, the setting of the spreader is taken into account by computing the coordinates of the landing points with the rotation angle a set around the disc axle.
The spread pattern produced by the right disc is computed by (1) using anew the same process to generate a second spread pattern for another set of granules; (2) changing the sign of the x coordinates of this second spread pattern.
The global spread pattern resulting from the twin-disc virtual spreader is finally deduced from the left and right spread patterns after translating the coordinates of the granules by half the disc spacing s disc in the left or right direction. This global spread pattern is defined by a set of granules for which each mass and each landing position is perfectly known.

Transverse distribution
The transverse distribution is deduced from the static spread pattern by considering a virtual row of collection trays placed continuously along a line perpendicular to the travel axis (along the xaxis) of the virtual spreader. For a given swath spacing L w , several static spread patterns were computed and translated on the right and the left at a multiple of L w of the central pass to reproduce the overlapping. The successive spread patterns were oriented to simulate overlaps resulting from adjacent swaths applied in alternate directions (i.e. back and forth mode).
Depending on the x-value of each granule of the spread patterns, the granules were affected to the corresponding collection trays, so that the sub-set of granules virtually collected by each tray was perfectly known in terms of granule masses and granule diameters. The transverse distribution of the fertiliser mass was obtained by summing the mass of all the granules virtually collected by each collection tray. Then, the transverse coefficient of variation CV is deduced from this mass distribution by dividing the standard deviation by the mean. In this paper the CV is expressed in percentage.
For a target application rate q t , the total mass m tot of fertiliser used to compute the static spread pattern was determined as the product of the application rate by the collection surface on the swath spacing. Calling l tray the length of the collection trays (measured in the travel direction), the total mass m tot is as follows: 3. Results and discussion

Comparison of the measured and simulated spread pattern
Although the objective of the HCSM was not the perfect description of the physical phenomena of the spreading process or the perfect reproduction of the experimental spread pattern, it is important to ensure the simulated spread pattern was in accordance with experimental results. Thus, the first simulation consisted in computing the spread pattern for a single disc for the same spreading conditions than those carried out with the experimental spreading device.
The 2D-representation of the spread pattern deposition was obtained by considering a grid with a sampling interval Dw grid = 0.25 m and Dl grid = 0.25 m in the transverse and longitudinal directions. Depending on the granule coordinates, the granules located in each grid cell were identified and the sum of their masses was affected to the corresponding cell. This matrix was the raw Cartesian representation of the spread pattern deposition.
Since the experimental measurement device (CEMIB) was a rotating system, and since the processing of the polar data included some interpolation and regularization steps, it was difficult to compare visually the raw Cartesian representation of the simulated spread pattern with the interpolated experimental measurement (especially for low fertiliser amount).
Thus, the raw Cartesian representations of the simulated spread pattern had been sampled in a polar coordinate system, regularized by a Gaussian filter and then re-interpolated into a Cartesian coordinate system to simulate the effect of the CEMIB acquisition and data processing. For three different amounts of fertiliser, Fig. 7 shows the raw Cartesian representation of the simulated spread pattern deposition and its representation when interpolations are applied in an intermediate polar system. Fig. 7 shows that the local variability increases inside the spread pattern when the fertiliser amount decreases. This corresponds to the well-known random variability observed in CV measurement from run to run. Moreover, when the fertiliser amount is 0.45 kg the representation is in good accordance with the reference spread pattern (Fig. 2) measured with the rotating test bench (when the same fertiliser amount was collected). A better comparison would have been obtained by modelling the rotating acquisition system, but this was out of scope of this paper.
When the fertiliser amount is very high (i.e. 50 kg), the relative local variability inside the spread pattern is reduced. Then, the representations of the spread patterns are similar whatever the use of intermediate steps in a polar system or not.

Setting map of the virtual spreader
A set of simulations were performed to draw the setting map of the virtual spreader. Thus, the spread pattern was computed by considering 10 6 particles (i.e. a total mass of approximately 26.8 kg) ejected by each disc. The value of the setting angle was from À30°to À60°, and for each value the transverse CV was computed for a set of swath spacing (from 15 to 45 m). The size of the virtual collection trays was 0.5 Â 0.5 m each. Fig. 8 shows the setting map deduced from the simulations. This map represents the CV value obtained with the virtual spreader with respect to setting angle and swath spacing. Since all the CV values were computed for a very high number of particles, these values reflect the geometry qualities or defects of the spread patterns related to the swath spacing and do not take into account other transverse variabilities that occur for lower application rates. This CV value is an estimation of CV geom which is only due to the geometrical shape of the spread pattern (for a specified swath spacing) regardless of the application rate. CV geom is the value of CV that would be reached if the application rate tended to infinity.

Influence of the application rate and collection tray size on the CV value
Simulations have been performed to investigate the effects of the application rate q and the collection tray size on CV measurements. For the simulations, the sizes of the collection trays (length Â width) were: 1 Â 1 m, 1 Â 0.5 m, 1 Â 0.25 m, 0.5 Â 1 m, 0.5 Â 0.5 m, 0.5 Â 0.25 m, 0.25 Â 1 m, and 0.25 Â 0.5 m. The application rates were: 50, 100, 200, 400, 800, 3000 kg/ha and for each rate the number of replication runs was respectively 2000, 1000, 500, 500, 400, and 200. The simulations have been computed for six spreading situations (SS1-SS6) located on the setting map (Fig. 8). Conditions SS1-SS4 correspond to a setting angle of À41°a nd swath spacing of respectively 16, 26, 32 and 42 m. Conditions SS5 and SS6 correspond to a setting angle of À55°and swath spacing of respectively 36 and 42 m. These situations have been chosen to illustrate various setting conditions: optimal settings and swath spacing (SS2 and SS6) and inadequate settings or swath spacing.
In the case of the situation SS2, Fig. 9 demonstrates that the mean value l CV and the standard deviation r CV of the CV increase when the application rate decreases. The curves also show that l CV and r CV depend on the size of the collection trays. The mean value l CV increases when the surface of the collection trays decreases, while the standard deviation r CV increases when the length of the collection trays decreases.
Concerning l CV , as shown in Fig. 9, it appears that the values are very similar when they are deduced from simulations computed with the same tray surface. Comparing the results obtained for the six spreading situations (SS1-SS6) at the six application rates, the maximum difference observed on the CV values is 0.36% for the following tray dimensions 1 Â 0.25 m, 0.5 Â 0.5 m, and 0.25 Â 1 m (this maximum difference is obtained for SS4 at 200 kg/ha). For the dimensions 1 Â 0.5 m and 0.5 Â 1 m, the maximum difference is 0.24% (obtained for SS4 at 3000 kg/ha). For the dimensions 0.5 Â 0.25 m and 0.25 Â 0.5 m, the maximum difference is 0.23% (for SS1 at 50 kg/ha). These differences are very low regarding the traditional range of CV values encountered in practice or regarding l CV values encountered here for spreading situations and all application rates (from 1.7% to 46.3%). Furthermore, considering one spreading situation, the mean values of the CV obtained for all tray size tend to converge.
All these observations lead to consider that the mean value of the CV depends on two components. One component reflects the spreading situation depending on the setting and on the swath spacing. This component is expressed by the constant value to which l CV tends for high application rates. The second component reflects the influence of the application rate and, more precisely, the influence of the mass collected in the trays (depending on the rate and the tray surface). Consequently, the mean value of the CV is expressed as a function of q and S tray as follows: where q is expressed in kg/ha, S tray is expressed in m 2 , a and b are two coefficients. This relationship is used as a regression model to fit the data obtained for the six studied spreading situations (SS1-SS6).  Table 1 shows the values of the parameters a, b and the correlation coefficient resulting from the use of Eq. (21) to fit the data. The regression curves are drawn on Fig. 10. The values of the correlation coefficient r demonstrate that Eq. (21) accurately describes the relationship between l CV , q and S tray . The lowest value of correlation coefficients is obtained for SS4 but is still higher than 0.996.
Considering Eq. (21), the value of l CV tends to b when the application rate approaches infinity. Thus, the parameter b corresponds to the CV geom , which depends on the geometrical shape of the spread pattern (for a specified swath spacing) regardless of the application rate. Regarding Table 1, the values of the parameter a are very similar whatever the spreading situation. Some additional simulations shows that this value depends on the particle size distribution of the fertiliser.
Another interesting aspect of the finding expressed in Eq. (21) is that, whatever the setting of the machine, the mean value of the CV is higher than a limit defined by the application rate and the collection tray surface. This limit is as follows: This means that trying to set the machine to obtain a CV below this limit does not make sense. Conversely, measuring CV values significantly higher than this value indicates that the setting or the design of the machine could be optimized to improve the spreading quality. Nevertheless, the practical use of this threshold remains dependent on the accuracy of the CV estimation regarding the measurement variability.
Concerning the standard deviation r CV of the CV, Fig. 9 shows that the values are close when they are deduced from simulations computed with the same collection tray length.
Comparing the results obtained for the six spreading situations (SS1-SS6) at the six application rates, the maximum difference observed on the standard deviation is 0.25% for the following tray dimensions 1 Â 1 m, 1 Â 0.25 m, 1 Â 0.5 m. For the dimensions 0.5 Â 1 m, 0.5 Â 0.5 m, 0.5 Â 0.25 m, the maximum difference is 0.19%. For the dimensions 0.25 Â 1 m, 0.25 Â 0.5 m, the maximum difference is 0.2%.
For each collection tray size and for each spreading situation, the standard deviation of the CV is well fitted by the following expression: Fig. 7. Simulated spread patterns for three fertiliser amounts: 0.1 kg, 0.45 kg, 50 kg (from top to bottom) and two kinds of representations: raw Cartesian representation of the simulated spread pattern (left) and after sampling and interpolating in a polar coordinate system (right) to reflect CEMIB data processing. The graduation of the colour scale reflects the fertiliser amount lying on the sampling area (0.25 Â 0.25 m) expressed in percentage of the total mass. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) The range of values of the parameter c and the range of value of correlation coefficients are presented in Table 2, when q is expressed in kg/ha.
In the literature, very few studies addressed the problem of the influence of the application rate or the collection tray size on the CV value. This is due to the difficulty in performing a high number of replications and in maintaining constant spreading conditions to establish a unbiased relationship when actual experiments have to be carried out.
Since the quality of the spreading results from numerous combined parameters, simulations afford the possibility of not only avoiding perturbations (e.g. humidity, wind or fertiliser property variations), but also analysing the studied parameter independently from the others. This is the case in this section: the effect of the application rate on the CV is studied without any change in the outlet angle distribution or in the angular mass flow distribution. This specific study is very difficult to carry out in practice, because the global shape of the spread pattern can be modified when the feeding flow rate is modified (Fulton et al., 2001) due    to the change in the feeding area on the spinning disc , or in the vane loading (Villette et al., 2012). Parish and de Visser (1989) suggested that given a collection tray width w 1 , and a resulting coefficient of variation CV 1 , a collection tray width of w 2 will result in a new coefficient of variation CV 2 , related to CV 1 as follows: The length of the collection tray was implicitly the same for the two kinds of trays. The authors underlined that, because of variability, this equation should be effective in dealing with averages of multiple tests.
Using Eq. (21), the ratio of the mean value of CV 2 on the mean value of CV 1 is as follows: where l is the length of the collection trays (which is the same for the collection trays used to measure CV 1 and CV 2 ). It appears that the relationship between CV 2 and CV 1 is more complex than the one suggested by Parish and de Visser (1989) and that Eq. (24) is not correct when the geometrical component of the CV (CV geom ) is not null. Nevertheless, in the case of a good quality spreading, CV geom1 and CV geom2 are low and are negligible compared to the component related to the effect of the application rate. In the case of a good quality spreading, Eq.(25) yields: This demonstrates that the relationship proposed by Parish and de Visser (1989) is a correct approximation only when the setting of the machine is very good, or when the CV value is widely due to the effect of the application rate (low application rate).
Simulation findings demonstrate the variability in CV measurement (i.e. r CV ) decreases with the application rate and the length of the collection trays used for the measurement. These results are in perfect accordance with recommendations of Jones et al. (2008) who concluded that multiple rows of trays, multiple passes and long trays can reduce experimental variability and improve the accuracy of bout width calculation. It is worth noting that doubling the number of passes corresponds to doubling the apparent application rate collected by the tray. Concerning the collection trays used for transverse tests, the European Standard EN 13739-2 (2011) recommends the size (length Â width) of 0.5 Â 0.5 m but also permits 1 Â 0.25 m. The simulation findings confirmed that the mean value expected for the CV is the same when these two collection tray sizes are used. Nevertheless, simulations demonstrated these two collection devices are not equivalent regarding the variability in CV measurement. Thus, the variability is reduced when the size of 1 Â 0.25 m is used and the confidence in bout width calculation is improved. This illustrates that simulations would be of practical interest when standard revision process are launched.

Influence of the test method on the CV value
Simulations had been carried out to compare the values of CV when they were obtained following two different measurement methods. The first method was a simple ''in-field" measurement consisting in measuring the CV when the machine was driven at the forward speed of 10 km/h and was set to apply the in-field target rate q f . The second method was a ''standard test" performed following the European Standard EN 13739-2 (2011). In this stan-dard test, the machine was driven at 4 km/h, the number of runs for each measurement was two, and the flow adjustment using 4 km/h was set to correspond to the flow rate obtained at a forward driving speed of 10 km/h. Thus, in practice, the application rate used in the simulation program for the virtual standard test was: Following the European Standard, two replications were done for each virtual standard test before computing the corresponding CV. Moreover, the mirror image of the transverse distribution of the central pass was used to compute the overlapped distribution with the adjacent passes on the swath spacing. In contrast, concerning the virtual ''in-field" measurement, no replication is done before computing the corresponding CV and the transverse distributions for each pass were completely independent (i.e. no mirror image was used). For both procedures, the same transverse line of collection trays was used. The size of these virtual trays was 0.5 Â 0.5 m.
For the spreading situation SS2, Fig. 11 presents the mean value and error bar (twice the standard deviation) of the CV with respect to the application rate. The mean value and the standard deviation of the CV decrease with the application rate. The figure also shows that the CV is lower when it is measured with the standard procedure than when it is measured with the simple in-field test. Moreover, the standard deviation is lower when it is measured following the standard procedure (the standard deviation is approximately divided by 2). These results illustrate that the use of the standard test reduce the variability on the measurement of the CV. Nevertheless, the standard test underestimates the value of the CV with respect to the value that would be obtained in field by considering the actual target application rate (i.e. with the actual application rate at the actual forward speed and without any replication). In this example, the ''in-field" CV is at least twice the ''standard" CV when the application rate is lower than 490 kg/ha.
The differences observed here only result from the measurement procedure since no additional perturbation (i.e. wind, ground topography, ground irregularity, guidance error. . .) was taken into account for the ''in field" measurements. Thus, the ''in-field" CV defined in this section should not be confused with the ''field CV" defined by Lawrence and Yule (2007). Fig. 11. Mean value and error bar (twice the standard deviation) of the CV with respect to the application rate for the spreading situation SS2. Simulations are performed with collection trays of 0.5 Â 0.5 m for virtual ''in-field" measurements (continuous line) and for virtual ''standard" measurements (dotted line).

Influence of granule size and drag coefficient
In this article, one particularity of the spread pattern simulations was that each fertiliser granule was tracked during the whole virtual spreading process. Consequently, at the end of the process, when all granules lied on the ground, the location and the diameter of each granule were perfectly known. Thus, simulations were used to study how fertiliser particles contribute to the transverse distribution in relation to their diameters. Fig. 12 presents the transverse distributions of the fertiliser sieve fractions for two spreading situations: SS2 and SS5. For each situation, simulations were performed with 10 6 fertiliser granules per disc. The size of the virtual collection trays was 0.5 Â 0.5 m. In the case of situation SS2, the proportions of each diameter class are approximately kept constant on all the working width (26 m) and correspond to the particle size analysis presented in Fig. 12. In this situation, the transverse distribution has a triangular shape before overlapping and adjacent passes overlap on a large part of the working width. In contrast, in the case of situation SS6, the proportions of diameter classes are modified at the extremities of the working width with small diameter vanishing. In this last situation, the swath spacing is 42 m, the transverse distribution has a trapezoidal shape before overlapping and adjacent passes overlap only on a low part of the working width. Thus, the overlapping area only concerns the external part of the spread pattern and the landing points of biggest granules.
In the case of situation SS6, the same rate (i.e. same mass per surface unit) is applied locally at 2.5 m and 21 m from the centreline of the virtual spreader, but the number of granules per surface unit is not the same. For the same application rate, the number of granules decreases when their sizes increase. Consequently, this affects the spatial variability of the fertiliser supply at very small scales. A further characterisation of this effect is worth of studying but is out of scope for this article.
The study of spreading segregation is also of particular interest for blended fertilisers. The HCSM is an interesting tool to investigate how ballistic segregation affects the spatial distribution of each fertiliser components. To illustrate this aspect, a simulation was performed by considering two fertiliser components and one spreading situation. The first fertiliser component corresponded to ammonium nitrate whose characteristics had been described and used in the previous sections. The second fertiliser component corresponded to a fictitious material, which only differ by the drag coefficient set at 0.60 (instead of 0.47 for the first component). This C d value was the one used by Grafton et al. (2015). Considering that the two fertiliser components were in the same weight proportion, the simulation was carried out for a setting angle of À54°and a swath spacing of 39 m. In these spreading conditions, Fig. 13 shows the transverse distribution of the blended fertilisers and of each component. The CV computed for the blended fertilisers was 6.1%. Nevertheless, the CV values were at least twice when each fertiliser component was considered independently. Thus, the CV reached 12.5% for the first component and 17.1% for the second one. This illustrates that the measurement of the CV for blended fertilisers (i.e. measured without differencing the components) does not systematically provide a good assessment of the spatial chemical variability.
In contrast with works of Antille et al. (2015) or Grafton et al. (2015), the benefit of this approach is that the results are not limited to comparisons of the ballistic lengths of individual particles. The present approach considers the two-dimensional spread pattern so that the transverse distribution is computed taking into account the overlapping of adjacent passes.
This section put the emphasis on the interest of HCSM to investigate the impact of fertiliser properties on the transverse distribution. The effect of the spreading segregation for blended fertilisers has been illustrated by considering a fictitious material so that the two components of the mixed fertilisers only differ by their drag coefficient. For actual mixed fertilisers, each material would usually differ in numerous parameters. Thus, the accurate study of a specific blended fertiliser would require characterising the mechanical behaviour of the mixture on the spinning disc to provide the horizontal and vertical outlet angle distributions and the angular mass flow distribution (as described in Section 2.2). It would also require characterising the properties of each fertiliser components in terms of specific density, diameter distribution and drag coefficient (as described in section 2.3). Using these input parameters, the HCSM would be an efficient and low cost strategy to study the behaviour of blended fertilisers and provide recommendations on swath spacing.

Conclusion
To simulate fertiliser spread pattern depositions, a hybrid model was proposed. This approach combined the use of theoretical motion equations and experimental results. The experimental data were used to adjust few constant parameters and provide the statistical distributions of other input parameters. This ensured the realistic nature of fertiliser mechanical behaviours and spread pattern simulations. The particularity of the hybrid model was the use of successive random selections to compute the spread pattern deposition of virtual particles whose size and motion parameters respected experimental statistical distributions. This Monte Carlo process took into account the variability of input parameters and made possible the use of simulation replications to access to statistical characteristics of the output variables.
Simulation results showed the Hybrid Centrifugal Spreading Model was worth of interest to study information that was difficult or impossible to access with actual experiments. In particular, results demonstrated the transverse CV not only depended on the spreader setting and the swath spacing but also increased when the application rate decreased. The CV value also increased when the collection tray surface decreased. A mathematical relationship had been derived from simulation results to describe these influences. The study also demonstrated the variability of CV measurements increased when the application rate or the collection tray length decreased. Differences observed in the CV value, when it was measured in field or following the standard specified in EN 13739-2, were highlighted.
An insight into the distribution of the fertiliser particles related to their diameter or their drag coefficient showed the Hybrid Centrifugal Spreading Model will be a powerful tool to analyse further the impact of fertiliser ballistic properties on the spread pattern.
More generally, the model and the associated Monte Carlo simulations open up the possibility of carrying out virtual and numerical experiments to avoid cumbersome experimental tasks for numerous research or development activities. For instance, this will be of particular interest in: studying the effects of perturbing factors such as wind; providing recommendations concerning the use of blended fertilisers for a selected swath spacing; comparing different test methods to assess the transverse distribution or the spread pattern (especially to design simplified tests); assessing the accuracy of test methods (especially in defining the optimal swath width); or defining the probability of obtaining a selected range of application rate for a selected spatial scale defined by agronomical criteria.