Introduction

Hungary’s climate is mild and continental with long, warm, often dry summers and fairly cold winters, according to the Köppen-Geiger classification (Cfb) (Kottek et al. 2006). During the 1976–2005 period, based on the Representative Concentration Pathway 8.5, significantly drier (~ 20%) summers are projected to occur in Hungary by the end of the twenty-first century (Torma et al. 2020). However, being located in Central Europe where the Mediterranean climate intersects with continental European air masses, the forecasted change in the amount and distribution of precipitation (PR) for Hungary is considerably uncertain (Gaál et al. 2014; IPCC 2013). A projected decline in summer PR together with elevated air temperatures (Ta) may have a serious adverse impact on the production of local rainfed soybean (Glycine max (L.) Merr.). Akhtar et al. (2019) concluded that changing patterns of climatic perturbation, including the amount and distribution of PR, together with elevated Ta, caused variation in physiological processes and yield components of soybean. Hatfield (2014) revealed that exposure to Ta above 23 °C reduced the growth rate of soybean. Water stress lowered seed yield in soybean by reducing photosynthesis (Brevedan and Egli 2003; Karam et al. 2005).

The distribution of solar radiation inside a plant’s canopy affects physiological processes, such as transpiration, growth and development (Gao et al. 2010). Actual interception of solar radiation depends on the overall geometric shape of the plant canopy and is affected by leaf area index (LAI) and leaf angles, the geometrical distribution of incoming radiation, which depends on the sun’s angle, and on the fraction of direct and diffuse radiation (Campbell and Norman 1998). The assimilatory surface area characterizes the potential of a canopy to intercept photosynthetically active radiation (PAR), which serves as an indicator of the effectiveness of photosynthesis, plant growth and development. Kross et al. (2015) suggested that maximum soybean seed yield requires an LAI of at least 3.5–4 during reproductive stages (R2–R5).

The resilience of soybean to water stress can also be expressed through its decline in assimilatory surface (Mondani et al. 2019). The magnitude of LAI is directly proportional to the crop’s ability to absorb solar radiation (Purcell 2000), and this is related to photosynthetic potential and crop productivity. Wijewardana et al. (2018) observed accelerated leaf senescence and leaf drop of mature leaves under moisture stress in two soybean cultivars (Asgrow and Progeny) at a Mississippi State University research facility. We hypothesized that leaf senescence, which is induced by water stress, and changes in canopy structure, which is a quickly identifiable physiological trait, can modify the amount of absorbed radiation and directly impact soybean yield.

Radiation use efficiency (RUE) is the ratio between biomass and solar radiation captured in plants. The gain in biomass by crops is strongly related to the amount of intercepted radiation. RUE in soybean depends on meteorological elements, such as vapour pressure-deficit (VPD), Ta and water supply (Muchow et al. 1993). Even though Muchow et al. concluded that RUE under well-watered conditions remained constant throughout most of the soybean growth cycle, it was unaffected by the ambient environment (incident and reflected PAR, CO2 and vapor fluxes). Nguy-Robertson et al. (2015) reported a positive linear correlation between biomass production and intercepted PAR just before physiological maturity (R7 stage), when translocation of photosynthates replace PAR interception, which is the main factor determining biomass accumulation (Black and Ong 2000).

Evapotranspiration rate (ET) in soybean can be estimated by direct measurements, such as the application of a lysimeter or evapotranspirometer (Karam et al. 2005), a flux tower (Lathuillière et al. 2018), empirical or statistical models (Lobell and Asseng 2017; Paredes et al. 2015), and crop coefficient (Allen et al. 1998). Although these techniques have both advantages and shortcomings, an evapotranspirometer is convenient because it can measure ET under field conditions. Soybean ET might be impacted by several factors, such as crop variety, fertility, water regime, climate, and plant protection (weed management), all of which directly affect seed yield (Mbangiwa et al. 2019). In this study, crop transpiration and soil evaporation, which contribute to ET, were investigated in soybean at a field scale.

The water stress tolerance of a crop depends not only on the physiological properties of a particular variety, but also on the duration of water stress and the stage at which soybean is exposed to water stress. Kumagai et al. (2020) reported that drought during flowering and at the pod-enlargement stage significantly reduced soybean yield by 29%. Payero and Irmak (2013) found that the best strategy was to start irrigation at the flowering stage if there was a full soil–water profile at planting. In their study, soybean seed yield was positively correlated with ET by a quadratic function (R2 = 0.75), independent of the timing of stress application. The response of soybean to water stress can also determine variety selection for cultivation, even under variable climate or geographic conditions. An increasing need for vegetable protein (oil) for humans and livestock, the use of new varieties, and the local impact of global warming are the main factors driving farmers’ decisions about whether to invest in irrigation, or not, under the variable continental climate of the Carpathian basin. In Hungary, soybean culture is currently expanding. Although many soybean studies involving the yield–water relationship and/or the link between yield and radiation have been conducted, analyses of soybean yield attributes with different watering and radiation properties under a continental climate are limited. Due to the expansion of negative impacts of global warming in this area (Kocsis et al. 2019), it is essential to identify a proper soybean variety that can be grown on-site and that is tolerant to water stress for higher crop production and to avoid or prevent reductions in yield.

The aims of this research were to evaluate three water supply levels in two soybean varieties for select crop characteristics (i.e., total aboveground dry matter (TDM), seed yield (and its attributes), PAR transmission, and RUE) and to test whether the differences in soybean yield traits were associated with water stress or a change in the attenuation of crops by PAR. Thus, this study aimed to provide new information to farmers and decision-makers at regional and local levels by (i) characterizing RUE, extinction coefficient (k) and PAR penetration into the soybean canopy grown under different water supplies using Thornthwaite–Matter type compensation evapotranspirometers and rain-fed conditions, and by (ii) quantifying variability in soybean ET and yield attributes under limited and non-limited water supplies, contributing to reliable field-scale estimates of yield. A specific hypothesis of this study was to quantify the different impacts of water stress during flowering on two soybean varieties with differing water demands. The novelty of this approach is that yield and ET data were not analysed separately for each year or variety. Rather, a uniform model was fitted for the entire period and the effect of year and variety were treated as random effects.

Materials and methods

Site description and agronomic practices

