CFD Modeling of the Microclimate in a Greenhouse Using a Rock Bed Thermal Storage Heating System

: The rock bed heating system is a more cost-effective concept for storing thermal energy use in greenhouses at night during the cold winter season. This system is considered an environmentally friendly solution compared to conventional heating systems that rely on fossil fuels. Despite the abundance of research on thermal energy-based heating systems, only limited work on climate modeling in greenhouses using rock bed heat storage systems has been reported. To ﬁll this research gap, this study aims to simulate the microclimate in a greenhouse equipped with a rock bed heating system using computational ﬂuid dynamics (CFD) models. User-deﬁned functions have been implemented to account for the interactions between the plants and the air within the greenhouse. Crop rows and rock bed blocks have been considered as porous media with their dynamic and thermal proprieties. The model’s accuracy was approved by comparing simulated and experimental climate parameter data from the greenhouse. The model’s ability to predict temperature, humidity, and air velocity ﬁelds in the greenhouse as well as in the rock bed system during both phases of energy storage and restitution was demonstrated. The thermal, dynamic, and hygric ﬁelds were accurately replicated with this numerical model. The growing zone had a vertical temperature gradient between the ground and the greenhouse roof, as well as high humidity. The distribution of temperature ﬁelds along the rock bed blocks showed a signiﬁcant temperature gradient between the air inlet and outlet in the blocks during the two phases of heat storage and restitution. As a result, the model could be useful for sensitivity studies to improve the performance of this thermal storage heating system.


Introduction
The good plant growth in greenhouses mainly depends on surrounding climate parameters such as CO 2 concentration, temperature, relative humidity, and solar radiation.A greenhouse provides the advantage of being able to better control these environmental conditions compared to the open field.In greenhouses, climate parameters can be maintained around set points using several climate control systems to provide optimal plant growth conditions.In winter, the optimum temperature can be maintained thanks to greenhouse heating systems [1].However, most conventional heating systems use fossil fuels, which have a negative impact on the environment.In view of the emergence of concerns and awareness on environment protection and sustainable development, the transition to more environmentally friendly heating systems is a necessity.
Among the eco-friendly heating systems using renewable energies are systems based on the storage of solar thermal energy.The principle is based on the use of materials with high thermal inertia that store energy from the sun during the day and release it into the greenhouse air at night.
Several thermal solar systems have been proposed by researchers to heat greenhouses.The most widespread are based on the storage of solar thermal energy in rock beds or in phase change material.
The rock bed heating system is based on sensible heat storage with air as the energy transport mechanism.The rock bed blocks are placed under the soil of the greenhouse.During the day, excess heat is transferred from inside the greenhouse to the underground storage with the help of fans.At night, the process is reversed.The cool air is moved through the storage, where heat is transferred from the rocks to the colder air and then returned to the greenhouse.
The phase change material (PCM) heating system is based on latent heat storage.The PCM absorbs solar energy during the day and releases that heat at night, keeping the greenhouse much warmer during cold winter nights.
A number of studies have described the use of these solar thermal storage systems for heating greenhouses and examined the relationship between the properties of the storage materials and their thermal performance [1][2][3][4][5][6].
Kooli et al. [7] experimentally tested a solar heating system with latent heat storage in a greenhouse in Tunisia.Following the same approach, Bazgaou et al. [8] analyzed the performance of a rock bed heating system in a conventional Canarian-type greenhouse in Morocco.For their side, Berroug et al. [9] investigated the impact of phase change materials (PCMs) on greenhouse temperature and humidity.
Guan et al. [10] experimentally studied the heat storage and release characteristics of a novel greenhouse wall with microheat pipe arrays (MHPAs) and phase change materials (PCMs).
The modeling approach has also been used as an effective tool to optimize solar thermal heating systems for greenhouses.In this vein, Ntinas et al. [11] developed a mathematical model to predict the thermal efficiency of a novel hybrid solar energy-saving system used to heat greenhouses.Mobtaker et al. [12] developed a dynamic model to predict the internal temperature in a solar greenhouse in northwestern Iran.
Recently, CFD models have been used to investigate the spatial and temporal distribution of climate parameters inside greenhouses equipped with heating systems in order to improve their performance [13,14].The advantage of using CFD models is that the climate in the greenhouse and heating system can be studied in more detail in contrast to energy balance models, where climatic variables in the greenhouse are considered to be uniform and homogeneous.
CFD modeling has been applied to several greenhouses equipped with different heating systems.Rouboa and Monteiro [15] simulated the effects on temperature and velocity fields by the introduction of hot water tubes into a greenhouse in night conditions.Chen et al. [16] proposed an optimized design of hot air piping in the greenhouse based on the CFD technique.Discrete ordinates (DO) and k-ε turbulence models were used.The discrete ordinates (DO) radiation model solves the radiative transfer equation (RTE) for a finite number of discrete solid angles.Mezrhab et al. [17] analyzed the effects of the heat exchange inside a greenhouse in winter according to the number of squared heating tubes used.Grigoriu et al. [18] proposed a CFD model to simulate the inside air temperature in a greenhouse equipped with a heating system based on the use of parabolic trough collectors (GHPTC).The proposed CFD model had two components: the greenhouse model and the collector's model.A proportional-integral-derivative (PID) controller was tuned and used to control the temperature and humidity inside the greenhouse based on the linearized model.This control method was suited for the system according to the obtained stationary and transient response performances.
For their part, Couto et al. [19] conducted a CFD analysis of the artificial heating tubes with natural ventilation.This CFD model was validated experimentally to obtain the temperature variation and air velocity inside the greenhouse.
Tadj et al. [20] conducted a CFD study to analyze the internal convective fluxes generated in a closed-tunnel greenhouse heated with a polyethylene blanket.
Morille et al. [21] performed a CFD study to predict the time evolution of the temperature and humidity distribution in the greenhouse, considering both radiative heat transfer and the interactions between plants with the indoor microclimate.Chen et al. [16] optimized the design of a hot air duct in a greenhouse through a CFD simulation of the greenhouse environment with a forced air heating system by the diameter of the hot air pipe (HAP) taking into account the shape, arrangement, and number of exhaust outlets in the pipe.
Despite the large volume of published studies on modeling heated greenhouses, the solar heating system based on sensible heat storage in a rock bed does not appear to have been studied much.Therefore, the aim of our work is to study the microclimate in greenhouses equipped with such a heating system using a CFD model in order to test the parameters affecting its performance.

