Climate Behaviour and Plant Heat Activity of a Citrus Tunnel Greenhouse: A Computational Fluid Dynamic Study

: Response to the expanding demand for high-quality citrus saplings plants requires optimisation and a deep understanding of production climate behaviour. In this context, greenhouse production is the most used technique because it allows farmers to effectively monitor plant growth through production condition control, especially climatic parameters. The current work presents an analysis of climate behaviour and plant heat activity of a citrus sapling tunnel greenhouse in the middle region of Morocco. In this regard, a computational fluid dynamic (CFD) model was developed and validated with respect to temperature and relative humidity measured values. The specificity of this model is the inclusion of a new non-grey radiative and heat transfers physical sub-models to couple the convective and radiative exchanges at the plastic roof cover and crop level. The findings showed that using a green shade net increased the greenhouse shadow, and the layering of plastic and shade net significantly reduced solar radiation inside the greenhouse by 50%. Also, the greenhouse airflow speed was deficient; it cannot exceed


Introduction
Most studies show the importance of orange fruits for overall health through nutrients and protective plant compounds, including vitamins, minerals, antioxidants, and essential fibre. So, according to USDA2021 [1], global orange production for 2021/22 is estimated to be 49 million Tons. In Morocco USDA2021 [2] report and Spreen et al. [3] study showed that an area of 59.600 Ha was cultivated with citrus fruit, with an average production of 2.55 million metric Tons. Oranges' productivity strongly depends on many factors, especially the saplings' quality, which needs strict control under specific climate conditions. Recently, with the increased demand for high citrus saplings quality in Mo-rocco, the production was oriented towards using greenhouses environment systems, especially tunnel type.
Several authors' research studies have recently focused on the microclimate in tunnel-type greenhouses. We can cite Boulard's [4] numerical and experimental study of crop transpiration heterogeneity based on CFD and a global solar radiation model. Nebbali [5] compares three CFD turbulence models for simulating an empty tunnel-type greenhouse with lateral vents in the same context. The same author [6] conducted a CFD simulation of a tomato tunnel greenhouse with a (discrete ordinates) DO radiation model. The results highlighted the combined influence of sun position, wind direction and intensity on crop evapotranspiration rate in a tunnel tomato greenhouse. For their side, Bartzanas et al. [7] investigated the screen influence on airflow and temperature patterns inside the greenhouse and the effect of different wind directions without considering plant evapotranspiration.
Fidaros [8] conducted a simulation study on a ventilated tomato tunnel greenhouse during a solar day, considering only the plants' shortwave radiation (PAR) band optical properties. They demonstrated the important effect of external temperature variation on greenhouse climate determination. As a result of two parametric CFD investigations, Baxevanou [9] emphasised the important role of buoyancy forces and roof cover material proprieties on the developed flow field and the greenhouse's temperature. The first study dealt with the variation of intensity and incident solar radiation angle, and the second with the optical properties differentiation of cover materials. on the other hand, most citrus studies were conducted on mature trees Consoli [10], Er-Raki [11], Villalobos [12], Rana [13], Yang [14].
More recently, Bekraoui et al. [15] characterised a citrus tunnel greenhouse microclimate and determined citrus plants' transpiration during winter. Results show the high variation of greenhouse air and plant leaves temperature value and the low plant transpiration.
Various authors carried out greenhouse ventilation studies. Majdoubi et al. [16] modelled a one-hectare tomato Canary-type greenhouse. The model shows the significant effect of insect screen mesh, ventilation openings locations, and crop row orientation on climate homogenisation. Also, the greenhouse condensation phenomenon was studied by several authors; Piscia [17] study results of roof cover condensation for a four-span plastic-covered greenhouse during night-time show that condensation rate can be modelled as a logistic function of time. Condensation inside a single-slope Chinese greenhouse was modelled by Liu [18] with a focus on leaf condensation on the cucumber canopy, concluding that the condensation appears first on the roof cover. Crop transpiration was studied by Boulard [19] using a validated model for a Venlo-type semi-closed glass greenhouse. The cited model includes radiation transfer, crop transpiration, heat transfer and 2 concentration due to photosynthesis with a variation of air conditioners arrangements. Mesmoud [20] compared different greenhouses geometries (tunnel, Venlo and plastic vertical wall greenhouse) with a variation of the roof cover material; with all greenhouses cropped with tomato, the study used a two-dimensional CFD model.
Among the numerical studies on greenhouse climate, too few consider roof cover and crop radiative and heat effect variability on climate determination; most studies imposed temperature or heat flux as a boundary condition for the greenhouse roof cover.
In this context, the present study focuses on the microclimate simulation of citrus saplings' tunnel greenhouse, and analyses plant citrus radiative interaction and heat activity. For these reasons, a CFD model was developed considering the mutual interaction between the roof cover and crop and greenhouse environment using a user-defined function (UDF). The developed model is validated against experimental data.

