Introduction

Using greenhouse agriculture, it is relatively easy to generate higher crop yields than outdoor-field agriculture. The enclosure makes it possible to improve the microclimate and photosynthetic rate via, for instance, CO2 enrichment, supplemental lighting, and ventilation. However, greenhouse yields can have high spatial heterogeneity, due to the heterogeneous microclimate and spatial variability of photosynthetic rates (Kimura et al., 2020b). Greenhouse structures interrupt sunlight, causing spatially heterogeneous incident light (measured as photosynthetic photon flux density, PPFD) (Cabrera-Bosquet et al., 2016; Kozai & Kimura, 1977; Stanhill et al., 1973). This spatial variability varies considerably diurnally and seasonally because of changes in sun position and external light intensity. Variation in PPFD thus causes spatiotemporal variation in photosynthesis, and in the resulting crop yield.

Assessing the spatiotemporal variability of PPFD in a greenhouse is highly challenging. Multipoint measurement using many photon sensors is costly, whereas mobile measurement using a single sensor is time-consuming and ignores fine-scale temporal changes in PPFD. As an alternative, numerical models of the effects of greenhouse architecture and solar trajectory on PPFD have been developed (e.g., Castellano et al., 2016; Cossu et al., 2017; Kozai & Kimura, 1977; Matsuda et al., 2020). Such models can accurately describe sunlight distribution in a greenhouse, but the simulated greenhouse geometry must be revised when applied to other greenhouses. Moreover, it is difficult to completely model the features of greenhouse architecture and instrumentation, which can include thin wall pipes, heaters, a CO2 generator, and circulation fans, and obstacles outside the greenhouse, such as mountains and tall trees.

However, multipoint analysis of hemispheric images (180° images from below of the greenhouse structure, produced using a fisheye camera), can be used to estimate PPFD. From these images, direct and diffuse PPFD can be calculated by analyzing the solar trajectory and the gap fraction (Steege, 1993, 2018). This technique can describe the structure of the greenhouse and external obstacles in detail; further, because it estimates the movement of the sun across the structure, it can estimate continuous changes in PPFD for a point of interest, from a single photograph (although if the greenhouse or external obstacles change, a new image is required). Hemispheric photography is widely applied in forest science to estimate incident light (Jonas et al., 2020; Schleppi & Paquette, 2017). Although this method was used in a glasshouse agriculture study (Cabrera-Bosquet et al., 2016), it has not yet been used in other indoor crop field studies.

It is more difficult to measure spatiotemporal variation in photosynthesis than in PPFD. Photosynthetic rates are often measured using open- or closed-chamber systems (Burkart et al., 2007; Davis et al., 1987; McDermitt et al., 1989; Song et al., 2016). However, these systems, which were developed for plot-sized experiments, are not suitable for multipoint measurement. Up to now, model simulation has been the sole approach for fine-resolution spatiotemporal analysis of photosynthesis. To estimate photosynthetic rate while accounting for microclimate, ecological and physiological studies routinely use models of C3 photosynthesis (Farquhar et al., 1980), stomatal conductance, CO2 diffusion into the leaf, and leaf energy budget.

Despite its importance in greenhouse agriculture, fine-scale spatiotemporal determination of variability of PPFD and photosynthesis is difficult to achieve. Further, the effects of this variability on yield remain unclear. Therefore, the objective of this study was to evaluate a novel approach for estimating spatiotemporal variability in PPFD and leaf photosynthetic rate in a greenhouse. This approach uses hemispheric image multipoint analysis to estimate PPFD, and a numerical model to estimate leaf photosynthetic rate. These fine-resolution spatiotemporal estimates are then assessed in relation to strawberry yield in a greenhouse, to determine how spatiotemporal variability in PPFD and leaf photosynthetic rate affects yield.

Materials and methods

Analysis of hemispheric images and PPFD

PPFD at the top of the plant canopy, at a point of interest in a greenhouse, is expressed as follows:

$$ {\text{PPFD}} = \tau (p_{{{\text{dir}}}} { } \times {\text{PPFD}}_{{{\text{ext}},{\text{dir}}}} + { }p_{{{\text{dif}}}} \times {\text{PPFD}}_{{{\text{ext}},{\text{dif}}}} ) $$
(1)

where τ is the transmissivity of the greenhouse covering material (here, τ = 0.78, for a polyvinyl-chloride sheet), assumed to be constant in space and time; and pdir and pdif are the proportions of incident to external PPFD, for direct (PPFDext,dir) and diffuse (PPFDext,dif) light, respectively (these depend on the architecture and obstacles in and around the greenhouse). PPFDext,dir and PPFDext,dif were separated from PPFDext, which was measured outside the greenhouse, using the method proposed by Spitters et al. (1986). This method is based on the ratio between measured solar radiation (Rs,ext) outside the greenhouse and calculated radiation outside the atmosphere.

To evaluate pdir and pdif, a hemispheric image (field of view of 180° or greater) (Fig. 1a) was analyzed (Cabrera-Bosquet et al., 2016). The image was binarized using the machine-learning-based software ‘ilastik’ (Berg et al., 2019), which learns obstacles and sky areas when these are manually drawn on the image, and realistically and flexibly classifies images (Fig. 1b). To calculate pdir, the position of the sun in the binarized image was determined for the time of interest (Fig. 1c) by matching the azimuth and elevational angles of the sun—calculated from astronomical formulae (Meeus, 1998; NOAA Solar Calculator, 2018)—with the angles in the image. Using solar position in the image, pdir was calculated:

$$ p_{{{\text{dir}}}} = \frac{{{\text{pixel}}_{{{\text{sun}}}} - {\text{pixel}}_{{{\text{obs}}}} }}{{{\text{pixel}}_{{{\text{sun}}}} }} $$
(2)
Fig. 1
figure 1

Schematic of the analysis of a hemispheric image of the experimental greenhouse, seen from below. A raw hemispheric image (a) is binarized using machine learning-based classification software (b). The image is then used to calculate the proportions, relative to external light, of incident direct (c) and diffuse (d) light. Section 2.1 describes this calculation

Full size image

where pixelsun is the number of pixels of the sun on the image, and pixelobs is the number of pixels of the obstacles overlapping the sun (that is, pdir is 1 when all pixels of the sun are on white areas on the binarized image, and 0 when the sun is behind black areas). The diameter of the sun in the image was calculated as follows:

$$ d_{{{\text{sun}}}} = N \times 0.53 $$
(3)

where dsun is the diameter (in pixels) of the sun, N is the number of pixels per 1° in the image, and 0.53 is the angular diameter (°) of the sun from the Earth. pdir varies spatiotemporally, because of spatial variation caused by obstacles inside and outside the greenhouse, and diurnal and seasonal variation in solar position.

To calculate pdif, 90 concentric rings were overlaid on the hemispheric area (indicating elevational angles of up to 90°) of the binarized image (Fig. 1d), generating 89 annuli per 1°. In each annulus, the ratio of its area without obstacles to its total area (fα, corrected by three terms, c1, c2, and c3) was calculated, and these ratios were summed to obtain pdif (Steege, 2018):

$$ p_{{{\text{dif}}}} = \mathop \sum \limits_{\alpha = 0.5}^{\alpha = 89.5} \left( { f_{\alpha } c_{\alpha ,1} c_{\alpha ,2} c_{\alpha ,3} } \right) $$
(4)
$$ f_{\alpha } = \frac{{{\text{pixel}}_{{{\text{tot}},\alpha }} - {\text{pixel}}_{{{\text{obs}},\alpha }} }}{{{\text{pixel}}_{{{\text{tot}},\alpha }} }} $$
(5)
$$ c_{\alpha ,1} = \sin \left( {\alpha + 0.5} \right) - \sin \left( {\alpha - 0.5} \right) $$
(6)
$$ c_{\alpha ,2} = \frac{1 + 2\sin \alpha }{3} $$
(7)
$$ c_{\alpha ,3} = \sin \alpha $$
(8)

where α is the elevational angle in the image (ranging from 0.5° to 89.5°, at 1° intervals along each annulus), pixeltot,α is the total number of pixels in each annulus, and pixelobs,α is the number of obstacle pixels in each annulus. cα,1 corrects the apparent area on the image to the actual area, using equidistant projection. cα,2 corrects the homogeneous distribution of diffuse light to a heterogeneous distribution, under the Standard Overcast Sky assumption that the sky is three times as bright at the zenith as near the horizon. cα,3 represents a cosine correction. pdif varies spatially, due to obstacles inside and outside the greenhouse, but not temporally.

Evaluation of leaf photosynthesis

Leaf photosynthetic rate (A) was estimated from the biochemical C3 photosynthesis model (Farquhar et al., 1980) integrated with leaf CO2 diffusion theory, stomatal conductance model (Medlyn et al., 2011), and leaf energy budget (Jones, 2013). In the photosynthesis model, A is limited by the rates of rubisco activity (Ac) or RuBP regeneration (Aj). The triose phosphate utilization (TPU)-limited rate was not considered here, because a TPU-limited phase was not observed via gas exchange measurements (described below) within the experimental CO2 concentrations in this study. Actual A is given by

$$ \theta_{{A}} A^{2} - A\left( {A_{{\text{c}}} + A_{{\text{j}}} } \right) + A_{{\text{c}}} A_{{\text{j}}} = 0 $$
(9)
$$ A_{{\text{c}}} = \frac{{V_{{{\text{c}}\max }} (C_{{\text{c}}} - \varGamma^{ *} )}}{{C_{{\text{c}}} + K_{{\text{c}}} (1 + O/K_{{\text{o}}} )}} - R_{{\text{d}}} $$
(10)
$$ A_{{\text{j}}} = \frac{{J (C_{{\text{c}}} - \varGamma^{*} )}}{{4C_{{\text{c}}} + 8\varGamma^{*} }} - R_{{\text{d}}} $$
(11)
$$ J = \frac{{\phi {\text{PPFD}} + J_{{{\text{high}}}} - \left\{ {\left( {\phi {\text{PPFD}} + J_{{{\text{high}}}} } \right)^{2} - 4\phi {\text{PPFD}}\theta_{J} J_{{{\text{high}}}} } \right\}^{1/2} }}{{2\theta_{{J}} }} $$
(12)

where θA is the curvature in the transition from one limitation to the other, Vcmax is the maximum carboxylation rate, Cc is the chloroplast CO2 concentration, Γ* is the CO2 compensation point in the absence of respiration, Kc and Ko are the Michaelis constants for carboxylation and oxygenation, respectively, O is the O2 concentration, Rd is the daytime respiration rate, J is the electron transport rate, Jhigh is the maximum rate of electron transport at high light intensity (Buckley & Diaz-Espejo, 2015), ϕ is the initial slope of the curve (the apparent quantum yield of electron transport under low light), and θJ is the convexity of the curve. Rd was estimated using its relationship to the dark respiration rate Rn (Villar et al., 1995):