Greenhouse Geometry and Description of the Heating System
The simulation study was designed and carried out in a Canarian-type greenhouse located in the Regional Center for Agronomic Research (CRRA) in Agadir, Morocco (latitude 30 • 13', longitude 9 • 23', altitude 80 m).The greenhouse area is 165 m 2 , i.e., 15 m long and 11 m wide and 4 m high at the gutter level and 5 m at the span (Figure 1).and used to control the temperature and humidity inside the greenhouse based on the linearized model.This control method was suited for the system according to the obtained stationary and transient response performances.
For their part, Couto et al. [19] conducted a CFD analysis of the artificial heating tubes with natural ventilation.This CFD model was validated experimentally to obtain the temperature variation and air velocity inside the greenhouse.Tadj et al. [20] conducted a CFD study to analyze the internal convective fluxes generated in a closed-tunnel greenhouse heated with a polyethylene blanket.
Morille et al. [21] performed a CFD study to predict the time evolution of the temperature and humidity distribution in the greenhouse, considering both radiative heat transfer and the interactions between plants with the indoor microclimate.Chen et al. [16] optimized the design of a hot air duct in a greenhouse through a CFD simulation of the greenhouse environment with a forced air heating system by the diameter of the hot air pipe (HAP) taking into account the shape, arrangement, and number of exhaust outlets in the pipe.
Despite the large volume of published studies on modeling heated greenhouses, the solar heating system based on sensible heat storage in a rock bed does not appear to have been studied much.Therefore, the aim of our work is to study the microclimate in greenhouses equipped with such a heating system using a CFD model in order to test the parameters affecting its performance.

Greenhouse Geometry and Description of the Heating System
The simulation study was designed and carried out in a Canarian-type greenhouse located in the Regional Center for Agronomic Research (CRRA) in Agadir, Morocco (latitude 30° 13', longitude 9° 23', altitude 80 m).The greenhouse area is 165 m 2 , i.e., 15 m long and 11 m wide and 4 m high at the gutter level and 5 m at the span (Figure 1).
The greenhouse was oriented perpendicular to the prevailing north-south wind direction.It was covered with polyethylene plastic with 75% light transmission and 200 µm thickness.The greenhouse was naturally ventilated by west and east roll-up side walls.The greenhouse was oriented perpendicular to the prevailing north-south wind direction.It was covered with polyethylene plastic with 75% light transmission and 200 µm thickness.The greenhouse was naturally ventilated by west and east roll-up side walls.
The greenhouse was equipped with a rock bed heating system.The system consisted of six blocks of rocks and 12 PVC pipes.Each block was connected at both ends with two pipes for supply and exhaust air.At the end of each pipe was a fan that was used to either draw in air or blow air out.The functional principle of the heating system is shown in Figure 2. The air is sucked through the air inlet pipes into the rock bed blocks, i.e., heat storage media, and passed through the spaces between the pebbles.Heat exchange takes place, and the pebbles store the heat contained in the air.By reversing the direction of the air flow during the night, the stored energy is released into the greenhouse via the air outlet pipes.
The greenhouse was equipped with a rock bed heating system.The system consisted of six blocks of rocks and 12 PVC pipes.Each block was connected at both ends with two pipes for supply and exhaust air.At the end of each pipe was a fan that was used to either draw in air or blow air out.The functional principle of the heating system is shown in Figure 2. The air is sucked through the air inlet pipes into the rock bed blocks, i.e., heat storage media, and passed through the spaces between the pebbles.Heat exchange takes place, and the pebbles store the heat contained in the air.By reversing the direction of the air flow during the night, the stored energy is released into the greenhouse via the air outlet pipes.

