—Prevention of frost in horticulture is important, but challenging, and its realization, especially using green-energy sources, has a huge societal impact. In this paper, ﬁrst we suggest an integrated solar photovoltaics (PV)-assisted framework where solar energy will be used as a secondary application for frost prevention (the primary application is electricity production for grid). Optimal design and operation of the suggested integrated system require detailed thermal modeling of air dynamics in the orchard, integrated system optimization and control tasks. Second, in this paper we address the ﬁrst task above: development of a novel, sophisticated parametric computational ﬂuid dynamics (CFD) model for orchard air thermal dynamics for different orchard parameters (such as fruit type, climate, the number of trees, their sizes, distance between them, etc.) and boundary conditions. Finally, the use of developed parametric CFD model is demonstrated through a case study to calculate the annual minimal thermal energy required to prevent frost in a test apricot orchard located in Malatya, Turkey, which is the world capital for dry apricot production. Index Terms —Horticulture, Frost prevention, Solar energy, Photovoltaics, CFD modeling. I. I NTRODUCTION High quality and diversity horticulture outputs (fruits) con- tribute to national health and wealth through providing a rich nutrition content and export incomes (for example,  15-20 billion efor EU member states [1]). However, the horti- culture activity and diversity level in each country is different, which depends on climate and land properties. For many fruit types which are vulnerable to low temperature levels, frost occurrence is a serious issue affecting both quantity and quality of yield. Although withering of fruit buds/ﬂowers due to frost is dependent not only on temperature levels below 0 °C, but also on humidity, wind, whether the period is budding, ﬂowering or early fruit period, cultivar type, the tree, etc., temperature levels below 0 °C are dangerous and are the main parameters for yield loss due to frost [2]-[3]. For example, as shown in Figure 1, on March 30, 2014 at Malatya, Turkey (a world capital for dry apricot production), the temperature level was -3.2 °C which caused 95 % of yield to be lost. The economic value of this yield loss was equivalent to 1.2 billion dollars. As observed, in some other years also the yield is very Ercan Atam is an independent researcher supported by T ¨UB˙ITAK. Address: Res ¸itpas ¸a Mahallesi, Katar caddesi, Teknokent ARI 1 Sitesi, No: 2/5/19, 34467, Sarıyer, Istanbul, Turkey. E-mail:ercan.atam@gmail.com Se-Woon Hong is with the Department of Rural and Biosystems Engineer- ing, Chonnam National University, Gwangju, Republic of Korea: Address: E-mail:hsewoon@chonnam.ac..krlow. These are years where frost has caused signiﬁcant yield loss. Fig. 1: Apricot production per year in Malatya, Turkey (the world capital for dry apricot production) showing 95 % yield loss on March 30, 2014 [4]. It is important to note that “frost” is deﬁned as withering of fruit buds/ﬂowers which are damaged and are not able to produce any yield. On the other hand, “freeze” is deﬁned as temperature levels equal to or below 0 °C. Freeze does not always mean frost but for some horticulture regions, for example in Malatya, Turkey, frost occurs every two to three years when freezing occurs [5]. Table I shows the number of freezing times over the last two decades for the Malatya region of Turkey for the spring season. From the presented data, one can see that during 1997-2016 at least one time (except 2001) freezing occurred during March-April months, which is potentially enough to cause frost. The currently used frost prevention measures in Europe and other countries in the world can be divided into two main categories: passive frost prevention methods and active frost prevention methods. The passive methods include site selec- tion, selection of species, cultivars and cultural techniques, and soil management [7]. Passive methods have very limited protection capabilities in general. Active protection methods include use of wind machines for air mixing, irrigation, man- made fog, ﬂooding, using heaters, using covers, and using chemicals [8]. Active methods have better frost protection capabilities compared to passive ones but their investment cost and/or operational cost is higher. As a result, they are not used in a large-scale extent. As a simple active frost protection example, in some regions distributed straw piles are burned to IEEE TRANSACTIONS ON SUSTAINABLE ENERGY 2 TABLE I: Temperature analysis for Malatya, Turkey over last 20 years during March & April which is the critical period for apricot frost phenomena [6]. year 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 # of freezing cases24 9 4 19 0 15 7 4 1 1 Daily minimum temperature (° C)-9 -5.2 -2.3 -7.3 2.6 -1.2 -11.7 -5.7 -3.2 -0.9 year 2007 2008 2009 2010 2011 2012 2013 2014 2015 2016 # of freezing cases1 1 4 1 3 17 5 2 4 2 Daily minimum temperature (° C)-0.9 -0.2 -1.8 0 -1.4 -6.5 -3.9 -3.4 -2.3 -2.7 Fig. 2: PV-integrated frost prevention system schematics. protect against frost in cold days, but this is not an effective, practical, and economically feasible solution in the large- scale. Recently, another method of active frost prevention was tried: researchers in ˙In¨on¨u University (a university in Malatya, Turkey) investigated growth of apricot trees in greenhouses. This project was successful, and during the year of 2014 when 95 % of apricot yield was lost due to frost, the apricot trees in the greenhouse gave a full crop. The researchers next decided to investigate the economic feasibility and wide- spread applicability of such a project. In this paper, we present a contribution to frost prevention in horticulture by proposing a new idea, which we believe it can have a huge societal impact for all nations. As illus- trated in Fig 2, we suggest an integrated solar photovoltaics (PV)-assisted framework where solar energy will be used as a secondary application for frost prevention (the primary application is electricity production for the grid). The PV- based system will mainly be used to produce electricity to be added to the grid. It will be necessary to store some of the produced electricity when the frost periods approach. The size of storage unit depends on the size of the fruit orchard and the predicted frost level. The electricity from the grid or storage device is used to heat air by an array of electricair heaters where, again, the number of heaters in the array depends on the size of the orchard. The heated air will be blown by fans on fruit trees during freezing days. The frost prevention control system will control the blowing time and the mass ﬂow rate of the heated air. PV-based frost prevention system is a good choice for regions with a rich solar radiation level. Optimal design of the integrated system is crucial for minimization of investment cost and payback time, and design of an automatic control system is essential for minimization of operational costs. However, to optimally design and control the integrated system, ﬁrst a detailed parametric computational ﬂuid dynamics (CFD) model is required to simulate the air thermal dynamics in the orchard under different parameters (such as fruit type, climate, the number of trees, their sizes, distance between them, etc) and boundary conditions. In this paper, we develop an advanced parametric CFD model for such purpose and demonstrate its use to calculate the annual minimal thermal energy required to prevent frost in a test apricot orchard located in Malatya, Turkey. In literature, there exist some parallel studies on climate simulations and determination of air ﬂow patterns, for ex- ample, CFD-based simulation of climate in greenhouses [9], [10], [11] and air ﬂow around and through trees [12], [13]. In IEEE TRANSACTIONS ON SUSTAINABLE ENERGY 3 the recent work [15], Hong et al. investigated airﬂow inside tree canopies discharged from air-assisted sprayers through a sophisticated CFD modeling. In this paper, both the considered research problem and the objectives are different than these previous studies, and hence the added contributions are new. Indeed, to our best knowledge, (i) the proposed integrated framework given in Figure 2 to prevent frost in horticulture applications is the ﬁrst one in the open-literature, (ii) the proposed idea is a “green” solution for the devastating frost problem in horticulture applications, (iii) the developed para- metric CFD model will be a backbone to develop an optimiza- tion algorithm for minimal sizing of the PV-based integrated system and to develop a frost preventing control system which together minimize the investment and operational costs. The rest of the paper is organized as follows. In Section II, the parametric CFD model development details are given. Section III presents the use of the developed parametric CFD model for annual minimal required frost preventing energy for a test apricot orchard located in Malatya, Turkey. Finally, the main ﬁndings of this study and conclusions along with some future research directions are given in Section IV. II. P ARAMETRIC CFD MODEL DEVELOPMENT A. Orchard parameters and computational domain A generic rectangular prism type fruit orchard with dimen- sions LmWmHm consisting of nsubzones is considered. We assume one hot air blower per zone to heat air around trees in each zone for frost prevention during frost periods. The distance between roots of trees is sm. The 2-D geometry of the desribed generic orchard is given in Figure 3. We assume the diameter (the approximately circular region covered by branches of a tree) and height of each tree to be dm and hm, respectively. The climate data (consisting of wind speed, wind direction, sky and ground temperature, etc.) is also taken as an orchard parameter. Fig. 3: Two dimensional domain for a generic rectangular prism type fruit orchard. Here, “x” denotes trees, “t” hot air blowers, vathe wind speed, and wind direction. A computational domain with a size of 2Lm2Wm 10Hm was designed by giving minimum L=2m and W=2m for buffer zones on longitudinal and transverse edges of the orchard, respectively. Additionally, these buffer zones should be greater than 10H(His the height of trees), whichis the minimum distance required between the objects and nearby boundaries according to Bjerg et al. [18] to ensure computational stability against the propagation of pressure between the objects and nearby boundaries. The height of the computational domain should be greater than 5Hsuggested by Bjerg et al [18] for computational stability. The fruit trees were modelled by porous media with a curved cone shape, which is similar to the shape of most of typical fruit trees (Figure 6). The trees were modelled in the computational domain by setting appropriate source terms for momentum loss and turbulence generation to grid cells. The source terms were calculated only for the grid cells located within the trees. This modeling method and the details of the source terms for momentum loss and turbulence generation due to tree canopies are described well by Hong et al. [15]. B. CFD for canopy air ﬂows Canopy air ﬂows within and around the orchard were solved by conservation equations for mass, momentum, and energy. For the mass and momentum equations, this study uses the Reynolds-averaged Navier-Stokes (RANS) equations which is a CFD approach to decompose the Navier-Stokes equations into the mean and ﬂuctuating components. The RANS equations are given as @ @t+@(ui) @xi= 0; (1a) @(ui) @t+@(uiuj) @xj=@p @xi+@ @xj[ (@ui @xj+@uj @xiu′ iu′ j)] +Su; (1b) where is the density of air, uis the mean air velocity, u′is the ﬂuctuating part of air velocity, pis the air pressure, is the air viscosity, and Suis the source term for the momentum loss due to tree canopies. The Reynolds stresses, u′ iu′ j, were modelled by the mean velocity gradients using the Boussinesq hypothesis [21] in order to close the RANS equations: u′ iu′ j=t(@ui @xj+@uj @xi) 2 3( k+t@uk @xk) ij; where tis called the turbulent air viscosity and it was computed as a function of k(turbulent kinetic energy) and w (speciﬁc dissipation rate), which were modelled by the kw turbulence transport model [22]. The source term for the momentum loss due to tree canopies were modelled by the pressure loss coefﬁcient, Cir, as Su={ Cir1 2juju;ifcg2Cg;tree 0; ifcg=2Cg;tree Here the pressure loss coefﬁcient Ciris inversely proportional to the density of the tree canopy; cgis the speciﬁc grid cell where the source term is solved; and Cg;tree is the group of grid cells located within trees [15]. Transport of energy in the canopy ﬂows is calculated by @ @t(E) +∇ (u(E+p)) =fkeff∇ThJg+Sh;(2) IEEE TRANSACTIONS ON SUSTAINABLE ENERGY 4 where Eis the total energy (sum of internal and kinetic energy) and calculated as E=hp +u2 2withhdenoting the sensible enthalpy (absolute enthalpy minus enthalpy of formation); keffis the effective conductivity, which is given by the thermal conductivity of the material and the product of the turbulent viscosity and the heat capacity divided by the turbulent Prandtl number; Tis the air temperature; Jis the diffusion ﬂux; and Shis the source term for the heat exchange between air and trees and the heat generation by heaters located at all zones. The source term for the energy equation can be decomposed into two terms: Sh=Sse+Ssh (3) where Sseis the sensible heat source from the plant to the air (W=m3); and Sshis the sensible heat generation from space heaters located in each zone ( W=m3). C. CFD-based thermophysics analysis of frost in fruit or- chards As shown in Figure 4, the thermophysics of fruit orchards involves complicated energy ﬂuxes between the trees and their surrounding. During the night time, trees in the orchard lose energy through radiation upward and at the same time gain the energy downward from the sky. During frost periods, the trees lose more energy than gained from the sky because at these times the sky temperature is always lower than the temperature of the tree or the ground. Under clear sky, more heat is lost than gained, while cloudy or foggy sky decreases the energy loss because the temperature of the clouds is higher than the sky temperature. Besides the radiation transfer, trees also exchange their energy with the ground through conductive heat transfer, and with the surrounding air through convective heat transfer. The combined heat transfer from and to the trees determines the variation of their canopy temperature. Fig. 4: Energy ﬂows around trees in a fruit orchard. Physical parameters associated with the canopy of a fruit orchard are described in Table II. While these parameters vary by individual trees and their growth conditions, this study assumed that all trees were identical with respect to the shape and canopy condition, which are, in general, reasonable assumptions for typical trees in an orchard. Next, the expressions for each of the energy transfers in Figure 4 using the canopy parameters in Table II are given. The outgoing radiation from the tree canopy can be calculated using the Stefan-Boltzmann law asTABLE II: Physical parameters associated with the canopy of a fruit orchard modelled in parametric CFD model. Variable description Variable notation Canopy width (m) wcanopy Tree height (m) htree Trunk diameter (m) dtrunk Trunk height (m) htrunk Leaf area density ( m2m3)LAD Leaf area index ( m2m3) LAI Leaf characteristic dimension (m) LCD Pus=1 LAIϵT4 c; (4) where Pusis the radiation heat ﬂux from the tree canopy to the sky per leaf area ( W=m2);LAIis the leaf area index of the tree ( m2=m2);ϵis the emissivity and approximately 0.74 for the night sky; is the Stefan-Boltzmann constant, 5:67108W=(m2K4); and Tcis the temperature of the tree canopy ( K). Similarly, the downward radiation ﬂux from the sky to the tree canopy can be calculated by applying the sky temperature to the Stefan-Boltzmann law. However, atmospheric and cloud conditions disturb the exact calculation of the downward radiation. In addition, water vapor in the atmosphere absorbs the radiation and creates another source of radiation. This complexity has been solved through empirical formulas, one of which is a modiﬁed version of the Swinbank model [23] of night time downward radiation Goforth et al. [19]: Pds=8:781012 LAI( 1 +KC2) T5:852R0:07195 H ;(5) where Pdsis the radiation heat ﬂux from the sky to the tree canopy per leaf area ( W=m2);Kis a constant, 0.34 for cloud height <2 km, 0.18 for 2 km <height <5 km, and 0.06 for height <5 km; Cis the cloud cover constant varying 0.0 for clear sky through 1.0 for totally overcast; Tis the air temperature (K); and RHis the relative humidity (in %). The convective heat exchange between the tree canopy and the surrounding air occurs on individual leaves, which is calculated based on the temperature difference air and canopy and the aerodynamic resistance of the boundary layer created on the leaf surfaces. In [24] the convective heat ﬂux from/to the leaves was calculated as Pa=aCp;aTcT ra; (6) where Pais the convective heat ﬂux from the tree canopy per leaf area ( W=m2);aCp;ais the air speciﬁc heat (J=(m3K)); andrais the aerodynamic resistance of the leaf boundary layer (s=m) calculated as ra= 305√ Lc juj+ju0jwith ju0j=√ 0:5gjTcTjLcpTc; IEEE TRANSACTIONS ON SUSTAINABLE ENERGY 5 where gis the gravitational acceleration, 9.81 m=s2;Lcis the characteristic length of leaves ( m);jujis the air speed ( m=s); andju0jis the buoyancy air speed ( m=s). The conductive heat exchange between the trees and the ground happens when the temperature difference between two bodies drives the heat ﬂow. The heat transfer through the tree stem is calculated by a general thermal conduction theory and can be adapted with respect to the tree canopy as follows Pg=UAt ACLAI(TcTg); (7) where Pgis the conductive heat ﬂux from the tree canopy to the ground per leaf area ( W=m2);Uis the conductance of the tree stem ( W=(m2K));Atis the area of the tree stem ( m2); Acis the ground area of the tree canopy ( m2); and Tgis the temperature of the ground ( K). The temperature of the tree canopy changes by gaining or losing energy from or to the surrounding as shown in Figure 4. Therefore, it can be calculated dynamically as cCp;cdTc dt=1 tc(Pus+Pds+Pa+Pg); (8) where cis the density of the tree canopy ( kg=m3);Cp;cis the speciﬁc heat of the leaves, 4186J=(kgK )for fresh as well as dry leaves (Jayalakshmy & Philip [20]); and tcis the thickness of the leaves, assuming 0:3103m. At the same time, the surrounding air will gain or lose energy by convective heat transfer from/to the tree canopy (which will affect the air temperature slightly) according to Sse={ LADPa;ifcg2Cg;tree 0; ifcg=2Cg;tree(9) where Sseis the sensible heat source from the plant to the air (W=m3) and LADis the leaf area density. Sseis included in the energy source term, Sh, of the energy equation for air ﬂow in the canopy. D. Boundary conditions As shown in Figure 5, the computational domain for the fruit orchard is taken as a rectangular prism. The temperature boundary conditions for the ground and the sky are speciﬁed as a constant assuming their surfaces are maintained at ﬁxed temperatures during the frost event. Four side boundaries, including left, right, front, and back, are set as an inlet or outlet depending on a wind direction. The temperature boundary condition for the inlet is speciﬁed as a logarithmic temperature proﬁle where the temperature at the bottom end (=the ground) and at the top end (=the sky) are set to Tground andTsky, respectively. The temperature boundary condition for the outlet is set as a no-ﬂux boundary. With reference to Figure 5, we have the following set of boundary conditions: TjEbottom =Tground ; TjEtop=Tsky; (10a) TjEinlet =f(z) =c1log(z+ 1) + c2; (10b) @T @n Eoutlet= (∇Tn)jEoutlet = 0; (10c) Fig. 5: Computational domain for the three dimensional or- chard. where ndenotes the wind direction; Einlet,Eoutlet are the edges associated with the inlet and outlet, respectively; f(z)is the temperature proﬁle (K) with respect to the altitude, z; and c1andc2are constants to satisfy the condition that temperature at the ground and the sky are Tground andTsky, respectively. E. CFD simulation procedure The simulation procedure of the parametric CFD model is as follows. First, we need to set the orchard parameters (fruit type, orchard size, etc) and boundary/initial conditions (which are considered as part of climate data). Second, to simulate the model we need to specify the operational principle of the control unit (the unit with number “6” in Figure 2). Here, the controller type we use is rule-based controller whose operational principle is as based on the following rules: For each zone in Figure 3, a separate controller will be used for controlling blow of heated air on trees. Let us denote the controller for the i-th zone as C i, i= 1;2;  ; n, where nis the number of zones. Controllers work independent of each other. The control-time step, ts(which is the interval at which a controller provides a ﬁx control input (heat ﬂux) to the orchard, but these control inputs can be vary between control time steps), is set as ts=30 seconds. In fact, ts can be considered as a “simulation parameter” and can be set to different values depending on the frost scenario. For example, in case of high wind speeds tscan be set lower values, such as 10 seconds. In principle, tscan be different for each C i, but better to be the same to simplify the simulation procedure. Let the heat ﬂux provided to each zone by CibeQi(t) with minimum and maximum values of QminandQmax, respectively. Qmincan be set to zero simply and Qmax depends on the capacity of air heaters in Figure 2. Let the control period over which to control the system be dived into a set of Nintervals where the length of each interval ists. Let us denote the time points at the beginning of each interval by t1= 0 < t 2< t 3<  < t N= (N1)ts.Qi(tk),i= 1;  ; n,k= 1;  ; N, is determined as follows. For zone iat time tkcalculate the IEEE TRANSACTIONS ON SUSTAINABLE ENERGY 6 TABLE III: Physical characteristics of apricot trees modelled in CFD simulation. Variables Values Canopy width (m) 6.0 Tree height (m) 9.0 Trunk diameter (m) 0.4 Trunk height (m) 1.5 Leaf area density ( m2m3) 2.0 Leaf area index ( m2m3) 3.0 Leaf characteristic dimension (m) 1.0 average temperature of each tree in the zone and let Tmin i be the minimum of all of these calculated average tree temperatures. Then, we determine Qi(tk)as Qi(tk) ={ 0; ifTmin i0°C Qmax;otherwise :(11) III. C ASE STUDY :USE OF PARAMETRIC CFD MODEL FOR FROST PREVENTING ENERGY CALCULATIONS In this section, we will demonstrate, among many, one useful application of the developed parametric CFD model: to calculate the required minimum frost prevention energy for each of three frost episodes (corresponding to mild level frost, medium level frost, and severe level frost) for a test apricot orchard located in Malatya, Turkey. With reference to Figure 3, the test orchard dimensions are given as L= 200m; W = 100 m; s = 10m;, and we assume that there are 10 zones ( n= 10 ). The parameters associated with apricot tree canopy are given in Table III. A typical apricot tree with the canopy parameters in Table III is given in Figure 6. Finally, Fig. 6: A typical ﬂowered apricot tree in the considered test orchard. the positions of 10 hot air blowers are speciﬁed in Table IV. For the simulations, the control time step is 1 minute and capacity of each hot air blower is Qmax = 100 W=m3. The ﬂuent model of the heated orchard, the variation of average tree temperatures for each zone, and control inputs for each zone are given throuh Figures 7-8.TABLE IV: Hot air blower identity numbers and their plane positions. There is one blower per each zone and in total 10 blowers. Blower number and position (m) Blower number and position (m) 1!(20, 25) 6!(100, 75) 2!(20, 75) 7!(140,25) 3!(60, 25) 8!(140, 75) 4!(60, 75) 9!(180, 25) 5!(100, 25) 10!(180, 75) Fig. 7: Fluent model of test orchard. The colors indicate the tree canopies. The vertical plane is located at the middle of orchard. The horizontal plane is located at 3m height. Fig. 8: Average temperature (per zone) variation of trees for each zone over a period of 15 minutes. Fig. 9: Control inputs for each zone over 900 seconds. The control time step is 1 minute and capacity of each hot air blower is 100W=m3. IEEE TRANSACTIONS ON SUSTAINABLE ENERGY 7 IV. C ONCLUSIONS AND SOME FUTURE RESEARCH DIRECTIONS In this paper, ﬁrst we proposed a novel and green energy- based active frost prevention method for horticulture applica- tions. The proposed integrated framework uses photovoltaics as a secondary application (the primary application being the classical electricity addition to grid) during frost periods for large-scale fruit orchards. Next, a parametric CFD model was developed for air thermal dynamics simulation of heat fruit orchards. The developed parametric CFD model is essential for (i) integrated system optimization to minimize investment cost, payback time and to maximize income from sold electricity to grid; (ii) to develop reduced-order control models for model- based control approaches (such as nonlinear model predictive control); and (iii) test the developed control systems on the CFD model under different parameters before testing the controller on the real system. Finally, in this study the use of developed parametric CFD model was demonstrated through simulation of the air thermal dynamics of a test apricot orchard located in Malatya, Turkey, which is the world capital for dry apricot production. This test case considered different frost levels and we calculated the required minimum thermal energy to prevent frost in each case, where a rule-based control system was used to control the operation of hot air blowers. There are many further research rooms on the presented subject. Some of them can be listed as follows. (i) Developing a control-oriented model and using this model to apply NMPC for frost prevention; (ii) Determining the optimal number and distribution of hot air blowers; (iii) Development of optimized rule-based controllers from NMPC-type controllers through machine learning; (iv) In addition to hot air blowers, inclusion of water sprayers around perimeter of the orchard during frost periods; (v) Partial block of wind via a mobile vertical covering system, (vi) Economic feasibility analysis of alternative energy storage choices instead of batteries (such as energy storage via compressed air) or green energy choices instead of PV (for example, ground-coupled heat pumps [16], [17]). All these approaches have potential in lowering the investment cost and payback time signiﬁcantly, and they will studied by the authors in the future. ACKNOWLEDGMENT The ﬁrst author would like to thank T ¨UB˙ITAK for funding (through T ¨UB˙ITAK 1512 program) his project entitled “De- velopment of tools for solar-energy assisted frost prevention systems in large-scale orchards”.