$$ R_{{\text{d}}} = R_{{\text{n}}} ({\text{if}}\,{\text{PPFD }} = 0\, \mu {\text{mol m}}^{ - 2} {\text{s}}^{ - 1} ) $$
(13)
$$ R_{{\text{d}}} = 0.4R_{{\text{n}}} ({\text{if}}\,{\text{PPFD}} > 0 \mu {\text{mol}}\,{\text{m}}^{ - 2} \,{\text{s}}^{ - 1} ) $$
(14)

Cc was described by diffusion theory (Fick’s law):

$$ C_{{\text{c}}} = C_{{\text{a}}} - A (1.27g_{{{\text{ah}}}}^{ - 1} + 1.56g_{{{\text{sw}}}}^{ - 1} + g_{{\text{m}}}^{ - 1} ) = C_{{\text{s}}} - A (1.56g_{{{\text{sw}}}}^{ - 1} + g_{{\text{m}}}^{ - 1} ) $$
(15)

where Ca and Cs are the CO2 concentrations in the ambient air and at the leaf surface, respectively, gah is the leaf boundary layer conductance for heat transfer, gsw is the stomatal conductance for water vapor transfer, and gm is the mesophyll conductance. gsw was estimated from the model proposed by Medlyn et al. (2011):

$$ g_{{{\text{sw}}}} = g_{0} + 1.6 \left( {1 + \frac{{g_{1} }}{{\sqrt {VPD} }}} \right)\frac{A}{{C_{{\text{s}}} }} $$
(16)

where g0 is the residual conductance, g1 is the empirical constant, and VPD is the leaf–air vapor pressure deficit. gah was estimated from wind speed u and leaf characteristic dimension dleaf using a semi-empirical relationship (Kimura et al., 2020b):

$$ g_{{{\text{ah}}}} = 0.21\left( {u/d_{{{\text{leaf}}}} } \right)^{0.5} $$
(17)

where the constant 0.21 was determined using the relationship between u, measured by anemometer, and gah, evaluated using the artificial leaf technique (Kimura et al., 2020a).

The model parameters (Γ*, Kc, Ko, Rn, Vcmax, Jhigh, and gm) vary with leaf temperature Tl, and the temperature dependence was described by the Arrhenius function or the peaked Arrhenius function:

$$ k_{{T_{l} }} = k_{25} \exp \left[ {\frac{{E_{{\text{a}}} \left( {T_{{{\text{l}},{\text{K}}}} - 298.15} \right)}}{{298.15 R T_{{{\text{l}},{\text{K}}}} }}} \right] \quad \text{for} \,\varGamma ^{*} , \,K_{{\text{c}}} , \,K_{{\text{o}}} , \,R_{{\text{n}}} $$
(18)
$$ k_{{T_{l} }} = k_{25} \exp \left[ {\frac{{E_{{\text{a}}} \left( {T_{{\text{l,K}}} - 298.15} \right)}}{{298.15 R T_{{\text{l,K}}} }}} \right]\frac{{1 + \exp \left( {\frac{{298.15\Delta S - H_{{\text{d}}} }}{298.15 R}} \right)}}{{1 + \exp \left( {\frac{{T_{{\text{l,K}}} \Delta S - H_{{\text{d}}} }}{{T_{{\text{l,K}}} R}}} \right)}}\quad \text{for}\, V_{{{\text{cmax}}}} , \,J_{{{\text{high}}}} , \,g_{{\text{m}}} $$
(19)

where \(k_{{T_{l} }}\) is the parameter value at a given leaf temperature, Tl,K is Tl in K, k25 is the parameter value at 25 °C, R is the universal gas constant, Ea is the activation energy, Hd is the deactivation energy, and ∆S is the entropy factor.

Tl was estimated from leaf energy budget assuming isothermal net radiation (Jones, 2013), modified to molar units (Buckley et al., 2014):

$$ T_{{\text{l}}} = T_{{\text{a}}} + \frac{{ \gamma R_{{\text{ni }}} / c_{{p}} - D_{{\text{a}}} (0.93g_{{{\text{ah}}}}^{ - 1} + g_{{{\text{sw}}}}^{ - 1} )^{ - 1} }}{{s(0.93g_{{{\text{ah}}}}^{ - 1} + g_{{{\text{sw}}}}^{ - 1} )^{ - 1} + \gamma (g_{{{\text{ah}}}} + 4\varepsilon_{{{\text{l}},{\text{leaf}}}} \sigma T_{{{\text{a}},{\text{K}}}}^{ 3} / c_{{\text{p}}} )}} $$
(20)

where Ta is the ambient air temperature, γ is the psychrometric constant, Rni is the isothermal net radiation, cp is the molar heat capacity of air, Da is the saturation vapor pressure deficit of air, s is the slope of the curve relating saturation water vapor pressure to temperature, εl,leaf is the leaf emissivity to longwave radiation, σ is the Stefan–Boltzmann constant, Ta,K is Ta in Kelvin. Rni was estimated from PPFD, and is given by:

$$ R_{{\text{ni }}} = \alpha_{{{\text{s}},{\text{leaf}}}} \times 0.495 \times {\text{PPFD}} + R_{{{\text{l}},{\text{abs}}}} - \varepsilon_{{{\text{l}},{\text{leaf}}}} \sigma T_{{{\text{a}},{\text{K}}}}^{ 4} $$
(21)

where αs,leaf is the leaf absorptivity of shortwave radiation, 0.495 is the ratio of total shortwave energy to photosynthetic photon flux (Mavi & Tupper, 2004), and Rl,abs is the absorbed longwave radiation. Rl,abs was described as the sum of three terms (De Boeck et al., 2012; Kimura et al., 2020b):