Experimental Measurements
To verify the accuracy and performance of the model, we performed experimental measurements of temperature and humidity in the greenhouse, which were used to compare experimental with simulated data from the numerical model.The temperature and relative humidity inside the greenhouse were measured using the HMP60 sensor at three heights in the greenhouse center: 1 m, 2 m, and 3 m from the ground.All measurements were sampled and recorded every 10 min using a CR3000 data logger (Campbell Scientific).For more details, see Bazgaou et al. [8].

CFD Governing Equations
Modeling (CFD) is widely used to study the behavior of mass, momentum, and heat transfer processes as well as to determine the distribution of air flows and climate param-

Experimental Measurements
To verify the accuracy and performance of the model, we performed experimental measurements of temperature and humidity in the greenhouse, which were used to compare experimental with simulated data from the numerical model.The temperature and relative humidity inside the greenhouse were measured using the HMP60 sensor at three heights in the greenhouse center: 1 m, 2 m, and 3 m from the ground.All measurements were sampled and recorded every 10 min using a CR3000 data logger (Campbell Scientific).For more details, see Bazgaou et al. [8].

CFD Governing Equations
Modeling (CFD) is widely used to study the behavior of mass, momentum, and heat transfer processes as well as to determine the distribution of air flows and climate parameters in a given geometry.
CFD modeling was used in this study to simulate the processes of heat transfer in the greenhouse, including solar radiation and heat convection, caused by natural ventilation and heating.This model's governing equations are as follows (Campen and Bot [22]): where ϕ is the transport concentration in dimensionless form, t (s) is the time, Γ m 2 /s is the diffusion coefficient, u j (m/s) represents each component of the velocity, X j is the coordinate in each direction, ρ kg/m 3 is the air density, and S ϕ is the generalized source term (per unit volume-time).
The Rayleigh number R a of a fluid in natural convection is given by Gholamalizadeh and Kim [20]: where µ is the dynamic viscosity, α is the thermal diffusivity, L is the characteristic length of the greenhouse, ∆T p is the temperature difference between the plants and the roof of the greenhouse, g is the gravity acceleration, and β is the coefficient of thermal expansion.
When the Rayleigh number is below 10 8 , the air flow is laminar.However, the air flow is turbulent when it is larger than 10 8 .When the heating system is operational, forced ventilation results in R a values greater than 10 10 , indicating that the air flow inside the greenhouse is turbulent.Therefore, a k − ε turbulence model is used to describe the flow in the heated greenhouse [23][24][25].
The k − ε turbulence model is the most common model used in CFD modeling to simulate airflow characteristics for turbulent flow conditions.This model is given in these two transport equations.
For turbulent kinetic energy k: For dissipation ε: where k is the turbulent kinetic energy, ε is the dissipation rate, σ k and σ ε are the Prandtl numbers of the turbulence kinetic energy and of the dissipation rate, respectively, P is the calculated pressure, and µ t is the turbulent viscosity.G B is the production of turbulent kinetic energy, which is defined as follows [26]: where T i is the inside air temperature, σ t is the turbulence Prandtl number, T ref is the reference temperature, β = 1 T ref is the voluminal coefficient of thermal expansion, and f 1 and f 2 are the auxiliary relations of viscous dissipation.The constants of the (k − ε) turbulence model are the following: 3, and σ t = 0.85.

Radiation Modeling
The discrete ordinates model (DO) is used to describe the radiative transfers inside the greenhouse.This model solves the radiative heat flow in a discrete number of directions in space.According to the multidirectional → s , the radiation intensity transfer equation can be solved.This model considers the radiative transfer equation in the → s direction as a field equation given in this formula [27]: where → r is the position vector, → s is the direction vector, I λ is the radiation intensity as a function of the position ( → r ) and of the direction ( → s ), ∇ is the divergence operator, → s is the diffusion direction vector, a λ is the absorption coefficient, σ is the Stefan-Boltzmann constant, σ s is the dispersion coefficient, n λ is the refractive index, T is the local temperature, Φ is the phase function, and Ω is the solid angle of the radiation.
Ansys Fluent 16 provides a solar loading model that can be used to calculate the effects of radiation from the sun entering a computational domain [25].In this work, the solar loading model ray tracing algorithm was used to predict the direct light energy source resulting from incidental solar radiation.The ray tracing approach is a very efficient and practical way to apply solar loads as a heat source term in the energy equation.Solar charge is only available in the three-dimensional solver and can be used to model stable and unstable fluxes.The resulting heat flux calculated with the solar ray tracing algorithm was coupled to the calculation through a source term in the energy equation.