Experimental Setup
The experiment was conducted on a tunnel plastic greenhouse (Figure 1) located in the middle region of Morocco (Latitude: 33°53′36″ N; Longitude: 5°32′50″ W; Elevation above sea level: 531 m). To reduce solar intensity inside the greenhouse, the roof is covered with a 300 µm polyethylene plastic film, superimposed with a green shade net. The greenhouse is cultivated in citrus saplings arranged in seven crop rows with a density of 32 plants per m −2 and a height of 0.9 m.
A weather station outside the greenhouse was used to acquire external climatic parameters, with an STH35 temperature and relative humidity sensor and an HYXC-HYGTR Pyranometer for global solar radiation measurements. SHT35 (A1, A2, A3, A4, A5, A7 and A9 Figure 1) sensors were used to measure greenhouse and crop rows' air temperature and relative humidity. Moreover, an ST-1307 solar power meter was used for net radiation measurements inside the greenhouse.

Computational Domain
The greenhouse was embedded in a computational domain ( Figure 2) (648 m, 648 m, 50 m). The computational mesh comprises 10.85 million total volume cells and 1.47 million cells inside the greenhouse. The mesh (Figure 3) was refined in regions of interest. The crop's lower part contains three layers, while the part with leaves is divided into five layers. To model conjugate heat transfers through the plastic cover, a boundary layer mesh was considered with a non-dimensional wall distance of 50 < + < 150 inside the greenhouse. The simulation was setup as provided in Table 1 and using boundary conditions from Appendix G.

Transport Equation
The governing equations for momentum, energy, species, turbulent kinetic energy, and turbulent dissipation rate can be written as transport equations: where can be any of the following quantities: mass, momentum, turbulent kinetic energy, turbulent kinetic energy dissipation rate, water vapour mass fraction and energy; the term represents the source term for each of the previously cited quantities.

•
Mass conservation equation • Momentum equation where is the Reynolds stress tensor • Turbulent kinetic energy equation where μ is the turbulent viscosity: Turbulent kinetic energy dissipation rate equation where the diffusion flux of species "i" is: • Energy equation , other symbols are listed in Appendix H.

Buoyancy Modelling
Given that the modelled greenhouse is almost closed, the buoyancy forces will play an essential role in convective transfers. To account for these forces, the variation of fluid density with temperature was considered through the ideal gas law:

Radiation Modelling
A two-band model with a cut-off wavelength = , = 4.2 was adopted for radiation as the covering material has maximal transmissivity in the shortwave (PAR) spectrum and minimal transmissivity in the longwave band [21]. The radiative transfers at position ⃗ in direction ⃗ are modelled considering the radiative transfer equation (RTE) for each wavelength band " ": where ∆F is the fraction of blackbody emissive power in the spectral wavelength band " ".

Heat Modelling inside Solid Cover
The cover of the greenhouse is subject to heat conduction and radiation absorption: To release constraints on meshing imposed by the very thin cover, we used a solid cover with a much larger thickness. The materials' proprieties were modified according to our model (Appendix F table A1).

Flow Through Plants
The crop was modelled as two parts; the lower one contains only trunks but no leaves; thus, having height porosity was modelled as air only. The upper one containing all the leaves was modelled as a porous media governed by a simplified version of Darcy-Forchheimer (Equation (5)) for the drag; giving a source term in momentum transport Equation (1): where is the air speed, and is the drag coefficient set to = 0.32 [22].

Canopy Water Vapour Transfers
The water vapour transfer from crop leaves to surrounding air must overcome two resistances in series. The first is the stomatal resistance , and the second is the aerodynamic resistance due to the boundary layer around the leaf. Therefore, a quantity of water vapour will be added as a source to the vapour transport equation in the crop's porous media: Using data from [23], We modelled the stomatal resistance as an exponential function of the net radiation: where = 576.64 s m −1 is the minimal stomatal resistance, a = 0.00848 and b = 5.6682 are two constants.
On the other hand, we adopted the model developed by Boulard et al. [24] For aerodynamic resistance: With "h" the convective transfers coefficient of the leaf: Canopy Heat Transfers The radiation absorbed throughout the canopy will be converted into latent heat , used for phase change of water, and sensible heat due to convective transfer between the leaf and the nearby air; another part ,1 will be emitted as radiative heat in the long wavelength band (in band 1): where is the net radiative flux in −2 , is the latent heat flux in −2 and is the sensible heat flux in −2 . The crop temperature can be computed as follows:

