Incorporating cultivar-specific stomatal traits into stomatal conductance models improves the estimation of evapotranspiration enhancing greenhouse climate management

Oliver K€ orner , Dimitrios Fanourakis , Michael Chung-Rung Hwang , Benita Hyldgaard , Georgios Tsaniklidis , Nikolaos Nikoloudakis , Dorthe Horn Larsen , Carl-Otto Ottosen , Eva Rosenqvist c a Leibniz-Institute of Vegetable and Ornamental Crops (IGZ), Großbeeren, Germany b Laboratory of Quality and Safety of Agricultural Products, Landscape and Environment, Department of Agriculture, School of Agricultural Sciences, Hellenic Mediterranean University, Estavromenos, 71004, Heraklion, Greece c Department of Plant and Environmental Sciences, University of Copenhagen, Hoejbakkegaard Alle 9, 2630 Taastrup, Denmark d Klasmann-Deilmann Asia Pacific Pte. Ltd, Singapore e Department of Food Science, Aarhus University, Agrofoodpark 48, 8200 Aarhus N, Denmark f Institute of Olive Tree, Subtropical Plants and Viticulture, Hellenic Agricultural Organization ‘Demeter’ (NAGREF), P.O. Box 2228, 71003 Heraklio, Greece g Cyprus University of Technology, Department of Agricultural Sciences, Biotechnology and Food Science, Cyprus h Horticulture and Product Physiology Group, Wageningen University, Wageningen, the Netherlands


Introduction
Climate-controlled greenhouses are high energy-consuming facilities, especially in hot or cold regions (Bot, 2001;Vanthoor, 2011, p. 307). Since the early 1990's, energy demand has been reduced by the introduction of crop directed climate control regimes with dynamic temperature boundaries (e.g. Aaslyng et al., 2003;K€ orner & Challa, 2003a). These exploit the naturally-occurring diurnal temperature fluctuations without the need for heating to a fixed temperature, achieving up to 18% savings in energy use for heating in e.g. Dutch climate conditions (K€ orner et al., 2004). Increasing greenhouse insulation (i.e., semi-closed and closed greenhouse systems) has further reduced energy use required for heating (Gelder et al., 2012;Opdam et al., 2005). However, temperature and relative air humidity (RH) control remain the two main challenges in these concepts (K€ orner & Challa, 2003b;Vadiee & Martin, 2012). Potential risks with accidental rise of RH above 85% include increased incidence of fungal diseases (e.g., powdery mildew and botrytis; Mortensen et al., 2007), induction of physiological disorders such as the blossom-end rot in fruit vegetables (Fanourakis, Aliniaeifard, et al., 2020;Ho et al., 1993) and, in the long-term, the development of stomata with reduced closing ability Giday et al., 2014). Although the occurrence of malfunctional stomata remains undetected during cultivation, it becomes apparent when plants rapidly wilt upon exposure to lower RH such as during display Fanourakis et al., 2015Fanourakis et al., , 2019a. However, the objective of keeping RH below 85% is counteracting the striving for low energy demand. The energy needed for dehumidification (by heating) increases when RH is adjusted to lower levels. Despite great efforts devoted to the development of more efficient dehumidification methods, reducing RH to a secure level ( 85%) remains rather costly, accounting for about 20% of annual energy demand in Northern European climates (K€ orner & Van Straten, 2008). As such, improved RH control will lead to a major decrease in both energy demand and CO 2 footprint, combined with a better (visual and inner) quality of the produce. In vegetable crops, inner quality refers to the nutritional value (Amitrano et al., 2021;Zhang et al., 2017), whereas in ornamental crops to the keeping quality Fanourakis et al., 2016bFanourakis et al., , 2020b. Since the crop itself is a major driver of RH via crop evapotranspiration (ET c ), its accurate estimation is essential for greenhouse climate control. The ET c depends on both environmental conditions [e.g. short-wave radiation, vapour pressure deficit (VPD) and air flow] and crop-related factors [e.g. crop/plant architecture and leaf stomatal conductance (g s )]. Despite its obvious importance, models estimating ET c mostly ignore some or all of Ball index (mol m À2 s À1 ) b empirical parameter Van Genuchten (À) BWB-model Ball-Woodrow-Berry model BWB-Liu-model BWB-model modified by Liu c measuring errors' correction (À) C s air CO 2 partial pressure at the leaf surface (mbar) D s stomatal density (number mm À2 ) DSS decision support system ET c crop evapotranspiration (L plant À1 ) g s stomatal conductance (mol m À2 s À1 ) g s,meas measured g s (mol m À2 s À1 ) g s,sim modelled g s (mol m À2 s À1 ) g 0 minimum (residual) g s as P nl reaches 0 (mol m À2 s À1 ) h s leaf surface RH fraction (À) I short-wave radiation (W m À2 ) lE latent heat of evaporation (J m À2 s À1 ) m s adjustment factor (À) m g empirical parameter Van Genuchten (À) n g empirical parameter Van Genuchten (À) n empirical parameter BWB-Liu model (À) Q h heat energy consumption (MJ m À2 ) p a Atmospheric pressure (kPa) P nl leaf net photosynthesis rate (mmol m À2 s À1 ) P nl,max maximum leaf net photosynthesis (mmol m À2 s À1 ) PPFD photosyntheticphotonfluxdensity(mmolm À2 s À1 ) J substrate water potential (kPa) RH relative air humidity (%) SLA specific leaf area (cm 2 g À1 ) q volumetric water content (v v À1 ) q/q s effective water content (v v À1 ) q r residual water content (v v À1 ) q s saturated water content (v v À1 ) VPD vapour pressure deficit (kPa) w pd dry pot mass (g) w pf fresh pot mass (g) z canopy levels these crop-specific factors, based on the lack of relevant parameters for the wide variety of crops and cultivars. Most investigators dealing with g s estimations use the concept of BalleWoodroweBerry (BWB) (Ball et al., 1987). The BWB concept describes the relationship between g s and net photosynthetic rate (P n ) using two factors, i.e. minimal stomatal conductance (g 0 ) and an empirical sensitivity coefficient for the slope (m s ). In that, the inputs to the combined P n e g s model are commonly measured climate or microclimate variables: short-wave radiation (I, W m À2 ), CO 2 (mmol mol À1 ), temperature ( o C) and RH (%). Although the model is widely used, it has been found to underestimate g s and transpiration in roses (Kim & Lieth, 2003) and to be less suitable at conditions with high temperature and elevated CO 2 (Janka et al., 2016). To improve the BWB model, soil water potential (J) as additional variable has been introduced adjusting the value of m s (Liu et al., 2009). The resultant BWB-Liu-model has shown improved accuracy of g s estimation under different irrigation regimes (Liu et al., 2009). However, two issues remain: 1) through the global parameterisation of the BWB type models applied to a wide range of crops and varieties, the crop and cultivar specific parameters induce a mismatch between model predictions and reality; 2) the newly introduced parameter J in the BWB-Liu-model is difficult to access in common greenhouses and climate control regimes, as it needs additional measurements of the root zone, while the remaining input variables of the BWB-model are commonly available climate variables in both greenhouses and simulators.
Despite the rather limited usability of the simple BWBmodel in commercial practise under dynamic climate conditions (e.g., high temperatures, elevated CO 2 , and highly fluctuating RH) and uncertain parameterisation, its practical advantage of not needing measured J as input could in some situations cover for the misfit. In this paper, we address the pros and cons of both models and present a possible improvement to them. We focus on the variability of ET c and g s in selected rose cultivars with similar external appearance, aiming to show the effect of cultivar specific BWB-model parameters on prediction quality and effect in model application to predictive greenhouse climate control.
We hypothesise that there are further possibilities for improvements to the BWB-model types by including cultivar specific traits into the model, i.e. the stomatal pore area (A p ) and stomatal density (i.e. stomatal number per unit area; D s ). These factors vary considerably not only between taxa, but also among cultivars of a given species . A p is dynamically adjusted by changes in pore aperture, since pore length is rather rigid during opening and closure of stomata Franks et al., 2009;Islam et al., 2019). Active pore aperture adjustments in response to internal (e.g. water status) and external (environmental) factors are physiologically regulated . Instead, pore length and D s are anatomical features, which are set during leaf elongation (Dow et al., 2014;Fanourakis et al., 2014). Therefore, by integrating A p or D s , as factors determining operating g s , the ET c estimation is expected to improve.
The current study presents the importance of cultivar specific parameterisation of a g s -model for application in predictive greenhouse climate control in a simulation study with an existing greenhouse simulator (Goddek & K€ orner, 2019;K€ orner & Hansen, 2012;K€ orner & Holst, 2017). To test the hypothesis that employing low-transpiring genotypes will improve energy saving (required for dehumidification), a comparative simulation study of energy demand at different climate setpoints with two cultivars was performed. Given that this hypothesis is validated, it is possible that 1) measures of g s may add value in selection indices of pot plant breeding programs, since these directly affect the energy demand which is related to RH control, and 2) parameterising the BWB-model types for each crop and cultivar in a monoculture enables energy saving options for model-based greenhouse climate controllers. Rose was employed as the model system since it is one of the most important ornamental crops, and a species in which genetic differences in terms of ET c have been observed.

