Spatiotemporally variable incident light, leaf photosynthesis, and yield across a greenhouse: fine-scale hemispherical photography and a photosynthesis model

Although greenhouse agriculture can generate high crop yields, they vary due to spatiotemporal differences in incident light and photosynthesis. To elucidate these dynamics, multipoint analysis of hemispheric images and a photosynthesis model were used to visualize the spatiotemporal distribution of photosynthetic photon flux density (PPFD) and leaf photosynthetic rate (A) and compared these with strawberry fruit yield in a greenhouse. This method enabled successful estimation of spatiotemporal variability in PPFD and A with relative root mean square errors of 4.4% and 11.0%, respectively. PPFD, captured at ca. 2 m resolution, varied diurnally and seasonally based on sun position and external light intensity. A showed less spatial variability, because it is reduced by physical and physiological mechanisms in the leaves at excessive leaf temperatures and becomes saturated at high PPFD. Yield spatial variability was better explained by A than by PPFD. The association between A and yield weakened over the cultivation period (R2 declined from 46% in winter to 12% in spring), thus suggesting that, over the cultivation period, factors such as photoassimilate availability replaced A as the primary limiting factor. The proposed method can be directly applied to other types of greenhouses, and the findings may facilitate spatiotemporal optimization in crop production, improving precision greenhouse agriculture.

1 3 k 25 (R n ) R N at 25 °C (μmol m −2 s −1 ) N The number of pixels per 1° in hemispheric image O O 2 concentration (mmol mol −1 ) P atm Atmospheric pressure (Pa) PPFD Photosynthetic photon flux density (μmol m −2 s −1 ) PPFD ext External PPFD (μmol m −2 s −1 ) PPFD ext,dir , PPFD ext,dif External PPFD for direct and diffuse light (μmol m −2 s −1 ) p dir , p dif Proportions of incident to external PPFD for direct and diffuse light pixel sun , pixel obs The number of pixels of the sun and obstacles overlapping the sun on hemispheric image pixel tot,α , pixel obs,α The total number of pixels and the number of obstacle pixels in each annulus generated on hemispheric image at the intervals of 1° elevation angle ϕ Initial slope of J-I curve (mol mol −1 ) R Universal gas constant (kJ mol −1 K −1 ) R d Day respiration rate (μmol m −2 s −1 ) R l,abs Absorbed longwave radiation (J m −2 s −1 ) R n Dark respiration rate (μmol m −2 s −1 ) R ni Net isothermal radiation (J m −2 s −1 )