Field experiments were conducted at the Keszthely Agrometeorological Research Station (ARS) (latitude: 46°44ʹ N, longitude: 17° 14ʹ E, elevation: 124 m above sea level) in 2017–2019. The experimental area is characterised by a continental climate with hot and mostly dry summers. The long-term (1971–2000) average Ta was 16.9 °C during the growing seasons. In the same period, the long-term seasonal PR sum was 384.3 mm with a significant inter- and intra-seasonal variabilities in the amount and distribution patterns. The soil in the experimental field is a clay loam, Haplic cambisol (FAO 2006) with a mean bulk density of 1150 kg m−3 in the top 1 m of the profile. Plant available water holding capacity to 1 m soil depth was 273 mm m−1. In the rainfed treatment, the water content profiles were measured gravimetrically in the top 1 m at 0.10-m intervals every 10 days. Soil samples were weighed immediately after sampling, and oven-dried at 105 °C for 72 h until constant weight. Available soil water (SM) was used, and water depth was measured in mm (Füleky 2004). The small size of the evapotranspirometer’s chambers did not allow the use of a gravimetric method in the controlled water treatments (WW and RO).

Two soybean varieties (Sinara, Sin: tolerant to water stress; Sigalia, Sig: variety bred for average weather conditions) with different water demands were grown under three conditions: unlimited water (WW), water stress applied only at the flowering stage (RO), and rainfed (P). Due to the fixed nature of the evapotranspirometer, the experimental design was a complete block with three replications. The rainfed experimental field, which was about 0.5 ha in size and lay adjacent to the growing pots of the evapotranspirometers, was divided into two sections with the two soybean varieties sown separately, in an area of 50 m (width) × 60 m (length) for each variety (Fig. 1). ET was measured using a Thornthwaite–Matter type compensation evapotranspirometer in which the surface area of pots was 4 m2 and the depth of pots was 1 m. Recording the changes in soil–water volumetric content in terms of a soil water budget (PR, irrigation and drainage), evapotranspiration rate can be expressed as a residual term on a daily basis (Anda 1986). In the RO treatment, water stress was simulated by closing the water supplier of the evapotranspirometer every second day during the flowering stage (R1) until the grain-filling stage (R4–R5). In RO, rain shelters were used to avoid rainfall (Anda et al. 2020) from the flowering stage (R1).

Fig. 1
figure 1

Layout of the field experiment at the Keszthely Agrometeorological Research Station (ARS). Legends in the field pots of the evapotranspirometer (left hand side) were; 1: Sinara with unlimited watering, 2: Sinara with water stress; 3: Sigalia with unlimited watering; 4: Sigalia with water stress. Rainfed plots of Sinara and Sigalia were placed in the field immediately adjacent to the evapotranspirometer pots. Evapotranspiration pots on the right-hand side with open quadrants were used only for destructive LAI measurements during flowering. The lower part of the sketch separated by a dashed line contains the meteorological garden with the measuring cellar of the evapotranspirometer, the wind mast, the automatic climate station of VAISALA (QLC-50 type; VAISALA, Helsinki, Finland), and a tipping bucket rain gauge of the QLC-50. Evaporation pans (standard A pan, GGI, WMO and U pans) were also included in the ARS

Full size image

Before planting, the experimental field and pots were tilled and fertilized with 300 kg of N, P, and K ha−1 (1:1:1, w/w) using commercial granular fertilizer (Kite Zrt, Nádudvar, Hungary). The planting dates were 10 May, 24 April and 24 May (cool spring) in 2017, 2018 and 2019, respectively. Plant density was 40 plant m−2 and inter-row spacing was 0.24 m. Cropping practices were those recommended locally by specialists of the University of Pannonia.

In accordance with the size of evapotranspirometer’s pots, total aboveground biomass and TDM samples of rainfed plots were randomly collected from five 4 m2 areas per variety at the end of the vegetation cycle from 18 to 30 September over the three-season study. Crop samples were oven-dried at 65 °C until constant weight, and yield components (seeds, pods) were separated by hand. Seed yield (g m−2) was normalized to 13% water content in seed. 1000-grain weight was also determined. An Inframatic 9200 NIR Grain Analyzer (PerkinElmer, Waltham, MA, USA) was used to measure the protein and oil content of seeds on a dry-matter basis.

Meteorological measurements, including PAR and LAI

Incoming solar radiation (Rs) was measured with a pyranometer (model CM-3; Kipp & Zonen Corp., Delft, the Netherlands) fitted to the climate station (model QLC-50, Vaisala, Helsinki, Finland) at the ARS, located at the border of Keszthely. The meteorological sensors, except for the anemometer, were placed at a standard height 2 m above the surface of the ground in the meteorological garden of the ARS. The height of the anemometer was 10.5 m (Anda et al. 2015). The standard meteorological station belongs to the national meteorological observation network operated and supervised by the Hungarian Meteorological Service under the provision of the World Meteorological Organization (WMO) (WMO 2008). The basic meteorological data of Ta, Tmax, PR, relative humidity (RH), VPD and Rs were used in the study.

Detailed estimates of PAR penetration were obtained by stratifying crop height into several (0.2 m thick) layers. PAR penetration into the canopy was measured by three sensors (model LI-190R; LI-COR Inc., Lincoln, NE, USA) positioned at three measuring points on sunny and shaded patches mounted on a mobile scaffold stubbed in the ground (to the north opposite an incoming light beam) in parallel to row directions (Fig. 2), during July (R5 stage). In July, PAR penetration samples were taken on two clear-sky and calm (no wind) days. For ease of presentation, only data of one selected day per season is presented (2017: 15th July; 2018: 20th July and 2019: 25th July). The vertical post of the scaffold was 150 cm high, and the movable horizontal arm was 50 cm long with a fixed levelled tray. All components of the scaffold were painted black to minimize multiple scattering of radiation that may have biased the measurements. The horizontal arm could be fixed at any height of the post, and calibrations in 20-cm intervals were marked on the vertical post. The lowest position was beneath the canopy on the surface of the ground. The horizontal arm maintained the tray with a red cardboard (ISO-A standard A4: 210 mm × 297 mm) on it, perpendicular to the plant rows, in an east–west orientation. PAR penetration measurements were made within 2 h of solar noon, at an incident angle of 55°–60°, over the entire period. This study followed the instructions of Purcell (2000), reporting only 2 h per day around the time of solar noon, under unobstructed sunlight for radiation interception measurements. Following the conclusion of Flénet et al. (1996), samples were taken at the R5 stage because between the R3 and R5 stages, soybean development has no effect on k.

Fig. 2
figure 2

