Preprint 3 1 INTRODUCTION Strawberries ( Fragaria x ananassa) are one of the most challenging fresh foods to preserve along the supply chain (Ramirez et al., 2020) . Postharvest losses can be as high as 40% (Porat et al., 2018; Siu, 2014; Trinetta et al., 2020) . Strawberries have a high respiration rate, due to which strawberries can only be stored for 1 –2 days at room temperature (Wang et al., 2019) . Even under optimal conditions, viz. 0 °C at 90 -95% relative humidity, the shelf life is only about 8 -12 days, depending on the cu ltivar (Nunes, 2009; Tano et al., 2009) . Due to the lack of a protective peel, strawberries are also prone to mechanical damage and moisture loss. Additionally, s ince strawberries are non -climacteric fruits, they are harvested almost fully ripe when they have a high water activity (~0.997), high sugar content (~7 -8%), and low pH (~3.84) (Caner et al., 2008; Chirife and Fontan, 1982; Manteau et al., 2003) . This makes strawberries particularly susceptible to spoilage by yeast and mold (Tournas and Katsoudas, 2005) . Botrytis cinerea (or grey mold) is the most common fungal pathogen infecting strawberries (Tournas and Katsoudas, 2005; Trinetta et al., 2020) . The conidia (or spores) of B. cinerea are omnipresent in fields and disperse d by wet splashes or air currents on the fruit surface (Ciliberti et al., 2016) . These spores remain dormant until the ripened berries encounter conditions feasible for spore germination, such as elevated humidity (Abd -Elkader et al., 2021) . Substantia l mycelium growth is often noticeable only in the later stages of the supply chain. A very low grey mold infection level (<0.5%) is acceptable for the retail and import of str awberries, making mold growth the most frequent reason for the rejection of straw berries at retail levels (Trinetta et al., 2020; Zoffoli et al., 1999) . A good packaging strategy, coupled with maintaining an optimal hygrothermal window around the fruit, slow s down the physiological and biological processes in the fruit, thus extend ing their shelf life. S trawberries are packed in a range of ventilated packag ing in several different designs, sizes, and materials. The most commonly employed packages are ventilated clamshells. These are transparent lidded plastic boxes with vents made of (recycled ) polyethylene terephthalate (PET) (Cagnon et al., 2013) . Other forms of packaging include unlidded (open) punnet trays made of paperboard or plastic that may be flow -wrapped with polyethylene films or top -sealed with re -closable polypropylene films (Matar et al., 2018) . These packages usually have vents on the top, bottom, and/or lateral surfaces to improve the cooling efficiency of the package and lower the risk of in -package condensation. However, the precise number, size, shape, and locations of these vents are still often determined arbitrary and designed by trial -and-error approaches (Cagnon et al., 2013) . Moreover, packaging materials such as PET have high water vapor barrier properties, therefore the relative humidity in the package headspace tends to attain high values (>94%) (Bai et al., 2019; Bovi et al., 2019) . High relative humidity is linked with th e risk of condensation of water on fruit surfaces or package walls, which can stimulate mold growth (Linke and Geyer, 2013) . In the last decades, there has been a vast amount of research on package design in postharvest supply chains. S everal studies ha ve employed computational fluid dynamics (CFD) modeling techniques to evaluate the impact of total vent area on the cooling efficiency of different fruit packages (Ferrua and Singh, 2011a, 2009a; Mercier et al., 2017; Nalbandi et al., 2016; Nalbandi and Seiiedlou, 2020; Pathare et al., 2012) . Most of these studies target one unit operation in the entire supply chain , for instance, pre - cooling (Ferrua and Singh, 2011b) , or quantify the cooling behavior but leave aside the impact to any fruit quality attr ibute. Few studies focus on the risk of condensation due to temperature fluctuations at high relative humidity or when a cooled fruit pa ckage is rewarmed to ambient temperatures (Gottschalk et al., 2007; Linke et al., 2021) . Additionally, w hile many stud ies have attempted the experimental characterization of condensation (Lufu et al., 2021) , the authors have reported the challenges in quantif ying the spatial distribution of the risk of condensation or residence time of the condensate for several fruit in a package (Bovi et al., 2019; Linke et al., 2021) . So far, very few studies have analyzed the spatio -temporal occurrences of condensation, its comparative risk amongst different packages, and the resulting impact on the microbiological quality of fruit. Moreover, every supply chain comprises a different combination of unit operations, and every shipment encounters unique logistic and environmental boundary conditions (Jalali et al., 20 20). In this sense, e very supply chain has different packaging requirement s for providing the optimal hygrothermal climate and ventilation to the product (Cagnon et al., 2013; Defraeye et al., 2015) . As such, every supply chain has its own optimal pack aging design, depending, among others, on the temperature profile over time throughout all unit operations. Current approaches for packaging design do not take into consideration the whole length of unit operations along the supply c hain, or the full bread th of relevant physical processes influencing fruit quality and marketability. In this study, we develop ed a physics -driven digital twin of strawberries in 21 different commercially available ventilated packag es, including ventilated clamshells, cardboard punnet trays, top -seal trays, and trays flow -wrapped in film packaging. A digital twin employs mechanistic simulations to mimic the hygrothermal, physiological, and microbiological evolution of the fruit in -silico (Defraeye et al., 2021, 2019; Verboven et al., 2020) . It is linked with sensor data of the hygrothermal environmental conditions in an actual supply chain. We link ed these twins with air temperature and humidity sensor data obtained from over 30 actual shipments of inland -transported straw berries from the farm to the retail store for a specific supply chain . With this digital twin, we predict ed actionable metrics of fruit quality and shelf life for fruit in each simulated package . These metrics include total mass loss, risk of putative mold infection due to B. cinerea , retention time of condensate, and the remaining shelf life based on respiration, transpiration, and mold growth . For every type of package, we pinpoint ed the critical unit operations in the supply chain where losses are aggravated. We explore d how the aerodynamic and thermodynamic metrics of the package , such as total vent area, degree of filling, pressure drop across the package, and seven -eighths cooling time, impact fruit quality attributes . We al so analyzed the spatial heterogeneity within these packages and highlighted the most vulnerable locations within selected packages. Ultimately, we compared the physical, biochemical, and microbiological shelf life of berries in the 21 packages to identify the optimal ventilated package for the supply chain from the farm to the retail store. To the best of our knowledge, such a comprehensive and holistic analysis of packaging strategies in postharvest supply chains for a wide range of package types is unique and is therefore the main merit of this research. Preprint 4 2 MATERIALS AND METHODS Simulated supply chain The hygrothermal profile of the strawberry postharvest supply chain from the farm to the retail store was reconstruct ed based on actual measure d data. We used a stochastic approach to compute air temperature and relative humidity based on the mean and variability in this segment of the measured values, using the 'rand ' function for a specific range, implemented in Microsoft Excel 2016. From pack -house to re tail stores, our values were derived directly from measured data for inland -transported berries from over 30 shipments from Spain to Switzerland , reported in (Shoji et al., 2022) . For other segments of the supply chain, we used data from literature (Coll ins and Perkins -Veazie, 1993; Ferrua and Singh, 2009a; Kuchi and Sharavani, 2019; Shoji et al., 2022) . The complete simulated supply chain is depicted in Figure 1. The typical air speed for every unit operation, based on the typical flow rates, are also indicated in this figure (Han et al., 2017; Li et al., 2021; Opara and Zou, 2007; Thompson et al., 2008a; Wu et al., 2019) . Figure 1: The simulated profile for the temperature, relative humidity , and superficial air speed of the air delivered to the packaged fruit in a typical supply chain of strawberries in -land transported from Spain to Switzerland (Collins and Perkins -Veazie, 1993; Ferrua and Singh, 2009a; Kuchi and Sharavani, 2019; Shoji et al., 2022) . These data were used as input for the physics -based digital twins . Image credits: This figure has been designed us ing resources from Flaticon.com . In a typical supply chain, straw berries are collected at harvest (at 20 °C) and field -packed in individual packages stacked in standard reusable p ackaging crate s (Ferrua and Singh, 2009a; IFCO, 2022) . These packed berries are then stored in a cold chamber (6-8 °C) for 1 -6 hours. Following this, they are palletized and pre -cooled in a forced -air cooling tunnel (2 -4 °C for 3 h) (Kuchi and Sharavani, 2019) . The berries are then loaded into refrigerated trailers and transported from the pre -cooling facility in Spain to the distribution center in Switzerland (~ 2-day journey). It is recommended that temperatures should be as low as possible (0-4 °C), but in a previous study, our analysis showed that the air temperature is in the range of 5 -6 °C during transport (Shoji et al., Preprint 5 2022) . One of the reasons for transporting at a higher temperature range is to reduce the risk of top freezing as the cold air enter ing the trailer is often not exactly at the set air temperature. Large retailers often hav e their own distribution centers (or warehouses) where the fruit are delivered and then redistributed to retail stores in the region. At the distribution center, the fruit spend anywhere between 12 and 18 hours (Shoji et al., 2022) . Air temperatures are relatively stable (6±1 °C) (Shoji et al., 2022) . The berries are then redistributed from the regional warehouse/ distribution center to the retail store, lasting about 18 hours. In order to lower the risk of condensation, the air temperature is often gradually increased ( or ramped up ) in this segment . From our measurements, we find that the air temperature increases linearly increase over time, for the ramp -up from 6 °C to 10 °C (Shoji et al., 2022) . Our measurements revealed that the ambient relative humidity during refri gerated transport and storage i s 80±3% from the pack house up to the retail store. On arriving at a retail store, the fruit are stored overnight (~8 h) in the cold storage room where air temperatures are usually 9-11 °C as the fruit are stored along with other fresh produce. Following this, we reconstructed a retail period of 8 h, when the fruit are displayed on retail store shelves at standard ambient temperatures ( 23 °C, 40% relative humidity ) (The British Standards Institution, 2017) . This is because strawberries are most often displayed o n retail store shelves with no refrigeration (Collins and Perkins -Veazie, 1993) . In total, a supply chain of 4.5 days is simulated with this digital twin. Since the strawberries are transferred from one unit operation to another, there is often an abrupt change in the environmental conditions. These rapid changes during handling are a result of different unit operations concatenating to compose the supply chain. Continuum multiphysics model for a package of strawberries 2.2.1 Computational system configuration The computational domain comprises a ventilated package filled with strawberries in an air channel without bypass, as shown in Figure 2. The fruit a re modeled explicitly since the package -to-product length ratio i s in the range of 2 -3 (so <10) , implying the fruit in a single package cannot be modeled as porous media (Dehghannya et al., 2010; Ferrua and Singh, 2009a; Nalbandi et al., 2016) . Each package is filled with equi-sized strawberries (equatorial fruit diameter = 3 0 mm, mass = 16.30 g, surface are a = 33.06 x 10-4 m2, volume = 16.9 x 10-6 m3, equivalent spherical diameter = 31.84 mm). The size is in accordance with the Class I commercial grade of berries imported to Switzerland (Schweizer Obstverband (Fruit -Union Suisse), 2005) . Every individual strawberry is modeled as a three -dimensional geometry created by axi -symmetricall y rotating a contour of a strawberry fruit (cv. Fortuna ). The number of strawberries in every package varie s depending on the intended weight of the stock -keeping unit , i.e., 250 g (16 berries), 300 g (20 berries), 400 g (25 berries), or 500 g (31 berries) . The packages are filled with fruit in 4 -7% excess of the intended weight to account for the reduction in saleable weight due to moisture loss along the su pply chain. The berries are stacked in a regular arrangement such that their longest edge rest s on the base of the package with a net -zero torque about the center of gravity. Horizontal airflow in a no bypass channel through a single package is modeled for the different packag ing configurations . The reason is that typically horizontal flow is encountered in a pre-cooling tunnel and also to some extent in a r efrigerated trailer (Ambaw et al., 2021; Ferrua and Singh, 2009b; Getahun et al., 2017; Wu et al., 2019) . Most packages also have vent -hole configurations that are designed for promoting horizontal airflow ventilation rather than vertical airflow ventilation. The sections upstream and downstream of the package are taken sufficiently large (~30 cm) to av oid an influence of the inlet and outlet boundary conditions on the airflow and heat transfer inside the package . All the selected packages are ventilated, i.e., these packages had vents or macro -perforations ( diameter 8- 12 mm) to enable the transfer of O 2 and CO 2 at rates much higher than the fruit's respiration. Therefore, modified atmosphere packaging or the O2 and CO 2 composition inside the package are not modeled (Paul and Clarke, 2002; Sousa-Gallagher and Mahajan, 2013) . Preprint 6 Figure 2: A schematic overview of the continuum multiphysics model depicting the computational system configuration, governing equations, and implemented boundary conditions for the air, package , and fruit domains for a representative ventil ated package. The equations corresponding to the air, package, and fruit domains have been described in the text and the nomenclature is elaborated at the end of this study. Image credits: This figure has been designed using resources from Flaticon.com and Adobe Stock (AdobeStock_49352459, AdobeStock_214803643). Preprint 7 2.2.2 Airflow A uniform upstream velocity ( Uinlet, m∙s-1) is defined at the inlet for airflow, based on the superficial air speed of the unit operation (Figure 1) . This value i s set to 0.01 m∙s-1 for retail and refrigerated storage, 0.1 m∙s-1 for refrigerated transport, and 1.0 m∙s-1 for pre-cooling (Han et al., 2017; Li et al., 2021; Opara and Zou, 2007; Thompson et al., 2008a; Wu et al., 2019) . Note that these ranges are indicative for the superficial air speeds in these unit operations. For our conceptual study, such indicative ranges of air speed are sufficient to draw relevant conclusions with respect to relative differences between packaging performances. The convective heat transfer coefficients and Biot numbers corresponding to these ranges of superficial air speed are presented in the Supplem entary mater ials, Section A.1. A turbulence intensity of 0.05% is assumed for the flow entering the computational domain (Berry et al., 2016) . An open boundary cond ition is defined for the outlet in the downstream section. As the contribution of moisture added due to evaporation from the fruit surface is small, a no -slip boundary condition i s applied to all solid surfaces, including all surfaces of the package and fr uit. We model incompressible flow, thus the continuity equation takes the form presented in Eq. 1. 0∇⋅ =u (1) where u is the velocity field vector and ∇ is the vector differential operator . The Navier Stokes equations take the form presented in Eq. 2. T( ) [( )( ( ) )]aT P ρ µµ ⋅∇ =−∇ +∇⋅ + ∇ + ∇uu u u (2) where ρa is the density of air (kg∙m-3), P is the atmospheric pressure (Pa), and (μ+ μT) represents the effective viscosity (Pa∙s). Buoyancy was not considered, as analysis of the Richardson number showed that buoyancy -driven flow is only significant at superficial air speeds below 0.05 m∙s-1. Note also that modeling buoyancy would significantly extend the computational cost of the simulations as, for changing thermal environments, this would require recalculating the flow field every time step. This is a re ason why it is rarely accounted for in CFD simulation studies in postharvest science. For a typical package, the flow regime is mo stly laminar around the package, but turbulent at the vents (Ferrua and Singh, 2009a) . Therefore, Eq. 2 i s solved in the turbulent regime using the Reynolds -averaged Navier– Stokes (RANS) equations, which account for turbulence effects on the mean flow fi eld through time averaging. To this end, we use a k-ε turbulence model with wall functions. This approach describes turbulence by an increase in effective fluid viscosity ( μ+ μT) (Norton et al., 2012) . The turbulent component of viscosity ( μT, Pa∙s) is c omputed using the turbulent kinetic energy ( kT), and the rate of dissipation of the turbulent kinetic energy (ε), expressed in Eq. 3. 2 aT TCkµρµε= (3) Here, Cμ is the turbulence model constant. The turbulence variables a re solved in the whole air domain. The no -slip boundary condition is implemented at th e walls. We used wall functions to model the near -wall regions and built a corresponding mesh with a desired y+ value of 1. As the thickness of the package is significantly smaller than the length scales of the adjacent domains, the walls of the package are modeled as interior boundaries, so a two-dimensional plane. This allows discontinuity in velocity, pressure, and turbulence variables across the package walls. We assume that the flow field is independent of the air temperature and moisture c ontent, and therefore calculate the flo w field in advance and then use it as input for the heat and moisture transport equations. 2.2.3 Heat transport Heat transfer in the air domain The heat transfer equation for the air domain is presented in Eq. 4. This equation describes how the air temperature ( Ta, K) changes with time, based on conduction and convection terms. ,, ( ) (k )a a pa a pa a a aTc cT Ttρρ∂+ ⋅∇ =∇⋅ ∇∂u (4) where cp,a is the specific heat capacity of air (c p,a = 1006 J∙kg-1∙K-1), and ka is the thermal conductivity of air (ka = 0.024 W ∙m-1∙K-1) (ASHRAE, 2010) . Heat transfer across the package walls The walls of the package are modelled as thin layers with a thermal resistance. Since the thermal conductivity of the packaging material is comparable to that of the surrounding air, a 'thermally thin approximation' is implemented. With this approximati on, we assume that the heat flux across the package is p roportional to the temperature difference between both sides of the package wall (Eq. 5). , ,, , ( )( ) ( )p in a a in a out a i n out a a out pkkT T T kTx− −∇ = − = −∇nn (5) Preprint 8 Here, the terms n in(-ka∇Ta,in) and n out(-ka∇Ta,out) correspond to the flux on either side of the package wall. The terms kp (W∙m-1∙K-1) and xp (m) refer to the thermal conductivity and thickness of the packaging material. The corresponding thermo-physical properties for the different packaging materials employed in this study are detailed in Table 1 . Table 1: Thermo -physical properties to quantify heat and moisture transport across the packaging materials corresponding to the different packages investigated in this study . Material Thickness (xp, mm) Thermal conductivity (kp, W∙m-1∙K-1) Specific heat (cp,p, J∙kg-1∙K-1) Density (ρp, kg∙m-3) Vapor permeability at 25 ºC , 90% RH (δp, kg∙m-1∙s-1∙Pa-1) Vapor resistance factor ( µp ,-)