Materials and methods
The current paper analyses the effect of genetic differences in stomatal conductance (a major driver for ET c ) on greenhouse humidity-related energy demand. The effect of climate variables is researched, and possible solutions for a combined strategy on genetic breeding and model-based climate control are suggested. In that, a step-by-step approach was employed, including physiological investigations, model development, calibration and validation, model selection and implementation, and eventually a simulation study (Fig. 1  b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 (Giday et al., 2013b(Giday et al., , 2015, six commercially available pot rose cultivars with uniform architecture were selected ('Alaska', 'Aloha', 'Apache', 'Felicitas', 'Flirt', and 'Pasadena'). The day/ night air temperatures were 22 ± 0.5 C/18 ± 0.4 C. Four treatments were applied, consisting of two irrigation regimes (70 and 100% of ET c , based on well-watered target weight, 330 g) combined with two RH levels (60 and 80% applied over time, see below). Water was supplied to the respective target mass, which was adjusted to account for plant growth [100% being 330 g (week 0e2), 340 g (week 3), 360 g (week 4) and 385 g (week 5e6)]. The RH was set to 60% (weeks 0e2 and 4e6) and 80% (week 3) to induce a different transpiration rate. In Exp. 2, conducted for 4 weeks (starting Jan. 12, 2016; Table 1), the effect of irrigation regime on plant transpiration was evaluated. Based on Exp. 1, four cultivars with expected large difference in transpiration were selected ('Alaska', 'Aloha', 'Apache' and 'Flirt'). Three treatments were applied based on irrigation regime (60, 80, and 100% of ET c ), adjusted as described in Exp. 1. RH was set to 80% throughout the fourweek experimental period, with remaining climate parameters set as in Exp. 1.
In Exp. 3, conducted for 6 weeks (starting Nov. 21, 2016; Table 1), the effect of air temperature on plant transpiration was determined on 'Alaska', 'Aloha', 'Pasadena' and 'Flirt' exposed to various watering regimes. Due to availability restrictions, 'Apache' used in Exp. 2 was replaced with 'Pasadena', which had similar ET c according to Exp. 1. Four treatments were applied, consisting of two irrigation regimes (60 and 100% of ET c , adjusted as in Exp. 1) combined with two air temperature levels [26/22 C (warm), 18/14 C (cold) for day/ night]. RH was adjusted to keep VPD constant among the two air temperature treatments (i.e., the day/night RH was 88/86% in the warm chamber, and 80/75% in the cool chamber), with remaining climate parameters were set as in Exp. 1.

Greenhouse setup and growth conditions (Exp. 4e6)
The greenhouse experiments were conducted in 50 m 2 greenhouse compartments fitted with pipe heating connected to a calorifier with ventilator under the tables to prevent climate gradients, and cooling by passive roof ventilation. The climate was controlled by a climate computer (LCC4, Senmatic A/S, Søndersø, Denmark). Air temperature set points were 21 C with ventilation at 30 C during daytime, and 17 C with ventilation at 20 C at night. RH was set to 70%, maintained by a high-pressure humidification system (Condair Systems, Pf€ affikon, Switzerland). The CO 2 concentration was set to 600 mmol mol À1 , which was realised by supplying pure CO 2 in the daytime when vents were closed. This carbon dioxide level is common in commercial rose crops. In addition to the climate sensors connected to the climate computer, air temperature and RH were also recorded with sensors (Sensohive Technologies ApS, Odense, Denmark) installed in the upper third of the plant canopy on each weighing table.
The plants were placed on six custom built (0.61 Â 1.08 m) ebb-flow tables (Staal & Plast A/S, Ringe Denmark) in a 2 Â 3 table pattern in each of two greenhouse compartments with automatic weight registration by single point load cells (1-PW12CC3, HBM, Darmstadt, Germany) connected to data loggers (DataTaker DT85, Thermo Fisher Scientific, Waltham, MA, USA). The irrigation system was constructed with commercially Table 2 e Parameters for various stomatal conductance models obtained from non-linear regression over the data of Experiment 6 (see Table 1 ); J e Substrate water potential (cm).
b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 available components for watering mobile tables (Anderup El, Odense, Denmark), having T-shaped pipes connected to the central watering/fertilisation computer (AMI Completa, Senmatic A/S, Søndersø, Denmark) pumping nutrient solution onto the tables without the pipes contacting the tables, and thus not affecting the weighing. The draining time was controlled by the 4 mm hole in the drain valve and the drain water was collected in gutters under the tables for recirculation via the central watering system. To control the water level on the tables during irrigation, a fully open 15 mm pipe was fitted through the table above the gutter, extended 15 mm above the table surface as overflow drain, when the nutrient solution was pumped up on the table. The plants were given a full nutrient solution automatically adjusted to pH of 6.0 and EC of 2.0 mS cm À1 (Fanourakis et al., 2008). Each weighing table was fitted with 35 cm high transparent green plastic filter (LEE Fern Green 122, LEE Filters, Andover, UK) around its perimeter to eliminate the difference in canopy transpiration between tables having one or two outer edges exposed to the surroundings. In the 6 weeks of Exp. 4 (starting June 14, 2016; Table 1), the effect of light intensity and plant spacing on plant transpiration was determined. Based on Exp. 1, the two cultivars with the largest difference in plant transpiration ('Alaska' and 'Flirt') were chosen. Cultivation took place during the summer period, where solar irradiation potentially reaches the highest annual values. Six treatments were applied, consisting of two light intensity regimes [(1) shading screens (50% shade) unfolded when outside radiation was >500 mmol m À2 s À1 PPFD; (2) no screens (full sunlight)] without supplementary light, combined by three plant spacing levels (12, 18 and 24 plants per weighing table, i.e., 25, 37.5, and 50 pots m À2 ). Plants were kept well-watered by using ebb-and-flow irrigation for a 6e9 min period every second day, resulting in water on the tables for ca 12e15 min, depending on plant size.
In Exp. 5 and 6, conducted for 6 and 9 weeks, respectively (starting Oct. 11, 2016 and Feb 17, 2017; Table 1), the effect of RH and plant spacing on plant transpiration was determined on 'Alaska' and 'Flirt', respectively. In both experiments, two RH regimes (30 and 70%) were combined with two plant spacing levels (12 and 24 plants per weighing table, i.e., 25 or 50 pots m À2 ). From 6:00 to 23:00, supplementary lighting (200 mmol m À2 s À1 PPFD) was provided at low solar radiation (PPFD < 100 mmol m À2 s À1 ) inside the greenhouse compartment, measured by a quantum sensor (LI-190SA, LI-COR, Lincoln, NE, USA) connected to the climate computer. In both experiments (5 and 6), the remaining climate parameters were set as in Exp. 4.

Biomass and leaf morphological components
In the greenhouse experiments, three pots per cultivar were harvested at day 0, and each week one plant per weighing table was randomly sampled for measurements of leaf area and biomass (see below). The daily plant leaf area was extrapolated from linear connection of the data points of the mean values per cultivar and treatment. Plants that were harvested from the weighing tables were replaced with spare plants grown on an adjacent table with the same plant densities and watering as the weighing tables. All harvest plants were randomly preselected, and all original plants were harvested before the substitute plants were used for harvest.
In all experiments, the final biomass production was assessed. Leaf number and area (leaf area meter LI-3100, LiCor, Lincoln, NE, USA), together with leaf and stem dry masses were recorded after drying the plant tissue for at least 48 h at 80 C in a drying oven (Seif et al., 2021). The specific leaf area (SLA; leaf area/leaf mass) was calculated. The measurements were conducted on three plants per weighing table (Exp. 3 and 4) or four plants per treatment (Exp. 1, 2, 5 and 6).
The effect of growth environment on stomatal pore length (major axis), pore aperture (minor axis), and D s was determined. Rosa x hybrida has compound leaves, where leaflets arise on both sides of the rachis (i.e. in pairs) besides the terminal leaflet (odd-pinnate arrangement). Measurements were carried out on one leaflet of the first pair of lateral leaflets from the first order pentafoliate leaf (counting from the apex). The silicon rubber impression technique was employed [elite HDþ, Zhermack, Badia Polesine, Italy ]. Imprints were made 2 h following the onset of the light period, since this time is required for plants exposed to prolonged darkness (i.e. night time period) to reach steady-state stomatal aperture. The sampling area (1 Â 1 cm) was located midway between the leaflet base and tip, and between the midrib and lateral margin Sørensen et al., 2020). The fields of view were located in interveinal areas, since veins lack stomata . Images were acquired using an optical microscope (Leitz Aristoplan; Ernst Leitz Wetzlar GmbH, Wetzlar, Germany) connected to a digital camera (Nikon DXM-1200; Nikon Corp., Tokyo, Japan). The rose is a hypostomatous species , so only the abaxial leaflet surface was assessed. Stomatal pore length and aperture were determined on ten randomly selected stomata (magnification Â 200), while D s was counted on five non-overlapping interveinal fields of view per leaflet (magnification Â 100). A p was calculated as the area of an ellipsis with major and minor axes being the pore length and pore aperture, Image processing was performed with ImageJ software (https:// imagej.nih.gov/, Koubouris et al., 2018;Fanourakis et al., 2021). Three (Exp. 5 and 6) or four (Exp. 2 and 3) leaflets were assessed per treatment.

Stomatal conductance (g s )
The effect of environmental conditions on g s was examined by porometry. In situ g s measurements were done on attached leaves of intact plants.

Plant transpiration
The effect of growth environment on whole plant transpiration was evaluated during growth. For the growth chamber b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 experiments (1e3) the water loss of individual pots was (automatically) determined by the mass difference in 5-min intervals using the two Drought Spotter units, each containing 24 single pot balances (Fanourakis et al., 2017). Growingmedia evaporation rate was assessed in three (Exp. 1) or five (Exp. 2, and 3) pots without plants, but with identical media water content as in pots with plants after a duration of four days with treatment. Direct evaporation from the media was subtracted from pot water loss to determine plant water loss (Fanourakis et al., 2017(Fanourakis et al., , 2019a(Fanourakis et al., , 2019b. Treatments were compared based on the average transpiration rate over the complete assessment period. For the greenhouse experiments (4e6), the canopy water loss of a 0.61 Â 1.08 m canopy was determined by the previously described weighing tables, where the mass was logged every 3 min. The plants were watered daily at 07:00, and the transpiration rate was calculated as the slope of the water loss over time from 11:00 to 16:00, since the transpiration curve was linear in this time interval and unaffected by the watering.
In all cases, transpiration rate was calculated as per plant for the last five days of the experimental periods.

Plant transpiration model
Evapotranspiration of a single plant was modelled as the sum of evapotranspiration of all leaves of that plant, and then it was up-scaled to the crop level for simulations. For vertical leaf distribution (within-plant or within-plant stand) and microclimate differences, the approach of separating the plant or crop into three horizontal layers was employed (described by K€ orner et al., 2007). Stomatal conductance was calculated for individual leaves with four different options: (1) the BWB-model (Ball et al., 1987), (2) the BWB-Liu-model (Liu et al., 2009), or (3, 4) the two adjusted versions of BWB-Liumodel (by employing A p and D s , respectively).
The empirical BWB-model as used by Kim and Lieth (2003) uses leaf net photosynthesis rate (P nl ), air CO 2 partial pressure (C s ), and h s (leaf surface RH fraction) as driving variables: where g 0 is the leaf minimum (residual) g s as P nl reaches 0, while m s is an empirical coefficient for the sensitivity of g s to P nl ; m s includes the empirical scaling to mol. Practically, the BWB-model describes the correlation of g s with physiological CO 2 assimilation and vice versa. P nl is input for the g s calculation, while inversely g s is input for the P nl calculation. This clustering of two sub-models (with interdependent input and output) demands start values and iteration in simulation.
The BWB-Liu-model uses substrate water potential (J) to explain the fluctuation of the slope (m s ). The average pot mass (obtained from the weighing tables) was converted into J, by using the Van Genuchten equation (Van Genuchten, 1980) with saturated water content (q s ), the residual water content (q r ; assumed to be 0), and the two dimensionless constants m g and n g . The effective water content (q/q s ) was obtained by using Equation (2). The volumetric water content (q) was obtained from linear least squares' regression by the empirical Equation (3) (using the fresh and dry pot mass data; w pf , w pd , respectively; Exp. 6), where c is the correction of measuring errors.
P nl was obtained by fitting the negative exponential light response curve (Thornley, 1976, p. 318) with leaf photochemical efficiency (a l , [mol CO 2 ] [mol photons] À1 ) and maximum leaf net photosynthesis (P nl,max, mmol m À2 s À1 ) for the different canopy levels (z) according to the PPFD (mmol m À2 s À1 ).
Absorbed PPFD values were calculated separately for diffuse and direct radiation and as function of z (Van Kraalingen & Rappoldt, 1989). Parameters a l and P nl,max were calculated based on a combined method of biochemical leaf photosynthesis models (Farquhar et al., 1980;Farquhar & Von Caemmerer, 1982) and approaches of Gijzen (1994), as shown by K€ orner (2004).
We have investigated the extension of the BWB and the BWB-Liu models with A p or D s . All models [(1) BWB, (2) BWB-Liu, (3, 4) BWB and BWB-Liu extended with A p , and (5, 6) BWB and BWB-Liu extended with D s ] were parameterised using data sets from both experimental measurements (leaf g s , leaf temperature, w pd ), environmental variables (air temperature, h s , and PPFD) of Exp. 4e6 or model calculations (P nl ).

Model parameterisation and implementation
For comparison of the two cultivars with extreme differences in potential ET c ('Alaska' and 'Flirt'), the two parameters of the BWB-model (m s and g 0 ) were parameterised with the fit procedure with non-linear least squares of the mathematical software environment MATLAB (ver. R2020a, The MathWorks, Inc., Nattics, USA). The parameterised BWB g s -model was implemented as sub-model in a crop microclimate and evaporation model (K€ orner et al., 2007) and validated as major driver for ET c . As such, the model was validated with measured ET c data obtained in Exp. 4 for the cultivars 'Alaska' and 'Flirt' with three plant spacing levels (12, 18 and 24 plants m À2 ) for an 11-d period each. For simulations, the crop specific parameters for pot rose and measured greenhouse climate data (3-min interval; temperature, RH, CO 2 concentration, and PPFD) were used as input to the model. Climate data were recorded with a combination of sensors (Sensohive Technologies ApS, Odense, Denmark; Parrot, Paris, France) installed close to the plant canopy. The cultivar specific parameters of the g s model were the only difference between the two cultivars as implemented in the ET c model.
For the data of Exp. 6, the CO 2 concentration was set constant at 600 ppm and the BWB model, the BWB-Liu, and the two extended BWB-Liu models were further parametrised for 'Flirt' by fitting the observed data with the mathematical software environment R (The R Project, r-project.org) using non-linear least squares method [nls()] and model prediction performance with linear regression analysis [lm()].

Simulation study
After validation of the g s model (parameterised with two opposite reacting rose cultivars; 'Alaska' and 'Flirt') as a submodel in an ET c model, the ET c model was further implemented in an existing greenhouse simulator (Goddek & K€ orner, 2019;K€ orner & Hansen, 2012) as implemented in HORTSIM v.0 (K€ orner & Holst, 2017). The simulator consists of a complete set of physical, biological and crop specific submodels, including a replica of a commercially-used climate controller. A standard multi-span greenhouse (1-ha floor area) with standard equipment for heating (pipes), ventilation (roof vents), screening (energy and shading), and supplementary lighting (max. 60 W m À2 high pressure sodium lamps) as well as set points for climate control of pot roses was used. Copenhagen (Denmark) was selected as location, and a data set of hourly values for a standard climate year in Denmark (Wang et al., 2013) was used. For climate control, set points for all actuators were calculated based on basic set points and climate inside and outside the greenhouse (as in commercial practise) e.g. heating and passive roof ventilation for temperature and RH control. Individual simulations with two cultivars ('Flirt' and 'Alaska'; all other parameters and conditions remained the same) were undertaken for a 365 d period with a 5 min time step, and integrated for each 60 min. Eight climate scenarios with different RH set points (ranging between 65% and 100%) were performed. Output from the simulations included hourly values of g s , greenhouse macro and microclimate variables (e.g. leaf temperature in three layers), ET c , as well as overall greenhouse energy household and demand.

Data analysis and statistics
All data analyses were performed with the mathematical software package Matlab (ver. R2020a; The MathWorks Inc.) incl. Curve Fitting Toolbox and Statistics and Machine Learning Toolbox. Statistical analyses were performed with analysis of variance (ANOVA) followed by a post-hoc test applying the Tukey's HSD test (TukeyeKramer) on the alpha level 0.05 using the procedures anovan and multcompare. Parameter estimation of the BWB, BWB-Liu and the extended versions of the BWB-Liu model was done with non-linear least square estimation (nlinfit and fit procedures).
Increasing RH (from 60 to 80%) and/or irrigation (from 70 to 100%) led to enhanced g s in all studied cultivars (Table 3).

Effect of irrigation level on g s and ET c
In Exp. 2, cultivar differences in leaf area were similar to Exp. 1 (Fig. 3A). 'Alaska' and 'Flirt' had the thickest and thinnest leaves with 218 and 287 cm 2 g ¡1 , respectively (as indicated by SLA, Fig. 3B). Similarly to Exp. 1, 'Flirt' and 'Aloha' had lower g s , as compared to 'Alaska' and 'Apache' (Fig. 3C).
A slight increase, owing to irrigation regime, was noted for leaf area (Fig. 4A) and SLA (Fig. 4B) in all four cultivars. While g s generally decreased with increased irrigation (100 as compared to 60%) in 'Aloha' and 'Flirt', a minor effect was observed for 'Apache' and 'Alaska' (100 as compared to 60%; Fig. 4C). However, the increased irrigation generally led to a higher ET c in all cultivars but 'Flirt' (Fig. 4D).

Effect of temperature (combined with RH) on g s and ET c
Based on the screening experiments (1 and 2), two cultivars ('Flirt' and 'Alaska') with extreme differences in ET c were selected for further analyses. The cultivation temperature (combined with RH to maintain a stable VPD) had only a minor effect on plant leaf area (Fig. 5A), while SLA decreased with increasing temperature in both cultivars (Fig. 5B). Instead, higher temperature enhanced both g s (Fig. 5C) and ET c (Fig. 5D). The effect of temperature on g s was more pronounced in 'Flirt' than in 'Alaska' (Fig. 5C). In contrast to the expected from increased g s and ET c with rising temperature, a higher temperature led to decreased stomatal pore aperture (Fig. 6A), pore length (Fig. 6B), and density in both cultivars (Fig. 6C).

Effect of plant density on g s and ET c
Plants were harvested at 7 or 14 d after the start of experiment, when plants were just after second cut. For either harvest date, the effect of plant density on plant leaf area was generally limited and not consistent among cultivars (Fig. 7A). Instead, increasing plant density clearly led to enhanced SLA (Fig. 7B). As the plants grew and the canopy closed, increasing plant density also led to a consistent increase of g s from cultivation day 6 and onwards (Fig. 8).

Effect of RH on g s and ET c
Increasing the set point of RH (from 30 to 70%) resulted in significantly different greenhouse climate (data not presented). The high RH (i.e. low VPD) led to both increased g s (Fig. 9A) and decreased ET c (Fig. 9B). Leaves developed at high RH had significant larger A p (Fig. 10A), while D s was not significantly different (Fig. 10B), as compared to leaves developed at low RH.

3.2.
Model and simulations

BWB-model parameterisation and validation
The BWB model was parametrised with data from Exp. 5 for the two cultivars 'Flirt' and 'Alaska' with and without additional term for CO 2 concentration (Table 4). Modelled g s (g s,sim ) was then implemented in an ET c model (K€ orner et al. 2007) and a validation study comparing g s,sim with measured g s (g s,meas ) on an independent data-set was performed. For that, ET c of the two cultivars with three plant densities (12, 18, 24 plants m À2 ) was used (Exp. 4). Generally, a very good agreement between simulated and measured ET c was noted, while the effect of plant density on ET c was also very well approximated by using the model (Fig. 11). It was clear that a higher plant density lowered plant ET c . At the lowest plant density (12 plants m À2 ), a slight underestimation of ET c was observed, and the model slightly underestimated the cultivar differences in ET c (Fig. 11).

Model extension
As first step, the BWB-model was extended with A p for the cultivar 'Flirt' in order to analyse the effect on improved fitting quality by adding an extra dynamic parameter. For that, this extended BWB-model was compared to measurements of Exp. 5 with two RH set-point levels (30, 70%) independently. The extension resulted in only a small difference to the original model (Fig. 12). In both model cases, a poor fit to measured g s was observed at low humidity levels (R 2 of 0.275 or 0.268), while at higher RH levels the fit for the two model versions was rather good (R 2 of 0.765 or 0.766; Fig. 12C,D). It was concluded that the current model was not capable of predicting g s at low RH with and without the extension of a dynamic A p , as it does not adequately take the greenhouse climate dynamics into account: in addition to the climate variables CO 2 and RH, it employs only one single variable (P nl ), incorporating the plant physiological responses to all climate variables. To improve the prediction quality at low RH levels and under dynamic climate conditions, the BWB-Liu model that includes substrate J as additional variable was therefore used in the next step, extended and parameterised with data of Exp. 6 for 'Flirt'. Four versions of the BWB-model were compared: (1) BWB-model, (2) BWB-Liu-model, (3) modified BWB-Liu-model by employing A p , and (4) modified BWB-Liumodel by employing D s (Fig. 13). In this parameterisation procedure, the complete climate data set with all RH values was used. As expected from Exp. 5 fittings, the model quality predicting g s was poor when the basic BWB-model was used (R 2 ¼ 0.12, Fig. 13A, Table 2). The BWB-Liu-model, which was extended with J, performed better, though the R 2 was still low (0.31; Fig. 13B; Table 2). Instead, by modifying the BWB-Liumodel further by adding stomatal features allowed much higher accuracy in estimating g s . When adding A p , R 2 increased to 0.52 (Fig. 13C, Table 2), while adding D s improved it to 0.59 ( Fig. 13D; Table 2). The extensions are thus capable of overcoming the weak model behaviour at low RH levels. At  Table 1). Data of the four treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). Different letters indicate statistically significant difference (alpha 0.05) with Tukey's HSD test (TukeyeKramer). Table 3 e Stomatal conductance of six pot rose cultivars in Experiment 1 (see Table 1). Plants were grown under two watering levels (70, 100%) and at moderate (60%) or high (80%) relative air humidity (RH). RH was set to 60% (0e2 and 4e6 weeks) or 80% (3rd week), depending on the week of growth. Measurements were conducted in the growing environment and on intact plants 2 h following the onset of the light period (n ¼ 3). Standard deviation is indicated in brackets.

Cultivar
Stomatal conductance (mmol m À2 s À1 )  (61) b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 higher RH levels (approximately >60%) the BWB-model, however, is a useful and simple predictor of g s without the necessity of additional input of stomatal traits or J, which are accessed with difficulty.

Greenhouse climate control simulations
Despite the improvements in model predictions with the extended BWB-Liu model versions, the parameterisation of the BWB-Liu model and the extensions is difficult in practise.  Table 1). Data of the three treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). Different letters indicate statistically significant difference (alpha 0.05) with Tukey's HSD test (TukeyeKramer).  Table 1). Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). Different letters indicate statistically significant difference between cultivars in each irrigation level (alpha 0.05) with Tukey's HSD test (Tukey-Kramer). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Fig. 6 e Effect of temperature (combined with relative air humidity) on stomatal pore aperture (A), stomatal pore length (B), and stomatal density (C) of the two most contrasting pot rose cultivars ('Alaska', magenta; 'Flirt', cyan) in Exp. 3 (see Table 1). Data of two irrigation treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)  Table 1). Data of two irrigation treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Fig. 8 e Effect of planting density on stomatal conductance (g s ) at 2 (A), 6 (B), 10 (C), 12 (D) days after the start of experiment (after second cut) of two most contrasting pot rose cultivars ('Alaska', magenta; 'Flirt', cyan) in Exp. 4 (see Table 1). Data of two light intensity treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Fig. 7 e Effect of planting density on leaf area (A), and specific leaf area (SLA; B) at 7 (lower boxplots) and 14 (upper boxplots) days after the start of experiment (after second cut) of two most contrasting pot rose cultivars ('Alaska', magenta; 'Flirt', cyan) in Exp. 4 (see Table 1). Data of two light intensity treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) As the simple BWB-model performs well under high RH levels, the implications of cultivar differences in g s for the greenhouse energy household was therefore quantified with the simple BWB-based ET c model using a range of elevated RH levels (with 65% as lowest). The comparative simulation study between the two cultivars 'Alaska' and 'Flirt' is illustrated in Fig. 14. Simulations showed that 'Flirt' has a consistently lower g s as compared to 'Alaska' (Fig. 14A) and a general decrease in latent heat production (measure for ET c ) was observed when increasing the RH set point (Fig. 14B). The relative difference in ET c between 'Flirt' and 'Alaska' was highest at lowest RH set point, while the difference decreased with increasing the RH set point ( Fig. 14C; absolute difference in Fig. 14D). The relative difference in energy demand between 'Flirt' and 'Alaska' was approximately 1% for 65%e75% RH, and between 2.35 and 5.75% for 85%e95% RH (Fig. 14E; absolute difference in Fig. 14F), peaking at the RH set point of 90%. As  Table 1). Data of two spacing treatments were pooled. Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). Different letters indicate statistically significant difference (alpha 0.05) with Tukey's HSD test (TukeyeKramer). Fig. 9 e Effect of relative air humidity on stomatal conductance (g s ; A) and evapotranspiration rate (ET c ; B) of two most contrasting pot rose cultivars ('Alaska', magenta; 'Flirt', cyan) in Exp. 1 (see Table 1). Whisker (W) was set to 2.5, and points were drawn as outliers if they were larger than Q 25 þW*(Q 75 -Q 25 ) or smaller than Q 25 -W*(Q 75 -Q 25 ). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) expected, due to the lack of RH control at a set point of 100%, no difference in heating energy consumption between the two cultivars was observed (Fig. 14E, F).