$$ R_{{{\text{l}},{\text{abs}}}} = \alpha_{{{\text{l}},{\text{leaf }}}} \tau_{{{\text{l}},{\text{gh }}}} \varepsilon_{{{\text{l}},{\text{atm}}}} \sigma T_{{{\text{a}},{\text{K}}}}^{ 4} + \alpha_{{{\text{l}},{\text{leaf }}}} \varepsilon_{{{\text{l}},{\text{gh}}}} \sigma T_{{{\text{gh}},{\text{K }}}}^{ 4} + \alpha_{{{\text{l}},{\text{leaf }}}} \rho_{{{\text{l}},{\text{gh }}}} \varepsilon_{{{\text{l}},{\text{leaf}}}} \sigma T_{{{\text{a}},{\text{K}}}}^{ 4} $$
(22)

where αl,leaf is the leaf absorptivity of longwave radiation; εl,atm is the atmospheric emissivity of longwave radiation; Tgh,K is the greenhouse temperature in Kelvin, assumed to equal Ta,K (De Boeck et al., 2012); and τl,gh, εl,gh, and ρl,gh are the longwave radiation transmissivity, emissivity, and reflectivity of the greenhouse material, respectively, using typical values for a polyvinyl-chloride greenhouse (De Boeck et al., 2012). εl,atm was given by Leuning et al., (1995):

$$ \varepsilon_{{{\text{l}},{\text{atm}}}} = 0.642\left( {0.001P_{{{\text{atm}}}} W_{{\text{a}}} /T_{{{\text{a}},{\text{K}}}} } \right)^{1/7} $$
(23)

where Patm is the total atmospheric pressure, and Wa is the water vapor concentration in the ambient air.

A can be obtained, together with other unknown variables (Cc, Cs, gsw, and Tl), by substituting the model parameters (Table 1) into Eqs. (9) to (23), and by iteratively solving these equations using a binary search (Gutschick, 2016; Kimura et al., 2020b).

Table 1 Model parameters for determining leaf photosynthetic rate
Full size table

Experimental greenhouse and plant materials

These methods were applied to a polyvinyl-chloride greenhouse, oriented NW–SE (32° westerly declination), in Fukuoka, Japan (33° 36′ 43″ N, 130° 13′ 56″ E, 9 m amsl). The greenhouse was 36 m wide, and 51 m and 32 m long on the western and eastern sides, respectively, 2.3 m high at the gutter, and 4.0 m high at the top, having six roofs (Fig. 2). In the greenhouse, 30 soil ridges were constructed, and strawberry plants (Fragaria × ananassa ‘Fukuoka S6’) were planted on October 14, 2018, in two rows on each ridge. The plants were cultivated until May 29, 2019, and the number of leaves on each plant was maintained at ca. 10 during the cultivation period, wherever possible.

Fig. 2
figure 2

Schematic of the structure of the experimental six-span greenhouse in which strawberry plants were cultivated. Various environmental controls, providing roof ventilation, sidewall ventilation, CO2 enrichment, heating, and air circulation, were installed in the greenhouse

Full size image

Several environmental control methods were used in the greenhouse (Fig. 2). Roof vents were installed on both sides of each gutter, with maximum opening areas of 80.0 m2 (1.6 m × 50.0 m) in the three vents on the western side and 49.6 m2 (1.6 m × 31.0 m) in the two vents on the eastern side. Sidewall vents were installed on both sides of the greenhouse, with maximum opening areas of 90.0 m2 (1.8 m × 50.0 m) on the western side and 55.8 m2 (1.8 m × 31.0 m) on the eastern side. These ventilation systems were operated to maintain greenhouse air temperature below 25 °C. A CO2 generator (CG-554T2, Nepon Inc., Tokyo, Japan), supplying CO2 at 8 kg h−1 and airflow at 33 m3 min−1, was situated near the greenhouse entrance at the southern side. The generator was operated in the daytime to maintain CO2 concentration in the greenhouse above 400 μmol mol−1. An oil-fueled heater (HK5027TFV, Nepon Inc.), with a thermal output of 145 kW, was situated near the CO2 generator. The heater was operated in the nighttime to maintain greenhouse air temperature above 8 °C. Circulation fans (Furaibou II, Nichinou Industrial Co., Ltd, Fukuoka, Japan), supplying airflow at 85 m3 min−1, were installed in each span of the greenhouse, and were operated to evenly distribute the greenhouse microclimate.

Measurement of hemispheric images, microclimate, model parameters, and fruit yield

To assess the fine-resolution spatiotemporal distribution of PPFD, hemispheric images were taken at 393 points at regular intervals in the greenhouse, before transplanting, using a fisheye camera (EX-FR200, Casio, Tokyo, Japan) with a 185° field of view. The photography was conducted just after sunrise or just before sunset, requiring six days to photograph all points. The images were taken 400 mm above the soil ridges, corresponding to the top of the plant canopy; the camera was leveled, and geographical north was checked for every photography. PPFDext and Rs,ext were measured during the cultivation period outside the greenhouse at 4.3 m above ground, using a quantum sensor (PQS1, Kipp & Zonen, Delft, Netherlands) and a pyranometer (CMP6, Kipp & Zonen), respectively. PPFD was calculated at the 393 points at intervals of 1 min, using a self-made program in R. To validate the PPFD estimates, PPFD was measured at five points in the greenhouse using quantum sensors (PAR-02DS, Prede, Tokyo, Japan).

