Digital twins to quantify the impact of growing variability on the harvest quality of orange

The quality of citrus fruit is influenced by various growing conditions, including weather. However, the impact of weather differences between growing regions on citrus quality at harvest is not well understood. This study employs a mechanistic-driven digital replica of the growth process of a Valencia orange from fruit set until harvest to quantify this impact. Measured daily temperature, humidity, rainfall, solar radiation, and wind speed data from various orange growing regions of South Africa, including Citrusdal, Nelspruit, Letsitele, and Sundays River Valley (SRV), inform and drive this digital model. The results suggest that the differences in weather conditions between growing regions significantly affect fruit diameter (FD), fruit weight (FW), rind thickness (RT), rind weight (RW), total soluble solids (TSS), titratable acidity (TA), and ◦ Brix/acid ratio of oranges at harvest. The differences between growing regions led to variations of up to threefold for FD, twofold for FW, RT, RW, TSS, ◦ Brix/acid ratio and up to fourfold for TA upon harvest. Particularly, oranges produced from warmer Letsitele and Nelspruit regions are found to be larger, sweeter and less acidic compared to those from coastal SRV upon harvest. The study also reveals the impact of the fruit growth process on the temperature gradient within the fruit, which varies across growing regions. The maximal temperature difference between the fruit core and surface during the growth process ranges from 4 ◦ C to ~7 ◦ C. These variations in fruit temperature gradient could lead to variations in temperature-driven quality decay of fruit from different climatic regions at the start of their postharvest journey. These findings provide valuable insights for the citrus industry, aiding to optimize practices, harvest planning, and postharvest logistics. The output of this digital twin will help identify areas needing extra precooling to extend postharvest shelf life and minimize quality decay. Real-world use allows growers to schedule harvests based on regional weather conditions.


Introduction
Citrus is one of the most cultivated fruit crops in the world and is grown in more than 80 countries, mostly under tropical and sub-tropical climates (M.S. Ladaniya, 2008).The fruit is known for its thirst-quenching ability, refreshing fragrance, and abundant vitamin C content (Vashisth and Kadyampakeni, 2019;M.S. Ladaniya, 2008).However, differences in their growing conditions can affect the quality of oranges at harvest (Davies, 1997;Lado et al., 2014).
The fruit growth and physiological quality properties of oranges often depend on several pre-harvest factors, such as horticultural practices, and climatic variability (Ashebre, 2011).However, weather conditions have the most predominant effect on fruit growth and quality (Ladaniya, 2008).A suitable growth temperature is a critical factor for successful citrus production, with optimal growth conditions between 13 and 32 • C. On the one hand, a very high temperature reduces fruit growth by decreasing total tree CO 2 assimilation and increasing fruit drop (Sarkar et al., 2021).These high temperatures also affect internal fruit quality parameters, such as sugar and acid contents (Balfagón et al., 2022).On the other hand, too low temperature results in frosts and freezing damage to fruit tissue (Sarkar et al., 2021).Additionally, moisture stress during fruit formation and development affects fruit quality (Ali et al., 2021).This means that weather conditions during crop growth in different production regions will affect crop growth and can result in variations in fruit quality (e.g., fruit size) at harvest.However, the actual impact of weather variation on fruit quality at harvest remains unclear.One approach to better understand and quantify this impact is to predict the evolution of fruit quality during growth and development under distinct growing conditions.
Crop growth models have been explored to describe the growth and development of plants under different growing conditions.For this purpose, mechanistic and empirical models such as LINTUL (van Ittersum et al., 2003), SUCROS (van Ittersum et al., 2003), WOFOST (BAI et al., 2020), SIMPLE (Zhao et al., 2019), von Bertalanffy (Tijskens et al., 2016;Ribeiro et al., 2018), and EPIC (Yang et al., 2019) have largely been deployed to predict the growth of fruits and vegetables in the past decades.However, relatively few studies explicitly modeled the impact of growing conditions, such as weather variability, on very few citrus fruit quality at harvest (Nawaz et al., 2020).Rarely is the effect of the weather variability on multiple physiochemical properties of citrus upon harvest accounted for.
To bridge this gap, we have pioneered a novel approach by creating a virtual model of oranges that encompasses their growth cycle from fruit