Crop Porous Medium
The crop was assimilated to an equivalent homogeneous porous medium of a solid matrix with interconnected pores with a porosity of 70% [28].The source term for momentum in the basic transport equation was caused by the drag action of plants.This effect was achieved through the environmental parameters and the transformation relationship between the equivalent diameters of the crop volume and the porosity rate [29].Plant transpiration and sensitive heat exchanges depend on local climate and convective exchange, which control their aerodynamic and stomatic resistance.
As the crop was considered a porous medium, the Darcy-Forchheimer equation can be used to describe the dynamic effect of the plant on the airflow inside the medium.The source term in the momentum equation S ϕ1 includes the viscous loss and the internal resistance, and is expressed as follows [26]: where a is the permeability of the porous medium, Y is the nonlinear coefficient of momentum loss, and u is the air velocity.
During the heating process of the greenhouse, an exchange of mass and energy takes place between the indoor air and the crop.The energy exchange includes the sensible heat Q sen caused by heat exchange between the leaves and the surrounding air as well as the latent heat Q lat caused by transpiration.The source terms in the energy equation S ϕ2 and the sensible heat equation Q sen are given by Kichah et al. [30]: where L ai is the leaf area index, C a is the air specific heat at atmospheric pressure, T f and T i are the temperature of the plant leaves and that of the inside air, respectively, and r a is the aerodynamic resistance of the canopy.According to air velocity, the aerodynamic resistance of the canopy in the greenhouse studied by Boulard and Wang [28] is as follows: where d is the characteristic length of the leaf.
The stomatal resistance of the canopy r s is given by Boulard and Wang [28]: where R n is the net solar radiation, D i is the pressure difference of saturated water vapor, and u is the air velocity.
The vegetation cover transpiration rate model is based on the Penman-Monteith equation.According to Impron et al. [31], the transpiration rate of the canopy is expressed as follows: where e s is the vapor pressure of saturated water at T, e a is the actual air vapor pressure, γ is the psychrometric constant, and ∆ is the slope of the vapor pressure curve for saturated water at the temperature T ( • C).This slope is defined by Boulard and Wang [28] as follows:

CFD Numerical Model CFD Model Settings
The widely used commercial CFD code Fluent module from ANSYS [27] was used to simulate the microclimate in the greenhouse equipped with a rock bed heating system.This code discretizes space and time using the finite volume method together with the pressure-velocity coupling algorithm SIMPLE (semi-implicit method for pressure-related equations).To obtain accurate results, second-order upwind discretization schemes were used for momentum and turbulence equations.The convergence criterion for all variables was 1 × 10 −6 .
The CFD software was customized in the C++ coding language to account for energy and water vapor transfers that occur at the plant level.Practically, the canopy was divided into several layers of adjacent cells.Each cell consisted of a solid matrix (i.e., leaves and stems) and air.It was assumed that the cell size was small enough compared to the volume of the canopy that the temperature of the solid matrix could be considered homogeneous within a given cell.
The physical and thermal parameters of the various materials that constituted the greenhouse and the crop are listed in Table 1.For the ground around the rock bed blocks, we assumed that its temperature remained constant.

Mesh Details
The geometry of the greenhouse with the rock bed heating system, the rows of plants, and the surrounding domain was created using the Design Modeler tool under Ansys Fluent.An unstructured mesh was created with 1,014,811 elements and 224,659 nodes.A finer mesh was generated near the heating system and crop rows, and a relatively coarser mesh was generated in areas remote from these zones, as shown in Figure 3a-c.The mesh quality was verified, and it satisfied the skewness and orthogonality parameters.For the ground around the rock bed blocks, we assumed that its temperature re mained constant.

Mesh Details
The geometry of the greenhouse with the rock bed heating system, the rows of plants and the surrounding domain was created using the Design Modeler tool under Ansy Fluent.An unstructured mesh was created with 1,014,811 elements and 224,659 nodes.A finer mesh was generated near the heating system and crop rows, and a relatively coarse mesh was generated in areas remote from these zones, as shown in Figure 3a-c.The mes quality was verified, and it satisfied the skewness and orthogonality parameters.

Boundary Conditions
The boundary conditions used are listed in Table 2 and correspond to the climat parameters measured on 15 February 2017 at noon and midnight.The selected day corre sponds to a typical winter day when the greenhouse has to be heated in order to achiev optimal climatic conditions for the plants in the greenhouse [32].

Boundary Conditions
The boundary conditions used are listed in Table 2 and correspond to the climate parameters measured on 15 February 2017 at noon and midnight.The selected day corresponds to a typical winter day when the greenhouse has to be heated in order to achieve optimal climatic conditions for the plants in the greenhouse [32].
Table 2 illustrates the average values of the outdoor temperature, relative humidity, wind speed, and solar radiation applied as boundary conditions to the computational domain.

Results and Discussion
This section is divided by subheadings.There should be a concise and accurate description of the experimental results, their interpretation, and the experimental conclusions that can be drawn.

Model Validation
The accuracy of the developed model was evaluated by comparing model simulations with the experimental data of air temperature and humidity inside the greenhouse.
Figures 4 and 5 show simulated and measured air temperature values at three different heights in the greenhouse center at noon (in Figure 4) and midnight (in Figure 5).This section is divided by subheadings.There should be a concise and accurate de scription of the experimental results, their interpretation, and the experimental conclu sions that can be drawn.