The microclimatic parameters (excepting PPFD) required for calculating A were measured 60 cm above ground, at the center of the greenhouse, during the cultivation period. Ta and relative humidity were measured by a thin-film capacitive sensor (HMP110, Vaisala, Helsinki, Finland) with a forced-ventilation system. Ca and u were measured by a non-dispersive infrared sensor (GMP343, Vaisala) and an omnidirectional anemometer (0965-00, Kanomax, Tokyo, Japan), respectively. Microclimatic data were recorded using a data logger (CR1000, Campbell Scientific, Logan, USA) every 1 min. Under the assumption that the microclimate (excepting PPFD) is spatially uniform across the greenhouse, the spatial distribution of A is affected by PPFD, via its direct and indirect effects on biochemical processes and on the leaf energy budget.

To determine Vcmax, Jhigh, gm, Rn, and their temperature dependencies, CO2 response curves of leaf photosynthesis and chlorophyll fluorescence were measured on 16 days during the experiment, using a portable open gas exchange system (LI-6400XT, LI-COR, Nebraska, USA), using a leaf chamber fluorometer (6400-40, LI-COR, Nebraska, USA). On each day, fully expanded leaves (n = 3) were randomly chosen, and leaf photosynthesis CO2 response curves were created under a PPFD of 1500 μmol m−2 s−1, using CO2 concentrations that were changed in a stepwise fashion (400, 300, 200, 100, 50, 400, 600, 800, 1000, 1200, and 1600 μmol mol−1), at different leaf temperatures (20, 25, 30, and 35 °C). During the measurements, the leaf–air vapor pressure deficit in the chamber was kept below 2.3 kPa. For each CO2 response curve, the curve-fitting program in the “plantecophys” R package (Duursma, 2015) was applied to estimate Vcmax and Jhigh. Chlorophyll fluorescence was measured at the same time as the response curves were generated. To estimate gm, the variable J method (Harley et al., 1992) was applied to these data to keep the intercellular CO2 concentration (Ci) of the curves within the range 240 to 600 μmol mol−1 (Xue et al., 2016). After each CO2 response curve measurement, Rn was measured at 0 μmol m−2 s−1 PPFD and 400 μmol mol−1 CO2 concentration. The temperature dependences of Vcmax, Jhigh, gm, and Rn was determined using a self-made fitting program in R. To determine g0 and g1, diurnal changes in A, gsw, VPD, and Cs were measured on 5 days during the experiment, using a LI-6400 device with a transparent top chamber (standard leaf chamber; LI-COR, Nebraska, USA) under ambient light, temperature, and humidity. On each day, fully expanded leaves (two or three) were randomly chosen, and the gas exchange measurements were conducted while changing leaves at intervals of 20 to 30 min. The fitting program in the “plantecophys” R package (Duursma, 2015) was applied to estimate g0 and g1.

Fruit yield was measured for 30 randomly selected plots, in the greenhouse at intervals of 2–3 days from December 4, 2018, to May 29, 2019. The fruits from 10 plants in each plot were harvested, then fresh weight was measured and averaged across the 10 plants at each plot for subsequent analysis.

Spatiotemporal and statistical analysis

Spatial variability in PPFD, A, and fruit yield was defined as

$$ {\text{Spatial variability}} = \frac{Max - Min}{{Mean}} \times 100 $$
(24)

where Max, Min, and Mean are the maximum, minimum, and mean values, respectively, at 393 points (PPFD and A) or 30 plots (fruit yield).

Spatial bias in PPFD and A was investigated via spatial autocorrelation analysis, using Moran’s I test. Moran’s I ranges from − 1 to 1, and negative and positive values tend to indicate negative and positive spatial autocorrelation (i.e., spatial dispersion and clustering), respectively. Values near 0 tend to show no spatial autocorrelation (i.e., spatial randomness).

To evaluate how the spatiotemporal variability of fruit yield was explained by PPFD and A, time courses for the coefficients of determination (R2) of the relationships between accumulated fruit yield, accumulated PPFD, and A until each harvest day were determined for the 30 plots where the fruit yield was measured.

Results

Incident light and leaf photosynthesis

The PPFD estimates based on hemispheric image analysis were validated at five points in the greenhouse (Fig. 3). The hemispheric images were successfully binarized via machine learning, on both cloudy and sunny days (rows 1 and 2, Fig. 3). Hemispheric image analysis accurately predicted PPFD, relative to the photon sensor measurements (row 3, Fig. 3), including capturing sudden reductions in direct light caused by thick beams (row 3, Fig. 3b, d). Owing to the thick beam, PPFD was lower under the gutter (Fig. 3b, d) than under the apex of the greenhouse (Fig. 3a, c, e). Hemispheric image analysis accurately predicted PPFD on most cloudy and clear days (row 4, Fig. 3), but underestimated it on a few days (Fig. 3c, e). The root mean square error (RMSE) between the predicted and observed 30 min averaged values was 47.1 μmol m−2 s−1, and relative RMSE was 4.4% of the mean of the five points.

Fig. 3
figure 3

Validation of incident light estimates derived from hemispheric images of the greenhouse (from below), at (a) the westernmost part of the six-span greenhouse, (b) a gutter near the heater and CO2 generator, (c) the roof of the greenhouse, (d) a gutter away from instruments, and (e) below the apex of the greenhouse. Rows 1 and 2: raw and binarized images, respectively. Row 3: representative time courses (30 min averaged values) of photosynthetic photon flux density (PPFD) measured (red line) and predicted (yellow line) in the greenhouse and measured outside the greenhouse (dashed line). Row 4: predicted versus observed PPFD (30 min averaged values) in the greenhouse

Full size image

The numerical model predictions of A were validated against changes in A measured in the greenhouse (Fig. 4). The photosynthesis model accurately predicted A, with RMSE of 1.57 μmol m−2 s−1 and relative RMSE of 11%, although with some scatter.