Introduction
Using greenhouse agriculture, it is relatively easy to generate higher crop yields than outdoorfield agriculture. The enclosure makes it possible to improve the microclimate and photosynthetic rate via, for instance, CO 2 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 1 3 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 CO 2 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(Steege, , 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 C 3 photosynthesis (Farquhar et al., 1980), stomatal conductance, CO 2 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.

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:

3
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 p dir and p dif are the proportions of incident to external PPFD, for direct (PPFD ext,dir ) and diffuse (PPFD ext,dif ) light, respectively (these depend on the architecture and obstacles in and around the greenhouse). PPFD ext,dir and PPFD ext,dif were separated from PPFD ext , 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 (R s,ext ) outside the greenhouse and calculated radiation outside the atmosphere.
To evaluate p dir and p dif , 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 machinelearning-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 p dir , 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, p dir was calculated: where pixel sun is the number of pixels of the sun on the image, and pixel obs is the number of pixels of the obstacles overlapping the sun (that is, p dir 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: Fraction of sky area in total area, weighted by the elevation angle in each annulus 2.
2. Fraction of sun area without obstacles at each sun position along the path Sun path Fig. 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 where d sun 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. p dir varies spatiotemporally, because of spatial variation caused by obstacles inside and outside the greenhouse, and diurnal and seasonal variation in solar position. To calculate p dif , 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, c 1 , c 2 , and c 3 ) was calculated, and these ratios were summed to obtain p dif (Steege, 2018): where α is the elevational angle in the image (ranging from 0.5° to 89.5°, at 1° intervals along each annulus), pixel tot,α is the total number of pixels in each annulus, and pixel obs,α 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. p dif 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 C 3 photosynthesis model (Farquhar et al., 1980) integrated with leaf CO 2 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 (A c ) or RuBP regeneration (A j ). 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 CO 2 concentrations in this study. Actual A is given by where θ A is the curvature in the transition from one limitation to the other, V cmax is the maximum carboxylation rate, C c is the chloroplast CO 2 concentration, Γ * is the CO 2 compensation point in the absence of respiration, K c and K o are the Michaelis constants for carboxylation and oxygenation, respectively, O is the O 2 concentration, R d is the daytime respiration rate, J is the electron transport rate, J high 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. R d was estimated using its relationship to the dark respiration rate R n (Villar et al., 1995): C c was described by diffusion theory (Fick's law): where C a and C s are the CO 2 concentrations in the ambient air and at the leaf surface, respectively, g ah is the leaf boundary layer conductance for heat transfer, g sw is the stomatal conductance for water vapor transfer, and g m is the mesophyll conductance. g sw was estimated from the model proposed by Medlyn et al. (2011): where g 0 is the residual conductance, g 1 is the empirical constant, and VPD is the leaf-air vapor pressure deficit. g ah was estimated from wind speed u and leaf characteristic dimension d leaf using a semi-empirical relationship (Kimura et al., 2020b): where the constant 0.21 was determined using the relationship between u, measured by anemometer, and g ah , evaluated using the artificial leaf technique (Kimura et al., 2020a). The model parameters (Γ * , K c , K o , R n , V cmax , J high , and g m ) vary with leaf temperature T l , and the temperature dependence was described by the Arrhenius function or the peaked Arrhenius function: where k T l is the parameter value at a given leaf temperature, T l,K is T l in K, k 25 is the parameter value at 25 °C, R is the universal gas constant, E a is the activation energy, H d is the deactivation energy, and ∆S is the entropy factor. T l was estimated from leaf energy budget assuming isothermal net radiation (Jones, 2013), modified to molar units (Buckley et al., 2014): where T a is the ambient air temperature, γ is the psychrometric constant, R ni is the isothermal net radiation, c p is the molar heat capacity of air, D a 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, T a,K is T a in Kelvin. R ni was estimated from PPFD, and is given by: 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 R l,abs is the absorbed longwave radiation. R l,abs was described as the sum of three terms (De Boeck et al., 2012;Kimura et al., 2020b): where α l,leaf is the leaf absorptivity of longwave radiation; ε l,atm is the atmospheric emissivity of longwave radiation; T gh,K is the greenhouse temperature in Kelvin, assumed to equal T a,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): where P atm is the total atmospheric pressure, and W a is the water vapor concentration in the ambient air.
A can be obtained, together with other unknown variables (C c , C s , g sw , and T l ), 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).

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 R l,abs = l,leaf l,gh l,atm T 4 a,K + l,leaf l,gh T 4 gh,K + l,leaf l,gh l,leaf T 4 a,K 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 m 2 (1.6 m × 50.0 m) in the three vents on the western side and 49.6 m 2 (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 m 2 (1.8 m × 50.0 m) on the western side and 55.8 m 2 (1.8 m × 31.0 m) on the eastern side. These ventilation systems were operated to maintain greenhouse air temperature below 25 °C. A CO 2 generator (CG-554T2, Nepon Inc., Tokyo, Japan), supplying CO 2 at 8 kg h −1 and airflow at 33 m 3 min −1 , was situated near the greenhouse entrance at the southern side. The generator was operated in the daytime to maintain CO 2 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 CO 2 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 m 3 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. PPFD ext and R s,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. T a and relative humidity were measured by a thin-film capacitive sensor (HMP110, Vaisala, Helsinki, Finland) with a forced-ventilation system. C a 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 V cmax , J high , g m , R n , and their temperature dependencies, CO 2 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 CO 2 response curves were created under a PPFD of 1500 μmol m −2 s −1 , using CO 2 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 CO 2 response curve, the curve-fitting program in the "plantecophys" R package (Duursma, 2015) was applied to estimate V cmax and J high . Chlorophyll fluorescence was measured at the same time as the response curves were generated. To estimate g m , the variable J method (Harley et al., 1992) was applied to these data to keep the intercellular CO 2 concentration (C i ) of the curves within the range 240 to 600 μmol mol −1 (Xue et al., 2016). After each CO 2 response curve measurement, R n was measured at 0 μmol m −2 s −1 PPFD and 400 μmol mol −1 CO 2 concentration. The temperature dependences of V cmax , J high , g m , and R n was determined using a self-made fitting program in R. To determine g 0 and g 1 , diurnal changes in A, g sw , VPD, and C s 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 g 0 and g 1 .
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 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 (24) Spatial variability = Max − Min Mean × 100 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 (R 2 ) 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.

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.
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.

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.

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 T l (exceeding ca. 30 °C), which reduces photosynthetic capacity and stomatal conductance. Spatial variability in daily-accumulated PPFD increased with daily-accumulated PPFD ext , 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%.
In contrast to daily-accumulated PPFD, spatial variability of daily accumulated A increased with daily accumulated PPFD ext 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 T l , 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 PPFD ext < 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 PPFD ext (Fig. 7c). Similarly, for daily accumulated A, Moran's I showed large scatter, but it decreased suddenly at higher daily accumulated PPFD ext (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 CO 2 generator in the greenhouse. Accumulated PPFD changed by up to 22% over distances < 2 m. 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 R 2 values of the relationships between accumulated PPFD and A and accumulated daily fruit yield show three peaks (Fig. 10a-c). R 2 values were the highest in winter (the first truss) and declined markedly over time (Fig. 10a-c). Interpreting these R 2 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).

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 m 2 (2 × 2 m) was achieved via multipoint photography before transplanting and using one-point measurement of external light conditions (PPFD ext and R s,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   Moran's I 0.50 *** 0.06 0.11 *** 0.30 *** 0.78 *** 0.47 *** 0.07 0.13 *** 0.18 *** 0.77 *** 1 3 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 (V cmax and J high ) and stomatal characteristics (g 0 and g 1 ) used in the model. V cmax and J high vary seasonally due to leaf senescence and environmental acclimation (Wilson et al., 2001), as do g 0 and g 1 (Ono et al., 2013), and these parameters vary between leaves. Although V cmax and J high 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., . 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 m 2 ; 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 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 T l caused by strong light intensity also affected its spatial variability. Here, the spatial Month variability of A increased in spring when T l exceeded ca. 30 °C, which reduced photosynthetic capacity (V cmax and J high ) 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 (Y p ) 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): where i is the given time within the cultivation period n. The accumulated A in this study is similar to n ∑ i PPFD i × RUE i in Eq. (25), and the physical and physiological processes of A can be assumed to reflect variation in RUE i . 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), f i , 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 (25) Y p = HI × n ∑ i PPFD i × RUE i × f i 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.