Model Validation
The accuracy of the developed model was evaluated by comparing model simula tions with the experimental data of air temperature and humidity inside the greenhouse Figures 4 and 5 show simulated and measured air temperature values at three differ ent heights in the greenhouse center at noon (in Figure 4) and midnight (in Figure 5).
When these values are compared, they show that the simulated and calculated values agree well.The average difference between observed and simulated air temperature is 0.65 °C during the day, and 0.25 °C at night.When these values are compared, they show that the simulated and calculated values agree well.The average difference between observed and simulated air temperature is 0.65 • C during the day, and 0.25 • C at night.To further confirm the validity of the model, we also compared the air humidity data from the measurements and simulations.Figure 6 shows the simulated and measured values of the inside relative air humidity variations as a function of greenhouse height at noon.Accordingly, relative humidity correlates inversely with air temperature.It decreases from the bottom to the top of the greenhouse.The measured and simulated values are almost identical with minor differences of 0.42% to 1.59%.These results show that there is no significant difference between the experimental measurements and those from the simulation, thus giving credibility to the CFD model for simulating the thermal performance of a rock bed heating system and confirming the accuracy of the simulation results.To further confirm the validity of the model, we also compared the air humidity data from the measurements and simulations.Figure 6 shows the simulated and measured values of the inside relative air humidity variations as a function of greenhouse height at noon.Accordingly, relative humidity correlates inversely with air temperature.It decreases from the bottom to the top of the greenhouse.The measured and simulated values are almost identical with minor differences of 0.42% to 1.59%.To further confirm the validity of the model, we also compared the air humidity data from the measurements and simulations.Figure 6 shows the simulated and measured values of the inside relative air humidity variations as a function of greenhouse height at noon.Accordingly, relative humidity correlates inversely with air temperature.It decreases from the bottom to the top of the greenhouse.The measured and simulated values are almost identical with minor differences of 0.42% to 1.59%.These results show that there is no significant difference between the experimental measurements and those from the simulation, thus giving credibility to the CFD model for simulating the thermal performance of a rock bed heating system and confirming the accuracy of the simulation results.These results show that there is no significant difference between the experimental measurements and those from the simulation, thus giving credibility to the CFD model for simulating the thermal performance of a rock bed heating system and confirming the accuracy of the simulation results.

Solar Radiation
Using this model with boundary conditions corresponding to noon allowed us to simulate the solar radiation received on the sidewalls and roof of the greenhouse.The results are shown in Figure 7. Accordingly, the solar flux on the roof was greater than on the side walls.Its value was around 850 W/m 2 on the roof and from 708 to 786 W/m 2 from bottom to top on the side walls.

Solar Radiation
Using this model with boundary conditions corresponding to noon allowed us to simulate the solar radiation received on the sidewalls and roof of the greenhouse.The results are shown in Figure 7. Accordingly, the solar flux on the roof was greater than on the side walls.Its value was around 850 W/m 2 on the roof and from 708 to 786 W/m 2 from bottom to top on the side walls.Radiation received inside the greenhouse is important for the thermal heating process.It warms the air in the greenhouse by turning it into heat.The excess heat is fed to the storage medium through the pipes from inside the greenhouse during the day and released again at night.

Solar Radiation
Using this model with boundary conditions corresponding to noon allowed us to simulate the solar radiation received on the sidewalls and roof of the greenhouse.The results are shown in Figure 7. Accordingly, the solar flux on the roof was greater than on the side walls.Its value was around 850 W/m 2 on the roof and from 708 to 786 W/m 2 from bottom to top on the side walls.Radiation received inside the greenhouse is important for the thermal heating process.It warms the air in the greenhouse by turning it into heat.The excess heat is fed to the storage medium through the pipes from inside the greenhouse during the day and released again at night.Radiation received inside the greenhouse is important for the thermal heating process.It warms the air in the greenhouse by turning it into heat.The excess heat is fed to the storage medium through the pipes from inside the greenhouse during the day and released again at night.
In order to improve the efficiency of the thermal heating system, the greenhouse cover must be as transparent as possible for solar radiation to reach almost 100% efficiency.