Fig. 4
figure 4

Predicted versus observed leaf photosynthetic rate (A) in strawberry leaves

Full size image

Spatiotemporal variability in incident light and leaf photosynthesis

Diurnal variability

PPFD was spatially clustered at low solar elevations (Fig. 5a, e; Moran’s I, Table 2): PPFD was higher on the eastern side at 8:00 and on the western side at 17:00, corresponding to the direction of the sun. As the solar elevational angle approached its maximum in the day, the spatial distribution of PPFD showed a clearly striped pattern, due to shading by the gutter (Fig. 5c). The stripe position varied with the solar position, with the shade moving from west to east, contrary to the solar trajectory (Fig. 5b–d). The spatial distribution of A was similar to that of PPFD, but with less variability, particularly at high solar elevational angles (Fig. 5b–d and Table 2) because A becomes saturated at high PPFD values.

Fig. 5
figure 5

Diurnal changes in the spatial distributions of photosynthetic photon flux density (PPFD) and leaf photosynthetic rate (A) in the greenhouse. Columns (a–e): relative distributions averaged for the previous 30 min, at the respective solar elevational angles (α), on a representative clear day (March 8, 2019)

Full size image
Table 2 Maximum, minimum, and mean values of photosynthetic photon flux density (PPFD) and leaf photosynthetic rate (A), at the 393 points in the greenhouse, at the respective solar elevational angles (α), on a representative clear day (March 8, 2019, as in Fig. 5)
Full size table

Seasonal variability

Daily accumulated PPFD showed typical seasonal changes, decreasing with the approach of winter, and increasing with that of spring (Fig. 6a). Daily accumulated A showed similar seasonal changes, although it decreased in April and May (Fig. 6b) because of their excessive Tl (exceeding ca. 30 °C), which reduces photosynthetic capacity and stomatal conductance.

Fig. 6
figure 6

Time course of (a) daily accumulated photosynthetic photon flux density (PPFD) and (b) leaf photosynthetic rate (A) at all 393 measurement points in the hemispheric images, during the strawberry cultivation period in the greenhouse

Full size image

Spatial variability in daily-accumulated PPFD increased with daily-accumulated PPFDext, being high on sunny days and low on cloudy days (Fig. 7a). The rate of the increase depended on the season: spatial variability depended on external light intensity more in winter and less in spring (Fig. 7a). Spatial variability of daily-accumulated PPFD ranged from 25 and 56%.

Fig. 7
figure 7

Relationships between spatial variability and Moran’s I of (a, c) daily accumulated photosynthetic photon flux density (PPFD) and (b, d) daily accumulated leaf photosynthetic rate (A) in the greenhouse, with daily accumulated external PPFD (PPFDext) measured outside the greenhouse, from October to November, December to January, February to March, and April to May. The Moran’s I represents significance at P < 0.001

Full size image

In contrast to daily-accumulated PPFD, spatial variability of daily accumulated A increased with daily accumulated PPFDext only in winter (December to January) (Fig. 7b). Again, this is because A becomes saturated at high PPFD. Spatial variability of daily accumulated A increased suddenly in spring (April to May) (Fig. 7b) because the strong light intensity resulted in excessive Tl, causing A to decline. Nevertheless, spatial variability of daily accumulated A remained at 33% (Fig. 7b), markedly below that of PPFD.

For daily accumulated PPFD, Moran’s I showed large scatter at daily accumulated PPFDext < 30 mol m−2 day−1 (Fig. 7c), indicating that clustering of daily accumulated PPFD in the greenhouse was independent of external light intensity and season; in contrast, it converged at higher daily accumulated PPFDext (Fig. 7c). Similarly, for daily accumulated A, Moran’s I showed large scatter, but it decreased suddenly at higher daily accumulated PPFDext (Fig. 7d), indicating that the spatial distribution of daily accumulated A tended to be random under high external light intensity.

Variability during the cultivation period

The spatial distribution of accumulated PPFD over the cultivation period showed a striped pattern along the length of the greenhouse: accumulated PPFD was higher beneath the apex of the greenhouse, and lower below the gutter (Fig. 8a). Accumulated PPFD was higher on the western side, due to the lack of obstacles around the greenhouse on this side, and lower around the heater and CO2 generator in the greenhouse. Accumulated PPFD changed by up to 22% over distances < 2 m.

Fig. 8
figure 8

Spatial distributions of (a) photosynthetic photon flux density (PPFD) and (b) leaf photosynthetic rate (A) accumulated over the cultivation period (from October 14, 2018, to May 28, 2019)

Full size image

The spatial distribution of accumulated A over the cultivation period was similar to that of accumulated PPFD, whereas its spatial variability was half that of accumulated PPFD (Fig. 8b). Accumulated A varied by up to 6% over distances < 2 m.

Effects of spatiotemporal variability of PPFD and A on fruit yield

Daily fruit yield showed three peaks, in the first (December), second (March), and third (May) trusses, peaking in the second truss (Fig. 9a). Total accumulated fruit yield was 721 g/plant (mean for all harvest points) (Fig. 9b). The spatial variability of accumulated fruit yield decreased from more than 100% at the start of the harvest to ca. 50% at the end of the harvest (Fig. 9b). The time course of accumulated daily fruit yield is shown in Fig. 9a. Matching this time course, the R2 values of the relationships between accumulated PPFD and A and accumulated daily fruit yield show three peaks (Fig. 10a–c). R2 values were the highest in winter (the first truss) and declined markedly over time (Fig. 10a–c). Interpreting these R2 values, fruit yield was explained better by spatial variability in accumulated A than in accumulated PPFD over the cultivation period: 46%, 16%, and 12% of the spatial variability in fruit yield was explained by the variability in accumulated A at these peaks; in contrast, the corresponding values for accumulated PPFD are 41%, 7%, and 3% (Fig. 10a–c).