Symbols ρ i
Density (kg m − 3 ) C pi Specific heat capacity (J kg Volumetric heat of respiration (W m − 3 ) ∇ Spatial derivative operator (-) T i Fruit temperature (K) i Fruit pulp and rind domain n Unit vector normal to the surface (-) Q o Convective heat flux at the air-fruit interface (W Convective mass transfer coefficient (s m − 1 ) P v,rind Vapor pressure on the fruit surface (Pa) P v,air Ambient vapor pressure (Pa) P sat Saturated pressure (Pa) a w Water activity below the fruit surface (-)  set until harvest.This virtual modeling approach is commonly referred to as a digital twin or digital shadow.Our digital twin of oranges, from fruit set until harvest, is a virtual model that replicates the key characteristics of oranges in a computer-based environment.This model incorporates essential fruit properties and simulates a wide range of physiological, metabolic, and thermal processes that occur during the growth and development of oranges.It is also linked to real-world processes through sensor data (Onwude et al., 2022;Defraeye et al., 2021).Using this approach, the impact of weather variability between growing regions on the quality of oranges upon harvest was evaluated.
Our study includes metrics such as fruit weight (FW), fruit diameter (FD), rind thickness (RT), rind weight (RW), total soluble solids (TSS), and total titratable acidity (TA).We also simulate the physics-based thermal stress experienced by fruit as it grows.This allows us to ascertain if there is a large temperature gradient of the fruit at the time of harvest, a crucial factor influencing the postharvest longevity of the fruits.

Study area
This study considered 'Turkey' Valencia oranges produced over the 2018/2019 and 2019/2020 consecutive seasons during the generic fruit growth stages phase I (1 September -31 December), II (1 January -15 May) and III (16 May -31 July/August).Four major citrus growing regions of South Africa (SA) with distinctly different climate zones (Fig. 1) were used for data collection.These areas are Letsitele in Limpopo, Nelspruit in Mpumalanga, Sundays River Valley (SRV) in Eastern Cape, and Citrusdal in Western Cape.More details of the production information of the experimental sites can be found in Table S1 in the supplementary material.

Growing regions and fruit sampling
This study was carried out for 'Turkey' Valencia oranges grown in the four regions, as Valencia is the most produced orange cultivar in South Africa, accounting for 44 % of citrus exports from southern Africa (Citrus Academy 2017).'Turkey' is one of the first Valencia cultivars to mature and is harvested between the last week of May and June in warmer climates.For this study, we harvested orange fruit from five commercial orchards in each region.Ten healthy and uniformly sized adjacent trees were specifically chosen and tagged within each orchard for consistent harvesting across two production seasons (Table S1).The decision to sample from adjacent trees within the orchard was guided by the intention to minimize variability within the selected sample.This approach ensured a more uniform representation of the prevailing orchard conditions.Fifteen fruit per tree were harvested from the outside up to 30 cm into the canopy to obtain 150 fruit per orchard.All fruit were harvested at commercial harvesting maturity as determined by the citrus producers per region; therefore, the harvesting times differed between regions.More details are available in (North-Dewing, 2023).

External and internal fruit quality evaluation
The external fruit quality, such as fruit total weight (g) and rind weight (g) of each orange fruit at harvest, was determined using an electronic scale NBK-30 (Model NWH 10,422, UWE South Africa).The fruit size (mm) and rind thickness (mm) were measured using a caliper (CD-6′' C, Mitutoyo Corp, Tokyo, Japan).
The internal quality of the fruit at harvest was assessed by cutting along the longitudinal plane of the fruit for juice extraction using a citrus juicer (8-SA10, Sunkist®, Chicago, USA).Pulp particles were removed from the juice by straining through a muslin cloth to produce a pure citrus juice sample for each replicate.Titratable acidity (TA) was determined by titrating 20 ml of juice against 0.16 N sodium hydroxide.Phenolphthalein was used as an indicator and titration was complete when the liquid turned pink in color.Data was expressed as a percentage citric acid content.The total soluble solids (TSS) of the fruit pulp (measured as • Brix and expressed as% TSS in the pulp) was determined using a digital refractometer (PR-32 Palette, ATAGO CO, Tokyo, Japan) (North-Dewing, 2023).The TSS/TA ratio was calculated from those data ( • Brix ÷%Acid) (Jayasena and Cameron, 2008).

Weather data
Weather data per region from September 2018 to August 2020, covering the entire fruit development period for two consecutive seasons, were obtained using the Agricultural Research Council (ARC) automated weather stations.Daily minimum and maximum temperature ( • C), minimum and maximum relative humidity (% RH), and rainfall (mm) data were obtained.Wind speed (U speed , m s − 1 ), short-wave (W m − 2 ), and long-wave (W m − 2 ) radiation data during the growth process for these locations were collected from NASA, via https://power.larc.nasa.gov/data-access-viewer/.Tables S2 and S3 present the selected weather differences between the four major citrus-growing areas of South Africa.
The vapor pressure deficit (VPD), which is a measure of the difference between the amount of moisture in the air and the maximum amount of moisture the air can hold when it is saturated at a given temperature, was also estimated.VPD (kPa) was derived from temperature and humidity parameters using the following calculation (Eq.( 1)) (Koverda): where T is the average temperature [((Maximum and RH is the average relative humidity [((Maximum + Minimum)/2) (%)].
Heat units are frequently used to describe the timing of biological processes, assuming no other limiting factors such as water and nutrition stress exist (Hardy and Khurshid, 2021).According to (Hardy and Khurshid, 2021), HU ( • C) is a useful tool for assessing the suitability of a region for growing citrus, as well as for estimating the duration of different growth stages and predicting the time of fruit maturation.
Heat Units (HU) were calculated for each day (i) per region using the following formula (Eq.( 2)) (Hardy and Khurshid, 2021): where T max is the daily maximum temperature, T min is the minimum daily temperature, and T b is the minimum temperature to sustain growth (for citrus, T b = 13 ∘ C).

Digital twin configuration of a fruit
Digital replica of 'Turkey' Valencia oranges were created to quantify the impact of weather conditions of citrus growing regions on the fruit quality during growth.We modeled a single fruit as a two-dimensional axisymmetric geometry of a sphere (Fig. 2).The domain was divided into two sections of the fruitthe rind (initial thickness ((rind thick,0 ) = 1.8 mm) and the fruit pulp (initial pulp radius (r pulp,0 ) = 2.5 mm).The initial values are based on the experimental data of Valencia orange from the literature (Bain, 1957;Abu-Goukh et al., 2010).Using the EPIC crop growth model, we first simulated the daily fruit weight over time based on heat unit (HU) and maximum fruit weight value from fruit set until harvest.The EPIC model was calibrated with real data and linked to daily weather data during the growing season obtained using automated weather stations (see Section 2.2.3).The EPIC model can simulate several fruit using unique input parameter values for each fruit.This mechanistic model can also simulate the growth process for both annual and perennial crops under varying environmental conditions (Jame and Cutforth, 1996;Williams et al., 1989).Other fruit quality parameters, including fruit diameter, rind thickness, rind weight, total soluble solids, and titratable acidity, were modeled using empirical regression models based on daily potential fruit weight (see Section 2.4).The finite element method was then used to simulate the physics-driven thermal stress of fruit during the growth process to determine the fruit temperature.Taken together, we developed a semi-mechanistic digital twin of Valencia oranges.Using this digital twin, we quantified the impact of growth conditions between the four citrus growing regions on important external and internal fruit quality properties at harvest.

Mechanistic model for pre-harvest fruit growth process
A mechanistic model was developed to predict the fruit growth process, internal and external fruit quality at harvest, and the fruit temperature evolution.This model couples several empirical functions with physics-based models.Additionally, empirical regression models were incorporated to predict FD, RW, RT, TSS, and TA based on the predicted FW using the EPIC model.Fig. 3 illustrates a flow chart outlining the methodology for creating a digital twin of oranges grown in various regions of South Africa, spanning from fruit set to harvest.

EPIC crop growth model
The EPIC crop model is based mainly on a collection of empirical functions for various plant growth processes (Williams et al., 1989).These processes are mainly modeled based on potential biomass growth, nutrient absorption, water use, crop yield, and growth stress.This study used a modified EPIC model to simulate the fruit weight of a single Valencia orange on a daily scale (Yang et al., 2019).The leaves, flower bud, and trees are not taken into account.Uniform expansion on all fruit sides during growth was also assumed.
The potential fruit growth depends on weather variables such as maximum and minimum daily temperature, and heat units of the previous day (Fig. 3A).The fruit growth process is based on the daily heat unit calculated from Eq. ( 2) (Yang et al., 2019;Hardy and Khurshid, 2021;Williams et al., 1989).
Based on the effective heat unit values, we then calculate the heat unit index (HUI) with values ranging from 0 at the start of fruit set to 1 at physiological maturity (time to harvest) (Yang et al., 2019;Williams et al., 1989): where HUI k (0-1) is the heat unit index on day k and PHU is the potential heat unit for crop maturation.PHU ( With Eq. ( 3), we compute the heat unit factor (HUF) (Yang et al., 2019;Williams et al., 1989): where a 1 and a 2 are two crop parameters calculated by solving Eq. ( 4) simultaneously using two distinct pairs of arbitrary values for HUI (0.7941 & 0.6790) and HUF (0.9964 & 0.9901).We solved for the two crop parameters using the numerical iteration method to approximate the values of a 1 and a 2 that satisfy both equations simultaneously.The values for a 1 and a 2 in this study are 0.04 and 7.4, respectively.Note that a 1 governs how fast the crop biomass accumulates in the early stages after planting while a 2 influences how the growth rate slows down as the crop approaches its maximum potential biomass (Bain, 1957).
To predict the growth process of a single Valencia orange, we use a relation between heat unit, maximum fruit weight, and the minimum crop stress factor (REG) (Yang et al., 2019): where FW k is the fruit weight on day k (g), and FW max is the maximum fruit weight (g).Using the experimental data of FW max estimated over different growing seasons per region, we modeled the fruit growth process of 150 oranges from 10 trees in five orchards within a region to capture microclimate variability.A graph of predicted fruit size with experimental fruit size data from the literature is shown in the supplementary data (Figure S1)

Empirical regression model for external and internal fruit parameters
In this study, the empirical regression functions used to model fruit diameter (mm), rind weight (g), rind thickness (mm), total soluble solids ( • Brix TSS), and titratable acidity (% TA) upon harvest were selected based on the following features: 1.They are monotonic, i.e., they generate varying straight lines or curves that are always decreasing or increasing, consistent with the general knowledge of physiochemical fruit properties during growth and development.2. The model consists of as few independent parameters as possible to describe the given data sets with acceptable accuracy.3. The model-predicted data agree reasonably well with the experimental data.
Based on these features, polynomial, linear, power, logarithmic, piecewise and exponential functions were identified as likely to predict orange's external and internal qualities accurately during the growth process.
All the above models were fitted to data obtained from the literature using ORIGINPRO 2022 (64-bit) SR1 (Government) (OriginLab, Northampton, Massachusetts, USA).The performance of the models was evaluated by comparing the coefficient of determination (R 2 ) and the root mean square error (RMSE).Acceptable predictions were determined based on a criterion of a coefficient of determination (R 2 ) greater than or equal to 0.5, along with a minimized root mean square error (RMSE) (Yang et al., 2019;Roy et al., 2009;Roy and Roy, 2008).The power function gave a reasonable prediction of fruit diameter (FD, mm) and rind weight (RW, g) during plant growth and are represented as: The models were calibrated using growth data of 1067 Valencia oranges produced in Gosford district of New South Wales (Bain, 1957).The calibration plots are shown in the supplementary material (Figure S2-S3).Similar relationship between single fruit weight (FW) and fruit diameters (FD) or single fruit weight and rind weight (RW) for fruits and vegetables has also been reported in the literature (Yang et al., 2019;Turrell et al., 1969;Sasikumar et al., 2021).
The rind thickness (RT, mm) for Valencia oranges sometimes has an unimodal relationship with fruit size during the plant growth process (Reuther et al., 1969;Ladaniya and Mahalle, 2011).To capture this relationship, a piecewise empirical model based on fruit size (mm) was developed and calibrated with data from (see supplementary material for details-Figure S4) (Reuther et al., 1969).This model provided the most accurate representation of rind thickness compared to polynomial, power, logarithmic, and exponential functions: Similar to literature reports on other fruits (Parra-Coronado et al., 2017;Colauto Stenzel et al., 2006), polynomial models best explained the TSS quantified as • Brix and TA (%) of orange pulp during growth as expressed below: where RP is the heat of respiration expressed in W kg − 1 (see Section 2 in the supplementary material for details).The TSS model was calibrated using fruit weight and TSS growth data of Valencia orange from (Quaggio et al., 2006).The model for TA was calibrated using citrus growth data of grapefruit (A.B.Abu-Goukh et al., 2010).The calibration plots are shown in the supplementary material (Figure S5 -S7).

Thermal model
Following the approach put forward by (Saudreau et al., 2011;Saudreau et al., 2007), we solved the spatial and temporal heat load (thermal stress) variations on the fruit surface due to its interaction with the environment using the equation below.To reduce the complexity of the model, we simplify the system by assuming thermal equilibrium between all components and phases inside the fruit.In addition, we assumed heat loss from only respiration activity inside the fruit.
ρ i is the density (kg m − 3 ), C p,i is the specific heat capacity (J kg − 1 K − 1 ), T i is the fruit temperature (K), k i is the thermal conductivity of the material (W m − 1 K − 1 ), ∇ is the spatial derivative operator, and i represents the fruit pulp and rind domain.The material and thermal properties used in this study are given in Table 1.Q resp,i (W m − 3 ) represents the volumetric heat of respiration during plant growth.We calculated Q resp,i as a product of heat of respiration (RP, W kg − 1 ) and fruit pulp density (ρ pulp, kg m − 3 ).The energy transferred through the boundary layer by conduction is carried away by convection through the air.We modeled the energy exchanges as boundary conditions by assuming that the normal heat flux at any point of the fruit surface was equal to the sum of sensible energy loss or gain by convection (Q o ), losss of energy by transpirational cooling (Q E ) and radiation (Q R ), taking into account the direction of heat flow along the surface normal: where n is the surface normal vector, i.e. the vector that is perpendicular (normal) to the fruit surface at a specific point.
We indirectly modeled the convective heat flux (Q o ) between the fruit surface and the atmospheric air using the heat transfer coefficient as given below: where Q o is the convective heat flux at the air-fruit interface (W m − 2 ), h c is the convective heat transfer coefficient (W m − 2 K − 1 ), T air,ave is the average daily air temperature during plant growth process (K), and T rind is the temperature at the fruit surface (K).The entire system (air and fruit domain) was assumed to be at an initial temperature for growth, T ini,r (26.0 • C) for Nelspruit, (24.0 • C) for Letsitele, (22.0 • C) for Citrusdal, and (23.0 • C) for SRV.These values correspond to the average ambient temperature recorded at each location on the first day of fruit set for consecutive growing seasons.
To estimate the spatially-varying convective heat transfer coefficient (h c ) over the fruit surface, we used an empirical relation of Nusselt number (Nu) for forced convection of a spherical shape as shown below (Cengel and Ghajar, 2011;Whitaker, 1972): .4Re 0.5 + 0.06Re 0.667 ) Pr 0.4 ( μ air μ air,rind ) 0.25 (16) where k air is the thermal conductivity of air (= 0.02 W m − 1 K − 1 ), Re (=FD U speed /υ air ) is the dimensionless Reynolds number as a function of the airspeed (-), Pr is the Prandtl number for air (= 0.73).The values for µ air and µ air, rind correspond respectively to the absolute viscosity of air (kg m − 1 s − 1 ) and the viscosity of air at fruit surface (kg m − 1 s − 1 ), which were considered to be equal, and υ air is the kinematic viscosity of air (= 1.43E − 5 m 2 s − 1 ).Buoyancy forces were neglected as forced convective heat transfer was dominant (see Section 2 in the supplementary material for details).
Transpiration involves the evaporation of water from the surfaces of leaves and other aerial tree parts.This process leads to the cooling of the tree and its surroundings, including the fruits.As water evaporates from the surface of the fruit, it carries away heat, leading to cooling of the fruit.The transpiration heat flux (Q E , W m − 2 ) was modeled using the formula (Saudreau et al., 2007): where k t is the convective mass transfer coefficient (s m − 1 ), P v,rind , is the vapor pressure on the fruit surface (Pa) and P v,air is the ambient vapor pressure (Pa) (see supplementary material for more details).
We determined k t from the contribution of the resistance due to moisture migration through the rind (k rind , s m − 1 ) and the resistance to mass transfer due to the air boundary layer (k air s m − 1 ): The air film mass transfer coefficient (k air ) was estimated based on the airspeed (u speed , m s − 1 ) using the Sherwood correlation for a sphere, as presented below (Becker et al., 1996;ASHRAE 2018).
where Sc (=υ air • δ wv,air − 1 ) is the Schmidt number (-), δ wv,air is the diffusion coefficient of water vapor in the air (m 2 s − 1 ).We modeled the radiative heat flux (Q R , W m − 2 ) with the following equation (Saudreau et al., 2011): To simplify the model, we assumed no radiation -radiation interaction between fruit due to direct exposure to sunlight.We also assumed isotropic radiation fields on the fruit surface during growth.
Here, a sw is the reflectance of orange rind to short-wave radiation (= 0.63) (Sharpe and Barber, 1972;Gaffney, 1973), R sw is daily short-wave radiation received by the fruit (W m − 2 ) from fruit set until harvest per regions, ε is the emissivity of infrared radiation for orange (= 0.98) (Caselles et al., 1988), R lw is the daily long-wave radiation received by the fruit (W m − 2 ) from fruit set until harvest per regions, and σ is Stefan-Boltzmann radiation constant (= 5.67E − 8 W m − 2 K − 4 ) (Onwude et al., 2018).

Numerical simulation
The pre-harvest semi-mechanistic model was implemented in COMSOL Multiphysics (version 6.1), a commercial finite element-based software.Using the linear interpolation function, we solved the single fruit growth model and external and internal fruit parameters (Eqs.(2)-12) (Fig. 3B-G).To capture the looping of weather data as a function of time in the model, we used the "primitive relation function" under interpolation function for each day after fruit set.We simulated the transient conduction, radiative, and convective heat fluxes in the fruit during the growth process using the 'Bioheat Transfer' physics.As the boundary of the computational domain translates and deforms during the growth process, we modeled the uniform deformation using the "Deformed Geometry" physics and coupled it with the heat transfer physics.We considered the fruit diameter change by integrating over the normal fruit displacement.The coupling was performed in a way such that the rind (FD/t) and pulp ((FD-RT)/t) domain velocities due to deformation are applied as the normal mesh velocities.A normal mesh was adopted for the simulation after mesh-sensitivity analysis (see supplementary material -Figure S12).We used parametric sweeps of transient studies to simulate 150 oranges from 10 trees in five orchards within a growing area (Fig. 3F).Quadratic Lagrange geometry shape function, with automatic remeshing selection and a MUMPS (MUltifrontal Massively Parallel sparse direct Solver) fully-coupled direct solver was used for the simulation.The remeshing condition was set to "distortion" to reduce mesh distortion and enhance convergence.The relative tolerance was set to 10 − 7 , while 0.05 was the scaling factor.The time-stepping for the simulation was set to 86,400 s, corresponding to the time interval of weather data recorded.

Statistical data analysis
To investigate the impact of climatic variability on the quality of oranges at harvest, a one-way analysis of variance (ANOVA) was performed.Prior to this, Levene's test was used to test the homogeneity of variance.Digital twin output for different growing regions was presented as median, 75th upper and 25th lower quartiles (box limits), and 1.5 × the interquartile range (IQR, whiskers) with a 0.95 confidence D. Onwude et al. level.These statistical indicators were used to show the impact of climatic differences between growing areas on FW, FD, RW, RT, TSS and TA.Fisher LSD test at p ≤ 0.05 significant level and 95 % confidence interval was used to assess if there were significant differences in the fruit quality attributes within groups.We used the Tukey test to compare the simulated external and internal fruit quality data at harvest with the measured experimental data of the 2019 and 2020 harvest years.Pearson's correlation coefficients (Pearson's r) at p = 0.05, 0.01, and 0.001 significant levels were used to statistically determine the correlation between the fruit quality attributes and climatic parameters.We also performed a sensitivity analysis to assess the contribution of each model parameter to the variability of the different quality attributes at harvest.Eq. ( 21) (Jorgensen, 1995) was used for this sensitivity analysis.The result of the significance of specific model parameters such as T 0 , FW max , T max , T min , T g , and REG in accurately predicting the influence of growing conditions on fruit growth and quality of oranges is presented in Section 3 of the supplementary material.All analyses were conducted using ORIGINPRO 2022 (64-bit) SR1 (Government) (OriginLab, Northampton, Massachusetts, USA).
where S is the percentage scaled relative sensitivity, ∂x denotes the change in fruit quality parameter values, x represents the fruit quality parameter, p represents the model input parameters, and ∂p signifies the change in value for the model input parameters.These changes are considered at a ± 20 % level on a temporal scale.

Model verification for fruit quality parameters
The predicted fruit weight (Fig. 4A) (using Eqn.(2)-5) and diameter (Fig. 4B) (using Eqn. ( 6)) of Valencia orange for all growing regions agree reasonably well with that of our average experimental data (max difference ≤ 0.5 g for FW; max difference ≤ 3 mm for FD) for consecutive growing seasons.The accuracy of the models is relatively high for fruit weight, with an average difference of less than 0.5 % between the predicted and experimental data across all growing regions.The average   D. Onwude et al. difference between predicted and experimental data is less than 5 % for fruit size, indicating a reasonably accurate performance across all production regions.Our model also performed well in predicting the fruit weight and size of oranges for individual growing seasons (refer to Figure S13 in the supplementary material).This finding suggest that the virtual modeling approach effectively captures the influence of climatic variability on fruit weight and size at harvest.Nevertheless, we acknowledge potential uncertainties, given that the model verification is based on measured data of two consecutive growing data.
Similarly, the consistency between our predicted values and average experimental data over two seasons is evident in Fig. 5 for various other quality parameters related to Valencia orange.The predicted values for rind weight (Fig. 5A) (Eqn.( 7)), rind thickness (Fig. 5B) (Eqn.8-10), TA (Fig. 5C) (Eqn.( 12)), and TSS (Fig. 5D) (Eqn.( 11)) of Valencia orange are reasonably consistent with our average experimental data for two consecutive growing seasons across all regions.The maximum variance of ≤ 7 g for rind weight, ≤ 0.2 mm for rind thickness, ≤ 0.3 % for TA, and ≤ 0.5 • Brix for TSS are obtained across all production regions.Our model also performed well in predicting these quality parameters of oranges for individual growing seasons (refer to Figure S14 in the supplementary material).These findings also demonstrate the effectiveness of our modeling approach in capturing the influence of weather variability on multiple fruit parameters, including rind weight, rind thickness, acidity, and total soluble solids at harvest.This reinforces its utility for understanding and managing variability in orange production across diverse geographic areas.

What is the impact of weather variation between growing regions on fruit quality attributes at harvest?
The various fruit quality outputs of our digital twin at harvest for the 2018/2019 and 2019/2020 production season per growing region based on the coupled mechanistic and empirical models (Eqn.(3)-12) were assessed (see Section 2.6).Fig. 6 illustrates the impact of the differences between growing regions on the external (so FW, FD, and RT) and internal (so TA, TSS, and • Brix/acidity ratio) fruit quality parameters at harvest.There is a significant effect of weather differences between growing regions on fruit weight (FW) (Fig. 6A), fruit diameter (FD) (Fig. 6B), and rind thickness (RT) (Fig. 6C) of oranges at harvest.Oranges cultivated in the coastal SRV region exhibit significantly lower fruit weight and smaller fruit size, compared to oranges grown in the Mediterranean Citrusdal region and in the warmer climates of Letsitele and Nelspruit.Oranges from SRV have thicker fruit rinds than those from warmer Letsitele and Nelspruit regions.The fruit thickness of oranges from Nelspruit does not show a significant difference when compared to that from Letsitele over two consecutive growing seasons.These findings align with our experimental results (results not shown).Generally, the warmer regions of Nelspruit and Letsitele, consistently yielded larger fruits with thinner rinds.This aligns with the observations made by Goldschmidt and Reuther (Goldschmidt, 1997;Reuther, 1988).For internal quality attributes, there is a significant impact of variations in growing regions on acidity (TA) (Fig. 6D), total soluble solids (TSS) (Fig. 6E), and • Brix/acidity ratio (Fig. 6F) of oranges at harvest.Specifically, the acidity levels in oranges from SRV were significantly higher compared to those from Citrusdal, Letsitele, and Nelspruit.No significant differences in TA were observed between oranges from Nelspruit and Letsitele over two consecutive growing seasons (Fig. 6D).Regarding TSS, oranges from Nelspruit and Letsitele significantly differed from those grown in Citrusdal, with Citrusdal displaying the lowest TSS (Fig. 6E).The balance between sweetness and acidity, represented by the • Brix/acidity ratio, in oranges from Letsitele did not differ significantly from those in Nelspruit but showed significant distinctions from those originating in Citrusdal and SRV (Fig. 6F).Likewise, oranges from Citrusdal demonstrated a statistically equivalent balance of sweetness and acidity compared to those from SRV.These findings align consistently with our experimental results.These findings show that the variations in internal fruit quality at harvest are significantly influenced by the contrasting growth conditions among the orange growing regions.
Taking together, the oranges produced from Letsitele and Nelspruit are larger with thinner rinds and have a good balance of sweetness and acidity, contributing to a better taste compared to those from SRV and Citrusdal.This could be because the warmer climates of Letsitele and Nelspruit have higher average growing temperatures and heat units compared to those of cooler climates such as SRV.Oranges thrive in moderately hot, sunny, dry and humid conditions, and as such, they tend to grow better in climates with higher average temperatures and effective heat units (Nawaz et al., 2020).In addition, during stages I and II, which encompass early fruit growth and cell division, the size of the fruit is primarily determined.In warmer climates, the conducive optimum temperatures (see Figure S10A and B in supplementary material) during these stages promote greater cell division and enlargement, leading to larger-sized oranges than those grown in cooler coastal regions.Furthermore, the low relative humidity in stage III (see Figure S10C and D in supplementary material) of citrus production in dry semitropical and subtropical regions could also produce citrus with excellent qualities (Davies, 1997;Goldschmidt, 1997).These qualities include rind color, sugar content, and acidity level.Compared to Nelspruit and Letsitele, the cooler climate of SRV is characterized by lower minimum temperatures.Low temperatures also have depressive effects on fruitlet growth and can result in smaller fruit (Agustí and Primo-Millo, 2020).
Citrusdal and SRV exhibit relatively similar average yearly temperatures and heat units (HU).However, the harvest quality attributes of oranges from Citrusdal are interestingly significantly different from those grown in SRV, and similar to those produced in warmer climates.This may be attributed to variations in temperature, specifically the lower minimum and maximum temperatures observed in Citrusdal, particularly during stage III (Figure S10A-B).Moreover, Citrusdal experiences higher maximum and minimum relative humidity during stage III (Figure S10C and D) in comparison to SRV.Another contributing factor to the disparity in fruit quality between oranges from SRV and Citrusdal is the difference in vapor pressure deficit (VPD) (Figure S10F).In Citrusdal, there is a pronounced decrease in VPD from stage II to III, significantly lower than that of SRV.Lower air temperatures and higher humidity reduce VPD, indicating a lower water loss from the fruit rind.This, in turn, facilitates better growth and fruit quality (Ladaniya and Mahalle, 2011;Parra-Coronado et al., 2017).Therefore, regions with optimal VPD and minimal seasonal variation, particularly during the later stages of fruit growth, as seen in Citrusdal, tend to produce high-quality citrus fruit (refer to Section 4 in the supplementary material, Figure S9).These complex interactions between temperature, humidity, and VPD can explain the large differences observed between the oranges grown in Citrusdal and SRV.However, the influence of weather variability on rind thickness in humid climates

Using digital twins to unveil the spatial-temporal temperature variations of fruit in different growing regions
The digital twin outputs the variations in fruit temperature at both spatial and temporal scales, accounting for changes in physical and thermal parameters during fruit growth and development and between fruit tissues.Fig. 7 compares the digital twin output of temperature variability within oranges produced in Letsitele, Nelspruit, Citrusdal, and SRV (using Eqn.(2)-10; 13-20).The daily changes in surface and core fruit temperatures showed similar patterns for all citrus growing regions (Fig. 7A).We observe differences between atmospheric air temperature and fruit temperature for all growing regions.Similar to atmospheric air temperature, fruit temperature was highest for oranges produced in the hotter and sunnier climates of Letsitele and Nelspruit compared to the cooler climates of Citrusdal and SRV.The core or center temperature for all growing regions exceeded the fruit surface temperature (Fig. 7A).Experimental results similar to those reported were also observed for naval oranges harvested in Citrsudal in 2017 (Goedhals--Gerber and Khumalo, 2020).This means that the accumulation of heat generated by the metabolic activities of the fruit and surrounding tissues during growth and development is greater than the effect of convective and radiative heat transfer at the fruit surface.In addition, there is also the cooling effect of evaporation, which significantly reduces the fruit surface temperature.In most cases, the temperature gradients (difference between core and surface temperature) increased with the time of fruit growth and development for most regions (Fig. 7B).This means that the temperature gradients in the fruit at stages II and III were higher than those at stage I.These significant temperature differences within the fruit can even reach up to 1.5 • C within a very short time (within 5 weeks after fruit set).The maximum temperature gradients during fruit growth varied across different regions.Letsitele and Citrusdal exhibited a maximum gradient of up to ~5.0 • C, Nelspruit had a similar gradient of 4.0 • C, and SRV had the highest maximum gradient of ~7 • C (Fig. 7B).The cooler climate of Citrusdal had the highest average fruit temperature gradient during plant growth, while the fruit from the warmer climate of Nelspruit displayed the lowest average gradient (Fig. 7B).The higher temperature gradients observed in cooler climates than warmer climates can be attributed to the enhanced radiative cooling effect in those regions.Cooler climates generally have lower average temperatures, especially at night, which promotes more significant radiative cooling.Consequently, this temperature differential results in a larger gradient between the core and surface of the fruit.Another contributing factor could be the wider diurnal temperature range commonly experienced in cooler climates.During the day, the fruit surface may be exposed to warmer temperatures, elevating the surface temperature (as seen in Fig. 7A).However, during the night, the cooler temperatures impede heat transfer from the fruit core, leading to a larger temperature gradient.Oranges grown in SRV displayed a more pronounced radiative cooling effect compared to Citrusdal despite Citrusdal having a lower growing temperature (Figure S10A&B).This could be attributed to a greater long-wave radiation cooling effect in SRV compared to Citrusdal (refer to Supplementary material, Figure S11), emphasizing the complexity of factors influencing fruit temperature dynamics.
This finding is significant because the temperature gradient affects enzymatic reactions during storage, transport and processing.However, these temperature heterogeneities within the fruit can become even lower in reality due to the small variability in local air movement within different orchards for different growing regions.This effect can be captured by explicitly modeling both free (see Section 2.4.3) and forced convection, as well as the fruit's transpiration.As the fruit's surface loses water through transpiration, the air inside the fruit becomes cooler and more humid, thereby reducing core temperature.Note that other factors, such as solar intensity, vapor pressure, canopy structures and microclimates, and orchard management practices, may, in reality, also contribute to temperature differences within the oranges in these regions.

Limitations and perspective for the citrus industry
The digital twin of pre-harvest citrus value chains can simulate and predict the impact of variability on the quality evolution of fruit at harvest.The quality parameters at harvest can then serve as an input parameter for a postharvest digital twin from farm to fork.This pioneering study also simulated the impact of a realistic environment (unsteady and non-homogeneous heat fluxes) on the temperature dynamics during the growth of fruit.The final volumetric fruit temperature at harvest serves as an initial condition for the hygrothermal modeling process during postharvest quality evolution of citrus.Current technologies used during citrus production rarely consider measuring fruit temperature during growth.This is because it is challenging to measure directly in the field due to the physical barriers presented by the plant, such as leaves, branches, and other fruit.Obtaining accurate measurements without damaging the fruit or disturbing its growth can also be difficult.
With this digital twin concept, stakeholders in the citrus industry can spatially measure and assess how changes in weather, such as temperature, humidity, wind and VPD patterns, could affect citrus postharvest citrus quality.This can help them optimize their operations, improve decision-making, and identify potential risks and opportunities related to climate change.For example, growers can use the digital twin to evaluate the impact of changing climate patterns of a particular growing region on fruit quality.Ultimately, a digital twin of the citrus plant growth process can support more sustainable and resilient citrus production in the face of increasing climatic and microclimate variability.Nonetheless, as a pioneering preliminary study, this study has some limitations due to the following concerns: • Weather data from meteorological stations may not accurately represent the microclimate conditions in specific orchards.• Model calibration data on fruit quality during growth were obtained from literature sources and may differ from actual growth data of the assed locations.• Model simplification: for instance, the EPIC model used in the study did not account for all physiological and structural changes during plant growth, such as tree size and canopy volume, which can vary between warmer and cooler regions.• The empirical models used for predicting fruit quality evolution may need recalibration for different citrus cultivars.• Certain model assumptions and factors, such as mass transfer during growth, radiation isotropy, wind direction, fruit shape, and energy exchanges between the plant, flower, and fruit, require further investigation and refinement.

Conclusion
This study quantified the impact of weather differences between Letsitele, Nelspruit, Citrusdal and Sunday River Valley (SRV) growing regions of South Africa, on citrus external and internal fruit quality at harvest.We did this using the concept of a digital replica of oranges from fruit set to harvest.We developed semi-mechanistic pre-harvest digital twins of 50 Valencia oranges representing an average of 15 fruits per tree from 10 trees in 5 different orchards per growing region for two consecutive seasons.These digital twins were coupled with measured real-world pre-harvest climate data such as daily average temperature, relative humidity, heat unit, wind speed, and solar radiation, and predicted fruit quality at harvest.Our study demonstrated that climate variability between citrus growing regions significantly impacts fruit quality at harvest.The key conclusions derived from this study are as follows: • Our models gave a reasonable prediction for fruit diameter (FD), fruit weight (FW), rind thickness (RT), rind weight (RW), and total soluble solids (TSS) of oranges at harvest, with an error of less than 5 % for consecutive growing seasons across all regions.• The model showed a maximum difference of ≤ 0.5 g for fruit weight and ≤ 3 mm for fruit size across all producing regions using the average data of two consecutive growing seasons.Given that citrus is typically sold by weight, and the fruit sorting process is significantly influenced by size, the practical implications of this model could be of great interest to the citrus industry.Farmers and exporters can leverage this model to enhance their sorting processes, ensuring an enhanced categorization of fruits.Additionally, retailers stand to benefit by proactively maximizing profits through improved inventory management and strategic pricing based on the predicted size and weight data provided by the model.
• Due to differences in climatic growing conditions, the fruit diameter (FD), fruit weight (FW), rind thickness (RT), rind weight (RW), titratable acidity (TA), total soluble solids (TSS) and • Brix/acidity ratio of oranges from different regions are significantly different at harvest.Therefore, stakeholders in the citrus industry should expect statistically significant differences in the harvest quality of citrus grown in hot semi-arid and humid sub-tropical climates with those of Mediterranean and temperate oceanic climates of South Africa.• Specifically, oranges originating from warmer climates of Letsitele and Nelspruit are anticipated to have higher weight and larger size than those from Coastal SRV.Fruit from these regions also have a more balanced taste than those from cooler climes of SRV.This is because the former locations exhibit a • Brix/acidity ratio of approximately 16, signifying a harmonious balance between sweetness and acidity, while the latter regions show a ratio of around 8. In addition, cooler climates are expected to produce oranges with thicker skin compared to more tropical climates.• Significant differences in fruit temperatures and temperature gradients were found for the different growing regions.The maximal temperature difference within the fruit during the growth and development of oranges produced is up to ~5 • C for those from Letsitele and Citrusdal, 4 • C for Nelspruit, and that of SRV is up to ~7 • C. Understanding these temperature variations will help the stakeholders determine the optimal postharvest storage and transport conditions for oranges from different regions.For instance, oranges from regions experiencing higher temperature differentials may require more controlled storage environments to prevent accelerated quality deterioration.• Overall, the growth process of Valencia oranges is complex, involving dynamic changes in FW, FD, RW, RT, TSS, TA and respiration rate that follow a non-linear pattern.
This study has provided an initial understanding into how varying growing conditions affect citrus quality at harvest.Our approach of digital twins of oranges from fruit set to harvest can potentially help growers predict the impact of weather conditions and evolving cultural practices on fruit quality.At the same time, this concept has the potential to be seamlessly integrated into soft sensors or web and mobile applications.Doing so could empower various actors within the citrus supply chain, including suppliers and retailers, to decide when and from which growing regions to source their supply.Real-time quality predictions could enhance supply chain efficiency and enable stakeholders to respond promptly to fluctuating market demands, ultimately ensuring the consistent delivery of high-quality citrus to consumers.

Declaration of generative AI and AI-assisted technologies
During the preparation of this work, we utilized Grammarly and ChatGPT to enhance the spelling, grammar, and style of the text.We did not generate any original content using these AI-assisted technologies.Following the use of these tools and services, we diligently reviewed and edited the content as necessary.As a result, we take full responsibility for the content presented in this publication.

CRediT authorship contribution statement
ave Average relative humidity (%)h c Convective heat transfer coefficient (W m − 2 K − 1 ) T air Air temperature ( • C) T air,aveAverage daily air temperature during plant growth process( • C) k air Thermal conductivity of air (W m − 1 K − 1 ) U speed Wind speed (m s − 1 ) G rGrashof number (-) ΔT Temperature difference between fruit surface and air (• C) β Air volume expansion coefficient (K − 1 ) a swReflectance of orange rind to short-wave radiation R sw Short-wave radiation received by the fruit (W m − 2 ) ε Emissivity of infrared radiation for orange R lw Long-wave radiation received by the fruit (W m − 2 ) σ Stefan-Boltzmann radiation constant C p,pulp Specific heat capacity of orange pulp (J kg − 1 K − 1 ) C p,rind Specific heat capacity of orange rind (J kg − 1 K − 1 ) ρ pulp Density of orange pulp (kg m − 3 ) ρ rind Density of orange rind (kg m − 3 ) k pulp Thermal conductivity of the orange pulp (W m − 1 K − 1 ) k rind Thermal conductivity of the orange rind (W m − 1 K −

Fig. 1 .
Fig. 1.Map of several climate zones and the four major citrus-growing regions within South Africa, with indicated annual maximum and minimum temperatures (Alexander).

Fig. 2 .
Fig. 2. Geometry and boundary conditions of orange during the plant growth process (figure not to scale).

Fig. 3 .
Fig. 3.A flowchart of the overall methodology implemented from weather data collection to the fruit quality at harvest (figure not to scale).

Fig. 4 .
Fig. 4. Average experimental and digital twin simulated data of [A] individual fruit weight (FW, g) and [B] individual fruit diameter (FD, mm) at harvest for different regions over the 2019 and 2020 harvest years.The error bar represents the mean ± 1.5 standard error.The percentage differences (Δ) between the experimental and simulated data for FW and FD are shown with the red-colored text.

Fig. 5 .
Fig. 5. Average experimental and digital twin simulated data of [A] rind weight (g), [B] rind thickness (mm), [C] titratable acidity (%), and [D] total soluble solids ( • Brix) at harvest for different regions over 2019 and 2020 harvest years.The error bar represents the mean ± 1.5 standard error.The percentage differences (Δ) between the experimental and simulated data for RW, RT, TA and TSS are shown in the red-colored text.

Fig. 7 .
Fig. 7. Digital twin output for Letsitele, Nelspruit, Citrusdal and SRV growing regions as a function of time showing [A] Average fruit surface temperature, core temperature and air temperature; [B] Temperature difference between the fruit surface and core.The boxplots in [B] represent the mean, median (centre line), 75th upper and 25th lower quartiles (box limits) and 1.5 × the interquartile range (IQR, whiskers) of the temperature gradient for different farming regions.
For 2018/2019 production season, the total duration from fruit set to harvest for oranges from Letsitele is 200 days, 193 days for Nelspruit, 235 days for Citrusdal and 214 days for SRV.These days to harvest were estimated from experimental data.Based on this and daily heat unit data calibration, the 2018/2019 PHU values for Valencia orange produced in Letsitele is 2559, 2557 for Nelspruit, that of Citrusdal is 1623 and SRV is 1918.For the 2019/2020 production season, the duration from fruit set to harvest also varied across different location: 183 days for Letsitele, 179 days for Nelspruit, 235 days for Citrusdal, and 200 days for SRV.The PHU values for Valencia orange produced in this season is 2451 for Letsitele, 2418 for Nelspruit, that of Citrusdal is 2223 and SRV is 1918.

Table 1
Material and thermal properties of orange fruit.