Air Temperature Profile
The developed model was also used to simulate the air temperature inside the greenhouse due to its importance in plant photosynthesis and metabolic reactions as well as in water and mineral intake and transpiration.
Figure 9 shows the air temperature distribution in the frontal and longitudinal planes drawn in the center of the greenhouse.According to this figure, the air temperature varied from 302 K as the minimum value in the lower part of the greenhouse to 313 K in the upper part of the greenhouse.This temperature gradient can be explained by the difference in density between warm and cold air: hot air rises from the bottom to the top, forming an air pocket the rows of plants.The relatively lower temperature at ground level was due to the intrusion of fresh cooling air from the rock bed system during the storage phase.
The developed model was also used to simulate the air temperature inside the greenhouse due to its importance in plant photosynthesis and metabolic reactions as well as in water and mineral intake and transpiration.
Figure 9 shows the air temperature distribution in the frontal and longitudinal planes drawn in the center of the greenhouse.According to this figure, the air temperature varied from 302 K as the minimum value in the lower part of the greenhouse to 313 K in the upper part of the greenhouse.This temperature gradient can be explained by the difference in density between warm and cold air: hot air rises from the bottom to the top, forming an air pocket over the rows of plants.The relatively lower temperature at ground level was due to the intrusion of fresh cooling air from the rock bed system during the storage phase.
This thermal stratification is common to all greenhouses.Warm and humid air rises to the top of the greenhouse, whereas cool air remains trapped in the lower part at plantlevel.Saberian et al. [33] observed the same phenomenon in a gabled greenhouse covered with semitransparent materials.
This temperature stratification on the vertical plane, with hot air rising towards the top of the greenhouse, explains the positioning of the air intake ducts at the level of the greenhouse vault in order to better recover heat from the warm air during the day.The thermal energy is drawn through the pipes and stored in the rock bed.The analysis of Figure 10 clearly shows this heat storage.This figure illustrates the temperature gradient in the rock bed blocks during the thermal storage phase.This gradient is symmetrical on both sides of the rock blocks.The supplied air was warmer than the extract air, with input and output temperatures of 34 °C and 17 °C, respectively.This temperature stratification on the vertical plane, with hot air rising towards the top of the greenhouse, explains the positioning of the air intake ducts at the level of the greenhouse vault in order to better recover heat from the warm air during the day.The thermal energy is drawn through the pipes and stored in the rock bed.The analysis of Figure 10 clearly shows this heat storage.This figure illustrates the temperature gradient in the rock bed blocks during the thermal storage phase.This gradient is symmetrical on both sides of the rock blocks.The supplied air was warmer than the extract air, with input and output temperatures of 34

Relative Air Humidity Profile
The contours of relative humidity in vertical and longitudinal sections in the center of the greenhouse at noon are shown in Figure 12.This figure shows that the relative humidity correlated inversely with the temperature inside the greenhouse.The highest humidity values were observed in the lower part of the greenhouse, and the lowest values were recorded in the upper part of the greenhouse.These values were between 20% and 36% from the top to the bottom of the greenhouse.The observed high values of humidity in the crop rows can be explained by the transpiration of the plants during the day.The heat release phase at night is visualized in Figure 11.This figure shows the crosssection of the simulated air temperature gradient at night in the vertical median plane of the greenhouse and the rock bed row.The temperature at the extremities of the rock bed blocks dropped to 10.5 °C due to the cold air ingested by the inlet fans.The outlet fans in the center released the stored heat, raising the greenhouse's internal air temperature to 13.5 °C.

Relative Air Humidity Profile
The contours of relative humidity in vertical and longitudinal sections in the center of the greenhouse at noon are shown in Figure 12.This figure shows that the relative humidity correlated inversely with the temperature inside the greenhouse.The highest humidity values were observed in the lower part of the greenhouse, and the lowest values were recorded in the upper part of the greenhouse.These values were between 20% and 36% from the top to the bottom of the greenhouse.The observed high values of humidity in the crop rows can be explained by the transpiration of the plants during the day.

Relative Air Humidity Profile
The contours of relative humidity in vertical and longitudinal sections in the center of the greenhouse at noon are shown in Figure 12.This figure shows that the relative humidity correlated inversely with the temperature inside the greenhouse.The highest humidity values were observed in the lower part of the greenhouse, and the lowest values were recorded in the upper part of the greenhouse.These values were between 20% and 36% from the top to the bottom of the greenhouse.The observed high values of humidity in the crop rows can be explained by the transpiration of the plants during the day.
This finding was also reported by Kim et al. [34], who have shown that the evolution of relative humidity in the greenhouse was strongly related to temperature variation.
The humidity stratification observed in the greenhouse and particularly in the plant canopies was mainly due to the moisture generated by transpiration.The air movement caused by the heating process of the rock bed can play an important role in controlling the humidity in the growing area.It can also help prevent condensation on the surface of the leaves.Moving air keeps the leaves dry, which helps prevent the development of cryptogamic diseases, especially at night when the temperature drops.This finding was also reported by Kim et al. [34], who have shown that the evolution of relative humidity in the greenhouse was strongly related to temperature variation.

Air Velocity Profile Inside the Greenhouse
The humidity stratification observed in the greenhouse and particularly in the plant canopies was mainly due to the moisture generated by transpiration.The air movement caused by the heating process of the rock bed can play an important role in controlling the humidity in the growing area.It can also help prevent condensation on the surface of the leaves.Moving air keeps the leaves dry, which helps prevent the development of cryptogamic diseases, especially at night when the temperature drops.During the day, warm air from the greenhouse moves through the rock bed during this storage phase.This passage causes heat to be transferred from the air to the rocks and stored for later use.Remember that this air circulation is forced convection caused by fans.

