Abstract
In order to optimize the energy requirements (heating/cooling) of a multi-zone greenhouse, and investigate its heat recovery potential, a mathematical, dynamic energy model, coded in the MATLAB/Simulink platform, is developed. For validation, a case study in cold climate conditions is evaluated. This dynamic model, based on both energy and water vapor mass balances, was able to calculate the year-round monthly energy demand for the case study. The model calculations were compared with actual energy consumption data and were shown to have an accuracy between 6% and 15.5% for different months. The results highlighted the potential of applying a heat recovery strategy, whether with a Phase Change Material (PCM) or a Heat Recovery Ventilator (HRV). It is shown that using a HRV can reduce the energy demand of the greenhouse by 5% for January and 4% for December. Regarding the greenhouse radiation performance, the south roof contributes the most to solar heat gain in winter and summer, while the north wall makes the minimum contribution. Consequently, it is proposed to increase the area of the south roof and insulate the north wall. Thus, an asymmetrical roof configuration can receive 6% more solar radiation. Calculations show that an east-west greenhouse orientation lowers energy demand by 3%.
Introduction
Currently, most countries are suffering from a huge shortage in sustainable food and energy resources. This situation is due in part to COVID-19, extreme weather conditions, and the ravages of war. Greenhouses are the most effective solution for the ever-increasing demand for sustainable agriculture (Ahamed et al., 2019). However, greenhouse energy consumption, which can be significantly higher than that of conventional cultivation, can affect the benefits of industrial farming (Barbosa et al., 2015). Beshada et al. (2006) focused on the energy consumption of greenhouses, and confirmed the vast amount of energy that they onsum. Generally, greenhouses are considered as a solar collector system, which means that by optimizing their properties (dimensions, orientation, materials, etc.) a net zero energy greenhouse in achievable (Vadiee and Martin, 2014; van ’t Ooster et al., 2007). The design parameters of a greenhouse are highly dependent on climate conditions. This topic was extensively reviewed by Ghani et al. (2019). Zaragoza et al. (2007) proposed a closed concept solar greenhouse that in addition to reducing energy consumption, has advantages such as lower operating cost, recovery of low-quality water (gray water), and lower air pollution. To have the highest productivity in arid regions such as the Middle East or Southern Europe, where farmers face obstacles such as extreme air temperature, insufficient water supply resources, and poor soil quality, a seawater greenhouse system (Al-Ismaili and Jayasuriya, 2016) has been proposed. Based on a thermal storage concept, the Chinese style greenhouse configuration (Gao et al., 2010; Tong et al., 2013) allows removal or reduction of the external heating utility.
To find potential solutions for reducing the energy consumption of a solar greenhouse in northern regions, developing a mathematical energy model is a key step. However, developing this model could be a very sophisticated task. To calculate the energy demand, the contribution of the most important factors must be taken into the account, including the type of plants (tomato, potato, roses, poinsettia, etc.), soil properties, interior air conditions (temperature set point, relative humidity, etc.), greenhouse structure, outdoor climate (temperature, wind speed, solar radiation, etc.), and greenhouse orientation. In fact, the number of considered factors largely determines to what extent a model is accurate and functional. Such a mathematical model with numerous parameters and their dynamic behavior will undoubtedly be very complex. The principal purpose of this model is to calculate the heat and mass transfer rates within a greenhouse. Having an accurate energy model enables energy analysts and researchers to propose appropriate strategies to target the minimum energy consumption and operating cost.
Early studies regarding the development of energy models for greenhouses date back to the 1980s when, to calculate the energy demand of a greenhouse, Bot (1988) utilized an energy model that considered the greenhouse cover, air, crop canopy, and four layers of soil. Later, this study was improved by Takakura (1993). In another study, with the objective of optimizing the net profit of a tomato greenhouse, Trigui et al. (2001) utilized a dynamic energy model to compute the heating/cooling load of a greenhouse. They proposed a multi-step approach, in which the first step was defining the profit. They defined the profit as the crop yield value less the energy costs, for heating and dehumidification, and the CO2 injection cost. Then, in the second step, they proposed an algorithm with a Hamiltonian function to solve the objective function J, which maximizes the net profit during harvest time.
De Zwart (1996) and later Luo et al. (2005) proposed a mathematical model based on an energy and mass balance which was capable of predicting the energy demand, humidity, temperature, and CO2 concentration within a greenhouse. Also, Golzar et al. (2018) proposed an integrated mathematical model that calculates both the energy demand of a greenhouse as well as the crop growth, including models for both. They studied and accounted for the interaction between energy consumption and crop yield. Also, they reviewed 30 different energy models for greenhouses. Esmaeli and Roshandel (2020) proposed a new energy model, based on climate conditions of Iran, and optimized a greenhouse where cucumbers, melons, and watermelons are cultivated. However, they did not consider the effect of the type of product on energy demand. Chen et al. (2015) used hybrid particle swarm optimization and a genetic algorithm (HPSO-GA) technique to optimize the energy demand of a greenhouse. They concluded that PSO would bring more accurate and reliable results for dynamic parameters. The dynamic model developed by Mohammadi et al. (2018) was capable of predicting the inside air and soil temperature.
There are two other computer-based simulation approaches that have been applied by researchers in their efforts to predict the thermal load of greenhouses. The first one is Building Energy Modeling (BEM) commercial software like EnergyPlus and more importantly TRNSYS. However, their ability to predict the greenhouse’s parameters is questionable. Commercial software mostly suffers from its inability to consider specific greenhouse phenomena such as evaporation and plant transpiration (evapotranspiration), which barely occur in residential or commercial buildings. For example, water moisture production in a greenhouse is a function of solar radiation, plant type, total plant floor coverage, plant growth, watering schedule, etc. Also, most practical greenhouses are equipped with thermal and blackout screens. Building software packages are unable to accurately simulate these features or account for several of the named greenhouse parameters. Another issue with commercial software is that the user is very limited in the ability to manipulate the software’s mathematical formulas. Having said that, while commercial software products can consider multi-zone features, the mentioned limitations make them an unreliable tool for greenhouse modeling and design. Ahamed et al. (2020) used TRNSYS software to obtain the energy demand of a greenhouse. However, his study revealed that this software is completely unable to produce accurate data. As another example, Chen et al. (2018) used EnergyPlus software to optimize a Chinese style greenhouse’s north wall to have maximum solar radiation in order to target minimum energy demand. The second simulation approach, which has significant capability but is very case dependent, is computational fluid dynamics (CFD). Using CFD, Zhang et al. (2016) obtained the temperature and moisture distribution within a Chinese style greenhouse. They used an ANSYS package to solve the 2-D Navier-Stokes equations with dynamic boundary conditions to obtain operational data within a greenhouse.
Previous studies considered the greenhouse as a single micro-climate, whereas in practice greenhouses are typically subdivided into several zones, each with their own temperature and humidity setpoints. Hereafter, this kind of greenhouse will be referred to as “multi-zone.” Multi-zone greenhouses can potentially provide an important opportunity for energy saving and heat recovery. In this study, an in-house MATLAB-Simulink code, which was previously validated by this research team (Ghaderi et al., 2023), is selected as the starting point. The earlier code was only able to predict the heat load of a single micro-climate zone greenhouse, and did not consider the effect of artificial lighting, thermal screens, and fans. Hence, the current study has two main objectives. The first objective it to develop and present a mathematical energy model for calculating the energy demand (heating/cooling load as well as dehumidification load) of a multi-zone greenhouse in a northern growing region. This model will consider thermal screens, artificial lighting, and equipment loads. To the best knowledge of the authors, no previous studies presented an energy model that calculates the energy demand of a multiple micro-climate greenhouse. By developing this energy model, the opportunities for taking advantage of heat recovery or heat storage in different zones can be investigated. The second objective is to propose the best physical properties of a greenhouse in cold climate agricultural regions, such as Quebec, Canada.
Model development
In the following sections, after introducing the governing equations and the properties of a chosen case study, the model validation results will be presented. Then, in the results and discussion chapter, the heat recovery potential and solar radiation factors will be investigated.
In this section the formulas that cover the energy and mass conservation in a multi-zone greenhouse are presented. A lumped model was found to be sufficiently accurate to capture the transient behavior of a single-zone greenhouse micro-climate (Van Beveren et al., 2015). A lumped energy model is a simplified representation of a complex physical system, where the system’s behavior is approximated using discrete components characterized by parameters such as mass, capacitance, and resistance. In this model, the system is divided into distinct elements, and interactions between them are defined by simple mathematical equations. The lumped energy model is particularly valuable in engineering and physics for analyzing and predicting the behavior of systems, such as electrical circuits or thermal systems, by focusing on energy storage and dissipation while disregarding spatial variations. Though it offers a streamlined approach to system analysis, it may not account for all real-world complexities in highly distributed or intricately interacting systems.
Energy balance
The overall greenhouse heat balance equation is presented in equation (1). Figure 1 schematically illustrates the different mechanisms that can affect the greenhouse energy consumption. The formulation of each mechanism is described in section 2.1 and 2.2.