Canopy Radiation Interaction
Several authors ( [16,[25][26][27][28]) used Beer-Lamber's law to describe the canopy interaction with photosynthetic radiation. The previous method suggests that the direction of the radiation is known, which necessitates modelling of the sun's direction. Other authors ( [19,29]) used a more straightforward model setting an absorption coefficient equivalent to the radiation extinction in the crop. In the current study, the latter approach is used; but with further developments to account for absorption, scattering and emission in each direction of the discrete ordinate model (as developed in Appendices A, B, C and D). the emission from the crop was also accounted for in a modified RTE (Equation (27)): Band 0: Band 1: The radiative transfer equation then becomes: Appropriate source terms were used to couple energy and radiation as described in the previous subsection (Canopy heat transfers, section 2.3.5.3). And the emission term ,1 of equation 19 was developed in Appendix E.

Model Validation
The numerical model validation with respect to experimental data shows, firstly ( Figure 4), that the simulated temperature values along the greenhouse length at 1.34 height above the ground ranged between 24.94 °C and 32.85 °C with an average value of 29.87 °C. The comparison between the simulated and measured temperature values demonstrates that the differences are weak, the average difference reaches 0.18 °C, and its maximum does not exceed 2.43 °C. On the other hand, relative humidity analysis ( Figure  5) at the same locations as temperature illustrates its variation between 59.5% and 100%, with an average of 74.5%. According to relative humidity experimental values, the difference varies between 1.53% and 10.5%, with an average difference value of 1.57%. The weak deference between the simulated and the experimental values shows a good agreement between simulation and experiment.

Greenhouse Radiative and Heat Transfers Analysis
Examination of the greenhouse incident short wavelength radiation field ( Figure 6) shows that only about 53% of the external solar intensity was transmitted to the greenhouse, which means that the shading reduces global solar radiation transmission by more than 30%. This important quantity increases the plastic roof cover energy storage and temperature value and significantly increases convective heat transfers and long wavelength radiation effect on climate parameters variation. The greenhouse's short wavelength radiation distribution is heterogeneous and degrading horizontally from the South-East side to North-West one (SE at x = 0.6 m: 205 (W m −2 ), centre at x = 4.3 m: 306 (W m −2 ), NW at x = 8.0 m: 521 (W m −2 ) at height 1.34 m from the ground) and vertically from the roof cover to the soil level (near the roof cover: 396 (W m −2 ), just above the canopy 300 (W m −2 ), the area between the canopy and substrate and at soil level: lower than 37 (W m −2 )).

Greenhouse Microclimate Details
The sensible heat lost by the crop (especially the top parts) (Figure 7e ) and (Figures  8a,b) will heat air and reduce its density and leading to the creation of natural convection loops ( Figure 9). So, convective transfers between air and hotter cover (Figure 10a) accentuate these convective loops.  Despite the high shadow and the low solar radiation inside the greenhouse, it was observed from Figure 10a,d that greenhouse air temperature and relative humidity fields are heterogeneous. The observed gradients were caused by natural convection heat transfer currents transporting air with a higher water mass fraction from the bottom of the canopy to the top part of the greenhouse (Figure 10b). As discussed in the previous subsection, the North-West area produces more water vapour. The created air loops (Figures 9 and 10c) generate two distinct regions with distinctly different mass fractions; the two regions' climate parameters will vary along the "z" axis depending on loops recirculation in each section.
An increased relative humidity (exceeding 100%) (Figure 10d) was observed at the bottom of the greenhouse, especially in the area near the substrate. This increase is due to the accumulation of water vapour in these regions caused by the weaker greenhouse ventilation rate, which may lead to the development of fungus and condensation of water on the greenhouse surface and the leaves in the lower parts of the canopy. Therefore, condensation model development is necessary to improve the greenhouse simulation climate study (Liu [18]).
Investigation of climate parameters across the length of the greenhouse (Figure 10) shows considerable heterogeneity in temperature field and relative humidity at section L = 32.5 m; this is due to the absence of crops in this region. On the other hand, the crop in the other regions (L = 19 and L = 46) forms a barrier to airflow free circulation and contributes to air homogenisation in these areas. Therefore, crops rows orientation has an essential effect on climate homogenisation and optimisation.

Conclusions
The main contributions of the present research study can be summarised in analyse of a tunnel citrus greenhouse climate through the development of a 3D CFD model coupled with two sub-models to simulate the roof cover and crop activity contribution to greenhouse climate determination. The developed local model for non-grey radiative transfers inside the crop considers absorption, scattering and emission for each band and radiation direction. Moreover, a modified radiative transfer equation was proposed to account for the radiative transfers inside the porous medium consisting of air and crop. Results confirm the heterogeneity of temperature and relative humidity fields and indicate the significant role of the natural convection heat transfer mode in closed tunnel greenhouse climate distribution. Also, simulation results reveal the high mass fraction of water vapour inside the greenhouse, especially in areas at the bottom of the canopy, where water vapour exceeds the saturation level, and condensation may occur. The previous conditions and the plant heat activity simulated results, in particular, low transpiration value, allow us to declare the malfunction of the citrus plants' activity under this climate. As a result, we advise farmers to increase inside sun radiation across the greenhouse by using a more transmissive shading net with a smaller roof-covered area, to choose a North-South greenhouse orientation for good solar radiation distribution, to optimize crop row orientation for a better airflow circulation and climate homogeneity, and to improve ventilation rate through greenhouse ventilation openings perpendicular to the prevailing wind direction (East and West sides).
Even if the developed crop model was used, its mathematical formulation does not limit its use for such cases. As a perspective, it could be used and validated for orchards and urban green areas. Also, the cover model could be used to relax meshing constrain for regions near thin solids (opaque or semi-transparent) in energy-related domains such as heat exchangers and solar collectors.

Acknowledgments:
We are grateful and indebted to many members of the domain teams who helped with the tremendous amount of practical work that enabled us to carry out this study.

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

Appendix A. Crop Radiation Absorption Coefficient for Long Wavelength (Band 1)
The radiation heat absorbed by a leaf in band 1: The radiation absorbed by all "N" leaves inside a volume "V" would be: ,1 = ,1 = � | ⃗. �⃗| 1, ∆ Thus, we can define a crop absorption coefficient for each radiation direction in band 1 as follows: With " ": the ray direction to the leaves normal. (Band 1) The reflected radiation by a leaf in the long wavelength band is:

Appendix B. Crop Radiation Scattering Coefficient for Long Wavelength
Using the same methods as for absorption, the reflected radiation by volume "V" is: A scattering coefficient for each radiation direction in band 1 would be:

Appendix C. Crop Radiation Absorption Coefficient for Short Wavelength (Band 0)
For short wavelength, we consider each leaf of the crop as a grey opaque wall; thus, it can only reflect and absorb radiation: Using the same method as for long wavelengths we This leads to defining a crop absorption coefficient for each radiation direction in band 0 as follows:

Appendix D. Crop Radiation Scattering Coefficient for Short Wavelength (Band 0)
To account for direct solar radiation scattering inside the crop, we consider the reflection: Integrating over a volume as for long wavelength: A scattering coefficient for each radiation direction in band 0 would be:

Appendix E. Crop Radiation Emission
The radiation emission in short wavelength can be neglected (because ∆F 0 , the fraction of blackbody emissive power in the spectral wavelength band 0, is negligible), so we consider that all the radiation emission is in the long wavelength band. For a leaf of the crop, the emission of radiation is: The radiation emitted by all "N" leaves inside a volume "V" would be: With blackbody emission at leaf temperature. Finally:

Appendix F. Cover Thin Wall Model
The heat conduction inside a volume element ( Figure A1) of the greenhouse cover, considered a thin wall, happens mainly in the thickness direction. Equation (4) then becomes one-dimensional: where x, y and z are locale coordinates, and x is normal to the cover in the thickness direction.
Nondimensionalizing "x" with: ′ = , we get: Without counting the dimensionless variables. There are ten variables in this equation ( = 10), that is (with their respective dimensions): [ ] Such that M, , , , T and are the respective six physical dimensions ( = 6) of mass, lengths, time, and temperature.
According to Buckingham's theorem, the previous heat equation can be expressed using = − = 4 dimensionless numbers: For the boundary condition: The variables in this equation are four ( = 4), i.e., (with their respective dimensions): Therefore, the dimensional matrix of this equation is: The number of physical dimensions involved in the equation is ( = 3). By Buckingham's theorem, the previous Neumann boundary condition can be expressed using one dimensionless number = − = 1: The same procedure makes it possible to obtain a dimensionless number for the boundary condition equation along the "z" axis: A thermally equivalent greenhouse cover must have the same values for these dimensionless numbers. Assuming that the mass m, the time t, the temperature ( ′), the thermal energy on the boundary Q, the thicknesses and , the area and the Boltzmann constant are the same for the model as for the reality. And using the index "m" for the cover model, we find the relations between the model and reality: With respective material proprieties used for this study: