Fluid Dynamic Approaches for Prediction of Spray Drift from Ground Pesticide Applications: A Review

Spray drifts have been studied by mathematical models and computer simulations as an essential complement to lab and field tests, among which are fluid dynamic approaches that help to understand the transport of spray droplets in a turbulent atmosphere and their potential impacts to the environment. From earlier fluid mechanical models to highly computational models, scientific advancement has led to a more realistic prediction of spray drift, but the current literature lacks reviews showing the trends and limitations of the existing approaches. This paper is to review the literature on fluid-mechanical-based modelling of spray drift resulting from ground spray applications. Consequently, it provides comprehensive understanding of the transition and development of fluid dynamic approaches and the future directions in this research field.


Introduction
The ultimate objective of any spraying system is to provide optimal deposition of spray materials on the targeted canopies to effectively control pests and diseases [1,2]. However, numerous reports have shown that a significant fraction of released chemicals drift to non-target areas during applications [3]. The amount of such losses has been estimated up to 50-60%, causing significant economic loss [3,4].
Pesticide spray drift is the airborne movement of spray droplets and particles to any site other than the area intended [5]. Comprehensive effects of the physicochemical properties of the spray formulation, sprayer design, crop characteristics, and weather conditions influence the complex phenomenon of spray drift. Depending on the spray applicator, the aerial application is conducted by aerial vehicles, such as aircraft, helicopters, and unmanned aerial vehicles, while the ground application is traditionally conducted by ground sprayers, such as ground vehicles, hand-held sprayers, and backpack sprayers. Many practical guidelines [6] indicate that pesticide applications directed upwards or released at a higher altitude are likely to cause more drifts. Opposite results are also found from measurements [7] where drift from the ground mist blower is greater than the aerial application due to smaller droplets and greater initial horizontal droplet velocity. Meanwhile, drift from the aerial applications is significantly influenced by the effect of wingtip vortices [8,9].
For most arable and vegetable crops, the air-jets tend to be located above the crop producing a downward flow, while an upward flow is applied to fruit and tree crops by locating the sprayer within the canopy [10]. The benefits of the downward sprayer are an improved control of the spray distribution within dense canopies and reduced spray drift over crops. In some applications, spray droplets injected from arrays of nozzles are carried within orientated air-jet streams. The strong air-jets carry droplets and reduce the flight time of droplets, which can reduce the interferences of wind and weather conditions and increase spray deposition within the crop canopy [11]. However, the use of air-jets can cause excessive environmental contamination by spray drift and ground deposition if the air streams are poorly matched to the intended crops [10]. The off-target deposition, as a result, will decrease the spray effectiveness and increase the cost for spraying materials.
While pesticide drift has raised world-wide concerns on dietary risks, human exposure, and environmental contamination to our ecosystems, many researchers have conducted laboratory and field experiments to assess and minimize drift losses from pesticide applications. However, such experiments are generally very expensive and time consuming. Major limitations of the experiments are uncontrollable meteorological conditions and challenges in measuring the multiple scales of the flow and airborne droplet mass and momentum (i.e., from the large scale 3D flow field all the way down to the small droplet boundary layer where evaporation takes place) [12].
In this regard, mathematical models and computer simulations based on fluid mechanics have been an important complement to field tests that help to understand physical processes taking place during the transport of spray droplets [8,13]. In the early stage, modelling spray droplet movements in the air were conducted using atmospheric dispersion models [14,15]. The most common model, the Gaussian plume model, which determines contaminant transport in medium or long-range distances (0.5-10 km), was applied to predict pesticide concentrations over a range of distances. Even though it has a very limited resolution near the contaminant source, it is effective in simulating the influence of meteorological factors, including atmospheric stability, on long-range drift. Another common model is the Lagrangian particle tracking model. It calculates the trajectory of droplets, each of which is exposed to several forces depending on its aerodynamic characteristics, velocity, and meteorological location. Due to its simplicity, it has mostly been used to assess downwind spray drifts from aerial applications [16].
Updates of those models have been made to predict spray drift from ground pesticide applications [15,17]. Two important airflow processes were considered in the models; the entrained airflow caused by momentum transfer between the spray droplets and the surrounding air, and the cross airflow due to the three-dimensional nature of the wind flow were considered as very important factors [8]. A random-walk model is one of the early stage attempts to adapt the effect of turbulent airflows based on the fluid mechanics. It calculates the velocity of each droplet based on its velocity at the previous time step and an additional random component, which represents the effect of turbulence. Detailed formulations for the entrained airflow were also suggested by several literature works [18,19].
One of the most recent fluid-mechanical approaches is a computational fluid dynamics (CFD) simulation, which solves complex airflow and turbulence patterns by solving the Navier-Stokes equations. Although the above analytical methods are available to describe some feature of turbulent air jets, they are mostly suitable for studying isolated free jets [11].
The key advantage of CFD is the potential of investigating large-scale three-dimensional flows involving the nature of turbulent airflow created by the sprayer fans, ambient wind conditions, plant canopy disturbances, and operational parameters [20]. However, the CFD simulation requires enormous computational resources and costs compared to other methods. It also needs careful consideration for a trade-off between the size of target areas, the desired accuracy of simulations, and the computational costs. Despite the complexity of the CFD approach, it has increasingly been used to investigate the transport of spray droplets within or beyond the application areas because it can increase the degree of realism by incorporating physical processes related to spray droplet dynamics.
In the past decades, many studies have contributed to the knowledge of the topic, but it is still challenging to apply all key factors influencing spray drift (see Section 2) realistically to various models due to the complexity of phenomena and limited technical aspects of the models. Now it is necessary to analyse the trends, strengths and limitations of existing fluid mechanical approaches. The objective of this paper is to review the past and recent research works that studied the spray droplet transport and its drift from ground applications using fluid mechanical approaches. This review does not include the spray drift process from aerial applications, which is quite different from that of ground applications. This paper also discusses significant recent studies and future directions toward the advanced numerical simulation and prediction.