Fig. 9
figure 9

Time course and spatial variability of (a) daily fruit yield (mean values, 30 harvest plots), (b) accumulated fruit yield (mean values, 30 harvest plots), from the start (December 4, 2018) to the end (May 29, 2019) of the harvest

Full size image
Fig. 10
figure 10

Relationships between accumulated photosynthetic photon flux density (PPFD) and leaf photosynthetic rate (A) and accumulated fruit yield, over the cultivation period. Row 1: time courses of the coefficients of determination (R.2) of the relationships between accumulated PPFD and A, and accumulated daily fruit yield. Rows 2 and 3: relationships between accumulated PPFD and A and accumulated fruit yield up to (a) December 19, 2018, (b) March 24, 2019, and (c) May 29, 2019, for 30 harvest plots

Full size image

Discussion

Validation of hemispheric image analysis and numerical model predictions

Hemispheric photography can generate fine-spatial-resolution maps of incident light in greenhouses. In this study, a spatial resolution of approximately 4 m2 (2 × 2 m) was achieved via multipoint photography before transplanting and using one-point measurement of external light conditions (PPFDext and Rs,ext). This is more efficient than using multiple photon sensors, which would require ca. 400 sensors.

Further, hemispheric photograph makes it possible to assess the shading effects of internal and external obstacles. Structural shading causes high spatial variability of light intensity in indoor agriculture (Cabrera-Bosquet et al., 2016; Kozai & Kimura, 1977; Stanhill et al., 1973). Shading by mountains strongly affects spatiotemporal changes in light intensity (Oliphant et al., 2003), particularly around dawn and dusk. Further, tall buildings can reduce PPFD in a greenhouse (Hidaka et al., 2017).

Although numerical simulation is often used to estimate greenhouse light conditions, it requires many parameters to estimate the effects of shading, and they are costly to obtain. In contrast, hemispheric photography requires only the time to take the images (up to a few minutes per image). Further, its accuracy compares favorably with that of numerical simulation, with relative RMSE of 3.7% to 5.0% at the five measurement points used in this study, and 1.0% to 16.9% at the 10 measurement points in the simulation conducted by Cossu et al. (2017).

Although hemispheric photography successfully captured spatiotemporal changes in the incident light in the greenhouse, it has some drawbacks. Repeat photography is necessary when new instruments and buildings are installed inside or outside the greenhouse. Automatization of photography can increase the frequency of image collection and overcome this problem. For example, Teitel et al. (2012) moved a photoradiometer along the greenhouse span to evaluate light distributions in multi-span greenhouses. Such an approach using mobile systems can also be applied to automatize the hemispherical photography. Further, opening and closing the roof vent introduces errors into incident light predictions. When the roof vent is opened, sun directly penetrates the greenhouse (transmissivity τ = 1; Eq. 1), whereas the beam is intercepted by the roof film when the roof vent is closed (τ = film transmissivity). Thus, τ varies depending on the opening or closing of the roof vent. In this study, τ was assumed to be constant (0.78, assuming closed vents). Therefore, on some sunny days when the roof vent was opened, PPFD was underestimated, particularly below the gutter where the area of the roof vent accounts for large proportion of the hemispheric image. To avoid this problem, the roof vent must be detected in the hemispheric images, and τ must be synchronized with the opening and closing of the roof vent. The underestimation was not substantial in this study. Nonetheless, applying the correction would enable more accurate estimation of the variability of incident light.

The use of models to estimate A provided accurate but slightly scattered estimates. This scatter is attributable to the parameters for photosynthetic capacity (Vcmax and Jhigh) and stomatal characteristics (g0 and g1) used in the model. Vcmax and Jhigh vary seasonally due to leaf senescence and environmental acclimation (Wilson et al., 2001), as do g0 and g1 (Ono et al., 2013), and these parameters vary between leaves. Although Vcmax and Jhigh for 16 days from November to April were determined in this study, no remarkable seasonal variation was seen. Thus, the average values over the period for the experiment were used. To measure or model spatiotemporal changes in these model parameters, numerous data points must be obtained, requiring a higher-throughput method such as leaf reflectance spectroscopy (Serbin et al., 2012), as well as gas exchange measurements.

Spatiotemporal variability of incident light, leaf photosynthesis, and fruit yield

The spatial variability of PPFD (up to 46%) was larger than that reported previously for a glasshouse (up to 24%, based on weekly accumulated values; Cabrera-Bosquet et al., 2016). This is mainly because the greenhouse used here had fewer opaque structural components than the glasshouse, based on comparison of the binarized images. This enables greater direct sunlight transmission into the greenhouse, leading to steeper incident light gradients. The 22% change in accumulated PPFD at distances of < 2 m, throughout the cultivation period, indicates that the spatial resolution of photography in this study (ca. 4 m2; 2 m × 2 m) was not excessive for capturing spatial variability of PPFD in the greenhouse.

The spatial variability of PPFD varied with both sun position (elevational and azimuth angle) and external light intensity. Under the same external light intensity, increasing the solar elevational angle reduces the incident angle of sunlight entering the greenhouse, reducing the shading area and consequently reducing spatial variability. In contrast, for the same solar position, strong external light intensity generates steep gradients in incident light between sunlit and shaded areas, increasing spatial variability. In reality, however, solar position and external light intensity change simultaneously on both diurnal and seasonal scales. The diurnal and seasonal complexity of the interactions between solar position and light intensity—with elevational angle increasing from morning to noon, but decreasing from winter to spring (Cabrera-Bosquet et al., 2016)—make it difficult to estimate diurnal and seasonal spatial variability of PPFD.