Sketch (a: plan view; b: side view) of the construction used to measure penetration of photosynthetically active radiation (PAR, W m−2) into the soybean canopy at 20-cm intervals until the top plant height (below 100 cm). Three PAR sensors (1) were included in the study: two of them in the shade (grey), and one in a sunny spot of a mobile wooden tray of 20 × 35 cm (2) moved layer by layer. Each mounted flange of PAR sensors incorporated a bubble/spirit level and three adjustment screws for keeping the sensor horizontal. Two metal support rods (3) with mounted consoles (5), stubbed in the ground stabilised by baskets (4), were applied to hold the tray. 6 denotes the stems (plan view) or crops (side view). All parts of the tool were painted black. Leaves below the tray were cut to determine the leaf area every 20 cm from the bottom to the top of the canopy

Full size image

When sampling the height of plants, PAR values were recorded over a 2-min period for each arm-position by a light sensor logger (LI-1500, LI-COR Inc.) at 1-s intervals. Samples were collected, and photos were taken, within 1.5 h.

Three sensors were positioned at one sunny and two shaded patches on the red cardboard. The weighted average PAR (In) [W m−2] was calculated with measured sunny spot’s PAR (Is) [W m−2], sunny area (f), and shaded spot’s PAR (Ih) [W m−2]. For each layer, 2-min integrated radiation values were used:

$${I}_{n}=f {I}_{s}+(1-f){I}_{h}$$
(1)

The sizes of sunny and shaded areas were computed on the basis of digital photographs taken of various plant heights from the bottom to the top of each plant in 20-cm intervals using flat pieces of evenly coloured red cardboard and a 17.9-megapixel digital camera (Canon EOS 7D, Canon Inc., Tokyo, Japan). Images were stored in a raw format then transferred to a PC for analysis. Due to the irregular location of the shaded and sunny patches, segmentation was analysed using the region-growing method applied in an image processing application SGMW 0.1. This program was developed in-house by our team (Anda and Lőke 2005) as follows.

A new raster point of the growing region of interest was defined as:

$$g\left(x,y\right)=\left\{\begin{array}{c}yes\, if \,d(x,y)\le 50\\ no\, if\, d\left(x,y\right)>50\end{array}\right.$$
(2)
$$d\left(x,y\right)=\left(R\left(x,y\right)-{R}_{0}\right) \quad \mathrm{and}\quad \left(G\left(x,y\right)-{G}_{0}\right) \quad \mathrm{and}\quad \left(B\left(x,y\right)-{B}_{0}\right)$$
(3)

where g(x, y) is the region map after thresholding and d(x, y) is the thresholding distance between the initial pixel (R0, G0, B0) and the next pixel (R(x, y), G(x, y), B(x, y)). Different types of pixels were distinguished, but these were masked using the Magic Wand technique, and then counted (Hafsi et al. 2000).

In addition to weekly non-destructive leaf area determination of five selected plants per season (Anda et al. 2019), destructive leaf area measurements were also carried out on the same section of the row in which PAR penetration was detected. These individual evapotranspirometer pots were excluded from the water supply observation. To measure radiation penetration into the canopy, trifoliates were cut layer by layer from the bottom to the top height (every 20 cm), in accordance with levels of PAR measurements. Each trifoliate was placed on a flat surface of an evenly coloured red cardboard with a scale on the side, directly under a vertically mounted Canon EOS 7D digital camera. A histogram-based threshold segmentation was applied in the SGDIP 0.1 image processing program to count the area of the region of interest (Anda et al. 2019). Five fully developed plants from each treatment (WW, RO, and P) were included.

The fraction of incident beam radiation that penetrates the canopy was described by Beer’s law (Monsi and Saeki 1953):

$${I}_{n}= {I}_{0 }\mathrm{exp}(-kLAI)$$
(4)

where I0 and In are the flux densities of incoming solar radiation on a horizontal surface above and below the canopy, at a specific height. k is the extinction coefficient for the plant stand. The on-site measured k was calculated using Eq. 4.

To compute the theoretical extinction coefficient (ksim), the ellipsoidal leaf angle distribution function of Campbell and Norman (1998) was applied:

$${k}_{sim}=\frac{\sqrt{{x}^{2}+{tan}^{2}\varphi }}{{x}+1.774({x+1.182)}^{-0.733}}$$
(5)

The parameter, x for ellipsoidal, the ratio of vertical to horizontal projections of canopy elements (leaf angle distribution) equalled 0.81 for soybean (Campbell and van Evert 1994). Solar zenith angle, which was ψ, calculated from latitude, longitude and time of estimation.

RUE was calculated as follows:

$$\mathrm{RUE}=\frac{TDM}{{I}_{0} \times F}$$
(6)

where TDM was in g m−2, and F was the fraction of PAR intercepted by the canopy:

$$F=1-{\mathrm{exp}}^{-k LAI}.$$
(7)

Statistics

The effects of LAI, year, variety and meteorological variables on ET were analysed by hierarchical regression. It was assumed that ET is mainly influenced by LAI and meteorological variables, but that this relationship may vary slightly from year to year and from variety to variety. Therefore, mixed-effect regression models (Pinheiro et al. 2019) were fit, where LAI and the meteorological variables were fixed effects while the year and variety were random effects. In the first model, LAI was the only independent variable. Year entered the second model as a random factor, variety entered the third model as a random factor, and finally meteorological variables (Ta, RH, Rs) entered the model as fixed factors. Since LAI was measured on a weekly basis, LAI of intermediate days was interpolated by cubic spline interpolation. The ET measured on consecutive days is not independent, therefore, the residuals were modelled as a first-order autoregressive (AR(1)) process in each model. To compare the models, they were fitted with the maximum-likelihood (ML) method to be comparable, but the final model was refitted with the restricted maximum-likelihood (REML) method to obtain unbiased estimations of variance and covariance parameters. In each step, the model with lowest Akaike Information Criterion (AIC) was selected. If the newly entered variable increased AIC, it was excluded from further model-building. All calculations were performed using the nonlinear mixed effects (nlme) package in R (Pinheiro et al. 2019). The figures created by the functions of the nlme package were only used to explore data analysis and are not included in the paper. The explanatory power of the models was assessed using mean absolute error (MAE) and mean absolute percentage error (MAPE). The explanatory power of the variables was quantified by standardized beta, partial and semi-partial correlations.

The effects of water treatment and variety on yield traits were analysed with the linear mixed effects (LME) model. Water treatment, variety and their interaction were the fixed effects in the models. The random effects were nested four-fold, as year, treatment within year, variety within treatment, and within-group error:

$${y}_{ijkl}= {\beta }_{0}+ {\beta }_{i}^{tr}+ {\beta }_{j}^{var}+ {\beta }_{ij}^{trxvar}+{b}_{k}^{y}+{b}_{ki}^{y/tr}+{b}_{kij}^{y/tr/var}+{\epsilon }_{ijkl}$$
(8)

i = 1, 2, 3; j = 1, 2; k = 1, 2, 3; l = 1, …, 5;

$${b}_{k}^{y}\in N\left(0,{\sigma }_{1}^{2}\right); {b}_{ki}^{y/tr}\in N\left(0,{\sigma }_{2}^{2}\right); {b}_{kij}^{y/tr/var}\in N\left(0,{\sigma }_{3}^{2}\right);{\epsilon }_{ijkl}\in N\left(0,{\sigma }^{2}\right);$$
(9)

where y is a yield trait, \({\beta }_{0}\) is the overall mean, \({\beta }_{i}^{tr}\) is the fixed effect of treatment i, \({\beta }_{j}^{var}\) is the fixed effect of variety j, \({\beta }_{ij}^{trxvar}\) is the interaction term of treatment i and variety j, \({b}_{k}^{y}\) is the random effect of year k, \({b}_{ki}^{y/tr}\) is the random effect of treatment i within year k, \({b}_{kij}^{y/tr/var}\) is the random effect of variety j within treatment i within year k and \({\epsilon }_{ijkl}\) is the within-group error. The random effects were assumed to be independent, and normally distributed variables. Models were fitted by REML.

Limited water supply is expected to decrease both oil and protein content but not to the same extent. To explore the relation between the quantity and composition of yield, two canonical variables were generated, as follows. Variables describing yield were divided into two groups. One group consisted of variables associated with the quantity of yield (seed yield, 1000-seed weight, TDM, oil mass and protein mass). The other group consisted of variables associated with the composition of seeds (percentage oil content and percentage protein content). Principal component analysis (PCA) was performed separately for both groups. The first principal component (PC) was retained from both groups and these canonical variables were used to analyze the treatment-variety-year combinations. Hereinafter, these derived variables are referred to as QuantityPC and CompositionPC.

The effects of the treatment, variety and season of seasonal mean LAI were analysed by three-way analysis of variance (ANOVA). A two-tailed one-sample t test was applied to compare the studied crop variables (LAImax, RUE, seed yield and its components). To facilitate the presentation of results, a 95% confidence interval was calculated.

Statistical computations were implemented using SPSS software version 17.0.

Results

Weather, SM, ET and its attributes

High seasonal mean Ta values were recorded by an automatic weather station installed at the research area of the ARS during the three-season study. A 1.3, 2.4 and 2.6 °C higher seasonal average Ta were measured in 2017, 2018 and 2019, respectively, compared to the long-term mean of 16.9 °C (1971–2000). Cooler springs in 2017 and 2019 forced soybean seeding to be postponed until May. At the study site, the growing season’s PR was 6.4% higher than the climatic norm during 2018 while the other two observation periods, 2017 and 2019, received 44.9% and 49.6% less PR than the long-term average (384.4 mm). Septembers were disregarded from the analysis because water is not utilized by soybean at the end of the growing season. It is important to note that there were acute differences in VPD between the two warm seasons (2017 and 2019). Extremely high VPD values were observed in the summer of 2017 (July: 0.86 kPa; August: 0.92 kPa) and in August of 2019 (0.85 kPa). At the same time, the long-term monthly climate norm was about 0.6–0.7 kPa. The lowest monthly average VPD was in June of 2019 (0.35 kPa).

Due to insignificant differences between the studied varieties, their SM was pooled in treatment P, irrespective of the season. The initial SM during sowing was 126.3 mm in 2017 (5 May), 225.5 mm in 2018 (24 April) and 290.3 mm in 2019 (24 May) (Fig. 3). At the end of the seasons (end of August), 114.0, 152.1 and 137.1 mm of SM indicated the end of the vegetative period in 2017, 2018 and 2019, respectively. At about half the PR in 2017 and 2019, seasonal average SM values were reduced by 45.9% and 17.9% in the drier seasons of 2017 and 2019, respectively, relative to wet 2018. The smaller reduction of SM in 2019 was probably due to extremely high soil moisture in spring. It is important to note that the large decline in SM during the reproductive stage in 2019 was more intense than in 2017.

Fig. 3
figure 3

Ten-day variation in the fraction of available soil moisture, SM (mm) and precipitation sums, PR (mm) from April to the end of August on rainfed plots during the soybean growing seasons in a 2017, b 2018 and c 2019. DOY: days of the year. Arrows indicate sowing date and the beginning of the reproductive stage (R1)

Full size image

In 2017, a cumulative ET of 694.9 mm and 657.4 mm was measured in Sin WW and Sig WW, respectively. In the 2018 growing season, total ET of 615.6 mm in Sin WW and 576.1 mm in Sig WW were observed. The seasonal ET sums in 2019 were 634.7 mm and 646.1 mm in Sin WW and Sig WW, respectively. Irrespective of variety, water stress reduced total ET by about half, ranging from 291.9 mm (Sig RO, 2017) to 332.1 mm (Sig RO, 2019).

In the WW treatment, seasonal mean LAI varied little [from 4.94 (Sig, 2018) to 5.27 (Sin, 2017)] over the entire period. Seasonal average LAI in RO was also similar in both varieties, ranging from 2.94 (Sig, 2019) to 3.99 (Sin, 2018) in the three-season study.

Regression statistics demonstrated that LAI by itself estimated ET with a MAPE of 31% (Table 1). When year was included, MAPE decreased by 2.5%. Variety did not change the accuracy of the model, so it was omitted from further models. Inclusion of the meteorological variables strongly reduced AIC and improved the accuracy of the model. The effect of daily mean Ta, RH and Rs sum was also significant. The mean absolute percentage error of the model, including LAI, year and meteorological variables was 17.2%. The model estimated ET with a MAE of 0.85 mm day−1 (Table 1).

Table 1 Goodness of fit statistics for evapotranspiration models. Abbreviations: leaf area index (LAI), Akaike Information Criterion (AIC), mean absolute error (MAE), mean absolute percentage error (MAPE)
Full size table

The relative explanatory power of the variables is presented in Table 2. LAI and RH had the strongest explanatory power with a standardized beta of 0.34 and − 0.41, respectively.

Table 2 Explanatory power of the fixed-effect variables in the model for evapotranspiration
Full size table

The unstandardized coefficients of the final model that were refitted with the REML method are shown in Table 3.

Table 3 Unstandardized coefficients of the fixed-effect variables in the model for evapotranspiration
Full size table

Daily mean Ta showed the highest year-to-year variability (Table 3). Its estimated value for all years was 0.16, meaning that a 1 °C increase in daily mean Ta increased ET by on average 0.16 mm day−1. This coefficient varied with a standard deviation of 0.07 among the three years.

The coefficient of the AR(1) process describing the residuals was 0.55, which indicates a fairly strong autocorrelation.

Some factors of canopy architecture influencing PAR penetration

Weekly LAI trends were similar to each treatment over the entire study period (Fig. 4): low LAI values from early June at about 0.1–0.6 peaked in the R5 stage (see also LAImax below) and then declined gradually until 1.0. However, the size of the assimilatory surface varied among different treatments. Across the three seasons, seasonal mean LAI varied between 4.94 (Sig, 2018) and 5.27 (Sin, 2017) in WW, between 2.94 (Sig, 2019) and 3.99 in RO, and between 3.67 (Sig, 2017) and 5.55 (Sin, 2018) in P.

Fig. 4
figure 4

a–c Weekly variation in leaf-area index (LAI) of two soybean varieties (Sinara: Sin; Sigalia: Sig) at three watering levels (unlimited: WW; water-stressed: RO; and rainfed: P) during the 2017–2019 observation period (a: 2017, b: 2018, c: 2019). DAS: days after sowing

Full size image

ANOVA showed a significant effect of season (S; p < 0.001), water (W; p < 0.001) and their interaction (S × W; p < 0.001) on seasonal mean LAI. The impact of variety on LAI and all the other interactions were insignificant. LAI and water supply were inversely related in each growing season. The highest seasonal mean LAI values were obtained in WW. A 29.77% (p < 0.001) decrease in Sin RO and a 33.8% (p < 0.001) decrease in Sig RO were observed in relation to the LAI of WW over the entire observation period. Over the same period, only a 10.63% (Sin; p < 0.001) and 15.34% (Sig; p < 0.001) decline in seasonal mean LAI of P compared to WW was observed. Except for WW, the seasonal mean LAIs increased significantly in 2018 compared to 2017 and 2019. Pronounced declines in seasonal mean LAI of 14.1% (Sig; p < 0.001) and 15.9% (Sig; p < 0.001) in Sin P, or 25.9% (Sig; p < 0.001) and 8.6% (Sig; p < 0.007) in Sig P, were observed in 2017 and 2019, respectively, compared to their 2018 counterparts. The regulation of water supply, in both WW and RO, likely smoothed the seasonal differences in the size of the assimilatory surface.

Differences in LAImax between Sin and Sig were not significant, thus results were pooled across both varieties. Irrespective of the season and treatment, soybean LAI peaked in the R5 stage. A second-order polynomial function was fitted to obtain a vertical distribution of the trifoliate areas (cm2) in different leaf layers for a fully developed canopy (Table 4). The maximum trifoliate area was between the 7th (WW, 2018) and 10th (P, 2018) leaf layers. This suggests that trifoliates above the 10th and below the 7th leaf layer utilized radiation ineffectively. There was no difference between WW and P for most architectural attributes (Table 4).

Table 4 Canopy attributes (LAImax, crop height (cm), number, and distribution of trifoliate areas) of soybean in response to differing water supply (WW, RO and P  = unlimited, water-stressed and rainfed crops, respectively). x =  the number of leaf layers; y = vertical trifoliate area distribution
Full size table

Considering the results of ANOVA when assessing LAImax, only crop evapotranspiration (p < 0.001) was significant. When averaged across the three seasons, 34.1% lower LAImax occurred in RO compared with WW. All the other responses and interactions related to LAImax were insignificant. Plant height also tended to be 17% less in RO than in WW.

k for fully developed soybean at about solar noon

The calculation performed by Eq. 5 yielded 0.48 for ksim between July 10 and 20, at Keszthely. At the same time, measured k values under a clear-sky at high solar radiation (solar elevation > 55°) around 13:00 h (LMT) ranged from 0.42 ± 0.03 (WW, 2017) to 0.47 ± 0.05 (RO, 2018) (Fig. 5). Since there was no difference in radiation properties between the two studied varieties, only the data of Sig was presented.

Fig. 5
figure 5

Change in the simulated extinction coefficient (ksim) and measured extinction coefficients of variety Sigalia under high solar angles in July over the entire study (2017–2019). Measured WW (unlimited), RO (water-stressed) and P (rainfed) values were obtained under clear-sky days over July between 2017 and 2019. Dotted line represents Ksim for 13:00 LMT from the middle of June to the end of July, assessed according to Campbell and Norman (1989, 1998). DOY: days of the year

Full size image

There were hardly any differences between measured and theoretically obtained k values (Fig. 5). Irrespective of the treatment, a moderate overestimation (4.4–13.3%) in estimated ksim, relative to measured k, was observed.

PAR transmission

Daily incoming solar radiation for each sample day during July was 25.3 MJ m−2 in 2017, 23.1 MJ m−2 in 2018 and 22.7 MJ m−2 in 2019. PAR penetration at different plant heights was modelled by an exponential function. The model was fit to data measured above the ground surface from 0 to 100 cm at every 20 cm (Fig. 6).

Fig. 6
figure 6

Vertical profile of photosynthetically active radiation (PAR, W m−2) flux density transmitted by soybean in 20-cm intervals from the bottom (0 cm) to the top of the canopy. Samples were taken under full-sun conditions around 13:00 h (LMT), on 15th July 2017, 20th July 2018 and 10th July 2019. Different water supplies were controlled (unlimited WW and water-stressed RO) and rainfed (P). RMSE: root mean square error, W m−2

Full size image

The shapes of radiation penetration patterns during the three seasons were essentially identical for every treatment (Fig. 6). Values of flux densities for the upper third of plant height were strongly scattered. PAR flux density at the bottom of the canopy ranged from 2.3% (WW) to 4.9% (RO) compared to above-canopy measurement over the three seasons. More radiation penetrated the two higher layers than the lower ones.

Radiation use efficiency

RUE showed slight variation across the three seasons (Table 5) ranging from 0.81 (Sig RO, 2019) to 1.48 (Sin ET, 2017).

Table 5 Radiation use efficiency (RUE; g MJ m−2) for soybean in three growing seasons between 2017 and 2019
Full size table

The combined analysis demonstrated that there were significant effects of treatment (p < 0.001) and season (p < 0.001) on RUE. All the other interactions were insignificant, indicating that the effect of water on RUE was similar in both varieties during the three seasons.

Quantitative components of soybean production (TDM, seed yield and 1000-seed weight)

Three quantitative traits of soybean production, namely TDM, seed yield and 1000-seed weight, all decreased with water stress, irrespective of the variety or season. TDM values ranged from 1.07 kg m−2 (Sin RO, 2017) to 1.95 kg m−2 (Sin ET, 2018) for all three growing seasons (Fig. 7a). During 2017–2019, additional water increased TDM by 37.3% in Sin WW (p<0.001), and by 25.4% in Sig WW (p < 0.001) compared to RO. TDM of both varieties was not consistently different between WW and P. The lowest TDM values were observed in 2019.

Fig. 7
figure 7

Boxplots of yield-related variables by year, variety (Sinara, Sin; Sigalia, Sig) and water treatment: a total dry matter (g m−2), b seed yield (g m−2), c 1000-seed weight (g), d oil mass (g m−2), e protein mass (g m−2). WW, RO and P denote unlimited, water stress and rainfed conditions, respectively

Full size image

During 2017–2019, mean grain yield of Sin was 536.03 ± 30.5 g m−2 in WW, 277.48 ± 13.9 g m−2 in RO, and 424.31 ± 17.5 g m−2 in P (Fig. 7b). During the same period, seed yield for Sig was 472.31 ± 8.3, 274.69 ± 21.5, and 334.49 ± 15.1 g m−2 in WW, RO and P, respectively. Variation in seed yield was strongly associated with changes in watering (p  < 0.001), but less with varieties. In WW, varieties had 48.2% (Sin, p < 0.001) and 41.8% (Sig, p < 0.001) higher seed yield than in RO, but in P, there were no differences between P and RO for Sig over the entire period. Grain yield was 11.9% (p < 0.001) and 26.9% (p < 0.001) higher in WW and P, respectively, of Sin relative to Sig values. There were no differences between varieties when RO treatments were compared.

Compared to WW, declines in 1000-seed weight of RO ranged from 17.0% (Sig, 2017, p < 0.001) to 26.8% (Sin, 2017, p < 0.001) (Fig. 7c). Somewhat lower decreases of 19.2% (Sin, 2017–2019, p < 0.001) and 8.7% (Sig, 2017–2019, p < 0.001) were observed in 1000-seed weight of Sin P compared to Sin WW.

Compositional traits of soybean production (oil and protein percentage)

Due to wet weather in 2018, the decline in seed oil mass ranged from 10.8% (Sin WW, p < 0.001) to 31.8% (Sig WW, p < 0.001) relative to the values for the dry 2017 (Fig. 7d). Seed oil mass of rainfed P was also sensitive to actual weather in 2018; 24.3% and 31.1% decreases were observed in Sin P (p < 0.001) and Sig P (p < 0.001), respectively, compared to 2017 values. With the exception of Sig RO in 2018, water stress significantly increased the oil mass percentage of seeds, ranging from 5.3% (Sig, 2018, p < 0.001) to 19.5% (Sin, 2017, p <  0.001) relative to corresponding WW values over the entire observation period. There were hardly any differences in oil mass between soybean varieties during the three-season study.

The protein mass in WW was 8.8% (p < 0.001) and 6.2% (p < 0.001) higher in Sin than in Sig, during 2017 and 2018, respectively. Except for 2017, Sig produced slightly higher protein mass in RO (2018: 9%, p < 0.001; 2019: 1.2%, p <  0.002) than Sin. Water stress during flowering decreased protein mass by 15.9% (p <  0.001) in Sin and by 10% (p<0.001) in Sig over the entire period (Fig. 7e).

Linear mixed effects model for yield traits

The LME models fitted for quantitative traits (Table 6) also revealed that the effect of water supply was significant for all traits. WW and RO produced the highest and lowest yield, respectively. However, the effect of variety and their interaction were significant only for seed yield. Sin produced higher seed yield in both P and WW, but the difference was less for RO.

Table 6 Coefficients of fixed effects and standard deviations of random effects from LME models fitted to yield traits. RO and P denote unlimited and water stress conditions, respectively. Sin = variety Sinara
Full size table

In the case of seed yield, two random effects were negligible, but the random effect of treatment within each year was about half of that of the within-group error (Table 6). The reason is that in 2017, WW, and in 2019, RO, resulted in lower seed yield compared to the other two seasons. In the case of 1000-seed weight, the random effect of year was about 60% of the within-group error, while that of the treatment within each year was negligible (Table 6). This is due to the fact that in 2019, 1000-seed weight was lower for each treatment. In the case of TDM, the random effect of year was only somewhat lower than the within-group error. Other random effects were negligible due to the lower TDM in 2019 for each treatment.

Regarding the qualitative traits, the effect of water treatment was significant on both oil mass and protein mass (Table 6). WW and RO resulted in 30 g m−2 higher and 40 g m−2 lower oil mass, respectively, than P. Similarly, protein mass was 45 g m−2 higher in WW and 81 g m−2 lower in RO. Even though the effect of variety and the treatment × variety interaction was significant only for protein mass, similar effects were observed for both oil and protein mass. Sin produced higher oil and protein mass, but the difference was lower for RO.

A small random effect of year was observed for both oil and protein mass (Table 6), where 2017 was different from the other two years. The random effect of year/treatment was considerable on protein mass. This was due to the relatively low value measured at RO in 2017.

Overview analysis with principal components, PCA

PCA was performed separately for both the quantity- and composition-related variables. The first PCs explained 73.4% and 89.7% of total variance of the quantity- and composition-related variables, respectively. The association between original variables and the first PCs, referred to as QuantityPC and CompositionPC, can be seen in Fig. 8. QuantityPC was positively correlated with all five original variables (Fig. 8a), while CompositionPC was positively correlated with percentage protein content and negatively correlated with percentage oil content (Fig. 8b). Observations were projected onto the plane spanned by QuantityPC and CompositionPC (Fig. 9).

Fig. 8
figure 8

Quantity-related (a) and composition-related (b) variables in the plane of principal components (PCs). Coordinates of variables are the correlations between PCs and the variables (total dry matter: TDM (g m−2), seed yield: seed (g m−2), 1000-seed weight (g), oil mass (g m−2), protein mass (g m−2))

Full size image
Fig. 9
figure 9

Observations described by the quantity-related and composition-related first principal components. Color of symbol refers to treatment, Sinara (Sin) is depicted by empty symbols, Sigalia (Sig) is represented by filled symbols. Higher values of CompositionPC refer to higher protein and lower oil content. Higher values of QuantityPC refer to higher yield. WW, RO and P denote unlimited, water stress and rainfed conditions, respectively

Full size image

Figure 9 reveals that RO remarkably decreased the quantity of yield and resulted in greater oil content. WW increased the percentage protein content. In 2018 and 2019, a similar pattern was observed, but quite different from the 2017 pattern. The percentage oil content in 2017 was relatively higher than in 2018 and 2019, but the ranking of RO < P < WW in terms of quantity was the same as in 2017. In P, Sin produced similar percentage oil content but higher yield than Sig. The difference was similar in all years. Comparing the three years’ yield component data of P, 2018 and 2019 were similar in magnitude to WW while 2017 was similar in magnitude to RO. In RO, the difference in yield components between varieties were small in 2017 and 2019. In 2018, the quantity was similar, but Sin produced a lower percentage protein content than Sig. In WW, the only difference between varieties was in 2017: Sig produced somewhat lower yield but a remarkably higher percentage oil content.

Discussion

One of the critical points of soybean cultivation is the timing of sowing and emergence based on actual SM (Egly and Cornelius 2009). The currently recommended sowing period for soybean in Hungary is late April to early May. However, variation in sowing times related to SM during the springs of the study differed (see also Fig. 3). Low SM in 2017 delayed sowing until 9th May. In 2019, due to extremely wet soil and about 2.0 °C cooler air temperatures relative to the climate norm (from the middle of April to the second decade of May) planting was postponed until 24th May. High and relatively stable SM contents were measured throughout 2018 due to high and evenly distributed PR. In 2017 and 2019, the SM curves showed slightly declining trends from the end of June, corresponding to the early reproductive stage of soybean (R1). In the same time period, a moderate decrease in SM was also observed throughout the wet 2018.

Leaf area, photosynthesis and biomass, the main crop features affected by water stress, are indicative of final soybean yield (Kross et al. 2015). Among varieties, there were no significant differences in seasonal mean LAI and LAImax. Seasonal mean LAI varied significantly among seasons. As the growing season advanced towards warmer and drier periods at later growth stages, weekly mean LAI among the treatments diverged. The wet growing season in 2018 exhibited similar LAI trajectories among treatments. Except for WW, the highest seasonal mean LAI was observed in humid 2018.

In Sin WW and Sig WW, LAImax was 35.8% and 43.5% higher than in their corresponding RO, on average, over the three seasons. A previous investigation found leaf cell division and enlargement to be more sensitive to water stress, reducing leaf size, causing the protoplasm to dehydrate, and reducing tissue expansion due to reduced cell division compared to leaf growth in irrigated soybean (Khan et al. 2014). Variability in the number of trifoliate layers and the location of the leaf layer with maximum assimilatory surface size were not pronounced. Similar to the results by Lémeur (1973) for soybean, lower leaf density was observed at low trifoliate inclination of the upper canopy height. Water-deficit exerts its effect on soybean growth by reducing LAImax and plant height. The number of trifoliate layers was hardly impacted by variety.

During water stress-induced senescence, the movement of nutrients in plants from their source tissues to reproductive organs is controlled by the regulation of physiological mechanisms (Sade et al. 2018). Two groups of hormones impact leaf senescence antagonistically, cytokinins and abscisic acid (Verslues 2016). In this study, early leaf senescence in water-stressed crops may be connected to a decline in endogenous cytokinin levels in leaves, although this would need to be assessed.

Inhibition of cell elongation caused by the interruption of water flow from xylem to cells (Wijewardana et al. 2018) decreases plant photosynthetic efficiency (Farooq et al. 2009) in water-stressed soybean. More details of this process may be found in Sade et al. (2018). A detailed eco-physiological analysis of the influence of water stress on soybean is beyond the scope of the study. Lower top assimilatory surface size in soybean with limited water might be due to decreased photosynthetic rate and worsened partitioning of dry matter to produce more leaves, as was shown in Iran (43°N) during 2016 (Mondani et al. 2019). Our results agreed well with those of Mondani et al. (2019), who also reported a decline in dry matter production of water-stressed soybean due to a reduction in maximum leaf expansion.

The results of this study did not support our hypothesis that change in canopy architecture and solar radiation penetration could be explained by differences among the two studied soybean varieties.

The magnitude of variation in PAR penetration between the upper and lower half of the canopy suggests that canopy structure is more erectophilic in the upper part, but more planophilic in the lower half of the plant. At the bottom of the canopy, RO consistently showed higher I/I0 ratios than WW or P for the entire three-year period. The greatest difference between RO and WW occurred during the warmest year (2019). Stand architecture affected the light environment of the canopy, decreasing PAR flux density in leaf layers of WW compared to RO. The fraction of incident beam radiation that penetrated the canopy (ground surface) in RO agreed well with an about 5% value by Sameshi (1995) in Tsukuba, Japan (36°N), at a plant density of 33 m−2 with a LAI of about 6. When averaged across three seasons at flowering, 97.7, 95.1 and 97.26% of incoming PAR were intercepted by the canopies of WW, RO and P, respectively. The close values in these treatments were not unexpected, because intercepted radiation in soybean generally reaches 95–98% of daily solar radiation when LAI >4 (Mondani et al. 2019). This PAR fraction is available for scattering, transpiration and photosynthesis (Campbell and Norman 1998).

Similar to the findings of Bajgain et al. (2015), RO in our study exhibited the highest k value of 0.46, versus only 0.42 and 0.44 in WW and P, respectively. More upright upper leaves resulted in higher radiation transmission to the upper third of the plant stand. It is important to note that k, when measured at high solar elevation, may result in large underestimations of daily radiation interception, because k, among other factors, is also strongly affected by solar angle (Campbell and Norman 1998).

The ratio of harvested seed yield to total aboveground biomass ranged from 20.8% (Sig RO) to 31.1% (Sin WW), values that were lower than the 40% published by Hussain et al. (2019) in Michigan, USA (42°N) for rainfed soybean. Water stress deteriorated the yield/TDM ratio to about 10%. Prasad et al. (2008), who revealed that drought desiccated the soybean endosperm resulting in limited embryo size, also highlighted that water stress can significantly influence yield-related processes when it occurs at a sensitive stage of crop growth, namely during flowering, as occurred in this study.

As expected, water treatment had a significant effect on all yield components: RO decreased yield whereas WW increased yield. Although the treatment × variety interaction was statistically significant only for seed yield and protein mass, the difference between varieties was also slightly different for each treatment in the case of protein mass. These interactions showed similar patterns. There were some differences between varieties in WW and P, but these differences diminished in RO. Since there was no considerable difference between varieties in 1000-seed weight and TDM, it can be stated that the difference between varieties was negligible during strong water stress.

A projection of the observations onto the plane spanned by QuantityPC and QualityPC provides a clear and straightforward overview of the year-treatment-variety combinations (see Fig. 8). When interpreting this figure, it is important to keep in mind that VPD was highest in 2017 and lowest in 2018. When comparing the treatments, they are markedly separated along the QuantityPC but strongly overlapped along the CompositionPC. Only RO shifted the composition towards a greater percentage oil content. When comparing the years, they completely overlapped along the QuantityPC and somewhat separated along the CompositionPC. The dry year, 2017, shifted toward a greater percentage oil content, while the most humid year, 2018, shifted toward a greater percentage protein content. Overall, drought stress, caused by either water deprivation or low air moisture, shifted the composition of soybean toward a higher percentage oil content. Comparing the two varieties, in the case of P, there was no difference along the CompositionPC and hardly any difference along the QuantityPC. Sin produced slightly higher yield. In the case of RO, there was no difference at all along the QuantityPC and along the CompositionPC only in 2018. Similarly, in the case of WW, differences among varieties were observable only in one year, 2017. Overall, in the case of water stress, higher air moisture accentuated differences among varieties, but in the case of unlimited water, higher VPD caused differences between the two varieties. In 2018, when there was low VPD, in RO, Sin produced more oil than Sig.

Similar to the results of Akhtar et al. (2019), one of the carbohydrate sinks, the seed’s oil content, was significantly greater in warm 2017 than in wet 2018 and dry 2019. The only difference in weather conditions between the two warm seasons (2017 and 2019) might be the higher VPD during 2017 relative to 2019. Lobell and Asner (2003) highlighted the extreme sensitivity of soybean to atmospheric moisture content during the reproductive stage. Kross et al. (2015) reported that mid-seasons with higher VPD (2011–2012) lowered the mean LAI and yield of rainfed soybean in Canada (45°N).

The three-season mean oil content of 20.4–20.5% in P (both varieties) was slightly higher than the 19.1% oil registered by Procházka et al. (2017) in a rainfed soybean variety, Merlin, grown from 2012 to 2015 in the Czech Republic.

The seed’s protein content followed an opposite trend to oil content, as was also observed by Aydinsakir (2018). The higher protein content of seeds was related to the lack of crop water stress. However, these slight changes in protein content (about 2–3%) in WW relative to P were not always significant throughout the entire three-year period in our study. A minimal effect of additional water on seed protein was reported by Sweeney et al. (2003) when irrigation was applied between 1991 and 1994 to an early maturing soybean variety at the reproductive stage. Protein content in our study was comparable to values reported by Candogan and Yazgan (2016) ranging from 30.7% to 35% in irrigated and rainfed soybean, variety Nova, in Bursa, Turkey (40°N), between 2015 and 2016.

RUE, which drives TDM production through photosynthesis per unit leaf area, is also affected by solar radiation. Water stress caused a decline in photosynthesis (Sinclair and Muchow 1999), which probably reduced leaf expansion (LAImax) also reducing the final TDM and RUE of RO, irrespective of the variety. RUE of 1.31 g MJ PAR−1 (Sig WW) and 1.32 g MJ PAR−1 (Sin WW) for the entire three-year period fell within the lower range of reviewed and reported soybean RUE values (1.32–2.52 g MJ PAR−1) by Sinclair and Muchow (1999) from six different studies. Singer et al. (2011) obtained a slightly higher soybean RUE of 1.44 g MJ PAR−1 at Ames, USA (42°N) under rainfed conditions during 2015.

In accordance with lowest TDM and increased seasonal incoming solar radiation during the extremely warm 2019, lowest RUE was obtained. With the exception of Sig WW and Sin P, significant increments in RUE ranged from 21.9% (Sin RO, p < 0.04) to 34.8% (Sin WW, p < 0.001) in 2017, compared to 2019 values. About 1 °C cooler Ta in the summer of 2017 relative to 2019, may have provided soybean with more favourable climatic conditions for biomass accumulation. The 34% higher PR total in 2017 decreased seasonal mean VPD by 0.8 kPa, which might also be desirable for soybean, which is of tropical origin. Due to higher VPD in 2017, ET totals were also slightly higher than in 2019.

This study also confirmed that RUE derived from TDM and accumulated PAR were not sensitive over a range of LAI<1.0 (Sands 1996). Suyker and Verma (2009) also reported no relative increase in intercepted PAR after full canopy cover in soybean, even when LAI exceeded 3–4. Variations in the literature related to RUE in soybean might also be attributed to geographic and varietal differences (Wang et al. 2001).

Conclusion

Based on our results, the water stress during flowering was a threat to both studied soybean varieties and had the greatest impact on yield and RUE. This is not surprising as the canopy structure was almost the same in both varieties for each water treatment. Mean LAI of the three seasons declined with decreasing water supply. In agreement with LAI values, the highest k was observed in RO, probably due to a stress-induced increase in I/I0 ratios at the bottom of the water-stressed canopies. This highest k might be the result of enhanced leaf wilting of lower leaf layers in RO compared to WW and P. The lower seed yield in water-stressed soybean was attributed to smaller sizes of seasonal mean LAI that contributed to less photosythates (production), irrespective of the variety. Under water stress, decreased yield potential was observed in both varieties. However, Sin still manifested favourable yield under unlimited watering and rainfed conditions. It is important to mention that only the water demands of the two soybean varieties used in this study, which are both of the same indeterminate type, differed (Karintia 2019).

Soybean that experienced water stress during flowering produced lower seed yield and TDM due to a reduced PAR transmission, a lower assimilatory surface area, and a decline in RUE, compared to soybean under no water stress. Rainfed soybean ameliorated, to a certain extent, the reduction in seed yield caused by RO, probably due to the even distribution in PR across the three-season study.

In the LME model of ET, the year-to-year variability of the coefficients was much lower than the absolute values of the coefficients, indicating that the same biophysical law controls ET. Our study had one limitation: given the limited number of observations, only the main effects of explanatory variables could be investigated. However, it would also be worth investigating the interaction between LAI and meteorological variables in a future study. An additional variable that describes the developmental state of the plant would also be useful in a more complex model. In addition, the high value of the coefficient in the residual AR(1) process indicates the need for such a variable.