Key Factors Influencing Spray Drift
Spray drift is dominated by many factors between which intricate relation makes drift estimation more challenging. Several papers were reviewed to summarize the factors that influence spray drift [8,13,[21][22][23][24][25]. The details of the influential factors are not necessary in this review, but this section shortly reviews the key factors from a modelling point of view.

Droplet Size
Spray droplet size has been considered as a primary factor affecting spray drift. Smaller droplets remain in the air longer and are carried farther away by winds. Larger droplets maintain their initial velocity longer than smaller ones and are more likely to be deposited on the intended targets. Droplets with a diameter of <100 µm were considered highly driftable by most studies, but a diameter of <50 µm, <150 µm, and <200 µm was also suggested as drift thresholds [21,23].
The size of droplets discharged from nozzles depends on the atomizer type, nozzle pressure, liquid flow rate, liquid properties, and atmospheric temperature and relative humidity [1]. While various nozzle types produce a wide range of droplet sizes, the effect of other factors is quite clear, as smaller droplets are generated by higher pressures and low liquid flow rate. Higher flow rates will increase the droplet size, and too high liquid flow rate may cause a poor atomization [26]. High air temperature and low relative humidity also decrease the droplet size during transportation and result in more drift because of droplet evaporation.
The droplet diameter can be decreased by the breakup of droplets. High discharge velocity of droplets or air stream by fans can bring droplets to the intended target and increase spray penetration into the plant canopy. However, whether or not air-blown fans are used, excessive relative velocity between the droplet and the surrounding air can result in breaking up of droplets [27].
Additives in the spray solution modify the physicochemical properties of the spray liquid, such as surface tension, viscosity, and evaporation rate, thereby influencing the droplet size [21] and consequently its drift [28].

Meteorological Conditions
While the droplet size is related to susceptibility to drift, wind speed and direction and atmospheric stability are meteorological factors that influence spray drift and deposition [29]. When droplets injected from the nozzles lose their momentum in the air, they are susceptible to be blown by the wind. As wind velocity increases, the drift distance increases and more spraying droplets are lost from the target area. When droplets are drifting, the wind profile is also of importance. In this regard, nozzle height has a great influence on drift because droplets sprayed at higher positions will be blown by higher wind speeds and will have a longer time for transport and evaporation before depositing on the ground [30]. Although atmospheric stability was not commonly measured in the field, it has been considered an important factor for spray drifts [22]. It was found that the spray droplets drifted along the mean wind direction in stable atmospheric conditions while the sprayed liquid plume was more extensive with more materials dispersed to higher levels in unstable conditions during air-assisted spray application [31]. Therefore, [32] did some early measurements, observing decreases in downwind spray deposition with the change of atmospheric stability from very stable to unstable. Most drift models did not consider atmospheric stability effect on spray drifts and sometimes underestimated tiny spray droplets and cloud displacement [13]. Especially, models based on accurate airflow calculation like CFD have limitations in modeling the atmospheric stability notwithstanding some attempts [33]. Several models consider the effect of atmospheric stability by adapting diffusion parameters of the Gaussian diffusion process that depend on atmospheric stability category [8]. As concluded from previous findings [34], the effects of the atmospheric stability are more significant at greater downwind distances. While wind speed dominates deposition in the near field where larger droplets are deposited by gravitational forces, the stability is more significant in the far field where smaller droplets deposit.

Plant Canopy
It is generally assumed that spray drift has an inverse relationship with canopy density. Chen et al. [35,36] measured spray retentions and off-target depositions, respectively, for tree canopies of three growth stages, which were leafing, half-foliage, and full-foliage stages. The increase of canopy density apparently decreased the amount of off-target drift. However, it was also reported that dense canopies of grapefruit and orange trees, on the contrary, increased spray drift because the dense foliage deflected the airflow over the top of trees [37]. This meant the spray retention resulted from the complex interaction between the canopy density and air penetration through the canopy [38,39].
Early numerical models have focused on spray drift and assumed that if a plant intercepted a droplet, it was always retained [16,20]. Schou et al. [40] suggested that this assumption might overestimate total spray retention because droplets that might fail to be retained on the leaf could rebound and continue their journey through the canopy. The droplet bounce and shatter have been adopted in recent models [40,41]. It was found that the amount of spray retained on a leaf surface was influenced by the size and velocity of incoming droplets, spray formulations, leaf surface characteristics and properties of any shatter or bounce droplets after impact [41]. Increasing droplet size and velocity decreased the spray retention within the canopy. The spray droplet retention by plant canopies is discussed in detail in Section 4.4, Spray droplet retention by plant canopies.

Early Models for Spray Drift Prediction
Numerical models of spray drifts have been developed and used mainly to predict fate and transport of spray droplets under specific conditions, which are very difficult and expensive to measure during field tests; aid understanding of the drift phenomena; and complement experiments with a limited number of point measurements. This chapter provides an overview of the early model approaches.

Plume Models and Droplet Trajectory Models
The Gaussian plume theory is the most common model applied to atmospheric pollutant dispersion. The primary factors that determine the extent of dispersion are the lateral and vertical standard deviation of the cloud, which are a function of a downwind distance. The resulting concentration decreases exponentially with distance from the source, and this is roughly reasonable in most steady-state dispersion phenomena. The simple plume models predict dispersion using a single mean wind speed, which gives a bad approximation near ground, but has the advantage of computational simplicity [42].
RTDrift model based on the Gaussian tilting plume model was developed to describe spray drifts from a boom sprayer [43]. It calculated the spray drift deposits up to 30 m from the boom edge by integrating the contribution of individual nozzles concerning time and droplet size classes. However, it has a limitation to consider the wind turbulence effect. Changes in wind direction during the experiment resulted in curved wind flows, which were not consistent with the assumption of the Gaussian plume theory, and caused large fluctuations in drift deposits.
Lagrangian models are commonly used for spray particles that move along trajectories in the atmosphere determined by the wind field, the buoyancy and the turbulence effects [8,13]. After calculating the trajectories of each particle, the final distribution of all particles gives a stochastic estimation of the spray concentration.
Puff models are an intersection or a hybrid between Gaussian and Lagrangian models assuming that the centre of plumes is moving along a transient Lagrangian trajectory, but the concentration pattern from the centre of each plume is calculated by a Gaussian distribution. Puff models are now most widely used in long-range atmospheric dispersion processes and include the Danish models, DERMA and RIMPUFF, in Europe and the CALPUFF in the US [44].

Empirical Models
Field measurements or wind tunnel experiments produced meaningful findings and sometimes established statistical relationships between spray drifts and other variables as shown in Table 1. In most regression studies, spray drifts were estimated as a function of the downwind distance from the sprayer or the downwind edge of the application area. The shape of decay curve was expressed as power, log, or exponential functions. Wind speed and air temperature as meteorological variables were also considered as essential variables in some studies [45][46][47][48]. For boom sprayers, nozzle height, nozzle pressure, or nozzle flowrate were introduced into the regression model [47,48].
For grazing and cereal crop fields, spray drift from a wide boom sprayer was investigated by measuring airborne, and fall-out drifts separately and modelled empirically [45]. The experiments found that the most decisive factors influencing the total spray drift were the wind speed and the boom height. However, when the wind speed and the boom height were consistent, the spray droplet sizes (especially the volume fraction of droplets ≤ 102 µm) was a significant variable to determine airborne or fall-out drifts.
While most regression models showed relatively high R2 values, they still had several limitations: the models were valid within limited ranges for each variable, and the drift distances were rarely measured over 30 m because of difficulties in conducting field experiments [49]. Another disadvantage of this approach was that it required a quantity of data for the large number of variables known to influence drift [50]. However, the most significant advantage of this approach is still that it is much less computationally expensive compared to other analytical models [51]. Table 1. Empirical models for spray drift prediction.

References Regression Models Conditions
Ganzelmeier and Rautmann (2000) [52] y is the spray drifts at a downwind distance of x or corrected distance x * ; S is the shape factor; ∅ is the optical porosity of plants; V or V k is the wind velocity at a height of k m (m s −1 ); f 102 is the volume fraction of droplets ≤ 102 µm (%); T is the air temperature ( • C); d vp is the vapor pressure deficit (mbar); A H is the absolute humidity (g kg −1 ); H n is the nozzle height (m); V id is the initial downward droplet speed (m s −1 ); p n is the nozzle pressure (kPa); F n is the nozzle flowrate (mL s −1 ); a is the max value (the intercept); b is the exponent; and c or c i is the constant.

Fluid Dynamic Considerations for Modelling Ground Pesticide Applications
According to the early reviews [8,13], the most commonly used models can be classified as plume models and droplet trajectory models, which were well summarized in Section 3.1, Plume models and droplet trajectory models. This paper complements models derived from fluid dynamic approaches including recent CFD simulations.

Air Jets and Entrained Air Currents
Most boom sprayers do not generate air jets, but spray droplets discharged from nozzles entrain the surrounding air and create air currents. Several studies emphasized that the entrained air affects the initial trajectories of small droplets and included the effect of the entrained air in drift models [50,54]. After droplets are discharged into the air, there is friction between the droplets and the surrounding air, resulting in momentum transfer between them, which causes an entrained airflow. The earliest model [55] described the velocity of the entrained air for flat fan nozzles. The entrained air velocity decreases with the distance from the nozzle outlet following a power law. The model was further developed into two dimensional jets [50] and three dimensional jets [54] by assuming that the entrained air velocity was maximum at the centre of the spray jet and decreased to zero at the lateral ends of the jet or their vicinity. The entrained air velocity for boom sprayers becomes a dominant air velocity acting on spay droplets in short distances, and fades away quickly out of the distances where the spray droplets are subject to ambient winds. Additionally, multiple nozzles may create some interaction between airflows around neighbouring nozzles [50].
Air-assisted sprayers generate strong air jets, and the airflow pattern and the air velocity distribution induced by the air jets are the most critical determining factor of carrying pesticide droplets to canopies [56]. The air jet model was initially developed from the turbulent jet theory [57] and then revised in various forms by further works [19,56,[58][59][60][61]. The jet model describes the geometric properties of an ideal turbulent jet by dividing the airflow into three regions (initial, transitional, and main regions) and applying the conservation of momentum equations. The fully developed flow of air jets in a stagnant ambient air exhibits a decay of the centreline air velocity with distance from the outlet. The turbulent jet theory gave reasonable predictions of the air jet velocities by distances in most cases, but some studies [58,62] indicated that the model prediction was not accurate at far distances or where the air velocity was quite low. Another study also showed that calculations with turbulent jet theory overestimated the measured velocities [60]. It may be because the jet models considered ideal air jets with no obstacles and ambient winds. Therefore, some studies considered more complex effects of canopies, the motion of sprayers, and ambient winds on turbulent air jets and reproduced the air jet velocity profiles at the nozzle outlets through measurements [56,63].

Droplet Release
Droplets have their maximum velocity immediately after being released from nozzles. The initial injection velocity of droplets can be simply estimated by Bernoulli's equation [64]. Most studies assumed that all droplets are discharged at a constant velocity, but an analytical model for variable droplet velocities was also proposed [65]. This model assumed that droplets decelerate during separation from a liquid sheet based on energy balance. In the model calculations as well as measurements, velocities of small droplets (≤70 µm) consistently decreased with the decrease in droplet size while the velocity of large droplets (>70 µm) remained fairly constant.
When droplets released from nozzles encounter an ambient airflow field, the relative velocity between the droplet velocity and the nearby air velocity cause an aerodynamic force on the droplet. The aerodynamic force deforms the shape of the droplet and, if it is excessive, breaks the droplet structure, which is called secondary atomization [66]. The critical relative velocity at which water droplets start to break up was derived from Weber numbers in many studies [40,64,67,68]. If the Weber number of the droplet is greater than the critical Weber number, then the droplet will break up. According to simple calculations [27,69], droplets smaller than 100 µm start to break up at over a relative velocity of 78.4 m s −1 while droplets within a range of 200 to 500 µm will break up at relative velocities of 35-55 m s −1 . This means larger droplets are more susceptible to breakup due to strong air jets of air-assisted sprayers than smaller droplets.

Turbulent Air Flows in a Plant Canopy and a Crop Field
In log-or power-law wind profiles, air velocity within a vegetative canopy is assumed to be zero, as the vegetative canopies are considered as a concept of surface roughness and zero plane displacement. This is practically feasible for short vegetation, but not sufficient for tall trees [70] because the spray drift mostly occurs within and above tree height. Observations made for tree canopies found that wind speeds within tree canopies were considerably lower than those over trees, but still effective enough to carry spray droplets and complicated due to high turbulence [70][71][72].
The effect of the plant canopy on air flows was considered as the loss of momentum of the air flows per unit ground area of the canopy in a manner similar to the definition of the drag coefficient [73]. For modelling in more than two dimensions, the ground area of the canopy was refined into the leaf area density [74], which was defined as total leaf area divided by the canopy volume. In recent studies [75][76][77], a plant canopy was assumed as a porous media, and its momentum loss was expressed using the Darcy-Forchheimer equation, which was composed of two parts, a viscous loss term and an inertial loss term (or called an aerodynamic pressure loss term). The viscous loss term is often smaller or negligible compared with the inertial loss term [39]. The leaf area density and pressure loss coefficient of the vegetation can be heterogeneous due to irregular distribution of leaves and branches, but most studies simplified the model with constant values for convenience. A few studies used vertical distribution of leaf area density that varied by tree height [70,78,79].
The above concepts should be coupled with air flow calculation because the velocity of wind field and its momentum loss by vegetation are interactive, and thus require an iterative computation to obtain a converged wind flow field within canopy height. However, there are other models that directly estimated the profile of wind speed within a canopy zone, such as an exponential wind profile [80,81] and a squared reciprocal wind profile [54]. Similarly, the profile of wind speed above vegetation could be estimated as a logarithmic profile assuming that the wind speed increased logarithmically with height above the top of the canopy [81,82] and more elaborated second-order profiles [83][84][85].
Regarding the turbulence due to a plant canopy, earlier study conducted in a wheat field reported that the turbulence intensity within crop canopies was constant [81]. Since then, however, most studies revealed that the turbulence intensity decreased gradually downward from the top of canopies in proportion to the wind speed [82]. Especially for tall crops and trees, vertical profiles of the turbulence quantities were influenced by the vertical distribution of canopy foliage, stems, and trunks [86]. For deciduous trees, the maximum turbulence intensity was observed at 0.8H t (H t is the tree height) where over 70% of the foliage was concentrated. The measurement also showed that the turbulence intensity decreased slightly at night compared to that during the day due to a stable atmospheric condition. Tree spacing also influenced the distribution of turbulence intensity. The measurement in Sitka spruce forest showed that the highest levels of wind speed and velocity fluctuation were observed between tree rows rather than within the tree rows because of wakes formed behind individual trees [87]. The turbulence intensity decreased as the tree spacing increased.
While the canopy turbulence has been studied based on extensive field and laboratory measurements, recent computational studies have directly solved relevant governing equations for spatial and temporal distribution of turbulence quantities, especially turbulent kinetic energy and its dissipation rate, caused by individual crop and, sometimes, individual leaves. When the wind penetrates through crop canopies, it loses part of its momentum by vegetative drag. Assuming that all energy extracted from the mean flow by vegetative drag changes to turbulence, the turbulence production associated with canopy wakes is calculated as [88]: where S k,c is the source of the turbulent kinetic energy; k, by canopy (m 2 s −3 ); L ad is the leaf area density of canopy (m −1 ); U is the mean wind speed (m s −1 ); and 1/2 can be omitted when the drag coefficient, C D , is expressed in order of 1/2. Tree canopy models in recent reports [79,[89][90][91][92] proposed that a sink term of the turbulence should be added to Equation (1) because the turbulence decreased with energy cascade from large to small scales due to the presence of plant foliage, that is called the spectral short-circuiting of the energy cascade. Therefore, the budget equation including the production and sink terms is currently used for the expression of turbulence due to a plant canopy. Most of all, measurements of wind velocity components are needed to accurately investigate vertical profiles of turbulence and wind speed because irregular shape and distribution of canopy structure make discrepancies between model assumption and field conditions especially for tall plants.

Spray Droplet Retention by Plant Canopies
While spray droplets travel in crop canopies by advection and diffusion, some droplets impact on a leaf surface. When a droplet hits a surface, it ends up in three possible states: adhesion, bounce, and shatter [93]. The bounced droplet or shattered tiny droplets can adhere on a surface through secondary impaction, and the final result of all adhered droplets is considered as the retention. The outcome of impaction depends on physical characteristics (surface tension and viscosity) and chemical formulation of the spray solution, aerodynamic characteristics (size and velocity) of the airborne droplets, and leaf surface structures (roughness, hydrophilicity, and angle to droplet impaction) [94].
Earlier studies mostly did not include the bounce and shatter, but considered the impaction mechanism as 'interception' [95]. The probabilistic analysis derived an empirical equation for the efficiency of droplet capture by inertial impaction, which was dominated by the droplet's Stokes number [95,96]. According to the equation, larger droplets have greater Stokes numbers and higher collection efficiency. Besides, small leaves like needle vegetation would be better at capturing spray droplets than large leaves. The overall efficiency of droplet capture by the whole plant was then calculated considering the Agronomy 2021, 11, 1182 9 of 20 optical porosity of the plant and the probability of droplets encountering any vegetation element [96,97].
The effect of droplet rebound was considered in an early study by simple assumption that the droplet retention efficiency was a product of the average collection efficiency by the canopy and the reduction ratio in collection due to particle rebound [98]. Recent studies have investigated physical processes including droplets' adhesion, bounce, and shatter, and their associated physiological parameters [41,93]. Key physical parameters that influenced the dynamics of a droplet impaction were the dynamic viscosity, surface tension and density of the droplet fluid, the droplet diameter, and the droplet approaching velocity. These parameters were used to characterize the droplet dynamics in forms of dimensionless numbers, such as the Weber number and Ohnesorge number. The energy equation was also used to describe the advancing and receding motion of a droplet upon impaction on the leaf surface [99]. The energy equation predicted the change in the size of a deformed droplet taking into account the changes in the kinetic energy and potential energy of the droplet and energy loss due to viscous dissipation. Using the energy equation, the ranges of the Weber number and the Ohnesorge number where a droplet became bounce or shatter were determined [93,100,101].

Wind Field and Downwind Spray Drift
Most studies, that used particle tracking, plume dispersion or other approaches to model the spray drift, require a description of appropriate wind fields on the downwind. Earlier studies that used the Gaussian plume model assumed a constant wind speed throughout the whole computational domain [43,102] while they considered the meteorological effects by adjusting the plume dispersion. In case of the Gaussian puff model, puffs are advected according to transient wind fields while being dispersed similarly to the plume of the Gaussian plume model. The transient wind field is usually obtained from the Navier-Stokes equations [24].
The horizontal change in wind speed and direction is considered only in studies that used wind flow models, such as the diagnostic wind model and Navier-Stokes equations, to calculate the wind field. The Navier-Stokes equations apply Newton's law to fluid motion in spatially discretized domains and solve the distribution of fluid characteristics, especially fluid mass and velocity, in a steady-state or unsteady-state regime. Many studies that used computational fluid dynamics simulations solved the Navier-Stokes equations to describe complex wind field formed by entrained air movement by spray injection and drag effect of plant structures [63,103,104]. However, most CFD studies have described three-dimensional, or sometimes two-dimensional, airflow patterns in the near-field of sprayers, and little simulated the wind field in the far-field. The most prominent reasons for this were a heavy burden on the computational resource required for a large-scale modelling, not even for spray drift but for all airborne droplet dispersion, and difficulty in modelling the thermal effect of different stability classes numerically [33]. However, some CFD studies [33,105] began to include stability effects by changing vertical profiles for thermal and turbulence variables depending on the atmospheric stability classes using the Monin-Obukhov universal function [106].
Wind speed profiles in most downwind regions conform to common logarithmic profiles [29,54,107] and are easy to evaluate, whereas wind speeds in the near-field of plants are much decreased by the shelter effect ( Figure 1). As the wind approaches plants, the plants adsorb momentum from the wind flow and result in a quite zone with considerably reduced wind speed and a mixing zone with strong wind shear and turbulence within a certain distance downwind [108]. The flow in the quiet zone can reverse direction, creating recirculation flows (or eddy flow) if the plant canopy is dense. An experimental study found that spray drift was mostly deposited in the quiet zone [109]. The strong wind shear in the mixing zone moves the air momentum downward and eventually re-establishes the logarithmic wind profile in the far-field. The wind profile of horizontal and vertical wind speed in the quiet and mixing zones is different from the typical wind profile, and the spray drift is affected to a certain extent by the sheltered areas in the way of strong mixing in the air and drift deposition on the ground [109,110]. The length of the quiet zone and mixing zone, known as recovery distance, depends on canopy density, plant height, and characteristics of the approach wind [108,109,[111][112][113][114]. Many studies did not suggest the recovery distance clearly, but provided a wide range from 9H p to 35H p (H p is the plant height). strong wind shear in the mixing zone moves the air momentum downward and eventually re-establishes the logarithmic wind profile in the far-field. The wind profile of horizontal and vertical wind speed in the quiet and mixing zones is different from the typical wind profile, and the spray drift is affected to a certain extent by the sheltered areas in the way of strong mixing in the air and drift deposition on the ground [109,110]. The length of the quiet zone and mixing zone, known as recovery distance, depends on canopy density, plant height, and characteristics of the approach wind [108,109,[111][112][113][114]. Many studies did not suggest the recovery distance clearly, but provided a wide range from 9Hp to 35Hp (Hp is the plant height).
Another important region is the roughness sublayer, which is associated with the inflected velocity profile in Figure 1C. This layer, being located around plant height up to 3Hp, dominates the transport of droplets above the plant into the canopy. Since the roughness sublayer is difficult to model by numerical and theoretical models, investigation based on a measurement will be important.

Computational Fluid Dynamics Approaches
The airflow patterns from sprayers are inherently three-dimensional (3D) according to spraying cloud, wind, and plant structure. Full-scale experiments have been carried out to understand the complex airflow patterns and resulting spray drift phenomena, but they are usually expensive and difficult to perform because of variations in meteorological conditions and plant structures [63]. Recent studies use CFD simulations to produce numerical solutions for complex 3D airflow patterns by solving physical conservation equations for mass, energy, and momentum [63,103]. It is still challenging to get a realistic solution for the spray drift because of difficulties in modelling complex atmospheric phenomena and plant structures. However, CFD can be a practical method to compare various sprayer Another important region is the roughness sublayer, which is associated with the inflected velocity profile in Figure 1C. This layer, being located around plant height up to 3H p , dominates the transport of droplets above the plant into the canopy. Since the roughness sublayer is difficult to model by numerical and theoretical models, investigation based on a measurement will be important.

Computational Fluid Dynamics Approaches
The airflow patterns from sprayers are inherently three-dimensional (3D) according to spraying cloud, wind, and plant structure. Full-scale experiments have been carried out to understand the complex airflow patterns and resulting spray drift phenomena, but they are usually expensive and difficult to perform because of variations in meteorological conditions and plant structures [63]. Recent studies use CFD simulations to produce numerical solutions for complex 3D airflow patterns by solving physical conservation equations for mass, energy, and momentum [63,103]. It is still challenging to get a realistic solution for the spray drift because of difficulties in modelling complex atmospheric phenomena and plant structures. However, CFD can be a practical method to compare various sprayer designs and operating conditions in terms of spraying efficiency and predict the drift of spray droplets under ideal meteorological conditions [20,25].

Modelling Considerations in CFD Prediction of Spray Drift
The most common approach to modelling spray drift is the Eulerian-Lagrangian method [9,20,115], which solves the problem in two steps. The first step is to solve the Navier-Stokes equations for the continuous phase, which is usually the air. In this step, airflow patterns are determined. The most common way to solve the Navier-Stokes equations is to transform the equations into the Reynolds Averaged Navier-Stokes (RANS) equations and additional turbulence models. The next step is to solve the Lagrangian particle tracking equation for the discrete phase, which is spraying materials. The discrete droplets released into the continuous phase are moved to new positions every time step according to mainly the drag and buoyancy forces. The continuous and discrete phase can be interactive by exchanging mass, energy, and momentum. The most common form of the Lagrangian particle motion equation is [11]: where u p is the particle velocity (m s −1 ); f D is the drag force factor per unit particle mass (s −1 ); u ∞ is the velocity of the continuous phase (m s −1 ); g is the gravitational acceleration (m s −2 ); ρ p is the particle density (kg m −3 ); ρ ∞ is the density of the continuous phase (kg m −3 ); and t is the time (s). The term on the left-hand side indicates the change of the particle velocity by time while the first and second terms on the right-hand side calculate the drag and buoyancy forces acting on each particle, respectively. Mass, momentum, and energy exchanges between the discrete and the continuous phases, mostly droplets and the air, are well described in the literature [2]. Droplets lose their masses from evaporation caused by these continuous exchanges in forms of convective and latent heat transfers. The simultaneous changes in droplet sizes because of evaporation in transport are included in computer simulations [12,104]. The initial size of droplets discharged from the sprayer can be determined by the atomization model or direct measurement. Most studies [2,56,116] used the atomization model developed by the liquid sheet atomization model [117]. The atomization model predicts the droplet size spectrum at the exit of the nozzle using the appropriate input parameters, such as spray angle, liquid flow rate and pressure, nozzle size, sheet constant, and ligament constants [2]. Some studies determined the droplet size spectrum by measuring droplets from the nozzle using size analysers [29,118,119]. In most studies, regardless of whether the droplet size was obtained by measurement or prediction, the obtained size distribution was compared with or fitted to a Rosin-Rammler size distribution.
Plants were not considered in early CFD works, but the latest studies included the effect of trees in their CFD models [20,63,116]. The geometric model of full-scale pear trees, especially all the branches, was constructed using the measurement and modelled as a 3D object. Porous sub-domains were added around the branches to model the effect of very thin branches, flowers and leaves by adding appropriate source terms in momentum and turbulence equations. Besides these works, plants were modelled as a porous media with cuboidal or spherical geometries [97,112]. Recent studies modelled a tree as a porous media with various porosities, i.e., branches and trunks as extremely low porosity and leaves as higher porosity [103,104,120]. This method set various porosities at the structured cells where imaginary trees are assumed to be located and enables easy modelling of trees without geometrically modelling complex 3D objects of trees.
In all studies reviewed in this paper, the geometrical feature of sprayers and tractors was not included in the CFD model due to its complicated shape, but simplified as a square inlet or simple geometries. In the case that the tractor or sprayer was moving during spraying, two methods were used to replace actual movement of the tractor or sprayer: a pulse function [56,103,104,116] and a moving coordinate system [2,29].
While the ground applications are targeted at orchards or crop fields, most CFD simulations have limited the computational domain to a part of the target areas, such as the last row of trees on the leeward side or the leeward edge of the field, in order to reduce the computational cost. Here, some recent studies emphasized that the wind profile for the inlet boundaries should be different from typical atmospheric wind profiles because the wind at the inlet of the computational domain already passed through the whole orchard or crop fields changing its vertical velocity distribution due to drag forces by plant canopies [20,63,104,116,121]. That is called the canopy wind profile. It was obtained by a series of steady-state simulations until the simulated wind profile matched the measurement [20,63,116,121] or additional simulation of sub-models containing windward part of the target orchard [104].

Experimental Validation
Validation of the CFD model was usually conducted by measuring distributions of air velocity and/or drift mass concentration. For air-assisted sprayers, the distribution of air velocities created by centrifugal fans or air ducts is a key factor determining the efficiency of spray application.
The homogeneity of the vertical distribution was examined by field measurements and, at the same time, used for the validation of the air flow simulations. Vertical velocity profile was assessed usually near the sprayer outlets and at the location of the spray target, while plants were not considered [60,122,123]. Some studies examined the velocity profile at the rear of the target plant [116] or inside the plant [103] in order to validate the effect of drag by the plant canopy on the air flow. Other studies also compared the distribution of air velocities created by jet nozzle at multiple positions at further distances [11,124].
For the validation of drift mass concentration, field experiments were conducted mostly using tracer materials, such as fluorescent dyes [29,104,125] and metal solution tracers [2,20,121], as a substitute of the pesticide due to the risk of pesticide toxicity. The spray drift was expressed as percentage of the applied mass [29,121], a normalized form based on the peak value [2,20], or deposit mass per unit area [104,125]. Some studies compared the vertical liquid distribution near the plant location between the model simulations and the measurements [2,20,104], while other studies evaluated the spray drift up to far distances, such as 20 m [29,125] and 40 m [121] from the target.
Controlled experiments using wind tunnels were also used for validation. The benefit of the wind tunnels is the controlled wind speed. However, the spray drift can only be explored within a limited distance, which is the length of the wind tunnel. Scaled test for spray drift has rarely been conducted because the similitude requirements were challenging. An early study conducted the wind tunnel test to validate a computer program by measuring spray drift distance within the tunnel length with respect to wind speed and droplet size [14]. Wind tunnels were also used for validating an airflow pattern behind the tree, which became decelerated and more turbulent due to the effect of the tree canopy [126].
Since drift is a complex phenomenon involving various scaled physical processes, high levels of understanding obtained by field and laboratory observations [127] can provide empirical and semi-empirical models and real data for the CFD validation and calibration.

Promising Uses of CFD Applications
Since the first time of using CFD simulations in pesticide spray applications [14], CFD has been actively used to predict the spray drift for ground applications and provided the most advanced solutions on the movement of spray droplets in the air. Especially, many studies investigated the trajectory of spray droplets injected from the nozzles until they reached the target plants or the distribution of droplets in the near field around the sprayer [122,123,128]. Since the nuisances arisen from spray drift are mostly far-field issues, the prediction of CFD simulations needs to be extended up to some distances, at least tens of meters, on the downwind side of the target area. Due to computational burdens, only a few studies considered such a large-sized computational domain [104,120,121]. Hong et al. [104] modelled a part of the apple orchard as 6 × 5 trees (30 trees) and predicted the spray deposition onto the target trees and adjacent trees, off-target deposition on the ground, and airborne drift beyond the orchard. Duga et al. [121] simulated the ground and airborne drift in a large computational domain with dimensions of 40 m in length and 50 m in width. Baetens et al. [29] predicted spray drift from a field boom sprayer at a distance up to 100 m. Hong et al. [120] also modelled a large computational domain with a downwind distance of 200 m long to predict the ground and airborne drift.
The benefit of CFD simulations with such a large domain is the production of the drift curve [120,121], which expresses the mass of spray drift with respect to the distance from the edge of the crop field. The drift curve is helpful to optimise the pesticide application or to prevent the risk to human health and neighbouring environment against spray drift.
However, it is difficult to obtain the drift curve by experiments because it varies with the type of sprayers and their setting, characteristics of target plants, driving direction and speed of the sprayers, and weather conditions at the time of the application [48]. In this regard, CFD can generate a huge dataset of drift curves under various conditions. Back to the 1990s, Zhu et al. [17] developed a computer program DRIFTSIM to predict spray drift potentials up to 200 m based on a large database established from a CFD program (FLUENT ® ). This user-friendly program has been widely used for boom sprayers to select critical components that minimize spray drift potentials. Hong et al. [120] also developed a CFD-based program SAAS similar to DRIFTSIM to predict spray drift and potential setback distances for vineyards and orchards with respect to crop conditions, spray operating conditions, and weather conditions. However, since CFD simulations require a huge amount of computational time and costs, the use of CFD for far-field drift prediction needs to consider a trade-off between the size of domain, grid resolution, and the desired accuracy of the result [120].
The fluid-plant structure interaction is also a significant feature of spray application. As most existing studies considered, the plant affected the turbulence in and around plant canopies (see the Section 4.3, Turbulent air flows in a plant canopy and a crop field) while the air flow influenced the deposition and retention of airborne droplets within the plant canopies (see the Section 4.4, Spray droplet retention by plant canopies). However, especially for air-assisted sprayers, strong air jets may shake tree branches or and deform plant canopies. The fluctuation of leaves will affect the droplet retention and rebound, resulting in variation in spray drift potential.
The dynamic actions, such as wave, flutter and vibrations, of plants due to turbulent air flows were reviewed, and the relevant mechanisms were described well in a recent review work [129]. An earlier study developed a two-way coupled CFD model that solved the turbulent air flow around the plant canopy using the flow solver and calculated the forces acting on the plant canopy and displacement variables due to the forces using the structural solver [130]. The plant vibration and the relevant vortex creation were well simulated. The impact behaviour of spray droplets on the leaf surface, even though the leaves were not moving, was studied using CFD simulations [131,132]. The impact behaviours, such as adherence, rebound, splash and shattering, were affected by droplet impact velocity, droplet diameter, formulation, surface topology and wettability, pesticide formulation, and distribution of pressure and velocity around the leaf surface. The effect of plant motion on the droplet retention in plant canopies or drift possibility was not investigated until now, but CFD approach is a promising tool to consider all related mechanisms.

Discussion and Future Research Trends
Many studies have made efforts to achieve a broader understanding of spray drift potentials based on the knowledge of fluid mechanics. From 1960s, physical and semiempirical models described turbulent air jets, airflows entrained by nozzles, canopy flows, and droplet atomization [55,71,80,81], which were not made for spray drift study. While the risk of spray drift was recognized in the 1960s, especially by measurements for aerial spray application, the fluid mechanical simulations came to the solution to predict spray drift in the 1980s [107]. Early models suggested the Lagrangian theory of droplet trajectory, which was also a key feature of current studies, but had limitations in modelling air flows around the target areas. With the technological advance, CFD became a promising tool since the 1990s due to its potential to investigate turbulent air flows and their interactions with airborne droplets [11].
Once the CFD achieved an accurate prediction of multi-scale flow regimes, the use of the former theories, such as turbulent air jet, entrained air current, and canopy airflows, had faded. Most of the flow-related phenomena were resolved by the CFD simulation with geometrical modelling of air jet ducts [122] and plant architecture [104], instead of the physical or semi-empirical theories. However, the CFD approach still used the former theories for some significant features, such as the droplet release and the droplet retention by plant canopy. Such phenomena not solved by CFD remain a formidable theoretical task.
Droplets initially released from nozzles have various sizes, but the range of the droplet sizes cannot be determined by CFD simulations due to limited computational resources. Some engineering fields modelled the droplet atomization process using CFD simulations [133], but they dealt with a very small computational domain with high-resolution meshes, which were not applicable to spray drift problems with a large domain. For the same reason, droplet retention by plant canopy was considered in the CFD simulations as probabilistic models [95][96][97] instead of modelling the actual impaction of droplets on leaves and branches. To the best of our knowledge, no study has modelled realistic branches and leaves for the spray drift issue. A simplified structure with appropriate porosity and spray retention models provided a good description of plants as a spray target [20,104,115].
Based on this literature review to date, a better solution requires understanding of the fate of pesticide droplets from their release to deposition. Pesticides, after being released from nozzles, adhere onto leaves of a target plant or other plants. The pesticides not captured by plants are suspended in the air, blown by winds, and deposited on the ground in the near vicinity of plants or very far-fields. There have rarely been models that best describe the downwind deposition of pesticide spray drift over all distances. Some models showed the best fit with field data closer to the sprayer while they overestimate or underestimate the spray deposits farther from the sprayer, and vice versa. For main factors that most influence the spray drift, some focused on characteristics of nozzles and droplet dynamics, while others studied the effect of meteorological variables, such as wind speed, air temperature, and humidity. Meanwhile, some of the variables are only important in limited stages of droplet transport processes. For example, meteorological conditions are not significant at locations very close to the sprayer nozzles because the initial velocity of droplets is already very high and in most cases overwhelms the ambient wind speed. It is more noticeable especially when an air blower assists the spray application. Air jets from the nozzle dominate the transport, collision, coalescence, and breakup of droplets. On the other hand, the effect of sprayer and nozzles recedes at very far from the sprayer, and droplets are transported by ambient wind and at the same time changing their size due to evaporation.
It is obvious that the transport of spray droplets can be divided into several processes, each of which is better described by different mechanistic models. However, most importantly, all processes evolve into a continuum model in the end. A process-based model is defined as the mathematical representation of one or several processes characterizing the functioning of delimited systems of interest [134]. In a process-based model, each process has its best fit models, and outputs of one process can serve as input to the subsequent processes. Models for each process should be developed in a mechanistic way because empirical analyses and curve fits restrict the use of data constraining the ability to accurately improve drift prediction [8]. More studies are highly recommended to integrate the existing knowledge into advanced technologies to predict the drift of pesticides at each destination. Such studies are of importance to optimize the efficiency of spray application as well as mitigate the impact of off-target spray drift to the environment.
Some future research recommendations based on our understanding to date are summarized as follows: • A process-based modelling, as one of the solutions to complex phenomena, which can offer more comprehensive understanding and easier interpretation of the spray drift than detailed numerical models. • Developing hybrid models using more than two approaches that can complement each other. For example, plume dispersion models and CFD models are appropriate for large-scale and small-scale simulations, respectively [8,50], and, thus, can complement each other.
• Improving spray retention models in relation to realistic features of leaf-airflow interactions. Leaves flutter or vibrate in wind and turbulence, resulting in spray retention and impact behaviours [130]. • Promoting validation of spray drift models along with enhanced measuring technologies. Most spray drift models were validated by data measured in limited environmental and operational conditions due to difficulties of the field experiment. • Using new computational technologies to import three-dimensional field images consisting of crops and terrains as boundary conditions and real-time local weather data as the initial inputs into the simulation models. • Including sprayer physical parameters and spray cloud patterns along with sprayer travel conditions in the computer simulations to demonstrate real-time droplet trajectories.

•
Assessing environmental risks and mitigation measures using the fluid-mechanical computer simulations for users and regulatory authorities [135].

Conclusions
Pesticide spray drift has been a worldwide concern on potential environment pollution and ecosystem damage. The scientific community has been expanding efforts to minimize the emission and drift of pesticides for spray applications. From fundamental fluid mechanical models to highly computational models, scientific achievements have been made in the development of a more realistic representation of pesticide movement in the vicinity of the sprayer, within plant canopies, and in the atmosphere. While current studies have prompted the confidence in numerical methodologies within limited scales of drift processes, this review emphasizes the need to integrate the existing knowledge and new technologies to optimize each sequential process from spray injection to drift mitigation. It is also important to combine efforts from experiments and numerical fluidmechanical modelling to produce more accurate predictions of the fate and transport of droplets under various atmospheric conditions.

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