This study demonstrated an estimation method to account for diurnal and seasonal greenhouse microclimatic variability over an entire season, substantially longer than previous photosynthetic rate studies, which were limited to a few days. Further, this study examined spatial variability in incident light via hemispheric image analysis, which improves on previous approaches, although previous studies provided good estimates for other microclimates. For instance, Boulard et al. (2017) used a computational fluid dynamics (CFD) model integrated with a photosynthesis model to evaluate daytime microclimate in a closed greenhouse. Zhang et al. (2021) used a functional–structural plant model (FSPM) and 3D model of a greenhouse to simulate the limitations to photosynthesis on sunny and cloudy days. Kimura et al. (2020b) sampled microclimates from a moving platform, and used a photosynthesis model to evaluate the effect of greenhouse environmental controls on daytime spatiotemporal variability of A.

The findings of this study show that spatial variability of A in the greenhouse was affected by physical and physiological leaf processes, as well as differences in PPFD due to solar position and external light intensity. Saturation of A at high PPFD reduced its spatial variability. In this study, due to the early saturation of A at low PPFD in strawberry plants (ca. 600–800 μmol m−2 s−1; Hidaka et al., 2013), the spatial variability of A did not increase with external light intensity, except in winter. The reduction in A under excessive Tl caused by strong light intensity also affected its spatial variability. Here, the spatial variability of A increased in spring when Tl exceeded ca. 30 °C, which reduced photosynthetic capacity (Vcmax and Jhigh) and stomatal closure as a result of the high VPD. This highlights the importance of leaf temperature, as well as light, in predicting greenhouse leaf photosynthesis (Zhang et al., 2021).

Based on the spatiotemporal variability of A reported here, supplemental lighting would be locally effective where A is not saturated. In winter, local supplemental lighting is effective even on sunny days, because spatial variability of A (and thus shading) increases with external light intensity. However, it is not effective in spring, because A becomes locally saturated or reduced in many parts of the greenhouse, and shading is not clustered. However, responses of A to light and temperature vary even within species (Hikosaka et al., 2016). Decisions regarding supplemental lighting should thus be made based on the physical and physiological traits of the target crop.

In terms of spatial variability, for the cumulative values of the cultivation period, fruit yield was better explained by A than PPFD. This is consistent with a simple model of potential yield (Yp) expressed as a function of PPFD, the fraction of light intercepted by the crop (f), radiation-use efficiency (RUE, Monteith, 1977), and harvest index (HI):

$$ Y_{{\text{p}}} = {\text{HI}} \times \mathop \sum \limits_{i}^{n} {\text{PPFD}}_{i} \times {\text{ RUE}}_{i} \times f_{i} $$
(25)

where i is the given time within the cultivation period n. The accumulated A in this study is similar to \(\mathop \sum \limits_{i}^{n} {\text{PPFD}}_{i}\) × RUEi in Eq. (25), and the physical and physiological processes of A can be assumed to reflect variation in RUEi. Therefore, accumulated A is a more accurate indicator of yield than the accumulated PPFD alone.

Canopy photosynthesis generally explains yield better than leaf photosynthesis. In Eq. (25), fi, which is modeled on the basis of plant parameters such as leaf area index and leaf inclination angle, makes it possible to move from leaf to canopy scale. Here, the spatial variability of leaf photosynthesis alone explained up to 46% of the spatial variability of fruit yield because factors that affect plant architecture (such as planting density, transplant date, and environmental controls) were kept equal or similar across the greenhouse (Long et al., 2006). Reliable assessment of the spatial variability of yield in other species or cultivars requires fine-resolution mapping of canopy photosynthesis. To do this, the approach proposed here should be integrated with one that detects the spatial distribution of plant architectures (for instance, using an unmanned aerial vehicle; Maes & Steppe, 2019; Zhang & Kovacs, 2012).

The spatial association between A and fruit yield weakened over the cultivation period, indicating that the yield-limiting factors varied from winter to spring. In winter, the reduced incident light limits the photosynthetic rate, which consequently cannot satisfy the photoassimilate demands of the fruit organs, thereby limiting yield. In contrast, in spring, light intensity is sufficient to achieve maximum photosynthetic rate, and the ability of the fruit organ to absorb the photoassimilate becomes main yield-limiting factor (Gifford & Evans, 1981; Okello et al., 2015). This is consistent with the aforementioned suggestion that the photosynthetic rate can be effectively enhanced by supplemental lighting in winter, but not in spring. Therefore, when optimizing the spatiotemporal variability of yield, it is important to balance the photoassimilate source and sink dynamics, as well as spatiotemporal variability of photosynthesis.

Conclusions

The proposed approach for estimating the drivers of yield, using multipoint analysis of hemispheric images and a numerical photosynthesis model, revealed considerable diurnal, seasonal, and spatial variability of incident light and leaf photosynthesis. Yield was more strongly associated with spatial variability of leaf photosynthesis than in PPFD, and this effect was markedly greater in winter. This indicates that spatial variability of photosynthesis becomes the main limiting factor for yield only when there is insufficient incident light. With minor modifications to the photosynthesis model, the approach proposed here can be directly applied to other crops and types of greenhouse. This study and approach provide a basis for improving the spatiotemporal stability of greenhouse crop yields.