Air Velocity Profile Inside the Greenhouse
The inlet pipes draw in hot air at a velocity of 4.56 m/s from the top of the greenhouse and circulate it through the rock bed blocks.At a velocity of 4.21 m/s, the outlet pipes release cool air into the greenhouse.It should be noted that air flows slowly through the porous medium formed by the rock bed, allowing the rocks to increase their heat storage capacity and thermal efficiency.
Transpiration creates higher humidity near plants than in the rest of the greenhouse.The air flow provided by the outlet pipes allows for the renewal of air near the plant leaves as well as the removal of humidity, reducing the likelihood of free water condensing on the leaves.Knowing that the presence of this free water on the leaves favors the development of cryptogamic diseases, the movement of air over the rock bed blocks would have a positive impact on plant health and, as a result, would reduce the use of chemical products.
The air circulation generated by air inlets and outlets through the rock bed also allows for uniform air temperature and the breakdown of air stratification at the crop level.This temperature uniformity is highly desired by growers who want to unify crop production and avoid heat accumulation at the top of the greenhouse at all costs.Heat losses are mainly caused by thermal stratification.The use of the rock bed heating system could help save energy in protected agriculture.

Conclusions
In this study, the microclimate in a greenhouse with a rock bed thermal storage heating system was simulated using a CFD model.The heat and mass exchanges between the different elements of the greenhouse, the rock bed blocks, and the surrounding environment have been considered, making it possible to reproduce the dynamic, thermal, and hygric fields in the Canarian greenhouse equipped with the rock bed heating system.
This model was successfully validated by comparing simulation results to experimental data.The comparison of measured and simulated temperature and relative humidity values in the greenhouse demonstrates the model's good accuracy and gives confidence in its use as a tool to optimize the system and improve its performance.Using this model, it is thus possible to conduct a global sensitivity study to test several adjustments that would increase the efficiency of the rock bed heating system.This sensitivity study will be based on the variation of (i) the parameters related to the thermal properties of storage media namely density, thermal conductivity, and specific heat; (ii) the porosity of the materials used as the energy storage medium; (iii) airflow rates to increase heat discharge by adjusting fan speed.
The model developed in this study is a useful tool for greenhouse climate simulation as well as a decision support system for solar greenhouse management.Furthermore, this model can be used to predict energy savings in other types of greenhouses and climate zones when using rock bed heating systems.
Continued modeling efforts are required to optimize the design parameters of the heating system in order to improve its thermal performance.

Figure 1 .
Figure 1.Greenhouse geometry and rock bed heating system.Three-dimensional coordinates represent x as east, y as zenith, and z as south.

Figure 1 .
Figure 1.Greenhouse geometry and rock bed heating system.Three-dimensional coordinates represent x as east, y as zenith, and z as south.

Figure 2 .
Figure 2. The operating principle of the rock bed heating system.

Figure 2 .
Figure 2. The operating principle of the rock bed heating system.

Figure 3 .
Figure 3. Mesh used to simulate the calculation domain (a), the greenhouse equipped with the roc bed heating system (b), and the crop (c).

Figure 3 .
Figure 3. Mesh used to simulate the calculation domain (a), the greenhouse equipped with the rock bed heating system (b), and the crop (c).

Figure 4 .
Figure 4. Simulated and measured air temperatures at three different heights (1 m, 2 m and 3 m) a noon on 15 February 2017.

Figure 4 .
Figure 4. Simulated and measured air temperatures at three different heights (1 m, 2 m and 3 m) at noon on 15 February 2017.

Figure 5 .
Figure 5. Evolution of simulated and measured air temperatures inside the greenhouse as a function of height (1 m, 2 m and 3 m) at midnight.

Figure 6 .
Figure 6.Evolution of the simulated and measured values of the relative humidity as a function of the greenhouse height (1 m, 2 m, and 3 m) at noon.

Figure 5 .
Figure 5. Evolution of simulated and measured air temperatures inside the greenhouse as a function of height (1 m, 2 m and 3 m) at midnight.

17 Figure 5 .
Figure 5. Evolution of simulated and measured air temperatures inside the greenhouse as a function of height (1 m, 2 m and 3 m) at midnight.

Figure 6 .
Figure 6.Evolution of the simulated and measured values of the relative humidity as a function of the greenhouse height (1 m, 2 m, and 3 m) at noon.

Figure 6 .
Figure 6.Evolution of the simulated and measured values of the relative humidity as a function of the greenhouse height (1 m, 2 m, and 3 m) at noon.

Figure 7 .
Figure 7. Simulated solar flux received on the side walls and roof of the greenhouse at solar noon.Part of this solar flux received by greenhouse walls was transmitted inward and received by the internal components (air, plants, soil, and internal structure).The solar flux received on the ground was approximately 630 W/m 2 , as shown in Figure 8.This value corresponds to approximately 75% of the global solar flux received by the greenhouse, which corresponds to the transmittance of the plastic cover.

Figure 8 .
Figure 8. Simulated solar flux received at the indoor ground level of the greenhouse at solar noon.

Figure 7 .
Figure 7. Simulated solar flux received on the side walls and roof of the greenhouse at solar noon.Part of this solar flux received by greenhouse walls was transmitted inward and received by the internal components (air, plants, soil, and internal structure).The solar flux received on the ground was approximately 630 W/m 2 , as shown in Figure 8.This value corresponds to approximately 75% of the global solar flux received by the greenhouse, which corresponds to the transmittance of the plastic cover.