Different heat loss/gain mechanisms in a conventional greenhouse.
In the above equation, Tin(t) is the time dependent temperature of indoor air, C
a
is the heat capacity of air,
In a greenhouse, the main source of heat gain comes from solar radiation. In other word, the greenhouse can be considered as solar collector. The solar radiation
A
i
is the total surface area of the greenhouse cover (m2), I
i
(t) is total solar irradiance
After some simplifications, the total amount of solar radiation on a sloped surface
In equation (3),
In equations (4) and (5), ω,

(a) Schematic diagram of different angles that affect solar radiation and (b) solar azimuth angle.
Heat transfer rates
where ηth is the thermal ratio of lamps,
Heat loss due to convection and conduction
In equations (8) and (9),
Since the air velocity inside greenhouses is negligible (less than 0.2
In equations (11)–(13),
Usually, greenhouses have a non-transparent wall portion on the east, south and west sides, extending upwards from ground level for 60 cm. To calculate the non-transparent vertical surfaces, the indoor convection coefficient can be obtained by equation 14.
Since the main heat transfer mechanism on the outside walls is forced convection, the convective heat transfer coefficient can be calculated by equation 15, noting that it mainly depends on the Reynolds and Prandtl numbers.
It should be noted that to obtain the Re value, a length-based formulation must be considered (equation 16).
In equation 16, ρ is air density,
One of the most challenging tasks in developing a mathematical energy model for a greenhouse is modeling mechanical ventilation in open-concept greenhouses. In this study we only considered the natural ventilation, which is responsible of up to 20% of heat lost in a greenhouse (Ahamed et al., 2018).
In the above equation,
Greenhouses have different mechanisms to meet the greenhouse thermal load in the winter season. One of these mechanisms is a floor heating system. The greenhouse floor can lose energy in two ways. First, heat conduction to the ground (equation (18)) and second, perimeter heat loss (equation (19)).
Thermal radiation ql(t) from the inside to the outside of the greenhouse, which is also known as long-wave radiation heat loss, is one of the most important heat losses that affects the greenhouse micro-climate. This term is computed by equation (20):
Where σ is the Stefan–Boltzmann constant (
In cold regions, greenhouses are equipped with thermal screens to control the mentioned term
Plant transpiration affects both the temperature and water content of the greenhouse indoor air. The model that is used in this work is based on the transpiration model of Van Beveren et al. (2015) and is described in equation (23):
Where
In equation (24), LAI is the Leaf Area Index (m2 leaf divided by m2 ground), ε is the ratio of latent to sensible heat content of saturated air,
Where χair,sat(t) is the indoor saturated vapor concentration (
Stomatal resistance rs(t) is an index showing the resistance to transpiration of water vapor through the stomata (leaf pores) on the plant leaves and is written follows by a semi-empirical model:
Where γ is a crop specific parameter (e.g. its value is 0.027 for tomato and 0.008 for rose (Van Beveren et al., 2015)). Also, in equation (28),
Where S v is the opening ventilation area in m2, U is the external wind speed, C is the wind pressure coefficient, and C d is the vent discharge coefficient (are assumed to be 0.7 and 0.11 respectively). For the mechanical ventilation and dehumidifier terms the following equations are used:
The values for
Before concluding the energy balance section, a few comments concerning the application of equations (1)−(39) for multi-zone greenhouses are in order. To prepare equation (7) for a multi-zone greenhouse, some adjustments are necessary. Like commercial software, for walls shared between adjacent zones, the boundary conditions need some modifications. As we will see in the case study, Figure 3, there are three temperature zones that will be considered. To clarify this concept, consider zone 3, which has one shared wall on the north side with zone 2. The north wall of zone 3 (south wall of zone 2) is not exposed to outside temperature, but instead is exposed to the thermal conditions of the adjacent portion of the greenhouse (zone 2). Also, the dominant heat transfer mechanism for the shared wall in multi-zone formulation is conduction, whereas it is convection for a single zone, where all walls are exterior walls. Hence, to calculate the value for conductive heat transfer on the north wall of zone 3, the air temperature of zone 2 should replace the outside temperature. A similar approach is necessary for solving the mass transfer equations. Although these modifications increase the code complexity, they will help increase the code accuracy. The necessity of this modification is greatest during the winter season when the outdoor temperature is negative.

Schematic of case study greenhouse, dimensions, location of fans, and orientation.
Vapor mass balance
The humidity balance for a greenhouse is described by:
E(t), C(t), D(t) and V(t) denote evapotranspiration (evaporation from plant’s leaves), condensation (on the greenhouse inside surfaces), dehumidification and ventilation respectively (per unit area of the greenhouse). It should be noted that since the relative humidity will be considered constant at 75%, the left-hand side of equation (32) is equal to 0.
Where
By using equations (1)−(39) and applying hourly climate data as well as the indoor humidity and temperature set points, the amount of hourly energy demand that is needed to maintain the desired temperature and relative humidity in the greenhouse is calculated. It is worthy to highlighted the point that since the selected greenhouse for this study is only using natural ventilation through roof ventilators, the values of
Case study
To validate the energy model accuracy, a flower greenhouse in the city of Sherbrooke, Quebec, Canada (45.4043° N, 71.8937° W) has been selected. This greenhouse, named “Les Serres St-Élie,” has three different micro-climate zones and a warehouse to support logistical requirements. Also, thermal screens, blackout screens and artificial lighting are employed in the greenhouse to reduce energy consumption and increase productivity. Each of the three zones has a different temperature and relative humidity setpoint for different months. Overall greenhouse dimensions and micro-climate conditions of the St. Elie greenhouse are illustrated in Tables 1 to 3, and Figure 3. These three zones are separated by vertical roll-up curtains. In this way, operators can merge or isolate the zones, as needed. To meet the energy demand, two natural gas boilers (Viessman CA3V- 3.0) with an output power of up to 879 kW each are installed in the greenhouse. To reduce temperature stratification and increase the air recirculation, 18 ceiling fans (model ECOfan+ CSA5400) are installed throughout the greenhouse (see Figure 3). To secure the maximum amount of productivity, artificial lighting equipment and horizontal blackout screens are installed in zone 1. All three zones have horizontal thermal screens. Table 2 lists the information concerning the artificial lighting for different months of the year.
Thermal properties of materials used in St. Elie greenhouse.
Interior micro-climate conditions of the St. Elie greenhouse.
Thermal properties of greenhouse material.
Computer code description and validation
The first part of the modeling task was accomplished using custom MATLAB-SIMULINK code. The inputs to this code are physical properties (dimensions, thermal properties, orientation, etc.) and interior micro-climate conditions (temperature and relative humidity set points, artificial lighting and thermal screen schedules, etc.) of the St. Elie greenhouse as design parameters, and Sherbrooke meteorological data (outside temperature, relative humidity, global horizontal irradiance, and wind speed). Figure 4 shows typical meteorological year (TMY) data for Sherbrooke.

Local weather data of St. Elie greenhouse.
To verify the accuracy of the mathematical energy model, the actual and calculated energy consumption values of the St. Elie greenhouse are compared and depicted in Figure 5(a). The peak heating requirement occurred on February 3rd (with minimum temperature of −31.7 °C), which was 2.4% higher than the highest value in January (with minimum temperature of −30.5 °C). The results of the MATLAB code (MATLAB R2021b) show good agreement with energy consumption from St. Elie greenhouse, as extracted from the billing records for this period. The error of the developed MATLAB code for cold seasons, when the ventilation value is 0, is between 7% and 15%, which is acceptable. The early static energy model which were developed for greenhouse (Chiapale and Kittas, 1981; Morris, 1956) had accuracy of up to 25%, but this value for the recent dynamic models (Golzar et al., 2018; Singh R and Tiwari G, 2000) is around 11%. Weather data, and controlling strategy in greenhouse seems to have the highest contributions in generating errors.

(a) Annual variation of the heating requirements in the St. Elie greenhouse versus heating load calculated by code and (b) the effect of multi-zone formulation on calculated heating load for January and February in zone 3.
Furthermore, it can be noticed that for the cold season, when the mechanical ventilation is off (i.e. roof vents and wall vents are closed) the model accuracy is maximum. In cold climate greenhouses, energy consumption is at its greatest during the cool/cold non summer months. Given that ventilation is principally required during the warm summer months, the simplifying assumption was made that the greenhouse has no ventilation system throughout the year. Although this assumption will introduce some error, particularly during the summer, it will provide valuable information concerning the potential for energy recovery within the greenhouse.
To present the importance of multi-zone formulation, Figure 5(b) compares the heating load in January and February for two formulations: single-zone and multi-zone, in zone 3. To do so, first the code ran for zone 3 as a single zone. It means four side of this zone exposed with outdoor temperature. The same process was repeated for multi-zone formulation but with one modification: north wall considered exposed with set point temperature of zone 2, and its relevant convection and convection heat transfer coefficient. As shown, in the multi-zone formulation the calculated heating demand for the cold season is less than that of the single-zone formulation. This demonstrates that a multi-zone formulation can bring higher accuracy than a simplified single-zone formulation in greenhouse simulation and design.
Results and discussion
As explained, the St Elie greenhouse has three distinct production zones and a warehouse, each with different micro-climate conditions. This fact causes each zone to have a different heating/cooling demand throughout the year, even in the summer season when the solar radiation is in its maximum value. Figure 6 presents the heating load of the three zones and the greenhouse warehouse.

Heating demand of different zones of St. Elie greenhouse.
All commercial greenhouses have set points for temperature and relative humidity. These setpoints, as well as the target photoperiod, are presented in Table 2. The most convenient method to maintain the RH inside a greenhouse is natural ventilation through roof vents. However, while this method is the cheapest method, is also entails some disadvantages, especially in a cold agricultural region like Canada. It is reported that the dehumidification process by ventilation can cause up to 50% undesired heating load (Campen and Bot, 2001). Various studies (Kacira et al., 2004) have shown that there exists an optimum ratio (0.015–0.25) of the vent opening area to the ground area of a greenhouse to maximize the efficiency of natural ventilation. This ratio for St. Elie greenhouse is 0.21 when the roof vents are fully open. The main mechanism for natural ventilation is the existence of a pressure difference between the inside and the outside of the greenhouse. Parameters such as temperature, size, the location of the openings and wind speed can affect ventilation efficiency.
As part of this study, the dehumidification load of St. Elie greenhouse, as well as the contribution of different mechanisms involved in the water vapor balance of the indoor air were calculated. Figure 7(a) shows the total dehumidification load required to maintain the RH level of 75%. Figure 7(b) presents the contribution of different mechanisms involved in the water vapor balance of the indoor air: evapotranspiration, condensation, and infiltration (equation (32)). Also, it should be highlighted that

Hourly distribution of (a) dehumidification load and (b) contribution of different mechanisms affecting the humidity level.
Energy saving strategies
Heat recovery
To reduce the greenhouse energy consumption, different strategies can be implemented. Figure 8 shows the cooling demand of zones 2 and 3 for the month of January in the St. Elie greenhouse. As shown in Table 2, the set point temperature and relative humidity of these zones are 5°C and 75%, respectively. The classic method of satisfying the cooling demand in a commercial greenhouse is by means of mechanical ventilation. This means that the roof ventilators, through the intervention of a control system, are sensitive to the greenhouse inside humidity and temperature. As soon as the measured parameter passes the set point, the control system opens the ventilators. However, based on Figure 8(b), during the coldest period of year (typical day Jan.16), it is possible to transfer hot air from zone 2 to zone 3, since zone 2 needs to remove hot air, while zone 3 requires heating. Our study on the St. Elie greenhouse shows that by installing a HRV (Heat Recovery Ventilator) between these two zones or by partially opening the roll-up separator between these two zones, for January and December, when the setpoint temperature is 5°C for both zones, there is an opportunity to cut the energy consumption up to 4.69% for January and 4.13% for December.

(a) variation of cooling demand
There is another opportunity to take advantage of heat recovery between zone 1 and zone 3. Based on Table 1 and Figure 6, during December and January the set point temperatures for zones 1 and 3 are 24°C and 5°C in the morning and 18°C and 5°C at night, respectively. Investigating the hourly heating/cooling load which is calculated by mathematical code shows that it is possible to transfer the hot air of zone 1 (
As can be seen in Figure 9, the black area related to zone 1 corresponds to warm air that is at or above the set point of zone 1 (

Heating and cooling demand for the St. Elie greenhouse zones 1, 2, and 3.
Having discussed heat recovery opportunities, it would also be worthwhile to consider how some parameters can affect the energy performance of a greenhouse: Orientation, geometrical properties, material, etc.
Solar radiation
It is proved that every parameter that can affect the greenhouse radiation efficiency (absorbed radiation by the greenhouse envelope/total radiation) can affect its energy performance (Esmaeli and Roshandel, 2020). If we consider a greenhouse as a solar collector, the orientation, area, roof angle, and construction material are the parameters that can affect its radiation performance. Figure 10 shows the hourly solar intensity on the different walls and roofs of the most southern portion (zone 3) of the St. Elie greenhouse for a typical day in summer and in winter. As can be seen, during winter there is less solar irradiance than during the summer. The south roof and south wall receive the highest amount of solar energy throughout the year. However, in the winter season, the south wall receives more solar energy (+800

Hourly variation of the total solar irradiance on different walls and roofs of the case study greenhouse: (a) in a typical day in summer and (b) in a typical day in winter.
As illustrated in Figure 10, the “South Roof,” meaning the sum of all south facing roof portions of the entire greenhouse, makes the greatest contribution to the greenhouse solar energy gain, even more than the south-facing wall. Based on this fact, different configurations can be proposed. For instance, an asymmetrical roof configuration is only one of the possible greenhouse configurations of interest. Esmaeli and Roshandel (2020) conducted a mathematical multi-functional optimization to design a Chinese style greenhouse with a minimum energy demand. Our study proposes a similar study in the future.
Figure 11, shows both symmetrical and asymmetrical roof configurations. The results and literature review show that increasing the area of the south wall, while the volume of greenhouse is constant, can affect the greenhouse energy consumption. However, obtaining the optimal geometrical parameters of such a configuration needs further study.

Isometric view of proposed configuration for St. Elie greenhouse (zone 3). (a) Symmetric roof and (b) asymmetric roof.
Obtaining the effect of an asymmetrical roof configuration on energy consumption of a greenhouse in a cold climate area requires a full CFD simulation, however, increasing the south roof area can affect the solar energy gain. In Figure 11, the values of X and β are 5.35 m and 31.2° for the actual symmetric roof, and 6.7 m and 24.4° for a hypothetical asymmetric roof, respectively. In this simplified evaluation, the greenhouse height, material (thermal properties), and micro-climate condition of greenhouse (temperature and relative humidity set point) remained unchanged. Figure 12 shows and compares the available total solar irradiance (

Variation in total solar irradiance availability in January and December in the St. Elie region (45°N) for two different roof configurations.
The orientation of a greenhouse can also affect its energy consumption. To evaluate the relative importance of greenhouse orientation, a simulation was conducted wherein the greenhouse was rotated 90 degrees to the east, while maintaining the warehouse on the north wall. Our results showed that the current east-west orientation has 2.67% less energy consumption, in comparison with a north-south orientation. Based on Figure 13, for the coldest months in a year, like January, February, and December, solar radiation in an east-west oriented greenhouse is higher than that of a north-south orientation, while for the hot months it is the reverse.

Comparison of annual variation in total solar irradiance availability for even-span greenhouse in E-W and N-S orientation for the St Elie greenhouse.
A thermal screen, also known as an energy screen, shading screen, or thermal curtain, is an important component in modern greenhouse technology. It is used to help control and manage the micro-climate within a greenhouse by regulating the amount of sunlight and heat that enters or exits the structure; it can help regulate the internal temperature of the greenhouse. Note that one says that a thermal screen, or curtain, is open, when one can see the sky; it is closed when the sky is hidden. At night during cold weather, the thermal screens can be closed to retain heat within the greenhouse by reducing thermal loses to the cold night sky. They can block excessive sunlight during hot days, reducing the greenhouse temperature, especially during the summer. Greenhouses adopt a control strategy to activate the artificial lighting and thermal screens. Usually, artificial lighting, blackout screens, and thermal screens are controlled by means of a combination of a timing schedule and the solar radiation level. In our study, when the radiation level is less than 50 W/m2, the thermal screens are activated (i.e. closed). Also in zone 1, when artificial lighting is on, the blackout screens are activated. The schedule of artificial lighting is provided in Table 2.
Figure 14 shows and compares the heating demand of zone 3 of the St. Elie greenhouse in January, both with and without the use of thermal screens. The calculations confirmed that using a thermal screen can reduce the energy demand by as much as 20% for zone 3 in January.

Comparison of heating demand of zone 3 in January with and without thermal screen.
Conclusion
This study presented a mathematical energy model capable of calculating the heating, cooling and dehumidification loads of a multi-zone (i.e. multiple micro-climate) greenhouse. The development of the model is based on simultaneous energy and vapor mass balances. The resulting MATLAB-Simulink code was then used to investigate potential solutions for reducing the energy consumption in multi-zone greenhouses in cold climate agricultural regions, such as found in Quebec, Canada. As a case study, a greenhouse in Sherbrooke, Quebec, named Les Serres St-Élie, was selected. It covers more than 5000
Footnotes
Appendix
Notation.
| area of cover, floor, and plant, | |
| area of non-transparent and transparent surfaces, | |
| CF | cloud cover factor, Octas |
| specific heat of air, | |
| condensation on the cover, | |
| Dehumidification load, | |
| motor efficiency | |
| crop transpiration, | |
| inside/outside , kJ | |
| cover view factor, and sky view factor, dimensionless | |
| perimeter heat loss factor, W | |
| heat conversion factor, and lighting allowance factor, dimensionless | |
| motor load factor, and motor use factor, dimensionless | |
| Grashof number, dimensionless | |
| , , | specific ventilation, condensation conductance, transpiration conductance, m |
| H | depth of underground soil for constant temperature, m |
| HHV | high heating value of fuel MJ |
| thermal air conductance, W | |
| convection coefficient for indoor and outdoor surfaces, W | |
| direct beam radiation, and diffuse radiation on horizontal surfaces, W | |
| global solar radiation on horizontal surface W | |
| , | thermal conductivity of air, cover, soil, and non-transparent wall, W |
| K | thermal conductivity of ith section in composite wall, |
| characteristic length of convective surfaces and plant leaves, m | |
| latent heat of water vaporization J | |
| LAI | leaf area index |
| N | number of air changes per hour |
| number of layers in covering | |
| Number of fans | |
| Nu | Nusselt number, dimensionless |
| P | perimeter of greenhouse, m |
| motor power rating, W | |
| Prandtl number, dimensionless | |
| P | atmospheric pressure, kPa |
| partial pressure of the water vapor, kPa | |
| partial pressure at saturation, kPa | |
| electric power of artificial lighting, | |
| Solar heat gain, W | |
| Heat gains due lighting, W | |
| Heat gains due to fan, W | |
| Conduction heat transfer, W | |
| Heat gains/lost due to transpiration, W | |
| Heat transfer to ground, W | |
| Long wave radiation heat transfer, W | |
| perimeter heat loss, W | |
| Heat loss due to roof ventilation, W | |
| Heat loss due to mechanical ventilation, W | |
| Dehumidification load, W | |
| Natural ventilation heat loss, W | |
| , , | aerodynamic resistance, stomatal resistance, and boundary layer resistance |
| Reynold number, dimensionless | |
| S | total solar irradiance entering the greenhouse, W |
| , | indoor temperature, and outdoor temperature, K |
| cover temperature, K | |
| underground soil temperature, K | |
| and sky temperature, K | |
| , | heat transfer coefficient for transparent and nontransparent surfaces, W |
| volume of greenhouse, | |
| , | indoor airspeed, and outdoor airspeed, |
| vapor exchange with outdoor air by natural ventilation, | |
| saturated humidity ratio of air at plant temperature, Kg | |
| humidity ratio of air at indoor temperature, Kg | |
| Greek symbols | |
| X | Roof rafter length, m |
|
|
Sky clearness factor |
|
|
factor for estimation of effective solar radiation, |
|
|
angle of inclined surface with horizontal, ° |
| ϒ | surface azimuth angle, ° |
| η th | Thermal ratio of lamps |
| emissivity of cover and indoor components, |
|
|
|
angle between two radiative surfaces, ° |
|
|
zenith angle of sun, ° |
|
|
angle of incidence of surfaces, ° |
|
|
air density, |
|
|
reflectivity of outdoor ground, dimensionless |
|
|
transmissivity of cover, dimensionless |
|
|
transmissivity of cover to long-wave radiation, |
| μ | dynamic viscosity of air, |
|
|
local latitude ° |
| ∂ | volumetric thermal expansion coefficient, |
| σ | Stefan-Boltzmann Constant, W |
| ω | hour angle, ° |
| ΔT | temperature difference, °C |
| absolute water vapor concentration, |
Credit author statement
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors acknowledge the financial support of the Natural Sciences and Engineering Research Council of Canada for this project (RGPIN-2019-05826).