Discussion
In the current study, g s was measured employing a set of six pot rose cultivars and a wide range of climate conditions (varying RH, irrigation, temperature, light level and planting density; Table 1), and results were utilised to improve the coupling of estimated and actual ET c (Fig. 1), and in turn contribute to the energy saving potential. This study reveals the necessity of cultivar specific g s models for improved estimations of ET c , and in this way enhanced climate management. By using a cultivar specific parameter set within an ET c model, the ET c of individual cultivars could be more accurately estimated during cultivation.

Low g s genotypes improve energy saving related to RH control
In greenhouse production, reduced energy use has become increasingly important. A more efficient RH control, is important since it represents a principal energy flow and has consequences for growth, production, and quality of the produce. Our study clearly indicates that manipulating g s directly affects energy demand related to RH control. Selecting low-g s cultivars reduces annual energy demand (2.5e5.75% depending on the RH set point; Fig. 14E, F). As economic benefits are considerable due to the high share of energy costs in moderate and cool climate greenhouse production, and since the proposed methodology can be implemented without additional investments, this approach certainly deserves further exploration.
Arguably g s is an unexploited trait in ongoing breeding programs. Breeding programs should consider overall steady state g s value as a parameter critical to cultivation management and costs (Fig. 14). Besides this critical benefit, a growing body of evidence also suggests that decreased g s also contributes to prolonged postharvest longevity Fanourakis et al., 2016bFanourakis et al., , 2020b. However, successful utilisation of a selection trait in breeding programs requires a robust non-invasive methodology, which is suitable for largescale plant analyses. Measurements of leaf and canopy temperature by thermal imaging have been successfully introduced to evaluate overall g s (Kohtaro & Olajumoke, 2020).
Potentially commercially important genetic variation in g s was found in a small set of commercial pot rose cultivars (Fig. 2) as it is evident in a range of other species (Medlyn et al., Table 1). The differences in the x-axes should be noted as different start dates for the different spacings. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) Table 4 e Fitted parameters of the BalleBerry stomatal conductance model (BWB-model) obtained from nonlinear regression with of data of Experiment 4 (see Table  1) with pot rose cultivars 'Flirt' and 'Alaska'. [g 0 , minimum (residual) g s as leaf net photosynthesis rate reaches 0; m s , adjustment factor at 400 mmol mol ¡1 CO 2 (m s ) and as function of variable CO 2 m s [CO2] ]. b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 2011; Merilo et al., 2014). This variation in g s was reasonably persistent across different experiments (Figs. 2,3,5 and 8).