Figure 7 .
Figure 7. Simulated solar flux received on the side walls and roof of the greenhouse at solar noon.Part of this solar flux received by greenhouse walls was transmitted inward and received by the internal components (air, plants, soil, and internal structure).The solar flux received on the ground was approximately 630 W/m 2 , as shown in Figure 8.This value corresponds to approximately 75% of the global solar flux received by the greenhouse, which corresponds to the transmittance of the plastic cover.

Figure 8 .
Figure 8. Simulated solar flux received at the indoor ground level of the greenhouse at solar noon.

Figure 8 .
Figure 8. Simulated solar flux received at the indoor ground level of the greenhouse at solar noon.

Figure 9 .
Figure 9. Temperature distribution (K) in the center of the greenhouse at noon during 15 February 2017: frontal (a) and longitudinal (b) view.

Figure 9 .
Figure 9. Temperature distribution (K) in the center of the greenhouse at noon during 15 February 2017: frontal (a) and longitudinal (b) view.This thermal stratification is common to all greenhouses.Warm and humid air rises to the top of the greenhouse, whereas cool air remains trapped in the lower part at plant-level.Saberian et al. [33] observed the same phenomenon in a gabled greenhouse covered with semitransparent materials.This temperature stratification on the vertical plane, with hot air rising towards the top of the greenhouse, explains the positioning of the air intake ducts at the level of the greenhouse vault in order to better recover heat from the warm air during the day.The thermal energy is drawn through the pipes and stored in the rock bed.The analysis of Figure10clearly shows this heat storage.This figure illustrates the temperature gradient in the rock bed blocks during the thermal storage phase.This gradient is symmetrical on both sides of the rock blocks.The supplied air was warmer than the extract air, with input and output temperatures of 34 • C and 17 • C, respectively.

Figure 10 .
Figure 10.Temperature distribution (K) in the blocks of the rock bed heating system during the storage phase at noon (273 K = 0 °C).

Figure 11 .
Figure 11.Air temperature distribution (in K) in the greenhouse and the rock bed blocks during the destocking phase at midnight (273 K = 0 °C).

Figure 10 . 17 Figure 10 .
Figure 10.Temperature distribution (K) in the blocks of the rock bed heating system during the storage phase at noon (273 K = 0 • C).The heat release phase at night is visualized in Figure 11.This figure shows the crosssection of the simulated air temperature gradient at night in the vertical median plane of the greenhouse and the rock bed row.The temperature at the extremities of the rock bed blocks dropped to 10.5 • C due to the cold air ingested by the inlet fans.The outlet fans in the center released the stored heat, raising the greenhouse's internal air temperature to 13.5 • C.

Figure 11 .
Figure 11.Air temperature distribution (in K) in the greenhouse and the rock bed blocks during the destocking phase at midnight (273 K = 0 °C).

Figure 11 .
Figure 11.Air temperature distribution (in K) in the greenhouse and the rock bed blocks during the destocking phase at midnight (273 K = 0 • C).

Figure 12 .
Figure 12.Relative air humidity contour (%) at the center of the greenhouse: vertical section (a) and horizontal section (b) at noon.

Figure 13
Figure 13 depicts the simulation of air circulation within the greenhouse and in the rock bed blocks.The general air pathlines in the greenhouse and in the heat storage medium are shown in this figure.

Figure 12 .
Figure 12.Relative air humidity contour (%) at the center of the greenhouse: vertical section (a) and horizontal section (b) at noon.

Figure 13
Figure 13  depicts the simulation of air circulation within the greenhouse and in the rock bed blocks.The general air pathlines in the greenhouse and in the heat storage medium are shown in this figure.During the day, warm air from the greenhouse moves through the rock bed during this storage phase.This passage causes heat to be transferred from the air to the rocks and stored for later use.Remember that this air circulation is forced convection caused by fans.The inlet pipes draw in hot air at a velocity of 4.56 m/s from the top of the greenhouse and circulate it through the rock bed blocks.At a velocity of 4.21 m/s, the outlet pipes release cool air into the greenhouse.It should be noted that air flows slowly through the porous medium formed by the rock bed, allowing the rocks to increase their heat storage capacity and thermal efficiency.Transpiration creates higher humidity near plants than in the rest of the greenhouse.The air flow provided by the outlet pipes allows for the renewal of air near the plant leaves as well as the removal of humidity, reducing the likelihood of free water condensing on the leaves.

Figure 12 .
Figure 12.Relative air humidity contour (%) at the center of the greenhouse: vertical section (a) and horizontal section (b) at noon.

3. 5 .
Figure 13 depicts the simulation of air circulation within the greenhouse and in the rock bed blocks.The general air pathlines in the greenhouse and in the heat storage medium are shown in this figure.

Figure 13 .
Figure 13.Pathlines of air circulation and its velocity (m/s) at the greenhouse center and the rock bed of the longitudinal plane during the day.

Table 1 .
Physical and thermal parameters of the different materials that constituted the greenhouse and the crop.

Table 2 .
Used boundary conditions.

Table 2 .
Used boundary conditions.