Cultivar
The range in ET c of these cultivars was also similar to that recorded in previous studies (Giday et al., 2013b(Giday et al., , 2015, with the difference between the lowest and highest ET c being about 20% (Fig. 2D). The rankings and the magnitude of cultivar differences in g s and ET c cited here may differ from those that could occur under persistent abiotic stress, such as water deficit  or saline conditions (Hassanvand et al., 2019). However, the cultivar effects on g s in the absence of stress are mostly relevant in protected cultivation. Evidence presented by other studies has also revealed a large genetic variation in g s of rose (Carvalho et al., 2015bFanourakis et al., 2013Fanourakis et al., , 2016bFanourakis et al., , 2020b. The findings of this study together with those of other studies might be taken to indicate that progress in breeding programs for improving this trait may be readily made if selection for low g s could be applied, without the necessity to resort to wild germplasm. Still, it is highly likely that much greater variation in g s could be found with more extensive surveys of germplasm (within and outside) of commercial breeding programs.
Mechanisms resulting in concurrent reduced g s may lead to net reductions in growth and yield at least in some environments, through its negative effect on carbon dioxide intake (Jackson et al., 2016). Thus, prolonged selection pressure for low g s alone may eventually hamper plant growth and biomass. In the experiments reported here or earlier studies (Giday et al., 2013a(Giday et al., , 2015, g s was not correlated with plant biomass.

4.2.
Extended g s prediction model At lower RH levels, prediction quality of the BWB-model decreased. An earlier reported model artefact (e.g. Amitrano et al., 2020) lies in the fact that the steady-state BWB-g s model has been characterised as a correlation rather than a mechanistic equation (Gutschick & Simonneau, 2002) and is not designed for all RH conditions. The simplicity of this model, however, increased its strength in suitability as submodels in larger model systems (K€ orner & Holst, 2017). Two cultivars with a consistently large difference in g s and ET c ['Flirt' (low values) vs 'Alaska' (high values), respectively] across experiments of this and earlier (Giday et al., 2013b(Giday et al., , 2015 studies were selected out of a group of six pot rose cultivars. By using these two cultivars, we have first extended the existing g s prediction model (BWB-model) for rose by incorporating A p or D s as QTL-related input variables and then focussed on the extension with the same traits on the low transpiring cultivar 'Flirt' with the more complex extension of the BWB-model, i.e. the BWB-Liu-model (Liu et al., 2009). The extended model with cultivar specific parameters, as presented here, is an improvement to the BWB-Liu-model in terms of prediction quality. The BWB-Liu model without further extensions has shown Fig. 12 e Observed versus predicted stomatal conductance (g s ) in the cultivar 'Flirt' in Exp. 5 at 30% (A, B) or 70% (C, D) relative air humidity setpoints (RH) by using the BWB model (light blue) and the BWB model with an additional function for individual stomatal pore area with R 2 of 0.275 or 0.268 at 30% RH; 0.765 or 0.766 at 70% RH (magenta) compared to measured g s (A, C; blue stars). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) b i o s y s t e m s e n g i n e e r i n g 2 0 8 ( 2 0 2 1 ) 1 3 1 e1 5 1 improvements of model fit, which is in agreement with earlier findings (e.g. Wei et al., 2018). The extended model introduces only one additional variable, beyond those required in the existing g s prediction model (BWB-Liu-model) (Liu et al., 2009), remarkably further improving g s prediction quality over a range of RH levels.
The ability of the extended model was examined by comparing the estimated values with experimental measurements (Fig. 13). In contrast to substantial deviations of the existing models (BWB-model and BWB-Liu-model) (Ball et al., 1987;Liu et al., 2009), the extended version quoted here accurately reproduced the variation of g s within one cultivar. This improvement of the extended model is achieved with the accommodation of the stomatal feature A p or D s . These two (cultivar-dependent) features were added separately as dynamic parameters, i.e. changing with time. However, while A p is a rapidly changing parameter, D s only changes in the longterm and is determined constant for each leaf. The choice of extension will have important implications for model application in climate regulation of stomata in general (Ball et al., 1987) and its use in modelling water use efficiency of crops in relation to water availability and irrigation (Liu et al., 2009). With the extended BWB-Liu model, ET c errors associated with g s estimations were reduced in comparison to previous models (Fig. 13). While this is at the cost of only a slight increase of model complexity, the additional input variables to the model as the fast dynamics of A p and J or the long-term changing D s could be an obstacle in commercial usage in decision support systems (DSS) or a climate controller.
We have shown that environmental variables influence the stomatal traits and those influences are variety specific. Thus, for practical application in DSS or model-based climate control, the dynamics of these traits need to be incorporated, i.e. the dynamic model of stomatal aperture. As the stimulation of stomatal responsiveness is proportional to the increase in ABA (Giday, 2014), setting it as a major determinant of stomatal regulation (Carvalho, Torre, et al., 2015), ABA-based deterministic models (i.e. without measurements) would undoubtedly improve both model prediction quality and usability. For the BWB-Liu model application, a model or sensor for substrate J is an additional need.
It is worth noting that rose is a hypostomatous species (Fanourakis et al., , 2019b. Next to dynamic modelling of A p and D s , further investigations are thus needed on how D s is parameterised in species with uneven distribution of stomata on the two sides of the leaf. The stomatal aperture (employed to compute A p ) is dynamic and dependent on the daily climate fluctuations, and may therefore be complicated to parameterise under some growth conditions . Contrary to this, D s is only affected by the long-term climate pattern (during plant growth), and with the current development of digital microscope scanners is fairly simple to be quantified (Song et al., 2020). Although extensions with A p result in better improvements of fitting quality, the BWB-Liu model extended by D s may have wider applications in modelling involving g s both in protected cultivation and open field agriculture. However, for the time being, we have chosen to use the BWB-model, which has a reasonably good prediction Fig. 13 e Observed versus predicted stomatal conductance (g s ) by using different models in 'Flirt' in Exp. 6: A, BWB model (R 2 ¼ 0.12); B, BWB-Liu model (R 2 ¼ 0.31); C, modified BWB-Liu model with individual stomatal pore area, A p , (R 2 ¼ 0.52); D, modified BWB-Liu model with stomatal density, D s , (R 2 ¼ 0.59).
quality at the higher RH levels that prevail in protected cultivation of ornamentals, to show the effect of cultivar selection on greenhouse energy household.

Energy saving with cultivar selection
Adjusted RH-based climate control can strongly reduce energy demand in greenhouses (K€ orner & Challa, 2004), while implementing RH-related algorithms in a greenhouse climate controller especially in modern highly insulated greenhouses is the preferred choice (K€ orner & Van Straten, 2008). In our investigations, we created the basis for a cultivar specific model-based climate control system. Simulations have shown that the system is robust for estimating ET c over diverse pot rose cultivars and environmental conditions (Fig. 11e13). The satisfactory agreement between actual data and estimated values indicated the strengths of the extended model as a framework and a reference for validating other approaches of calculating g s . However, we have also shown that model extension would strongly improve prediction of g s . In our validation study, the effect of a low prediction quality of the basic BWB-model was not evident, probably due to the higher RH levels during cultivation, as the low prediction quality comes mainly from situations at low humidity (<60%). For conditions with low RH, a model improvement needs to be done. A possible starting point is incorporating dynamic stomatal behaviour, as suggested by Vialet-Chabrand et al. (2013). While these dynamics were introduced with mathematical measures, to understand and apply future dynamic g s models in a wide range of varieties and crops, deterministic approaches as discussed above with deterministic modelling of the role of ABA in g s regulation would be preferable. However, for any type of model for energy saving and crop-based climate control, crop and cultivar specific parameterisation of the underlying models is of utmost importance.

Conclusions
The RH control can currently account for over 20% of greenhouse energy demand. Our aim was to reduce this share by considering cultivar g s differences in estimating ET c . Large differences in ET c (>20%) were noted among pot rose cultivars, which together with environmentally-induced variation resulted to a wide ET c range. This range was used to adapt and extend the BWB-Liu model, parameterise, calibrate and validate the underlying parameters and finally implement it in a greenhouse simulator. The modified g s -model version (considering D s ) largely improved g s estimation (R 2 of 0.59 vs 0.31) at RHs lower than 60%. A very good agreement between simulated and measured ET c was achieved. Depending on the RH set point, selecting low-transpiring cultivars may save between 2.5 and 5.75% of energy demand on annual basis. Incorporating the presented cultivar specific model in decision support tools and climate control would enable model-based controllers to adjust climate more efficiently (with real-time DSS), and create ground for planning tools in greenhouse construction (in association with greenhouse planning tools such as Virtual Greenhouse).