We build a validated physics -based model of a refrigerated container filled with fruit in ventilated packaging. This model of a virtual container i s the basis for simulations in an accompanying paper on citrus fruit shipped overseas from South Africa to Europe. The model is used to understand better how the cargo cools and when and where food quality is lost in these supply chains. We build a computational fluid dynamics model with a two -phase porous media approach that simulates the airflow in the container and the cooling process of e very fruit. This container can be considered aerodynamically to be a slot -ventilated enclosure. W e also model the fruit's thermally -driven quality loss. Using a two -phase porous media approach for the ventilated cargo and modeling temperature -driven fruit quality evolution are two steps forward compared to most existing physics -based refrigerated contai ner models. We validate the porous media model implementation. We define and apply actionable metrics for every fruit inside the cargo, such as remaining shelf life upon arrival and seven -eighths cooling time. Preprint 3 2 MATERIALS AND METHODS Container, packaging, palletization , and cold chain specifications 2.1.1 Packaging BOX A telescopic corrugated fiberb oard carton type ‘ A15C Supervent’ was used (0.4 x 0.3 x 0.27 m, Figure 1a). The carton contains ‘Navel’ oranges ( Citrus senensis L. Osb). The citrus fruit is spherical and ha s a diameter of 80 mm. The carton box has half -circular vent holes located at the top and bottom of each side. Note that newer designs of this carton type have been developed in the past years. This carton was already evaluated for forced airflow cooling (Thijs Defraeye et al., 2014, 2013a) . The carton's total open area (TO A) on the long vertical side, short vertical side, bottom , and top is 3.5%, 3.1%, 9.1%, and 12.9%, respectively . All fruit are placed in the carton according to a staggered pattern. The Supervent cartons hold 64 orange fruit and ha ve a mass of 16.5 kg each . The amount of fruit in the box is 53% of the total box volume . The corresponding porosity is 47% (εfr). This porosity includes the volume of air and the packag ing material as 'air space' and only the fruit as the solid material . The porosity, when including the package and the fruit as a solid material (ε), is lower and equals 38%. CARTON The Supervent cartons are stacked to assemble a high -cube pallet (1.2 m x 1.0 m x 2.16 m, Figure 1b). On a pallet, 8 layers of cartons are stacked, where each layer holds 10 cartons (Figure 1b). In total, 3 rows are present (Row1, Row2 , and Row3) along the horizontal flow path for precooling and cold storage, and 8 layers (L1- L8) are present along the vertical flow path for refrige rated container transport. Each layer contains 10 cartons (C01- C10). Row1 and Row2 have 3 cartons for each layer , whereas Row3 has 4 cartons ( Figure 1b). The pallet cont ains 80 cartons and holds 5, 120 fruit. The ‘Supervent’ package was designed to form vertical and horizontal ventilation pathways via the vent holes during precooling and transport after palletization. Though this carton design also enables vertical pathway s, these are less evident in practice due to the staggered stacking pattern of consecutive layers ( Figure 1c). Practitioners deem this necessary to ensure the stability of the pallet during handling and transport , but it obstructs several vent holes. In addition, the slats of the wooden pallet base fully or partially cover some of the vent holes of the bottom carton box (Figure 1c). PALLET BASE The wooden pallets (Figure 2) on which the cartons are stacked measure 1.21 m x 1.01 m x 0.16 m. This pallet is common in the South African citrus industry. The bottom part of the pallet base has an open area of 60% of the total surface area of the pallet (1.01 x 1.21 m). The top part of the pallet base, which is in contact with the fruit cart ons, has an open area of only 24% of the pallet base. Other pallets can have an even smaller open area (15% for (EPAL, 2021) ). The wooden blocks between the top and bottom of the pallet base have an open area (or also open volume) of about 90% of the pallet area. The pallet base provides a significant blockage to the airflow , especially the top part of the pallet base . Compared to the pallet of fruit in ventilated packaging , the resulting pressure loss over the pallet base is , however, rather small. Preprint 4 Figure 1. (a) Geometry and dimensions of the Supervent carton used in the present study and by (Thijs Defraeye et al., 2014, 2013a) . (b) Illustration of idealized stacking of the cartons for an 8 -layer, high -cube citrus pallet (10 Supervent cartons per layer, each packed with 16 .5 kg of fruit), together with typical airflow for forced- air precooling (horizontal flow) and in a refrigerated container (vertical airflow). (c) Illustration of actual stacking used in commercial operations with a pallet base. Note that a slightly different design of the Supervent carton is shown here, compared to the one shown in (a) , with fewer ventilation holes. Dimensions are in mm. Figure 2. Pallet base top view and bottom view with dimensions (in mm). 2.1.2 Container DIMENSIONS A state -of-the-art reefer container was used in the present study (Figure 3). Such refrigerated containers have a lifetime of about 12-18 years (Lawton, 2016) . The container details are specified in (Getahun et al., 2018) . This 40- foot HC (High -Cube) container had inner dimensions of 11.56 x 2.29 x 2.58 m and external dimensions of 12.19 x 2.43 x 2.90 m. These inner dimensions do not include the T -bar floor , which is installed afterward. This floor is composed a b Flow (refrigerated transport) c 2160C01 C02 C03C04 C05 C06C08 C09 C10 L1L2L3L4L5L6L7L8C07 Row1Row2Row3 75 75 100 100 100 100 100 75 7575 50 80 80 50 75 1210 101095 9595 95 95a b Top view Bottom view Preprint 5 of extruded T -bars that go fro m the side of the refrigeration unit to the container door side (3 5 full T-bars of 63.5 mm height and 63.5 mm width ). The container was filled with 20 pallets of fruit (1 x 1.2 m). The loading pattern is shown in Figure 3, but other loading patterns are also possible. In this study, the container contains 1,600 cartons and 102,400 fruit. Due to a n inherent mismatch between the dimensions of the pallets with those of the container, not the entire floor area can be used. Even when filled with standard pallets ( pallet base = 1.01 x 1.21 m2), 7.7% of the container’s floor area remains uncovered (2.03 m2). This area is located mostly at the door end of the container. However, on the lateral sides of the container, a gap of 70 mm is present for such pallets. This is usually required to provide some clearance when loading the pallets, as they are not perfectly stac ked to a 1.01 x 1.21 m size. Nevertheless, the open area on the lateral sides remains. Also, note that the pallets of fruit, so the packages on the pallet base, only occupy 1.0 x 1.2 m when stacked ideally . The space on the lateral sides of the container for the packages is 90 mm since the container is 2,290 mm wide and two pallets only occupy 2,200 mm. In practice, pallets tend to bulge and thus become a bit wider after handling. As such, the actual gap s ize will likely be less than 90 mm and i rregularly distributed. The gaps will mostly be on the center line between the two pallet rows because it is more convenient for the forklift driver during loading . The air enters the container via the bottom through an inlet slot of a height of 63.5 mm (Getahun et al., 2018) . This inlet slot has the same height as the T- bar floor since the air is directed from the refrigeration unit via a baffle plate into the T -bar floor. Figure 3. Side and top cross -sectional views of a high -cube refrigerated container with an indication of the stacking pattern of the pallets and the box layers. The inner dimensions of the container are shown. AIRFLOW GENERATION & AIRFLOW RATE Airflow in the container is generated by the two axial evaporator fans in the refrigeration unit (Wild et al., 2005) . These fans circulate air through the container along the bottom via a grated T -bar floor. Air is drawn vertically upwards through the cargo and returns to the refrigeration unit via the ceiling ( Figure 3). For fresh fruit and vegetables, the air is forced through the cargo , packed in ventilated packaging, to dissipate possible respiratory heat TOP VIEW 1.2 m 1 m 1 m 1.2 m 11.56 mT-bar floor 2.29 mSIDE VIEWCooling unitL 1L 2L 3L 4L 5L 6L 7L 8 Optional void plug (to block airflow)Return air Supply air0.2 m 0.0635 m2.32 m2.16 m 2.58 m Inlet slot (0.0635 m)Fresh air intake Baffle plateEvaporator HeaterFans Fresh air intake Preprint 6 and metabolic gasses ( e.g., carbon dioxide and ethylene). For frozen cargo, however , the air is directed around the cargo. Typically two fans are installed. These fans operate continuously but are briefly switched off during a defrosting cycle. Nevertheles s, the container's evaporator fans typically have two operating regimes for which a different RPM can be set. For the Carrier Transicold container in this study, the two evaporator axial fans can run at two speeds: high (max) speed, so an RPM of 3460/2850 driven at 62/50 Hz, and low (half) speed, so an RPM of 1760/1425 driven at 62/50 Hz motor, and no speed, when the fans are switched off. These fans are driven by a 3 -phase (380/460 V) 62/50 Hz motor. There are two operating modes of these fans for a container: chilled mode (> ~ - 5ºC) and frozen mode (< ~ - 5ºC (Hamburd -Sud, 2021) ). Typically, in chilled mode, the evaporator fans run continuously at high speed, and the compressor operates continuously. In frozen mode, the evaporator fans run at a low speed , and the refrigeratio n system works in on/off control (Lawton, 2016) . Both container fans typically deliver a total flow rate of 1800- 5500 m3 h-1 together, depending on whether the container is loaded and on the selected speed regime of the evaporator fans. The corresponding air exchange rate lies between about 25-80 air changes per hour of the internal volume of the container at low and high speeds , respectively. The airflow conditions in a loaded container that were used in this study are detailed in Table 1 and were measured in a loaded co ntainer at different operating regimes at 50 Hz and 62 Hz (CRT, 2013) . These conditions result in an average airspeed of 3.4- 7.9 m s-1 at the inlet slot (0.0635 x 2.29 m2) and 0.019 -0.044 m s-1 over the container floor area ( Table 1). Nowadays, however, the low -speed regime is used to save energy and reduce fan heat production and is used in frozen or chilled mode with non/low -respiring fruit and vegetables . With more advanced cooling protocols, th e controller often chooses the low -speed regime based on the measured temperature difference between supply and return air (Lukasse et al., 2012, 2011) . This reduces the energy consumption of the fans in chilled mode by optimizing evaporator fan speed based on the heat load and also avoids inefficient part -load compressor operation while still maintaining quality (Lukasse et al., 2012, 2011) . Table 1. Specifications of airflow in a loaded high -cube refrigerated container (for 2 evaporator fans). Airflow regime Airflow rate [m3 h-1] Fan pressure rise [Pa] Fan power [W] Average airspeed over the container floor area x 103 [m s-1] Flow rate per kg of fruit [L s-1kg-1] Airspeed at inlet slot [m s-1] Low speed @ 62Hz (50Hz) 2175 (1770) 55 (37) 270 (200) 22.8 (18.6) 0.023 (0.019) 4.2 (3.4) High speed @ 62Hz (50Hz) 4150 (3490) 198 (146) 1710 (940) 43.5 (36.6) 0.044 (0.037) 7.9 (6.7) EVAPORATOR FANS Fans are typically characterized by their fan performance curve. This curve displays several airflow rates and the static pressure rise that the fan can deliver. This curve is determined by measuring the air flow rate through a partially obstructed channel where the fan is placed . By incrementally increasing the aerodynamic resistance of the channel, so increasing the obstruction, different airflow rates are obtained at different static pressures. Typical fan performance curves for the two fans installed in a refrigerated container are depicted in Figure 4 (Berry et al., 2023; CRT, 2013; Tiamiyu, 2020) . Note that these fan performance curves are independent o f the airflow resistance of the system in which the fan is placed. Preprint 7 Figure 4. Fan performance curve (for 2 fans together ), so airflow rate through these fans as a function of the fans' static pressure rise from two sources (Berry et al., 2023; CRT, 2013) . The system curve with individual working points of the fans for a loaded container is also shown (CRT, 2013) . AIRFLOW RESISTANCE THAT THE FAN NEEDS TO OVERCOME The actual airflow rate and the pressure increase the fans provide are, to some extent , dependent on the pressure resistance of the cargo and the rest of the pressure losses in the closed circuit of the container. The aerodynamic (airflow) resistance of t he container circuit i s expressed as the relation between the total pressure drop over the container circuit from the fan outlet to the fan inlet ( ∆Pt [Pa]) and the volumetric air flow rate through the fans (G a [m3 s-1]): 2 ,, t t i tot a fan iP PG ξ ∆= ∆ =∑ (1) Here ξ is the pressure loss coefficient of the entire circuit ([Pa s2 m-6] or [kg m-7]). This pressure drop only accounts for inertial effects (Forchheimer term). The viscous effects (Darcy term, proportional to G a) only become important at very low flow speeds (Verboven et al., 2006) and are negligible in most parts o f the container due to the high speeds. It should be noted that the total pressure (P t) is the sum of the static pressure (P s) and the dynamic pressure (P d = 0.5 ρ U2, with ρ the air density [kg m-3] and U the airspeed [m s-1]). The airspeed upstream and downstream of the fan is the same, so P d remains constant by which ∆P t = ∆Ps in Eq. (1) (ASHRAE, 2010a) . However, the dynamic pressure will change throughout the airflow circuit of the refrigerated container as spee ds in the cargo hold are much lower. All system components that cause pressure losses in the closed circuit of the container determine the container's system curve. The intersection of the system curve with the fan performance curve determines the corresponding working point of the system (∆P t,wp(Ga,wp)). The system curve for a loaded container with individual working points at different speeds is shown in Figure 4. The contribution of palletized fruit to the system curve, so to the working point ∆Pt,wp(Ga,wp) and the total power consumed by the fan (P w,fan), can be substantial. This is often the case fo r forced airflow cooling (De Castro et al., 2005) , where the p roduct -packaging ensemble is one of the main resistances in the system. For refrigerated containers , though, the pressure loss coe fficient of the pallets is high , but the airspeed through them is much lower than in other components in the refrigeration unit . Thereby t he contribution of the palletized fruit to the resulting pressu re loss over the cargo hold is rather small (∼ 100-101 Pa, from our study and other work reported in (CRT, 2013) ) compared to the pressure rise provided by th e fan (∼ 50-200 Pa, Table 1, Figure 4). The cargo, and thus the fruit geometry, packaging density , and packaging are relevant but not -dominant determinants in the aerodynamic resistance the fans need to overcome. The main aerodynamic resistance lies in the refrigeration unit. However, still, w ith an empty container, the fans can run at a different speed then for a filled container, d ue to the contribution of the cargo to the pressure loss over the entire circuit . FAN POWER CONSUMPTION The total power consumed by the fan (P w,fan) can be calculated as: Pw,fan = ∆P t,wp Ga,wp/(η fan ηmotor), where wp indicates the working point of the fan . Pw,fan also depends on the fan and motor efficiencies (η fan and ηmotor (Baird et al., 1988) ) 0100200300400500600700800 0 1000 2000 3000 4000 5000 6000 7000Pressure [Pa] Airflow rate [m3h-1]CRT (2013) - Low speed (50 Hz) CRT (2013) - Low speed (62 Hz) CRT (2013) - High speed (50 Hz) CRT (2013) - High speed (62 Hz) Working point fan Berry et al. (2023) - Cont. A Berry et al. (2023) - Cont. C Preprint 8 and is thereby proportional to G a,wp3. The fan energy consu mption (E fan [J]) can easily be calculated as Efan = Pw,fan top, with top the operational time of the fan. The power required to push air through the ventilated cargo (Pw,cargo [W]) can be determined (Thijs D efraeye et al., 2014; Ferrua and Singh, 2011) : 3 , arg , arg arg w co t co a co aP PG G ξ = ∆= (2) Pw only contains here the contribution of the packaging . Therefore the aerodynamic resistance of other system components has to be accounted for as well to determine the system curve (similar to Eq. (9) (ASHRAE, 2012, 2010a) ). As mentioned, the pressure losses over the other components will likely dominate the pressure the fan s need to deliver. CONTAINER THERMAL CONTROL The airflow, at high or low speeds, is driven by the fans through the evaporator. This evaporator needs to cool away: (1) heat produced b y the evaporator fans and other equipment; (2) the heat loss due to air infiltration by fresh air intake; (3) the product heat load, namely the heat still stored in the product that needs to be cooled away in case the cargo was not entirely precooled to th e supply air temperature; (4) the respiration heat from the fruit or vegetables; (5) the heat transmission losses through the container walls. Note that also latent heat is extracted from the air when moisture evaporates from the products. After the evapor ator, a heater is in place. The heater has three functions, in order of importance: 1. Periodic defrosting of the evaporator coil at set points below approx imately 5 °C. 2. Avoid too cold supply air temperature during chilled transport. 3. Dehumidify the supply air in chilled mode shipments where a humidity set -point is used . Here, the heater is turned on and invokes a colder evaporator temperature, hence extra dehumidification. After passing the evaporator, the cold air is delivered into the cargo hold at the supply /delivery air temperature (T sup, TDAT). The entry is via an inlet slot at the bottom of the container. The air is driven from bottom to top through the cargo in a chilled model for respiring fresh fruit and vegetables, which are packed in ventilated packaging. For frozen cargo, the air is driven around the cargo. This cold air then absorbs the heat from the cargo and the heat entering the cargo from outside. This air is extracted via an outlet slot at the top of the container, back to the reefer unit , and the evaporator, where the return air temperature is measured (T ret). The difference in s upply and return air temperature depends on the container's different heat loads (depicted in Figure 5) and the rate at which the heat is removed from the fruit. T he heat exchange rate depends on the airflow rate (so on the pressure loss over the pallets) and the air access to the ventilated cargo. This access depends on the ventilated packaging design , container stowing strategy, and airflow bypass. Figure 5. Heat loads in a container (side view) . Red indicates a heat gain that needs to be cooled away, and blue indicates a heat loss or extraction from the container. One pallet is indicated in dark grey to illustrate which heat is ext racted and gained from a pallet. The nomination f() is used to indicate to which parameters each heat load is a function of. HUMIDITY CONTROL The air extracted from the container is sent to the evaporator unit. The reduction in temperature in this region increases the air's relative humidity, which typically can induce condensation or frost formation on the evaporator fins. This process remo ves water vapor from the air and lower s the supply air temperature even to subzero Infiltration heat gain Evaporator HeaterFans Exfiltration heat lossTransmission heat gain f(species , fruit temperature )f(species , surface temperature & RH) f(initial fruit temperature , mass fruit , specific heat capacity )f(open vent area , outside air temperature ) f(outside temperature, airspeed, inside air temperature, thermal conductivity container)f(airspeed , porosity , packaging , stowing pallets ) Preprint 9 temperatures. Cooling thereby dehumidifies. In a small number of shipments (typically garlic, onions, and flower bulbs) , a humidity set -point of around 70% is specified to evoke extra dehumidification. The heaters are turned on when the humidity set -point is specified , and humidity is above the set -point . The supply air temperature controller responds by adjusting the evaporator temperature . To further magnify the effect , the evap orator fan speed is often reduced. When the air enters the cargo hold again, water vapor in the air can incr ease due to the evaporation of fruit transpiration and previously -condensed liquid water. Also , the temperature of the air increases, by which the air can hold more water vapor as well. The higher the return air tempe rature, the more moisture can be stored in the air , so the more the evaporator can remove by cooling. T he refrigeration unit typically leads to net dehumidification of the air and cargo. With fresh air intake, more moisture is introduced in to the container . But almost all inlet moisture i s immediately captured on the cold evaporator before it enters the cargo space, as the fresh air first passes the cold evaporator before entering the cargo space. The humidity in refrigerated containers typically equilibrates to 85- 95%. Figure 6. Psychrometric chart with an indication of the operation of a refrigerated container (adapted from (Ogawa, 2021) ). FRESH AIR INTAKE Air intake is often set to remove metabolic gasses (CO 2, C2H4). Still, it also leads to the intake of additional humidity (Lawton, 2016) regulated by two openings in the container upstream and downstream of the evaporator fans (Figure 3). The pressure difference over the fans draws in the fresh air and thereby only works efficiently when these fans are switched on. This fresh air is typically warmer and contains more water vapor. This air intake introduces additional moisture, which increases the need for defrosting . Furthermore, the warmer air can complicate temperature contro l. Typical air intake rates can be set from zero to about hundred m3 h-1 and are chosen depending on the commodity (Hamburd -Sud, 2021; Lawton, 2016) . In the case of citrus, air intake rates are typically set to 15 m3 h-1. REFRIGERATION UNIT & POWER CONSUMPTION Evaporator HeaterFans1 2 3 Heat & moisture absorption inside container1 2 3Cooling in evaporator Heating byheater Preprint 10 The refrigeration unit in the container has a refrigeration capacity [kW] that depends on the set temperature inside the container and the outside ambient temperature. The net refrigeration capacity is given in Figure 7a. The resulting power consumption of the container is also shown . The higher the internal temperature in the containe r is set , the higher the required electrical power when running at its maximal refrigeration capacity, and the higher the refrigeration capacity (Wild et al., 2005) . Typically, a refrigerated container's required power is measured to be 2.7- 3.6 kW per TEU (Twenty Foot Equivalent Unit, (Fitzgerald et al., 2011; Wild et al., 2005) ), so about 5.4- 7.2 kW for a 40-foot container. This is representative of the transport of products in chilled mode . Tra ditionally, transport in frozen mode used to require less power due to inefficient operation in the chilled range. The ratio of this net refrigeration capacity, so the heat removed from the container (Q cool) to the power required for the container to do so (Pcont), is the coefficient of performance (COP (James, 2019) ). The COP equals: cool contQCOPP= (3) The COP is , in this case , below 1.4 (Figure 7b). Of this power consumption (P cont), the fans require about 1.7 kW at high speed and produce additional heat that needs to be cooled away (T. Defraeye et al., 2015b) . The required power is supplied by clip-on generator sets ('gen -sets') during transport by road or rail and by the ship during marine shipping. These gen -sets deliver about 15 kW (Lawton, 2016) . Note that the refrigeration unit needs to dissipate th is heat that is cooled away from the container and the additional heat produced by the fans and other equipment into the environment. This amount depends on the type of respiring cargo loaded, the fan operating state, the operating state (chilled/frozen), and the corresponding internal and external temperatures . Typically between 7- 15 kW (Wild et al., 2005) (Qcool) needs to be cooled away . This heat is then released into the environment surrounding the container . Such heat dissipation must be accounted for on ships, especially as many containers are packed i n a tight space. Note that the numbers specified here on the power and energy use are subject to change as new containers with lower energy use are being developed and deployed. Figure 7. Net refrigeration capacity (Q cool), power consumption (P cont) (a) , and resulting COP (b) of a refrigerated container for two types of compressors (Comp.) and refrigerant s (Refr.) (Wild et al., 2005) . DEFROSTING CYCLE The evaporator cooler unit needs to be defrosted periodically. The reason is that frost builds up if surface temperatures in the evaporator are below zero. The frost on the evaporator reduces heat transfer performance from the evaporator to the air and induces additional airflow blockage . This higher airflow resistance , thereby , can change the operating point of the fans since the airspeed s through the evaporator are large. Frost formation often occurs if the cargo’s set (delivery) air temperature is below 5 -10 ºC (Wild et al., 2005) . Such defrosting implies that the fans are stopped, and the evaporator unit and drip tray are heated electrically . Also , hot-gas defrosting of the evaporator is used by some manufacturers . Such defrosting cycles are noticed in the return air temperature data, where higher temperatures are logged due to the rising hot air from the heating. The control of the defrosting cycle is typically done at regular time intervals (e.g. , 6h, 8h, 12h, (Wild et al., 2005) ). Automatic defrosting control is also done at variable time intervals, 0246810121416 -30 -20 -10 0 10 20Cooling capacity & Power (kW) Internal container temperature ( °C)Q (type 1) P (type 1) Q (type 2) P (type 2) 00.20.40.60.811.21.4 -30 -20 -10 0 10 20Coefficient of Performance (kW) Internal container temperature ( °C)Refr. + Comp. Type 1 Refr. + Comp. Type 2a b Preprint 11 depending upon need (Hamburd -Sud, 2021) , which optimizes defrosting activity. C ontrollers can pick up the need for defrosting if the supply air temperature is not reached over an extended period or if the temperature difference between supply and return air is too high. The duration of the defrosting period, however , is variable . A defrost cycle is typically terminated when a defro st termination sensor, mounted just above the evaporator, reads a value above a certain treshold (Wild et al., 2005) . Once the threshold for end temperature is reached, the defrosting is stopped, and the fans are switched on again. REMOTE CONTAINER MONITORING & MANAGEMENT A recent development in container monitoring is called remote container monitoring and management (RCM (CP, 2020; Freitag, 2021; Hamburd -Sud, 2021; Jedermann et al., 2017; Maersk, 2020) ). Now , real- time information is available in the cloud on the set point, supply and return air temperature, humidity in the container , and machine operation ( power -off periods and fan settings) . For controlled atmosphere transport , also oxygen and carbon dioxide levels are monitored. The transfer of information from the container is done via a cellular network -enabled device on the refrigerated container that transmits data when connected and logs data when no connection is available. The key beneficiaries of RCM are owners of refrigerated containers or transport -service providers. They can mitigate with the use of RCM the risk of erroneous set -point adjustments and malfunctioning containers (Jedermann et al., 2017) . In addition, RCM increases the shipments' visibility and transparency and allows to replan shipment logistics in advance based on the hygrothermal history of the cargo. N otifications and alarms can be sent during transit , including unwanted door openings . For ambient loading, without precooling, the initial pulldown can be monitored. This is essential for banana transport and some citrus shipments (Thijs Defraeye et al., 2015c) . RCM enables , in pri nciple, in-transit adjustments of the set points for delivery air temperature and ventilation settings. However, despite this opportunity to adjust the set points in transit, tailoring such interventions by the stakeholders for each shipment requires the r ight metrics for decision -making. Stakeholders need to have detailed, real- time information on how these environmental conditions within the cargo (H 2O, CO 2, O2, C2H4) affect the fruit quality attributes, which is , to date, still challenging. The data that is logged by the container via RCM can be made available to stakeholders, for example , importers. Note that each large refrigeration unit manufacturer often has its own (cloud- based) software to access these RCM data. However, a second line of monitoring of the hygrothermal conditions of the cargo is often done in parallel. One or more sensors are placed into the cargo, typically the last pallet near the door end , during loading. This sensor is placed, for example, by the company packag ing and shipping the fruit. These sensor data are available via another (often cloud - based) platform and are used, for example, by retailers or importers , to verify the conditions during transport. One reason is that this is an independent system with which retailers can trac k all their shipments, independent of which refrigerated unit brand the container was equipped with . There are , thereby , for a container often two datasets that monitor the hygrothermal behavior of the ca rgo: o ne sensed by the refrigeration unit company an d one sensed by independent sensors placed by the exporte r. Interfaces that connect the two datasets, unfortunately, seem still missing or are scarce (Freitag, 2021) . COLD DISINFESTATION Cold treatments are one of several effective phytosanitary risk mitigation measures required by an importing country. More recently, industrie s have started to move towards s ystem approaches, which also include low -temperature shipping, but drastically reduce costs, carbon footprints , and food losses (Follett and Neven, 2006; Hattingh et al., 2020) . However, standalone in -transit c old disinfestation treatment s are applied to eradicate live pests during the export voyage, such as false codling moth larvae or fruit fly larvae . For example, cold treatment protocols for false codling moth larvae (FCM) ( Lepidoptera: Tortricidae )) and fruit flies ( Ceratitis capitata and C. dorsalis ) typically entail fruit first being precooled to target temperature before loading . Subsequently, supply air temperatures are set to ensure the fruit s are maintained below 0 ºC for 16 d during transport as a probit 9 disinfestation treatment (Manrakhan et al., 2022; Moore et al., 2022) . Expo rt to countries that do not require disinfestation treatments is normally done by precooling the fruit to about 4 to 7 ºC before loading and shipping at the same temperature . Cold treatments represent a niche component of exports, as the treatment process is costly and substantially increases the risk and incidence of chilling injury to the fruit. Preprint 12 Modeling airflow and heat transfer in porous media 2.2.1 Airflow and convective h eat transfer in porous media There are four main approaches for physics -based modeling of airflow and the resulting convective heat transfer during refrigerated postharvest operations for fruit and vegetables ( Figure 8). They are listed in increasing order of complexity : porous medium modeling with the one-phase approach (PM one -phase), the two-phase approach (PM two- phase), multiscale porous -medium modeling (PM Multiscale) , and the object -resolved approach . Figure 8. Different approaches to model heat and mass transfer in porous bulks of fresh horticultural produce (Ω is the domain type for porous medium, fluid or solid; CHTC is the convective heat transfer coefficient ; REV is representative elementary volume ). DISCRETE , OBJECT -RESOLVED APPROACH In the discrete, object -resolved approach, indiv idual products and the heat and mass transport are modeled explicitly based on computational fluid dynamics (CFD). The Navier -Stoke s equations for turbulent flow are solved in a computational domain , and every fruit is resolved explicitly. This detail captures the turbulent microclimate around each fruit and the thermal gradient s within the fruit. The loss of momentum due to pressure drag and skin -friction drag is explicitly accounted for by solving the flow field and boundary layer around each fruit. This object -resolved methodology has been used in postharvest science to analyze a selected number of fruit in ventilated packaging in precooling processes (T. Defraeye et al., 2013a; Delele et al., 2013; Gruyters et al., 2020; Han et al., 2015; O’Sullivan et al., 2016) , refrigerated containers (Berry et al., 2022; Thijs Defraeye et al., 2015b; Wentao Wu et al., 2019; Wu et al., 2018) and cold storage at producers, wholesalers , retailers (A. Ambaw et al., 2013b; Delele et al., 2019) , among others . The focus was, among others, on improving uniformity in cooling w ithin the shipmen t or optimiz ing ventilated packaging design. This high -resolution methodology is computationally costly, especially when modeling large quantities of products, for example, a refrigerated container containing over 100, 000 products. To our best knowledge, the maximal amount of fruit that was modeled discretely was a full pallet containing 5,000 fruit (Wentao Wu et al., 2019; Wu and Defraeye, 2018) . POROUS MEDIUM APPROACH Alternatively, the less computationally -expensive porous medium approach can be used (Ambaw et al., 2014; Hoang et al., 2015) . This lumped approach simulates the porous assembly of packed fruit and packaging materials without PMone-phase PM two -phase PM Multiscale Object -resolved approach * Fruit s are modeled homogeneized * Local thermal equilibrium between fluid and solid* Fruit s are modeled homogeneized * No local thermal equilibrium between fluid and solid * Temperature in solid uniform* Fruit s are modeled homogeneized * No thermal equilibrium between fluid and solid * Temperature gradients in solid* Each fruit modeled individually * Thermal gradients in theairand within thefruit Porous medium Porous medium Porous medium REV REV REV Preprint 13 modeling each individual fruit. Instead , the fruit -carton -air system is represented as a porous medium. The transport equations for momentum and heat transport are adjusted accordingly. For momentum transport, volumetric sink terms are inserted in the Navier -Stokes equations to account f or additional pressure losses from viscous, inertial, and near -wall confinement effects (section 2.3.2 ). Convective heat or moisture transfer in the porous medium is tackled via the single -phase a nd two -phase approaches . The Biot number is an important parameter to evaluate which method can be accurately applied . This number is show n in Figure 9 as a function o f fruit diameter and airspeed for fruit with the thermal properties of an orange. Figure 9. Biot number as a function of the fruit diameter and the approach flow airspeed for a single fruit for orange fruit. The convective heat transfer coefficient to calculate the Biot number is taken from Eq. (49) for a single sphere , so using a porosity of 1 in this equation . THE BIOT NUMBER The Biot number is defined as the ratio of the convective heat transfer coefficient h c [W m-2 K-1] and the characteristic length L H [m] to the thermal conductivity of the solid λ s [W m-1 K-1]. cH shLBiλ= (4) The characteristic length for heat transfer L H can be estimated as the ratio between the product's volume V [m3] and the surface area via which heat is exchanged with the environment A sf [m2]. Small temperature gradients inside thefruitLarge temperature gradients inside thefruit Preprint 14 H sfVLA= (5) For a spherical product, this characteristic length can be written as a function of the particle radius r p [m]: 3 24 3 43pp H prrLrπ π= = (6) The characteristic leng th in the Biot number is thereby different from that of the Reynolds number, which is the diameter of the fruit. VALIDITY of PM APPROACH and the REV In principle, the porous medium approach models the fruit -packaging system as a continuum . This implies that each macroscopic poin t is assumed to be the center of a representative elementary volume (REV , (Defraeye, 2014) ), where all these phases are superimposed on one another , implying that they cannot be distinguished separately . The porous material is thus modeled as a homogeneous material in which microscopic heterogeneities are ave raged over each REV. These heterogeneities are, for example, vent holes or holes between the fruit. Transport inside the porous material is modeled by (effective) material properties, which are obtained experimentally or via numerical simulation at smaller scales. Examples are air permeability or thermal properties. The complex transport pathways and processes at the microscale level, so through the ventilated packaging, are included in a lumped way in the material transport properties and transport equatio ns (Defraeye, 2014) . As such, the packaging and the fruit are n ot modeled discre tely. For this continuum representation of the porous medium to be theoretically valid, the R EV size (LREV) should be sufficiently larger than the microscopic heterogeneities (e.g. , the size of a pore or a vent hole, l MIC). This is essential for the averaging procedure of a certain quantity of interest, e.g. , the density or porosity , to be accurate over that REV. In principle, the homogenization approach is not valid if the materials or phases are not homogeneous within a REV. If a REV is too small, i t can fall fully within a pore or the solid material (e.g., fruit), by which the used homogenized properties are not valid there. In principle , the REV should also be smaller than the distance over which macroscopic variations in material properties (l MAC) are found. On the other hand, a REV should also be small enough to avoid variations of quantities within the REV (e.g. , temperature) due to macroscopic gradients and associated non - equilibrium conditions at t his microscale level. The latter is called the principle of separation of scales, namely l MIC << l REV << l MAC. Since the continuum transport equations are solved on a computational mesh, the size of the computational cells with which the porous medium is meshed should typically be larger than the REV . If the mesh is smaller , then one mesh cell could include, for example, only part of a single fruit, by which the material properties in such a computational cell cannot be assumed to be homogenous. In the case of palletized fruit in ventilated packaging, we face the problem that the micro and macro scale s lie quite close to one another. The microscale would be the size of the fruit and the pores, so about 0.10 m. The macroscale would be the width or height of the pallets, so 1 -2 m. The size of the REV should lie well between these two values: 0.1 m << l REV << 1 m, which is difficult to achieve . In our case, the REV size and the mesh are typically not much larger than the pore size or vent holes . In some regions, the mesh can even be smaller to resolve airflow and heat transfer in a numerically stable way . It is well possible that these mesh elements that are smaller than the REV have no effect on the solution. Nevertheless, it can induce some differe nces, but this was not investigated to our best knowledge. However, it is computationally very demanding often to model the entire cargo with a discrete approach . Note that for very small fruits, such as berries, the size of the REV could be taken such that it satisfies better the separation of scales. SINGLE -PHASE APPROACH In th e single -phase approach, th e air and solid phases are lumped into a single phase. This approach assumes that both materials are well- mixed entities (Berry et al., 2023; Getahun et al., 2017a) . The effective thermal properties are Preprint 15 derived via averaging , based on the relative amount of fruit versus air and the porosity. Local thermal equilibriu m (LTE) between the gas and the solid phase (fruit) exists. Thereby, all materials are always at the same temperature. Such an approach is, however, only realistic if the differences in thermal properties between the two materials are not very large and if the surface area for heat exchange is large, so if the diameter of the products is small (ASHRAE, 2010b) . In addition, local thermal equilibrium assumes tha t airspeed s are rather low so that equilibrium can set in. This situation is not the case for fruit stacked in a pallet and cooled with air, as there we have: (1) rather large particles compared to the flow path length, i.e. , 3.7% (= 0.08 m / 2.16 m); (2) non -negligible airspeed (e.g., 0.02- 0.06 m s-1 for refrigerated transport); (3) a large difference between the thermal conductivity of air vs. (o range) fruit (λ f/λa = 22). Furthermore, t he air phase will respond much faster to a temperature change than the fruit due to the large specific heat capacity of the fruit, compared to that of the air (c p,f/cp,a ≈ 103). Therefore, the single -phase approach is rarely used for this application and is also not viable . For example, it is more viable for hydro -cooling or when only natural convection is present. Finally, this approach also assumes that n o thermal gradients are present within the fruit. This implies low Biot numbers, s o low airspeed and/or small products. TWO -PHASE APPROACH The two -phase approach s eparately solves transport equations for the air and solid phases (section 2.3.2 ). A convective transfer coefficient (h c [W m-2 K-1]) accounts for the convective exchange at t he air -product interface. In addition, the specific surface area for heat exchange between the fruit and the fluid needs to be specified (Asf [m2 m-3]) (section 2.3.2). The two phases are not necessarily at the same temperature , and a local thermal non -equilibrium (LTNE) is present. However, the thermal gradients inside the fruit are still not accounted fo r with this approach since fruit size or shape is not explicitly included in t he conservation equations. This is valid for low Biot numbers, so low airspeed and/or small products. However, as a workaround, the convective heat transfer coefficient and speci fic surface area for heat exchange can be artificially adjusted to account indirectly for fruit size and shape and, thus, thermal gradients inside the fruit. Both porous medium approaches are powerful. However, their key limitation is that t hey do not account for thermal gradients inside the fruit. At high Biot numbers, such gradients occur during precooling or for very large fruit . In this case, a porous medium approach does not enable distinguish ing between pulp and surface temperature. Workarounds have included adding an internal resistance in addition to the convective heat transfer coefficient to calculate the average temperature inside the fruit (Hoang et al., 2015) . This resistance includes an equivalent penetration depth, dependent on the cooling speed. This limitation of not including the temperature gradients in the fruit hinders pinpointing critical events occurring at high Biot numbers . We do not know if the core and surface temperatures are different, by which quality differences can occur within the fruit . We cannot accurately identify sites where surface condensation on the fruit occurs. The pulp temperature is also not directly available since the kinetics of the cooling of individual fruit and the local microclimate is not captured , although being a key parameter to evaluate fruit quality . MULTISCALE APPROACH Multiscale modeling is an alternative to the object -resolved approach that captures each fruit's thermal gradients while simultaneously capturing the entire shipment's thermal behavior . Multiscale modeling implies then modeling from the single -fruit -scale to the shipment scale. This method is not yet widespread in postharvest science , even though several multiscale modeling approaches exist for porous media (Helmig et al., 2013) . Both one -way (sequential) and two-way (parallel, concurrent) coupling can be done between different scales (Ho et al., 2013) . Multiscale models in postharvest research have mainly focused on momentum transport, so fluid flow. Rarely heat and mass transport are considered . Fruit quality evolution within shipments is also not treated. Example of advanced multiscale model setup . In postharvest refrigerated transport, a multiscale model could be built up as follows, among others. First, a two -phase porous medium model of the fruit in the shipment is buil t based on existing porous medium models (e.g. (Alemayehu Ambaw et al., 2013; Ambaw et al., 2014; Hoang et al., 2015) ). The fruit and packaging act as the porous medium's solid matrix , and (moist) air acts as the pore space. This macroscale model is extended to include relevant transport processes , for example : (1) loss of momentum via Darcy - Forchheimer- Brinkman; (2) turbulence generation and dissipation by both volumetric source and sink terms in the Preprint 16 transport equations, namely for turbu lent dissipation rate and turbulent kinetic energy; (3) convective heat and moisture transfer at the interface of the air with the fruit; (4) radiation scattering and absorption; (5) phase change (condensation, evaporation) based on the fruit's surface tem perature. This macroscale model requires different model coefficients to calculate airflow, heat, and mass transport . Examples are viscous and inertial pressure loss coefficients, turbulence model coefficients, and convective transfer coefficients. These t ransport parameters can be separately determined empirically or by calculating the same transport processes on a representative elementary volume (REV) of a periodic microstructure, typically a small ensemble of fruit such as a carton . This REV model calculates these transport parameters based on packing density, open porosity, tortuosity, product size, etc. This information on the spatial variability of the environmental conditions throughout the entire cargo at the shipment scale is then used in a further downscaling step (localization) to calculate – concurrently – the hygrothermal state and result ing key quality attributes (KQA ) of every single fruit in each REV. This step is essential to capture internal gradients within each fruit . The macroscale and microscale simulations run simultaneously with their internal time steps. For example, the information exchange between the models can be done via a temporal coupling using an explicit (staggered) coupling procedur e. This coupling is d one at predefined points in time, namely well below the typical time scales for heat transfer or quality decay within the fruit. Smaller time steps are required with an explicit coupling compared to implicit coupling, but an explicit coupling is preferred for its enhanced stability. ALTERNATIVE APPROACHES Hybrid approaches have also been used, where ventilated packaging and its ventilation openings have been modeled discretely, but the fruit inside these packages is modeled as a porous medium (Berry et al., 2019) . As mentioned above, the separation of scales is not guaranteed in our case since the size of the REV in the porous medium, and so the mesh is of the same size as the vent holes and pores between the fruit. As such, using the hybrid approach can also induce some differences with a discrete approach, but this needs to be investigated. APPROACH used in this study This study focuses on cooling fruit on pallets inside a refrigerated container. The maximal possible airspeed s in a container are: (1) a superficial airspeed of 0.048 m/s, assuming that all air (4150 m3 h-1) flows through the 20 pallets and no airflow bypass is present; (2) a corresponding phy sical airspeed of 0.126 m s-1, assuming a po rosity of 3 8%. Note that locally higher airspeed s can occur near the vent holes. A t these airspeed s, the Biot number will be below about 0.12 (Figure 9) for the fruit of a diameter of 80 mm when calculating the CHTC according to Eq. (52). This Biot number implies the thermal gradients within the fruit will be rather limited . As such, we use the two -phase approach in this study. In a previous study (T. Defraeye et al., 2015a) , we already identified rather limited differences between fruit average and core temperatures for this airspeed range . As an illustration, we com pared the resulting SECT in Preprint 17 Figure 10, calculated using core and average fruit temperature. We get rather small differences below 20%, except for the cartons on the bottom of the pallet. These differences bec ome larger at higher airspeed s, so higher Biot numbers. Figure 10. The relative difference in the seven -eighths cooling time , calculated based on fruit core temperature and averag e fruit temperature at different airspeed s for different layers of cartons in a pallet for vertical airflow through a stack of palletized boxes (data from (Thijs Defraeye et al., 2015b) ). 2.2.2 Momentum loss in porous media We model the momentum loss in porous media by using the Darcy -Forchheimer- Brinkman equation . This equation relates the pressure gradient over the porous medium ( p∇[Pa m-1]) to the velocity (u [m s-1]). This equation originates from Darcy's law . Darcy's law is derived from the Navier -Stokes equations, which is simplified to the Stokes equation by assuming stationary, incompressible , and creeping flow without any body force. When it is assumed the viscous resistance increases linearly with the velocity and an isotropic porous medium is present, Darcy's law is obtained. QkqpAµ== −∇ (7) Here q is the volumetric flux [m s-1], Q is the volumetric flow rate [m3 s-1], A is the surface area [m2], μ is the dynamic viscosity o f the fluid [Pa s] [kg m-1 s-1] and k is the permeability [m2]. The Darcy -Forchheimer- Brinkman equation is (A. Ambaw et al., 2013a; Datta, 2007a; Verboven et al., 2006) : 22 , 1/2a Fa a eff Da Fo a Brcp CC Ckkµρµρ ∇= − − + ∇ = + + ∇ u uu u u uu u (8) 1/2 ,a Da F Fo Br a effCk cCk Cµ µ=− =− = (9) Here, p is the pressure [Pa], μ a is the dynamic viscosity of the fluid, which is air in this case [Pa s] or [kg m-1 s-1], k is the Darcy permeability [m2], c F is a dimensionless For chheimer coefficient [ -], ρ a is the fluid density [kg m-3], μ a,eff (Verboven et al., 2006) is the effective viscosity ( μa,eff = μ aτ/ε), and u is the average air velocity. Sometimes an alternative Preprint 18 effective viscosity is used in the equations, namely μa,e (Datta, 2007b) (μa,e = μ a τ = ε μ a,eff). CDa, C Fo, and C Br are the Darcy, Forchheimer, and Brinkman loss coefficients. The mentioned velocity is the average velocity within the porous medium volume . This includes both solid matrix and pore space. This air velocity and the corresponding airspeed are called superficial (discharge) velocity u and speed U [m s-1]. Th e physical (intrinsic, seepage) velocity v or speed V in the pore space is higher and equals v = u/ε (or u = v.ε), where ε is the porosity (Datta, 2007a; Verboven et al., 2006) . When substituting the superficial velocity u with the physical velocity v , we get: 22 , 1/2 22 , 1/2a Fa a eff a Fa aecpkk cpkkµρε ε εµ µρε εµ∇= − − + ∇ ∇= − − + ∇v vv v v vv v (10) Darcy -Forchheimer -Brinkman terms In this equation, the first term (Darcy) on the right -hand side represents pressure loss due to viscous forces, which increases linear ly with the velocity. The second term (Forchheimer) represents pressure loss due to inertial forces occurring in the flow, which increases quadratically with the velocity. The Forchheimer term does not imply turbulent flow, however. The third term ( Brinkman ) represents pressure loss due to viscous forces occurring in near -wall flow due to confinement, which increases linearly with the velocity. Each of these terms is only significant under certain conditions . Darcy dominates under low -speed conditions , whereas Forchheimer dominates at high speeds. The Brinkman term is only relevant close to the exterior walls or surfaces with which the porous medium is in contact (via the divergence of the velocity gradient), as only at these locations are velocity gra dients present. Particle Reynolds number The particle Reynolds number (Re p) can be determined based on the approach flow airspeed around the particle and a characteristic length for fluid flow L F [m], for which the particle diameter ( dp) is typically chosen. This characteristic length is representative of flow around the sphere, for example , the flow structures such as vortices that are generated downstream of a sphere. The approach flow airspeed in the pore space is the physical speed V in the porous medium , so around the particles . The Reynolds number can be written as: ..Repp pVd Udρµ εν= = (11) This definition is in line with the typical definition of the Reynolds number . Previous work on flow in packed beds also used this (Naghash et al., 2016; Wakao et al., 1979) . This Reynolds number increases if the porosity decreases for a certain approach flow (superficial) airspeed U. The reason is that the physical airspeed to which the particles are subjected to, increases. This characteristic length for fluid flow (dp) differs from that of the Biot number for heat transfer, which was dp/6 (Eq. (6)). Note that other au thors used an alternative definition, which was derived differentl y for packed beds (Alvarez et al., 2003; Verboven et al., 2006) : ()()..Re11pp pUd Udρµεν ε= =−− (12) Here the Reynolds number decreases if t he porosity decreases, which seems counter -intuitive as the physical airspeed increases with decreasing porosity . This Reynolds number is used to characterize the flow regime, where Re p lower than 1 -10 are required to have purely Darcy flow (Datta, 2007a) . Sphere packing density The highest sphere/solid packing density (θs=1-ε) that one can achieve for equally -sized spheres is 74% so a porosity of 26%. Irregular sphere packaging of equally -sized spheres leads to a maximal packing density of 64% so a porosity of 36%. The minimal porosity of a porous ensemble of equally sized sph eres is thus higher than 26%. A simple cubic lattice, where spheres of a diameter D are stacked on top of one another, has a porosity and solidity of: Preprint 19 3 3 3 34 41 1 321 1 48%32 6 1 52%tot sphere a tot tot sDDVV V VV Dπ φ ππ θφ−−  == = = − = −=  =−= (13) Here V tot is the total volume [m3], Vsphere is the volume taken up by the sphere, V a is the volume taken up by the air. Ergun equations The permeability k , Forchheimer coefficient c F and the effective viscosity μ a,eff depend on several factors, including the pore size, pore geometry, porosity , and tortuosity of the porous material. Here μ a,eff (μa,eff = μ aτ/ε) is often simplified to μa,e/ε, so neglecting any tortuosity effects. The permeability k and Forchheimer coefficient cF can be determined via analytical/empirical formulations. A typical ex ample is the Ergun equations for s pherical particles (Ergun, 1952; Verboven et al., 2006) : ()2 1 23 2 1/2 3/2 11 1 p FK kd KcKε ε ε−= = (14) K1 and K2 are constants equal to 150 and 1.75 for randomly stacked spheres, respectively, although also values of 180 and 1.8 have been reported (Verboven et al., 2006) . Here, dp [m] is the particle diameter. This diameter is straightforward to determine for spherical particles, namely twice the particle radius r p. For spherical -like particles , the equivalent particle diameter can be calculated as 6V p/Ap, where V p [m3] and A p [m2] are the particle volume and surface area. This equation reduces to 2r p for a spherical particle: 3 24 366 24pp p pprVrArπ π= = (15) When neglecting the Brinkmann term, t he pressure gradient and pressure loss over a certain distance L can be written as: ()() ()() () ()1/2 1/222 11 2 23 23 1 / 23 / 2 1 2 1 2 23 3 2 23 311 11 150 1 1.75 1a Fa aa pp aa pp aa pp Da Focpkk KK Kpd dK K Kpdd pL Ldd p CL CLµρ εεµρε εε εµ ερεε µε ρε εε∇= − − −−∇= − − −∇= − −− −−∆= + ∆= +u uu u uu u uu u uu u uu (16) In postharvest problems, where ventilated and non -ventilated packaging is included (Figure 1), such analytical formulations for momentum loss over stacks of spheres are less realistic . The reason is that in addition to the fruit size, shape , and stacking within the packaging , the pressure loss is also dependent on the packaging design (vent hole size, shape , and placement) and the stacking pattern of the packages on the pallet. The packaging often introduces a large pressure drop over the limited ventilated openings. Therefore, k and c F for the composite system of packed fruit or vegetables are typically estimated experimentally or via physics -based simulations by quantifying the pressure Preprint 20 drop over the porous medium at a certain flow rate or vice- versa (Getahun, 2017; Getahun et al., 2017b; Tiamiyu, 2020) . Darcy versus Forchheimer for porous media We quantify Darcy and Forchheimer's relative contribution to the pressure gradient (Eq.(8) or Δp/L). We calculate the pressure loss over a unit length (L) of 1 m. We quantified each contribution for equally sized spheres of different diameters with irregular packing . We relied on the Ergun equations (Eq. (14), (16)), so a stack of spheres without the cardboard packaging. We assumed the following: (1) a porosity of 36%, which is the minimal porosity that can be achieved for an irregular packing of spheres, (2) a fruit diameter of 80 mm and 10 mm (e.g., orange fruit vs. mandarin fruit); (3) The results are shown in Figure 11. The dependency of the pressure loss ( Δp) on the Darcy and Forchheimer terms is as follows: (1) The larger the particle size becomes, the more dominant the Forchheimer term ( proportional to 1/d p) becomes at a particular speed compared to the Darcy term ( proportional to 1/d p2); (2) The larger the porosit y becomes, the more dominant the Forchheimer term (proportional to (1-ε)) becomes at a particular speed compared to the Darcy term ( proportional to (1-ε) 2). Preprint 21 Figure 11. (a),(b) Pressure loss as a function of airspeed , for the Darcy and Forchheimer terms , and both for two fruit sizes. In (a) also , the results at a lower (26%) and higher porosity (52%) are depicted with a dotted and dashed line, respectively. (c) Both terms together are depicted for four fruit sizes. The superfic ial airspeed range in a refriger ated container is also depicted by the grey area. The main conclusions can be made for high packing densities, so low porosities . The Forchheimer term must always be incorporated for airflow in refrigerated containers ( Figure 12). For very small particles, the Darcy term only overshadows the Forchheimer term (i.e. , CFo < 0.1 CDa) at airspeed s of 0.01 m/s, which is alre ady rather low. Darcy flow is mainly dominant, compared to the inertial term , when the particle Reynolds number is below 1- 10 (Datta, 2007a) . Even then, C Fo still plays a role. The critical velocity below which the Forchheimer term becomes negligible (CFo < 0.01 C Da) is 0.001 m s-1 for small fruit, which is very low , and higher airspeed s are always found with natural convection. 0.0010.010.11101001000100001000001000000 0.001 0.01 0.1 1 10Pressure drop [Pa] Air speed [m/s]10 mm Da. + Fo. 20 mm Da. + Fo. 40 mm Da. + Fo. 80 mm Da. + Fo.0.000010.00010.0010.010.1110100100010000100000 0.001 0.01 0.1 1 10Pressure drop [Pa] Air speed [m/s]Fruit 80mm 80 mm Da. 80 mm For. 80 mm Da. + Fo. 0.0010.010.11101001000100001000001000000 0.001 0.01 0.1 1 10Pressure drop [Pa] Air speed [m/s]Fruit 10mm 10 mm Da. 10 mm For. 10 mm Da. + Fo.Refrigerated container range Refrigerated container range Refrigerated container rangea b c Preprint 22 However, the Darcy term cannot be neglected , especially for airspeed s in refrigerated containers (0.01- 0.1 m/s). For large fruit, the Forchheimer dominates the Darcy term (i.e. , CFo > 10 C Da) at airspeed s of 0.1 m s-1. The Darcy term can, however , only be neglected (C Fo/CDa > 100) at airspeed s above 1- 10 m/s and at particle Reynolds numbers above about 1 0,000. However, even for very den sely packed stacks of spherical products (with regular or irregular packing , so a porosity of 26%- 36%), the Darcy term is often less important for standard cooling operations (refrigerated container and precool ing). Even in the worst case, where the Darcy term has the most impact (low porosities 36%, small particle sizes 0.01 m), both terms are of equal size at a speed of about 0.1 m/s. Figure 12. The ratio of Forchheimer to Darcy coefficient as a function of the airspeed for different particle sizes at a certain porosity. Also , the results at a lower porosity (26%) for a sphere of 80 mm are depicted with a dotted line. 0.010.11101001000 0.001 0.01 0.1 1 10Ratio of CFo/CDa[-] Air speed [m/s]Ratio F/D 10mm Ratio F/D 20mm Ratio F/D 40mm Ratio F/D 80mm Treshold F/D = 10 Treshold F/D = 0.1 Ratio F/D 80mmForchheimer dominant Darcy dominantRefrigerated container Preprint 23 Multiphysics continuum model 2.3.1 Computational system configuration A 2D computational model of the refrigerated container was constructed to study the cooling of citrus fruit during refrigerated transport. The model mimics the cooling of citrus fruit by vertical airflow for fruit packed in Supervent cartons and stacked on 20 pallets in a refrigerated container. This (base) model is shown in Figure 13 with its dimensions and boundary conditions, where the latter is detailed in section 2.3.6. A two -dimensional model was chosen in this study to model airflow and heat transfer in the container instead of a three -dimensional model. A key r eason for choosing a two- dimensional model is that we can run the model in an acceptable time ( ∼ hours ) on a standard desktop computer or laptop. This feature is essential if we would like to widely disseminate and use this model . Our starting point is that currently , our stakeholders have only temperature sensor data. We can upcycle these sensor data with a physics -based virtual container by turning them into actionable metrics for the entire cargo. From this viewpoint, we transform the zero -(spatial) dimensional temperature data to two-dimensional spatial information inside the container , which is a huge step forward. Running our application is impossible for 3D simulations of millions of control volumes as, apart from computational time, the RAM is als o lacking on standard computers . Airflow in a container is inherently three -dimensional, and most previous studies looked at three- dimensional container geometries. However, a two- dimensional representation can realistically approximate the 3D container fo r the following reasons: - The inlet and outlet slot s span almost the full width of the container, creating a quasi- 2D inlet and outlet flow conditions. - The height of the zones where air flows, excluding the porous pallets, is small compared to the contain er's width (2.29 m) . These cavities include the bottom cavity ( T-bar floor and the pallet base height = 0.224 m), the void at the door end of the container (width = 0. 56-0.76 m, depending on the pallet row ), and the top cavity above the pallets (height = 0.197 m). This leads to container width (Wco) to cavity height (H cav) ratios of 10.2, 3- 4, and 11.6, respectively (W co/Hcav). The minimal aspect ratio between a channel's width and height to have a two-dimensional flow is 7:1 (Dean, 1978) . Thus , flow in the enclosures in a container can be considered two -dimensional, except for the void at the door end of the container. This im plies that airflow in the inlet and outlet cavities is two -dimensional. In the void at the door end of the container, the side edges could lead to a somewhat different 3D airflow field prediction. Simplifying the computational model with a 2D representation has a significant advantage. The simulations can be run much faster due to the much lower control volum es required. The mesh size is reduced by a n estimated order of magnitude 102. This reduction also proportionally reduces the computational time. Another advantage is that the stability of the simulation can be improved by highly refining the grid near the inlet. A fine grid is required since a large airspeed reduction into the T -bar floor occurs here , which can introduc e numerical instabilities. A s in previous studies, a 3D grid would require millions of finite elements to model the airflow realistically . The main disadvantage of the two -dimensional model is that lateral heat excha nge with the outside environment cannot be directly included in the model as the side walls ar e not modeled. In addition, the impact of asymmetric stowage patterns and chimneys between four pallets (windmill configuration) can also not be studied. These losses are therefore included equivalently in the 2D model via the top, bottom , and vertical door end walls . Preprint 24 Figure 13. Computational model and boundary conditions. The computational model contains the following components . We describe each component indicated in Figure 13 sequentially , based on the air's path . We discuss the model simplifications that are made to model the component: - The inl et section, where the airflow produced by the evaporator fans , enters the container via the inlet slot. Here , our model included the entire supply air duct, starting from the outlet of the fans. This supply air duct in a refrigerated container was modeled similarly to (Senguttuvan et al., 2020) . This duct stretches along the entire container width, so it can be represented in 2D as well . The dimensions of the inlet duct were extracted from (Senguttuvan et al., 2020) , images , and CAD drawings of refrigerated containers since the ac tual dimensions are typically proprietary information. In our study, the 2D-inlet duct includes a 400 mm wide section where the fans are located, a converging nozzle part that reduces the cross -section to 100 mm, a straight bulkhead section of 100 mm wide, and a corner section with a baffle plate that guides the airflow into the narrow channels of the container floor. Note that no compone nts, such as the evaporator or heater after the fans, are included in the duct, although they are present in reality and will induce a large pressure loss. As such, the total pressure loss over the model will be lower than in reality. - A baffle plat e at the refrigerator unit side. The pallet base of the first pallets is usually placed on the baffle plate. The dimensions and shape of the baffle plate were based on those of a commercially -available baffle plate ((BuffersUSA, 2022), model 1311- 8598). - The narrow open space of channels within the T -bar floor (H tb = 63.5 mm), where the air is directed horizontally and upwards. - The orifices/openings in the T-bar floor, which induce a speed -up of the flow when exiting the bottom floor channel and a pressure loss. These openings lead to an open area (or porosity) of the T -bar floor of 55%. The openings in the T -bar floor are explicitly modeled in the 2D compu tational model, which means we include each T -bar plate separately . However, the channels normally have an opening that stretches out over the full length of the container (x -direction). In our 2D model, we needed to model them as stretching out in the y - direction. Note that the thickness of the T- bars was not modeled to simplify the model and the computational mesh. - The base of the wooden pallet s, which rests on the T -bar floor. The wooden slats cover 39% of the surface area of each pallet and thereby block a significant part of the airflow from the T- bar floor. The T- bar floor blocks 45% of the floor area, and a significant overlap occurs with the wooden pallet bases. The pallet bas es are not modeled explicitly , as the blockage by the T -bar floor already accounts for some blockage. It is difficult to realistically estimate or model how this blockage is translated in the 2D model. - The narrow open space between the wooden pallets' base and the pallet's top wooden slats. Here only wooden connectors between the bottom and top of the base are present, which take up about 10% of the volume in this space. These blocks are not explicitly included in the model. - The top of the wooden pallet, wh ere the wooden slats cover 76% of th e surface area of each pallet and thereby block a significant part of the airflow coming out of the T -bar floor. The pallet top wooden slats are modeled explicitly , namely , 9 slats of 102 mm wide per pallet . Since we hav e a 2D model, equivalent slat width and gaps between the slats (36 mm) w ere used to attain the same open area as in 3D, namely 24%. Inlet slot (width 0.4 m) Velocity inlet , Ga= 4150 m3h-1, TinVoid plugOutlet slot Pressure outlet, Ps,sl= 0 Pa T-bar floor (height 0.0635 m) Wooden pallet basePorous zone with pallets 1.2 m wideGap between palletswall (no -slip) interior wallBulkhead (width 0.1 m) Baffle plate wall (no -slip)wall (no -slip) Preprint 25 - The fruit that is palletized in cardboard cartons on the wooden pallet. For airflow, the porous -medium approach is used to represent the fruit palletized in carton boxes. E ach of the 20 pallets is represented as a homogeneous entity containing 80 cartons and 5, 120 fruit. T he porous media model represents 102, 400 fruit packed in 1, 600 cartons. In total, 9 pallets are modeled in the 2D model. The width of the 9 palletized fruit (y - direction 1.2 m wide) is correct ed to 1.212 m to account fo r the other row of 11 pallets in the container (y- direction 1 m wide). That way, the total volume of fruit in a full container in the 2D model is equivalent to the 3D model. The cartons were not included in the model for the thermal calculations , and only the fruit was considered . The carton has a much lower volume and mass than the fruit and, thereby , also a lower thermal storage capacity than the fruit (m x c p [J K-1]), only 1.5% of that of the fruit. Furthermore, heat transfer from individual fruit was governed by convection to the air , not by conduction to other fruit and then to the carton. As such, the carton's insulating effect is less critical. - Small vertical gaps betwee n the pallets of 20 mm are accounted for in the base model , by modeling them discretely . This value is based on observations in prac tice. Still, large variations in actual gaps between the pallets can occur, depending on the level of detail given when stacking the cartons on the pallet and when loading. Note that other studies have modeled the gaps as a porous medium (Moureh et al., 2009) . Since a 2D model is used, these gaps , to some extent , also account for the lateral gaps between the two pallet rows. - The open space above the pallets was modeled as an air gap. - The return air duct, where the air from the cargo hold flows back to the refrigeration unit of the container. Similar to the inlet duct, the dimensions were also based on (Senguttuvan et al., 2020), drawings and images. The outlet of the computational domain is assumed to be the inlet of the fans. As such, the only component not included in the airflow model is the fans. - The container walls are not explicitly included in the model. The thermal inertia of the walls is not accounted for. The heat that can be stored in the container walls ( ∼ m x c p [J/K]) is only a small part (4%) of the heat stored in the fruit. The field heat in the product comprises the main heat load. As such, this assumption will not largely affect the cooling behavior of ambient -loaded cargo. - The transmission heat exchange of the container with the exterior environment is modeled . These transmission losses are about 1.5- 2.6 kW (Thijs Defraeye et al., 2015c; Wild et al., 2005) . The total he at flow that enters the container from the outside equals Qtrans = U.A co,tot(Tout-Tin), where U is the thermal transmittance [W m-2 K-1], Aco,tot is the total exposed surface area of the container to the outside [m2]. This surface includes the top and bottom walls , the side walls, and the vertical wall at the door end. The vertical wall at the refrigeration unit end is not accounted for , as the refrigeration unit is pla ced there. T out and T in are the outside a nd inside air temperatures , respectively. A refrigerated container is well insulated with a theoretical U ≈ 0.2-0.3 W m-2 K-1 or U.A co,tot ≈ 40 W K-1) (ASHRAE, 2010a; Thijs Defraeye et al., 2015c) . In practice, measured values are often a bit higher ( U ≈ 0.3-0.4 W m-2 K-1 or U.A co,tot ≈ 60 W K-1). We used a U.A co,tot value of 42 W K-1, and recalculated an equivalent 2D U -value to use in the simulations. In the 2D simulations, only transmission heat losses occur via the 11.56 m long roof, floor , and 2.58 m high wall at the door end over a width of 2.29 m . We obtain an equivalent U -value of 0.71 W m-2 K-1 for the 2D simulation. - A void plug is installed at the back of the container to avoid airflow bypass via the back opening of the container . Airflow bypass would induce lower flow rates through the pallets a nd slower c argo cooling . The simplifications above were made to keep the model complexity and computational cost within limits. These simplifications render the computational model in to a more idealized case . 2.3.2 Governing equations for porous media modeling 2.3.2.1 Conservation equation for momentum in air and porous media The conservation of momentum in the air and the porous medium are modeled by solving the Navier -Stokes equations. We assume incompressible dry air flow , with a constant density and viscosity , so qu asi-isothermal conditions or small temperature differences. Solving these equations in the porous medium ( ΩPM) requires that additional sink terms be accounted for in the porous medium and that the effect of the porosity on the airflow is accounted for. Preprint 26 2 22 , 1/2PM PMaa a a Fa a a a a effpt cpt kkρρ µχ µρρρ µχ µΩ Ω∂+ ⋅∇ =−∇ + ∇ +∂ ∂ + ⋅∇ =−∇ + ∇ + − − + ∇ ∂ m,PMwww w S ww w w u uu u (17) 1 in in 0 in in PMPM PM ffχΩΩΩ= = ΩΩvwu (18) Here, ρ is the air den sity (1.29 kg m-3 at 0ºC), t is the time [s], μ a is the dynamic viscosity of air , and u is the superficial velocity [m s-1]. The impact of the porous medium on the fluid flow is accounted for by including a volumetric sink term for momentum in the porous medium S m,PM ([Pa m-1] or [kg m-2 s-2]) in the momentum equation . This source term equals the Darcy -Forchheimer- Brinkman equation (Eq. (8) (Datta, 2007a) ). Ω𝑃𝑃𝑃𝑃 and Ω𝑓𝑓 stand respectively for the porous media and free fluid domains. 𝜒𝜒Ω𝑃𝑃𝑃𝑃 equals 1 in the porous media and 0 in the free domain. Th e variable w corresponds to the interstitial/physical velocity v (v = u/ε) in the porous medium Ω𝑃𝑃𝑃𝑃 and to the superficial velocity u in the fluid domain Ω𝑓𝑓. In the porous medium, we rewrite this equation as a function of the superficial velocity or the physical velocity : 22 1/2 22 , 2 1/211 1F a a a eff a Fa a a a a effcpt kk cpt kkρµρρ µ µ µρρρ µ µεε ε∂+ ⋅∇ =−∇ + ∇ − − + ∇∂ ∂+ ⋅∇ =−∇ + ∇ − − + ∇∂vv v v u uu u uu u u u uu u (19) When neglecting the Brinkman term, we get the equation used in this study: 2 1/2 2 2 1/211 1a Fa aa a a Fa aa acpt kk cpt kkµρρρ µ µρρρ µεε ε∂+ ⋅∇ =−∇ + ∇ − −∂ ∂+ ⋅∇ =−∇ + ∇ − −∂vv v v u uu uu u u u uu (20) Written only in terms of the physical velocity, this becomes: 22 1/2a Fa aa acpt kkµρρρ µε ε∂+ ⋅∇ =−∇ + ∇ − −∂vv v v v vv (21) In the free -air domain , this equation becomes: 2 aa a ptρρ µ∂+ ⋅∇ =−∇ + ∇∂uuu u (22) Here the porosity ε equals 38 %, and assumes that the fruit and the packaging are 'solid' materials through which no air can flow. The remaining volume is available for airflow. Buoyancy Note that buoyancy effects were not included in the model equations. Buoyancy implies gravity -driven airflow generated by temperature differences in the air, which in turn generate s density differences. Neglecting buoyancy implies that no flow due to temperature -driven density difference s is modeled. In this case, we assume that natural convec tion does not play a large role, so the airflow is largely pressure -driven, and only forced -convective airflow is accounted for. A key reason for this assumption i s that the temperature does not impact the airflow field. As such, heat can be modeled as a passive scalar , and the airflow does not need to be recalculated each ti me step during the transient cooling process . This assumption greatly simplifies the solution procedure and reduces the computational cost. This is the reason why it is very rarely tak en into account in CFD studies in postharvest engineering . Preprint 27 However, s ince airspeed s in the refrigerated container are rather low in the porous medium, we need to assess how important natural convective flow is in addition to forced convection. The validit y of this assumption can be estimated based on the magnitude of buoyancy forces on the airflow by quantifying the Richardson number. The Richardson number (Ri) is the ratio of the buoyancy forces to the inertial forces: 2ReGrRi= (23) The dimensionless Grashof (Gr) and Reyno lds number (Re) are defined as : ()3 , 2 Refr s ref ref ref refgT T L Gr ULβ ν ν− = = (24) Here, β is the thermal volumetric expansion coefficient [ K-1], ν is the kinematic viscosity of air (1.461 x 10-5 m² s-1), g is the gravitational acceleration (9.81 m s-2), U ref is a reference airspeed due to forced convective airflow which is in our case the physical airspeed in the porous medium [m s-1], L ref is a characteristic length [m], for example , the diameter of the fruit or height of the stack of fruit. The stated temperature difference represents the difference between the regions of warm air, which induce a low air density, and the regions of cold air, which have a higher air density. In our ca se, this is the difference between the temperature of the fruit surface (T fr,s) and the local air temperature around the fruit (T ref = Ta,loc). Note that the Richardson number is defined here based on the local air temperature around the fruit, as these lo cal air density differences will matter. The coefficient β can be taken as equal to 1/T for ideal gasses, so the average temperature of the air and the fruit surface. We can rewrite the Richardson number as: () ()() ()() ()3 , 2,, ,, ,, 2 22 ,, ,, 22 2fr s ref ref fr s a loc ref fr s a loc ref fr s a loc ref ref ref fr s a loc fr s a loc ref ref refgT T L g TTL g TTL TT gLRiUU TT TT ULUβ βν ν− −− − = = = = ++   (25) We see that the Richardson number is inversely proportional to the square of the local airspeed and proportional to the temperature difference and the length scale . For Richardson numbers considerably smaller than one, the influence of buoyancy on the flow field can be neglected, and forced convective flow is found. In that case, heat is considered a passive scalar in the flow field, indicating that it does not influence the flow field significantly. We use this analytical estimate to calculate the Richardson number as a function of the temperature difference between the local air temperature and the fruit surface temperature (Figure 14). This temperature difference will change during the cooling process. We assume an airspeed in the porous stack of the fruit that lies roughly between roug hly 0.02- 0.2 m s-1 (Uref), depending on th e location . The superficial airspeed in our study is 0.048 m s-1, and the corresponding physical airspeed , for a porosity of 38 %, is 0.1 26 m s-1. The fruit diameter of 80 mm is taken as the length scale (L ref). Based on the results and a chosen physical s peed of 0.1 4 m s-1 (Figure 14), the Rich ardson number implies that buoyancy forces are not negligible at the start of the cooling process in a refrigerated container . When the cargo is still warm at the start of the cooli ng process, they play an important role , or if the physical airspeed s get low . The surface temperature of the fruit rapidly drops within the first hour of cooling (Defraeye et al., 2017) . In most refrigerated containers, the temperature differenc es between warm and colder zones, such as fruit surface and local air temperature, are , however , rather limited. This analytical estimate suggests that for carg o that is not well ventilated (physical airspeed V <<<) due to a lot of airflow bypass, or cargo with high porosity (so V<<), and temperature differences above 1 K, buoyancy effects will play a role in addition to forced convection . Investigating the impact of buoyancy in such cases is a topic of further research , but will pose challenges to numerical stability, convergence , and computational cost. Preprint 28 Figure 14. Richardson number as a function of the temperature difference between the fruit surface and the local air temperature for different airspeed s. Turbulence modeling Airflow in the container and within the palletized cargo is turbulent. The Reynolds number at the inlet of the container is 36,300 when based on the inlet height (63.5 mm) and the airspeed (7.93 m s-1) entering the T -bar floor via the baffle plate. The Reynolds number within the porous medium is 730, based on a characteristic (physical) airspeed (0.126 m s-1) and a length scale that equals the diameter of the fruit (80 mm) . Airflow will also be turbulent w ithin the container since length scales are typically larger and airspeed s higher. The continuity equation and Reynolds -averaged Navier -Stokes (RANS) equations for momentum and energy are solved using computational fluid dy namics (CFD). The airflow within the container and cargo is assumed to be a stead y, turbulent , and incompressible fluid flow of dry air. These corresponding transport equations can be found in several fluid -mechanics textbooks . The conservation equation fo r momentum transfe r (Eq.(17)) need s to be adjusted . In principle, only the dynamic viscosity of air μ a is adjusted , so it equals the effective viscosity for turbulent flow (μa = μa,eff,T = μ a,g + μ T). Here, μ a,a is the molecular dynamic viscosity of air ( 1.79 x 10-5 kg m-1s-1), and μ T is the turbulent viscosity [kg m-1s-1]. The Navier -Stokes equations for turbulent flow are solved by applying RANS with th e standard k -ε turbulence model. This turbulence model is still the most commonly used in computational fluid dynamics (CFD) engineering and has been extensively tested (Casey and Wintergerste, 2000). This is an eddy -viscosity model , which assume s the isotropic behavior of turbulence . The model includes two addi tional turbulence transport equations for turbulent kinetic energy and turbulence dissipation. The turbulent viscosity (μT), used to determine the Reynolds stresses, is calculated by assuming a s pecific relation between the transported variables, such as the turbulent kinetic energy (k TM [m2 s-2]) and the turbulence dissipation rate (ε TM [m2 s-3]). Here μ T [kg m-1 s-1] equals : 2 TM T TMkCµµρε= (26) Where Cμ is an empirical constant, often taken at 0.09. Note that eddy -viscosity models have some limitations: - The turbulent viscosity is presumed to be an isotropic scalar quantity. This assumption implies that turbulence is isotropic as well. For many fl ow fields, this is not the case. The rapid expansion of the jet into the co ntainer enclosure in this study likely does not result in isotropic turbulence. - The turbulent viscos ity, used to calculate the Reynolds stresses , is determined by assuming a relation be tween the transported variables, namely the turbulent kinetic energy (k) and the turbulence dissipation rate (ε) . This relation was established (empirically or theoretical ly) for specific flow conditions. No turbulence model Start ofcooling processCooling almost finished 0.00010.0010.010.11101001000 0.1 1 10Richardson number [ -] Temperature difference [K]0.02 m/s 0.08 m/s 0.14 m/s 0.2 m/sMixed Forced conv .Natural conv . Preprint 29 is universally valid for all classes of flows (Casey and Wintergerste, 2000) . As such, a model will perform better under specific flows than another. - Most eddy viscosity models are developed to predict turbulence in the high -Reynolds nu mber turbulent core of the flow. The model parameters are less accurate in the near -wall region, i.e. , the low -Reynolds number region. Challenges when m odeling slot-ventilated enclosures Refrigerated containers can be considered slot -ventilated enclosures. Air enters from the bottom of a refrigerated container at high speed into the enclosure and is extracted at the top. Much research has been done on this specific configuration (Kang and van Hooff, 2022) , both on empty enclosures and encl osures filled with ventilated packaging, for example. Modeling this flow field is rather challenging for several reasons. The airflow enters the cargo hold at high speed and then expands into a much larger space, by which a fast reduction in the speed occurs. The Coanda effect can make this air jet to remain attached to the bottom wall. This rapid airflow expansion combined with the Coanda effect is challenging to accurately model by RANS turbulence models . The jet ca n be laminar, turbulent , or in the transition regime. The Reynolds number in our refrigerated container is 35, 000, based on the inlet speed at the slot inlet (7.9 m s-1) and the slot height (63.5 mm). At these Reynolds numbers, the jet is highly turbulent, so not in the transitional regime. In addition, the Coanda effect and flow expansion will also be affected by the fact that the container is filled with porous ventilated packaging and that the T -bar floor is present. The main characteristics and findings of previous studies (Kang and van Hooff, 2022) are that often uniform inlet boundary conditions have been used, and most slots spanned the entire width of the en closure. The Reynolds numbers of most studies were below 10,000, so slightly lower than in our refrigerated container. Clear differences between turbulence models have been found. Boundary -layer modeling Wall f unctions are used to model transport in the b oundary layer. Wall functions are often preferred since they provide increased computational economy and easier grid generation compared to low -Reynolds number modeling (LRNM) . Instead of resolving the boundary layer explicitly , these functions calculate the airflow in the boundary -layer region using semi -empirical functions. Despite being less accurate in some cases (Defraeye et al., 2010; Thijs Defraeye et al., 2013b) , wall functions are often the only viable option as LRNM grid generation and computation is challen ging for complex 3 -D configurations at high speeds (Defraeye et al., 2012; Kondjoyan, 2006) . In our model, v ery high airspeed s are found in some regions. Therefor e it is challenging to mesh the boundary -layer region fine enough to apply low - Reynolds number modeling . Resolving the transport in the boundary layer on the no -slip surfaces (cont ainer walls, T - bar floor, etc.) required a huge computational cost . As such, wall functions are the most viable option in our case . In addition, our airflow field predictions are less prone to how th e boundary layer is resolved compared to flow around aerodynamic bod ies, for example. The reason is that airflow separation, for example , at the outlet or inlet, is prescribed by the geometry . For example , the edges of the T -bar floor are fixed separation points , so airflow separation is not impacted by transport processes in the boundary layer. In this case, the use of wall functions will not necessarily compromise accuracy. In our study, the dimensionless wall distance y+ was above 11 on the majority of the container walls. In these zones, the first grid cell lies outside the laminar sublayer or the buffer layer, meaning that the entire boundary layer is approximated by the wall functions. The highest values were reached near the inlet of the cargo hold, where y+ values up to 60 were found. As such, all computational cells lie within the lower part of the logarithmic layer of the boundary layer, where wall functions ar e valid. 2.3.2.2 Conservation equation for energy The two -phase approach is used to solve heat transfer in the porous medium in the air and the fruit. This approach does not require thermal equilibrium between the air and the fruit. Fruit and air can have differen t temperatures. The fruit is, however, assumed to cool uniformly, by which core and surface temperatures are assumed equal. This approach is valid for low Biot numbers, where thermal gradients inside the fruit are not present . This implies sufficiently low airspeed s and/or small products to have a Biot number well below 1. For fruit stored in refrigerated Preprint 30 containers, such low Biot numbers are found (Figure 9). The conservation equations for a two- phase approach are in the fluid and solid domains : ()() ()()( )(), , , 1 1fr p f f f f fr p f f fr f f sf sf s f evap fr p s s s fr s s sf sf f s respcTc T T Ah T T St cT T Ah T T Stερ ρ ε ελ ε ρ ελ∂+ ⋅∇ =−∇⋅ − ∇ + − +∂ ∂− =−∇⋅ − − ∇ + − +∂u (27) Here u is the superficial air velocity. When assuming a constant density, specific heat capacity , and porosity, and neglecting the source term for evaporation, we get the equation used in this study: ()() ()()()(),, ,11f sf sf f pf f pf f f f s f fr sf sf resp s s ps s s f s fr frT Ahc c T T TTt Ah S Tc T TTtρρ λε ρλ εε∂+ ⋅∇ =−∇⋅ − ∇ + −∂ ∂=−∇⋅ − ∇ + − +∂ −−u (28) The subscript s f and s indicate the fluid and solid phase, ρ is the density [kg m-3], cp is the specific heat capacity [J kg- 1K-1] (for a ir, 1005 J kg-1K-1 at 20º C), T is the temperature [K], λ is the thermal conductivity [W m-1K-1] (for air 0.025 5 W m-1K-1 at 20 ºC). The heat exchange at the interface of the two materials is quantified by the convective heat transfer coefficient (h sf [W m-2 K-1]) and the specific surface area for heat exchange between the fruit and the fluid A sf [m2 m-3]. If required, additional source terms for respiratory heat generation or evaporative heat loss can be accounted for (Sresp [W m-3], S evap [W m-3]). Often the contribution of respiration and evaporation is rather limited. Note that evaporation will cool down both the air around the fruit and the fruit surface. It was , however , now included as a source term in the air domain only . Other studies included this as a source term at the air -fruit interface as well. However, we neglected the evaporation source term in this study. The specific surface area Asf [m2 m-3], can be calculated for spherical fruit with radius R as: () 1fr fr sf fr t frAAAVVε = = − (29) ()()2 343114 3sf fr frRARRπεε π= −=− (30) Here , Vfr is the volume of fruit in the carton [m3], Vt is the total volume of the carton and A fr is the surface area of the fruit [m2]. The porosity ε fr of the fruit inside the cartons in the palletized fruit is calculated as: ,, ,, ,1box e fr box fr fr a fr box e box e box eV V VN V VV Vε−= = = − (31) Here, Nfr is the number of fruit in a carton . The carton 's volume is calculated using the external carton box dimensions (Vbox,e). The volume of the carton is neglected here for calculating the fruit porosity . Thermally, the carton will have little impact, which is why it is now considered thermally as air and not inclu ded in this porosity. In this study, the resulting fruit's porosity εfr is 47%. The porosity of the carton , so accounting for both fruit and carton (ε), equals 3 8% and is calculated as: , ,, , ,, ,box i fr box i ai a box e box e box eVV V V VV Vε−= = = (32) Here Va is the volume taken up by the air [m3]. Here we neglect the air space taken up by the ventilation openings (so Va = V a,i). Note that this porosity is used in the momentum transport equation to calculate the actual physical speed that the fruit experiences around it. Preprint 31 In this study, the porosity εfr was used in the energy conservation equation (Eq. (28)), and ε was used in the momentum conservation equation (Eqs. (20)). This porosity assumption ε fr is essential to invoke a realistic thermal response since the carton would be considered fruit if a lower porosity ( ε) were used. However, the porosity ε is essential to have the right p hysical speed inside the porous cargo since the carton leaves less spa ce for the air to flow as it takes up a certain volume. For turbulent flow, the energy conservation transport equation is adjusted. The conductive heat flux is defined as: () ,,f eff f f f T f TTλ λλ= − ∇= − + ∇q (33) Where the turbulent thermal conductivity of the air is defined as: , ,PrT pa fT Tcµλ= (34) Here, PrT is the turbulent Prandtl number , μT is the turbulent viscosity , and c p,a is the specific heat capacity of air. The turbulent Prandtl number is typically about 0.85. In our study, it was between 0.85 and 1.4 and was calculated in the software with a so-called Kays -Crawford approximation as: 12 ,, ,0.3 0.3 1Pr 1 exp2Pr Pr 0.3 Prpa T pa T a T Ta a T pa T Tcc cµµ λ λ λµ− ∞ ∞∞    = + − −−       (35) Here, Pr T∞ is the Prandt l number at infinity and equals 0.85 and λ a is the thermal conductivity of air. Heat of respiration Since the h eat of respiration is quite low for oranges ( at 5ºC about 0.014 -0.019 W kg-1 (ASHRAE, 2010b) ), it is unlikely to have a significant impact on the cooling of fresh horticultural produce, as reported already for forced- air precool ing (Brosnan and Sun, 2001; Sadashive Gowda et al., 1997) . Nevertheless, it was included in the model to make it generally valid for other fruit with high respiratory heat production. Examples are banana fruit or mangoes, where sensors can even pick up r espiratory heat via an increase in temperature (Defraeye et al., 2019) . The heat of respiration can be included via a source term in the energy conservation equation for the heat of respiration in the fruit. The respiratory source term S r is defined in the fruit domain as [W m fr-3]. r fr rSQρ= (36) Here, the subscript fr indicates fruit, ρ fr is the density of the fruit [kg fr mfr-3], and Q r is the mass -based heat of respiration [W kg fr-1]. Oranges have respiration heats that are temperature dependent and range from 9.2 mW kg fr-1 at 0 °C up to 105 mW kg fr-1 at 25 ° C (ASHRAE, 2010b; Kader, 1997) . We took these literature data to determine respiration heat as a function of temperature by linear interpolation. A constant -value extrapolation was used for values above 25 °C and below 0 °C (Figure 15). Note that these values ar e per cubic meter of citrus fruit. In the computational model, the volumetric source term is only accounted for in the fruit , not the entire porous medium. The entire cargo (26.3 tons) produces respiration heat of 2351 W at 20º C but only 498 W at 5 ºC. The fruit is cooled down rapidly , so the respiration values are rather low. Preprint 32 Figure 15. The heat of respiration from the literature by an interpolated function for orange fruit (ASHRAE, 2010b; Kader, 1997) . Radiation exchange between the fruit Long -wave radiation exchange between the fruit and the surrounding surfaces is not included. Long -wave radiation can be neglected inside the porous medium since the temperature difference between adjacent citru s fruit is small (typically < 1 °C) during cooling, so limited radiation heat would be exchanged compared to convection heat transfer . The temperature differences between the pallets and the container walls are higher but also limited, and in these zones, the airspeed s are also high, so convective heat transfer will dominate. Latent heat due to evaporation (transpiration) or con densation from and to the fruit If the water evaporates from the fruit, h eat is extracted from the fruit and the surrounding air, and additional air cooling occurs. The latent heat of evaporation is 2.5 MJ per kg of evaporated water. The heat extracted from the air by evaporation due to mass loss is small for refrigerated container transport of most fruit . The mass loss over a 3- week day trip is typically about 1- 2% (Thijs Defraeye et al., 2016; Shrivastava et al., 2022) . Higher mass loss of about 7-10% typically render s fruit unacceptable to the importer (Ladaniya, 2008) but is common in shipments of flowers and plants . This amount im plies 0.02 kg of evap orated water per kg of fruit. This amounts to a latent heat absorption of 28 mW kg-1 of fruit when averaged over a trip of 21 days . The evaporative heat loss is similar to the heat of respiration, which is low for oranges. The latent heat of evaporation resulting from the mass loss of the fruit wa s not included in the model. The condensation and evaporation that occurs in the refrigeration unit w ere not included, as the refrigeration unit of the container was not contained in the model. 2.3.2.3 Fruit quality loss A kinetic rate law model was implemented to predict the evolution of fruit quality over time throughout the cold chain and its dependency on temperature . We do not target s pecific quality attributes, such as firmness, color, soluble solids content, color, vitamin C, flavor, or titratable acidity (Owoyemi et al., 2022) . Our kinetic rate law model predicts the change of a general index of fruit quality (Tijskens and Polderdijk, 1996) . This quality index quantifies the combined effect of the thermal history to which the fruit was subjected after harvest on the fruit quality. The quality index , therefore, accounts for thermal quality loss and not for transpiration or microbiologically -related quality decay. Next to this generic modeling of fruit quality and shelf life, individual quality attributes can be modeled using the same approach (Defraeye et al., 2019; Hertog et al., 2014) . Below w e present the modeling framework for a generic quality index and individual quality attributes . This study only target s modeling the thermally -driven general index of quality, which is later translated into the remaining shelf life (section 2.6.2). The kinetic rate law model quantifies the change of the index of fruit quality A i (Robertson, 2016; Van Boekel, 2008) , but is also applicable to an y other quality attribute : in i iidAkA dt−= (37) 020406080100120 -5 0 5 10 15 20 25 30Heat of respiration (mW/kg) Temperature ( °C) Preprint 33 Here, the subscript i indicates the quality index or specific attribute (e.g., firmness), k i is the corresponding rate constant [s-1], and ni is the order of the reaction, which dictates if and how the rate is dependent on the value of A. The order of the reaction is dependent on the attribute’s kinetics. Examples of zero- order reactions are browning due to the Maillard reaction, lipid oxidation, and enzymatic degradation (Robertson, 2016; Van Boekel, 2008) . A typical first-order reaction is vitamin loss. Many quality attributes, such as firmness and soluble solids content , actually depend on many biochemical reactions occurring in the fruit during maturation . The effect is lumped into a single attribute, serving as a measure of quality loss. This ordinary differential equation (Eq. (37)) is solved. For a constant value of k, that is, at a constant temperature (as detailed below), the quality index or quality attribute decreases linearly over time (for zero -order reactions), where the magnitude of the slope equals k . An e xponential decrease is found for first -order reactions at a constant temperature : ()()0, i iiAt A kTt= − (38) ()() 0,ik Tt ii iAt Ae C−= + (39) Here, A0,i is the quality index or attribute at the start of the cooling process (t = 0 s) , and C i is an integration constant. However, the rate constant k i is not constant, and so Eq. (37) needs to be explicitly solved over time. The t emperature dependency of the quality index or attribute was incorporated into the rate constant through an Arrhenius relationship (Robertson, 1993) : (), 0,aiE RT iikT ke− = (40) Here, k0,i is a constant [s-1], E a,i is the activation energy [J mol-1], R is the ideal gas constant (8.314 J mol-1 K-1), and T is the absolute temperature [K]. The constants k 0,i and E a,i were calibrated from quality index data as a function of time. The exact procedure on how the parameters k 0 and E a were determined in this study for citrus fruit is explained and is based on (Shrivastava et al., 2022) . The overall fruit quality index ( Asl), linked to the remaining storage life (SL), was modeled using a first -order model. The f ollowing assumptions were applie d based on information in the literature on 'Navel' orange fruit (Cantwell, 2001; Lado et al., 2019) . First, we assume orange fruit can be stored for approximately 70 d, so 10 weeks (tend) at 0 °C. This implies A sl (0 °C, 0 hours) = 100% = A0,sl and Aend,sl(0 °C, t end = 70 d so 10 weeks) = 20 %, so when the fruit is considered to be lost . This lower threshold can be chosen arbitrar ily to some extent , as the choice does not affect the stora ge time at a certain temperature. The threshold affects the shape of the quality loss curve, however. Lower thresholds result in a curve that is more steeply decreasing. If data at different points in time is available for the fruit of interest at a certain temperature, a more precise choice of the threshold can be made. Second, we assume a n increase in temperature of 10 °C doubles the rate constant or decreases storage life by a factor of 2. This result means that th e Q10 value = 2 so that orange fruit can be stored for approxima tely 35 d at 10 ºC. This Q10 value equals the ratio of the rate constants at temperatures T and T+10K (Q 10=kT+10/kT) and is typically about 2- 3 for degradation reactions in fruit (Robertson, 2016; Thompson, 2004) . These values depend on the cultivar, harvest location, and exact ripeness classification . However, th is study aimed to identify differences between treatments for orange fruit. As such, the conclusions will remain the same, even if slightly different values are taken for tend or Q 10. Calculation of parameters Based on the relations above, these parameters, the constant k 0,i, and activation energy E a,i, were calculated as follows. , 10 10 ,iT iTkQk+= (41) Preprint 34 ()(), , ,,10 10 0, 10 0,ai ai ai aiE E RT RT i EE RT RT ike eQ ke e− − + + −−= = (42) Out of this Q 10 value equation , the activation energy is extracted, assuming k0,i is constant . ()()() ()()(), , ,, , , 10 210 10 10 ln 10 10 10 5ai ai ai ai ai ai E E ET E T E E Q R T RT RT T RT T RT− − −+ + = −= = ≈ + ++ + (43) ()()()()22 , 10 1055ln ln 10 10cal aiRT RTEQ Q++≈= (44) This activation energy needed to be calibrated at a certain temperature Tcal, as slight variation will occur and the value is not entirely constant. In this study, T cal = 273.15 K = 0 ºC was chosen, so the temperature for which data was available. This leads to an exact match of the final quality attribute with the data at 0 °C, but some discrepancy with the data at other temperatures. Once this activation energy is known, the constant k0,i is extracted: () ,0,ai cali cal i E RTkTk e−= (45) The rate constant at a constant temperature (T cal), is calculated via Eq. (39), by known data on the remaining quality at this storage temperature ( tend, Ai(tend) = A end,sl) as: () ()() 0, 0,1( ) lni cal endkT t i end i i i end i i cal end iAt Ae C At CkTtA−= + −=− (46) Using the data (t end = 70 days, A 0,i = 1, A i(tend) = 0. 2 and C i = 0), we then get ki(Tcal = 0 ºC) = 0.000958 h-1 or 2.66 x 10- 7 s-1. This leads to E a,i = 44,585 J mol-1, k0,i = 321,947 h-1 or 89 s-1. We can use these data to calculate the quality index loss over time and the rate constant at any temperature . The model data at several temperatures are shown in Figure 16 for the quality index. Note that the choice of A end,sl, and T cal affect s the shape of the curve to some extent. Preprint 35 Figure 16. The overall quality index of orange fruit as a function of a constant storage temperature and time. The threshold value, where quality is considered to be fully lost ( 20%), is shown by the dashed grey line. The red dots indicate the experimental calibrat ion points where the quality index intersects with the threshold, namely at 70 d at 0 ºC and 35 d at 10 ºC. 2.3.3 Constitutive equations 2.3.3.1 Convective heat transfer coefficient at the fruit surface The convective heat exchange between the fruit and the airflow in th e porous medium is quantified by a conve ctive heat transfer coefficient (CHTC [W m-2 K-1]). This coefficient mainly depends on the fruit 's size, shape , and airspeed approaching the fruit. We evaluate existing correlation s for airflow in a packed bed of spherical particles and airflow around a single spherical fruit. We assume a single, spatially -constant CHTC over the entire fruit surface. A spatial CHTC variation can be present over the fruit surface (T. Defraeye et al., 2013b; Tagliavini et al., 2019) , but this was not within the scope of the present study since we do not discretely model the single products . We highlight the correlations for flow around a single sphere and flow in a packed bed of spheres . Single sphere For a single sphere , we get (Whitaker, 1972) : ( )0.25 0.5 0.667 0.4 a p pp a,wallNu 2 0.4 Re 0.06 Re Prµ= ++ µ (47) with ..Repp pa aaVd Udρµ εν= = ; .p p aCHTC dNuλ= ; ;ref p refUL dU Vε= = = ; ()()()2 2 ,, , 42s ref loc p s ref loc p s ref locQQ QCHTC A TT d TT dTTπ π= = = −− − Here, μa and μa,wall are the absolute viscosities of the air (subscript a) and the air at the wall, respectively, which were considered as equal in this study (μ a = 1.7 9 x 10-5 kg m-1 s-1) as the air properties were assumed constant; Pr is the Prandtl number (Pr = νa/αa), which was 0.7 in the present study. Here, νa is the kinematic viscosity of air [m2 s-1], which equals μ a/ρa and Tref,loc is the local approach flow temperature of the air that impinges on the sphere. Q is the heat flow rate [W] and A is the s phere's sur face [m2]. In this single -sphere correlation, the airspeed is the physical 0102030405060708090100 0 20 40 60Overall quality Asl[%] Time [days] @ 0°C @ 10°C @ 15°C @ 20°C @ 25°C Preprint 36 airspeed (V) since this is the airspeed around the sphere. Using these equations, the CHTC can directly be calculated as: ( )0.25 p 0.5 0.667 0.4 a p pp a a,wall 0.25 0.5 0.667 pp 0.4 a a a a,wallCHTC.d Nu 2 0.4Re 0.06 Re Pr U.d U.d 2 0.4 0.06 Prµ= = ++  λµ     µ=++    εν εν µ     (48) 0.25 0.5 0.667 pp 0.4 aa p a a a,wallU.d U.d CHTC 2 0.4 0.06 Prd    λµ= ++    εν εν µ     (49) This equation lets us calculate the CHTC directly from the superficial airspeed , porosity, parti cle diameter , and air properties. Note that the porosity here assumes that the fruit and the packaging obstruct airflow and are considered solid materials. Packed bed For a packed bed , we get (Whitaker, 1972) : ( )0.5 0.667 0.33 pb pb pb Nu 0.5Re 0.2 Re Pr= + 22 < Repb < 8,000 (50) with ,..;11ref pb p pb pb p aaCHTC L CHTC dNu Nu Nuεε λ λε ε= = =−− (),...1Re ;Re Re11pref pb p pb pb p a aaUdVL Udε ε εε ν ν νε ε−= = = =−− , ;1ref pb p refUL d UVε εε= = =−;3 246663243ppp pp pprVrdrArπ π= = = = (), 1s s ref af pb a pQCHTC AT T NuCHTCdλε ε= − −= Here T ref,af is the approach flow temperature of the air that enters the packed bed. As the energy conservation equation in the physics -based model relies on the local fluid temperature in the packed bed (section 2.3.2.2), this correlation is likely less accurate since it was derived for an entire packed bed using the approach flow temperature . Using these equations, t he CHTC can directly be calculated as : ( )0.5 0.667 ppp 0.5 0.667 0.33 0.33 pb pb pb a aaUU.d .dCHTC.d11Nu 0.5Re 0.2 Re Pr 0.5 0.2 Pr1 εε   ε ε −ε ε −ε = = += +   λ −ε ν ν     (51) 0.5 0.667 pp0.33 a pa aUU.d .d1 11CHTC 0.5 0.2 Prd εε   λ−εε −ε ε −ε= +  εν ν    (52) Preprint 37 This equation enables us to calculate the CHTC directly from the superficial airspeed , porosity, parti cle diameter , and air properties. Note that these empirical correlation s were derived for forced convection . Buoyancy effects were not explicitly accounted for, as these depend on the temperature difference between the air and the fruit. The packed bed correlation was derived for a Re pb range from 22 – 8,000 . The single sphere correlation is valid for a Re p range from 3.5 – 76,000. These Reynolds number ranges for both correlations cover the airspeed s in packed fruit in refrigerated containers. These CHTC correlations are shown as a function of the superficial airspeed in Figure 17. Also, other correlations for packed beds are shown here. Note that these packed bed correlations were originally defined as functions of Nu pb' and Rep'. Re p' is defined here bas ed on the superficial airspeed and not physical airspeed , as used above. All these correlations are independent of the porosity as they include the superficial airspeed around the particles and are based on a length scale equal to the particle diameter. 0.6 0.33 pb p Nu ' 2 1.1Re' Pr= + 15 < Rep' < 8,500 (Hoang et al., 2015; Wakao et al., 1979) (53) with .'p pb aCHTC dNuλ= ; ..Re'pp pa aaV d Udερµν= = 0.5 0.33 pb p Nu ' 2 1.8Re Pr= + 50 < Rep' < 5,000 (Imai et al., 1994) (54) with .'p pb aCHTC dNuλ= ; .Re'p p aUd ν= We conclude that there is a clear difference between the single- sphere CHTC and the CHTC predicted for packed beds. The packed bed correlations all agree well with one another. In the superficial airspeed range of a refrigerated container, the CHTCs for a single sphere and a packed bed var y between about 3-5 W m-2K-1, and 5-11 W m-2K-1, respectively . In this study, we used the correlation for a single sphere (Eq. (49)) since this gives us the local CHTC of every single particle as a function of the local superficial speed. Figure 17. CHTC as a function of the superficial airspeed in a packed bed of spherical fruit and for a single spherical fruit of 80 mm with a porosity of 38%. The superficial airspeed range used for refrigerated containers is indicated with the grey area. Different definitions for single sphere and packed bed For a single sphere, Re p is the particle Reynolds number (Eq. (11)) in these equations. This Reynolds number (Re p) was based on the air speed of the approach flow (U ref), so the physical airspeed (V) and length scale of the fruit (L ref). This length scale for spherically -shaped products is the (equivalent) diameter of the sphere (d p), as this is representative of the flow structures behind the sphere (e.g. , vortices) , making it equivalent to the particle Reynolds 0123456789101112131415 0 0.02 0.04 0.06 0.08 0.1CHTC (W [m-2K-1] Superficial air speed U [m s-1]PB (Whitaker 1972) PB (Imai 1994) PB (Wakao 1979) Sphere (Whitaker 1972) Refrigerated container Preprint 38 number. Nu p is the particle Nusselt number and equals CHTC.L ref/λa, where λ a is the thermal conductivity of air (0.0242 W m-1 K-1) and the length scale equals d p. ..Repp pa aaVd Udρµ εν= = (55) For a packed bed, Re pb is the particle Reynolds number in the packed bed. This Reynolds number was sometimes defined in a different way, namely where the characteristic length scale is a measure of the size of the void spaces (Whitaker, 1972) . Both Reynolds and Nusselt number s are a function of the porosity which implies that the impact of the void space and packing density on the airflow and corresponding heat transfer is included in the correlations. The correlation can be rewritten as: ( )0.5 0.667 0.33 pb pb pb 0.5 0.667 0.33 pp pNu 0.5Re 0.2 Re Pr 1Nu 0.5 Re 0.2 Re Pr 11= + ε ε −ε  = +  −ε −ε ε   (56) Note that both of these correlations were defined using different Reynolds and Nusselt numbers. The difference is the length scale used , where both use the local (physica l) airspeed. This should be accounted for when extracting the convective heat transfer coefficients from one of these correlations or comparing them . 2.3.4 Material properties The applied material properties are given in Table 2 and their sources. The fruit shape , size, and material properties were assumed to be identical for all cases. Note that all these can differ according to the harvest time and place and due to the biological variability between individual fruit in a batch. Thi s variability comes from variation s in the local growing conditions of each tree and within each tree, ripeness degree, harvest location, and harvest year. S mall differences in thermal properties between different literature sources are present (Bon et al., 2010) . However, the rationale of the present study was to target specific questions about the cooling and quality evolution, and our simplified configuration is sufficient to serve this purpose. The material properties of fruit were taken constant, thus independent of temperature, which was sufficient for the scope of this study. The thermal properties of air were slightly dependent on temperature but varied little in the range of temperatures con sidered in the study. Note that the carton cardboard was neglected in the thermal properties of the porous medium. Table 2. Thermal properties of the materials. Density [kg m-3] Specific heat capacity [J kg-1K-1] Thermal conductivity [W m-1K-1] Dynam ic viscosity [kg m-1 s-1] Source Fruit 955 3598 0.543 - (Shrivastava et al., 2022) Air (at 0°C & 20°C) 1.29 - 1.21 1004.8 - 1005.4 0.0242 - 0.0258 1.71 x 10-5 - 1.81 x 10-5 (COMSOL, 2022) 2.3.5 Porous -medium parameters & calibration We determine and calibrate the porous medium parameters for the palletized fruit. We model flow through 80 cartons of 80 mm orange fruit, with a porosity of 38 % (including the fruit and the cardboard boxes for airflow ) via the porous media approach. The po rous-medium parameters in the volumetric sink term in the Navier -Stokes equations for Darcy and Forchheimer are determined via a calibration procedure. In this study, we did no t account for the Brinkman term . For a refrigerated container, this term account s for the viscous losses at the container walls . They are expected to have a limited impact . The first step is to obtai n empirical data on the pressure loss (Δ p [Pa]) versus flow rate (G a [m3 s-1]) through the entity that we want to model. In our case, this entity is a pallet of fruit. These empirical data can come from wind tunnel experiments (Gruyters et al., 2019) and CFD simulations (Thijs Defraeye et al., 2013a) on a similar geometrical setup Preprint 39 using the same fruit and packaging. The fl ow rate G a equals A.U, where A [m2] is the surface area perpendicular to the flow direction , and U is the superficial airspeed entering the porous medium. Second, the volumetric sink term in the momentum equation (Eq. (17)) is rewritten to estimate this source term in the Navier -Stokes equations based on the empirical data. We rewrite it in terms of th e pressure drop over the porous medium Δp [Pa], the average superficial airspeed in the porous medium U and the path length that the flow needs to travel in the porous medium L PM, but not accounting for tortuosity ( Figure 18). When also substituting the average superficial airspeed into the sink term, we get from Eq. (8), when neglecting the Brinkman n term : 2 , 1/2 1/2 2 1/2a Fa a Fa m PM PM a Fa PM PMcc pS p UUL kk k k cp LU LUkkµρ µρ µρ∆= ∇= − = − − = − − ∆= +u uu (57) The second term in the equation can be rewritten into the well -known form with the drag coefficient C D, where C D [m-1] or C D'' [kg m-4] equals 2 cF/k1/2: 2 2 22 1/2 1/2 1/2 1/221 11''2 22 2 2''Fa F PM a PM D a PM D PM F D F Dac cp UL UL C UL C ULkk cCk cCkρρρ ρ∆= = = = = = (58) We fit the empirically obtained data points ((Δ p1, U1), (Δp2, U2), (Δp3, U3), …) with Eq. (57) of a second -order polynomial to obtain the polynomial parameters a and b. From these two parameters, we can inversely determine the coefficients cF and k that provide the best fit to the model: 22 1/2 1/2 1/2 1/2; ;a Fa PM PM a Fa PM PM a PM a PM F a PM a PMcp L U L U aU bUkk ca Lb Lkk bLbk ak Lca LLµρ µρ µ µ ρρ∆= + = + = =  = = = (59) In this study, the porous -medium parameters are determined for a pallet of Supervent packaging from previous CFD studies for vertical and horizontal flow through the pallet. For vertical flow , the pallet was simplified as a single stack of cartons (T. Defraeye et al., 2015a) . For horizontal flow, the pallet was simplified as a single layer of cartons (T. Defraeye et al., 2014) . We assumed here that the resistance in both horizontal directions is the same. The differences between both direct ions are relatively small compared to those with the vertical flow. These empirical data and the best model fits are shown in Figure 18. Also , the results from the Ergun equations (Eq.(14)) for a stack of spherical products are given. The porous medium coefficients for the Darcy and Forchheimer term s are given in Table 3 in the vertical and horizontal direction s. Note that the pressure -velocity data were fitted with a regression on the p oints at the lowest airspeed s (vertic al: 0.02354 – 0.265 m s-1, horizontal: 0.144- 0.457 m s-1). The reason is that the fit becomes inaccurate at low airspeed s when we also include fitting the data points at higher speeds . These coefficients and i mplementation w ere validated on a 2D model of a single vertical stack of cartons (section 3.1). We find that for a pallet of fruit , the airflow resistance in the horizontal direction is higher than in the vertical direction. This means that horizontal flow will be more suppressed than vertical flow, so air prefers to travel vertically through the pallet s. This corresponds to the much lower total open area of the vent holes (3%) compare d to the vertical Preprint 40 direction (1 3%). The airflow resistance for a stack of spheres without cartons , calculated via the Ergun equations, is lower than with cartons . As such, the Ergun equations are not representative of estimating pressure losses over palletized fruit packed in ventilated packaging , which is in accordance with previous findings (Lufu et al., 2020; Verboven et al., 2006) . Also f or this packaging type, the ventilated packaging seems to be a key contributor to the pressure resistance over the pallet . From the ratio of Forchheimer to Darcy terms for a refrigerated container (Figure 12), both terms are relevant to include for the present pallet of Supervent. The pressure resistance in the vertical and horizontal direction of the pallet was found to be different . In this study, the porous medium was model ed anisotropically for the Darcy and the Forchheimer terms . However, it was not possible to implement an anisotropic value of the Forchheimer coeff icient (cF) in the software used . Therefore, we used the same c F in both directions , namely the one for the vertical flow. We then refitted the value of the horizontal Darcy permeability (k) , so it agrees as well as possible with the horizontal fit (Table 3, ' Horizontal (best fit with data)') in the lower speed range. We get a slightly different permeability ( Table 3, ' Horizontal ( used in model )'), but the data agree well in the speed range in a refrigerated container ( ∼ 0.05 m s-1). These data ( Table 3, ' Horizontal ( used in model )') were implemented in the model for the anisotropic porous medium. The discrepancies of this fit will be higher at speeds that are higher than those found in a re frigerated container. Since the airflow rate is prescr ibed in the simulations, the airflow resistance will affect only the flow field and pressure distribution in side the container , not the airflow rate through the refrigerated container. In reality, the higher pressure resistance of the cargo hold will lead to a different operating point for the evaporator fans. However, given the low -pressure loss over the cargo hold (∼101 Pa), compared to the pressure delivered by the fans (∼102 Pa), the impact of the cargo on the operating point is expected to be rather low. Figure 18. Pressure loss over a pallet in the vertical and horizontal direction as a function of the flow rate through the pallet , as well as the re sults from the Ergun equations and vertical flow through a stack of spheres without packaging . The data points from the CFD simulations (dots) , as well as the second -order polynomial approximations, are shown. (a) zoomed in on the speed range relevant for refrigerated containers, (b) full range for fitting, (c) schematic of the speed in a porous medium. pin Airflow rate GaPorous material Superficial speed U LPMPhysical speed VSuperficial speed Upout Airflow rate GaΔp = pin-pout010002000300040005000600070008000900010000 0 1 2 3Pressure [Pa] Superficial air speed [m/s]Vertical - CFD sims Vertical - Fit Vertical - Ergun Horizontal - CFD sims Horizontal - Fit (best fit) Horizontal - Fit (used in model) Superficial speed in container 05101520253035404550 0 0.02 0.04 0.06 0.08 0.1Pressure [Pa] Superficial air speed [m/s]Vertical - CFD sims Vertical - Fit Vertical - Ergun Horizontal - CFD sims Horizontal - Fit (best fit) Horizontal - Fit (used in model) Superficial speed in containera b c Preprint 41 Table 3. Porous medium coefficients derived from fitting data. Configuration Direction Darcy Forchheimer k [m2] cF [-] CD [m-1] Reference Palletized orange fruit Vertical 3.750 x 10-6 1.014 1,048 (T. Defr aeye et al., 2015a) Palletized orange fruit Horizontal (used in model) 1.92 x 10-7 1.014 4,630 (T. Defraeye et al., 2014) Palletized orange fruit Horizon tal (best fit with data) 1.36 x 10-6 4.53 7,774 (T. Defraeye et al., 2014) 2.3.6 Boundary and initial conditions The boundary conditions for airflow and heat transfer in the container are depicted in Figure 13. The boundary conditions for the base case are detailed. Other variants t hat were simulated are discussed elsewhere. Airflow was solved in a steady state since the airflow rate is kept constant over time. We assume that the fruit was at ambient temperatures after harvest for solving heat transfer . The fruit is cooled down within the container. This case is representative of the ambient or 'warm' loading of refrigerated containers (Thijs Defraeye et al., 2016, 2015c) . This practice alleviates the limited throughput of precool ing facilitie s during the peak harvest season by using the containers as mobile precoolers . The initial temperature of the fruit and cardboard boxes was equal to the characteristic fruit temperature after being packed and palletized, typically between 10 a nd 20 °C for South Africa. In this study, we assumed an initial temperature of 16 °C. At the inlet of the domain, a constant airflow rate is imposed. A uniform airspeed is imposed . Its magnitude is based on the flow rate of evaporator fans in the refrigeration unit of the container. We evaluated the high -speed fl ow rate (4150 m3 h-1, Table 1), whi ch induced a speed of 7.9 m s-1 at the inlet slot of the cargo hold . The corresponding specific flow rate, i.e. , per fruit mass , is 0.0 437 L s−1kg−1. A co nstant air temperature was imposed at the inlet for the base case, the so -called delivery -air temperature (DAT) . This assumption implies that the refrigeration unit can cool the return air entirely back to the set temperature, which is often the case in pract ice. The measured delivery air temperature from sensor data was also used in other simulations, such as the validation simulations , where the set delivery air temperat ure was set at 1 °C. The inlet (DAT) air temperature was 2 °C. The reason for this low temperature is that this research was conducted for the South African citrus industry, which exported 2.4 million tons during t he 20 21 season, making it the second -largest exporter of fresh citrus worldwide (CGA, 2022) . Some export markets require a cold di sinfestation treatment as a risk mitigation protocol for pests associated with the fruit export pathway. This protocol entails a lower throughput in the precool ing facilities , and this temperature must be maintained during overseas transport. Ambient loading is particularly of interest to the South African citrus industry . Since the container h as sufficient refrigeration capacity to cool the fruit within 3-6 days (Berry et al., 2021; T. Defraeye et al., 2016; Thijs Defraeye et al., 2015c) , the use of ambient loading significantly reduces demand for precool ing facilities , which are currently in short supply. At the outlet , the ambient atmospheric pressure was imposed , and a thermal ou tflow condition, which assumes no thermal gradients in the flow direction. The container walls were model ed as no -slip walls with zero roughness. The container walls on the top, bottom , and door end exchanged heat with the environment. The heat flux q [W m-2] was calculated as: ()2,D eq ext s qU T T= − (60) Where U 2D,eq is the thermal transmittance (0.71 W m-2 K-1, (section 2.3.1 )), T ext is the external air temperature, which was taken equal to 20 °C. Note that for heat load calculations, often higher temperatures are used, such as 38 °C (ASHRAE, 2010b; Thijs Def raeye et al., 2015c) , and T s is the surface temperature of the interior container wall , which results from the simulation . At the interface between the porous medium (pallets) and the air, also internal boundary conditions needed to be specified for heat transfer and airflow. For airflow, the continuity of fluid flow and temperature was modeled. Preprint 42 The vertical wall at the refri geration unit end is modeled as an internal, thin stainless steel wall ( 1 mm). Heat can be transferred by conduction from the refrigeration unit end to the cargo hold. Due to the high conductivity of the wall, the temperature on both sides will be quite si milar . The T- bar floor is also modeled as an internal, thin aluminum wall (1 mm) . The void plug blockage was modeled as a 1 mm thick plastic coated paper material. 2.3.7 Time scales We include a discussion on the different time scales for cooling and food quality decay. More details can be found in (Defraeye et al., 2019) . The reason is that we typically observe that the temperature differences among different speeds or within the fruit are substantial . Yet, their impact on fruit quality attributes can appear relatively insignificant. This phenomenon can be elucidated by analyzing the time constants associated with cooling and quality decay . For an exponential function, the time constant of parameter X is defined as follows. ()/ 0tX t Xeτ−= (61) In the proposed first -order reaction kinetics model for overall quality, as well as some individual quality attributes described by a similar exponential function , the time constant can be derived as τ = 1/k i(T). Furthermore, the temperature -time evolution during the cooling process can be represented by the following exponential function (Thijs Defraeye et al., 2015a; Thompson, 2008) , where P denotes the cooling coefficient. ()PtYt e−= (62) Hence, the time constant can be determined as τ = 1/P. The time constants for overall quality decay, various quality attributes, and the cooling process (at each speed) were found to be substantially lower - typically more than three orders of magnitude - than those pertaining to quality decay, indicating a significantly faster cooling process (Defraeye et al., 2019) . Notably, the time constants for cooling are highest at low speeds, making them more akin to those of quality decay. Consequent ly, in typical forced convective cooling scenarios, the process often occurs rapidly, thereby limiting the impact of cooling rate differences on quality decay. However, in forced convective cooling within ventilated palletized packaging, both high speeds n ear vent holes and low speeds due to dead zones can be encountered. As a result, within a cargo, time constants for cooling can exhibit significant variations in practice, even within a single unit operation. The time constant for fruit quality decreases as the fruit ripens, underscoring the imperative of faster and more uniform cooling for ripe fruit. For the quality attributes, all of them demonstrate equal susceptibility to temperature, exhibiting similar time constants. In conclusion, employing high airflow rates during post -harvest cooling is crucial, particularly for ripe fruit. However, the variations in quality within the fruit remain relatively constant and independent of the air -cooling speed. Note that the time constant for turbulent airflow is m uch lower than that of cooling. Spatial and temporal discretization An appropriate gri d (Figure 19) was constructed for the 2D container based on a grid sensitivity analysis . The grid consists of 487,803 quadrilateral and triangular finite elements. A few layers of quadrilateral elements w ere used on the walls , and then a transition to triangular elements was implemented. A gradual refinement toward the air –porous medium and air -container wall interfaces was applied to enhance numerical accuracy and stability, as the largest gradients occur there, including at the start of the cooling process. Also, the grid at the inlet and bottom cavity (T -bar floor and pallet base) was refined since high airspeed s occur at these points. The transient simulations applied adaptive time -stepping, with a maximal time step o f 600 s. T his time step was determined from sensitivity analysis and ensured a sufficient temporal resolution for the output data. Note that refrigerated containers typically log temperature data at a 1 -hour interval, and sensors used in the cold chain during precooling or refrigerated transport log at a 10- 15 min interval. Preprint 43 Figure 19. Computational grid that was used in the simulations, with several details highlighted (pallets with fruits are depicted in orange) . Numerical implementation and simulation This model was implemented in COMSOL Multiphysics (version 6 .1). COMSOL is a finite element -based commercial software. Airflow in the air and inside the porous medium (containing air and frui t) was solved using the ' Free and Porous Media Flow ' physics interface. The steady -state a irflow was solved only once before the transient heat transfer calculation. Since the flow field was assumed steady and stable over time , as no buoyancy was included in the model, the flow field did not need to be recomputed anymore during the transient simulations . Thus the flow equations were switched off during transient calculations . Transient heat transfer in the air and the fruit duri ng convective air cooling was solved using the ‘ Heat Transfer in Solids and Fluids ' interface. Heat transfer in the fruit and the air were coupled using the 'Nonisothermal Flow ' interface . The interfacial conditions between fruit and air were defined via t he convective heat transfer coefficient. The kinetic rate law model for the fruit quality was implemented using the ‘Ordinary Differential Equation’ interface. The equations were solved for the dependent variables velocity u (u,v,w), air temperature T a, fruit temperature T fr, and fruit quality index A. The quality attributes do not affect temperature and could be solved separately after the thermal calculations. We used the following discretization s chemes: P2+P2 for airflow and quadratic Lagrange elements for heat transport in the fruit and the air. A segregated direct solver was used for steady -state airflow , relying on the GMRES solver for fluid flow and the PARDISO solver for turbulence . A fully -coupled direct solver was used for transient heat transfer airflow , relying on the MUMPS (MUltifrontal Massively Parallel sparse direct Solver) solver scheme. The tolerances for convergence and other solver settings were determined based on sensitivity analysis, so further inc reases in the tolerance did not furthe r alter the solution results. These tolerances were 10-4 for the steady -state calculation and 10-4 for the tran sient simulation. Before simulating the transient cooling process, steady -state simulations were performed to obtain the flow field and the initi al temperature conditions. During these simulations, the temperature of the fruit was fixed to its initial value (16 °C), and the inlet temperature was taken equal to the set point of the container air temperature ( 2 °C). After the steady -state simulations , transient simulations of the precool ing proc ess were performed. During the transient simulations, the airflow field did not need to be recomputed anymore . The computational cost was significantly Preprint 44 reduced since only the ener gy equation needed to be solved over time. The two -step approach is often applied for forced -air cooling applications. Metrics to evaluate cooling & fruit quality 2.6.1 Seven -eighths cooling time The output of the simulations is the history of the average fruit temperature at each point in the container. For the cooling kinetics, the cooling rate of each fruit was assessed from the fruit temperature profiles. For this, the fractional unaccomplished temperature change (Y ) was determined: , 0()fr DAT T r DA fTt TYTT−=− (63) Here, subscripts fr,0 and DAT represent the initial temperature of the fruit and the set point temperature of the delivery air in the associated cold chain unit operations, respectively. Tfr(t) represents the fruit temperature, which can be, for examp le, the core temperature of the fruit pulp or the volume -averaged fruit temperature. Only the average temperature i n this study is a vailable and equal to the surface or core temperatures . The reason is that the two -phase model cannot capture possible thermal gradients inside the fruit . Given the relatively low airspeed s and Biot numbers, the thermal gradients will be limited , so the av erage fruit temperature is representative for the entire fruit. From the definition of Y , the seven -eighths cooling/heating time (SECT, t 7/8) and the half cooling time (HCT, t 1/2) can be determined. The t 7/8 is the time required to reduce the temperature d ifference between the initial temperature of the fruit and that of the set point/delivery air by seven -eighths (Y = 0.125). The SECT is a useful parameter to characterize the cooling behavior of the fruit in each of the unit operations. The SECT is frequen tly used in commercial (pre)cooling operations because the fruit temperature at that value is acceptably close to the required storage temperature (Brosnan and Sun, 2001). 2.6.2 Shelf -life prediction The remaining shelf life can be predicted based on the remaining quality at the end of the supply chain (A i(tend)) and the quality where the food is considered to be lost (A end,sl = 20%). This remaining shelf life or keeping quality (Tijskens and Polderdijk, 1996) is the remaining time until a commodity becomes unacceptable at shelf life conditions. The remaining shelf life can be pre dicted b y taking the known final quality at the end of the supply chain, prior to shelf life storage, A(t end), and by keeping the fruit as of that point in time at constant conditions until A end,sl is rea ched. The shelf -life temperature used in this study is 23 °C. At a constant shelf -life temperature, the rate constant k i also remains constant , and the quality attribute shows an exponential decrease for the shelf- life prediction (for first -order reactions): ()() 0,ik Tt ii iAt Ae C−= + (64) We reset the time for the shelf life calculation t' to zero at the start of the calculation, so t' = t -tend. Combining this equation with the known points (t '=0, A i= Ai(tend)), (t' = tRSL, A i = Aend,sl = 20%), C = 0 (assume no shift) and ()0a SLE RT i SLk T T ke− = =, we get: () ()0, '0i i i end At A At= = = (65) ()() 0, 20%i SL RSLkT t i RSL iAt Ae−= = (66) Out of these equations, we calculate t RSL by: Preprint 45 ()() () ()() () ()()()() ( )()()()()20% 11ln ln() () 1 11ln ln ln ln() ()SL RSLkT t i RSL i end i RSL i end RSL SL i end SL i RSL RSL i RSL i end i end SL SL i RSLAt At e At AttkT A t kT A t t At At AtkT kT A t−= =  =−=     =− −= +  (67) Comparison of porous medium simulations with experiments Comparison with discrete CFD simulation of a pallet for validation of airflow and cooling We compare the porous -media simulation for cooling a stack of palletized cartons with that of a fully -resolved CFD simulation (Thijs Defraeye et al., 2015b) . This comparison allows accuracy verification of our porous media simulations . The simulatio ns evaluate the cooling of orange fruit. The fruit are initially at ambient temperature (20ºC) and cooled t o -0.5ºC, similar to the fully -resolved CFD simulations. The pressure of 1 Pa was imposed at the computational domain's inlet , leading to a superficial airspeed of 0.0235 m s-1, which is representative of airflow in a refrigerated container. The lateral boundaries were modeled as a symmetry boundary condition (slip wall), which assumes that the normal velocity component and the normal gradie nts at the boundary are zero. The average fruit temperatures for each of the 8 cartons are monitored for 100 hours , so roughly 4 days . The fully -resolved CFD simulations did not include the thickness of the carton boxes, and only the fr uit contributed to the porosity . Therefore, the porosity for airflow was adjusted, compared to the container simulations (ε = 38%), and equals 47% (= ε fr). Full-scale experiment We compare the cooling behavior of a simulated refrigerated container during ambient loading wit h a previous full- scale experiment we performed (Thijs Defraeye et al., 2016) . We aim to show that the simulations capture the cool ing process sufficiently accurate ly. The experimental details are detailed in (Thijs Defraeye et al., 2016) . The experiment evaluat es the cooling of fruit at ambient temperatures . The fruit was loaded in a refrigerate d container at around 16 °C and cooled there. Twenty pallets of 'Washington Navel' oranges fruit are loaded into Supervent packages. A void plug was installed at the bottom of the last pallets to reduce flow , bypassing the pallets. The delivery air temperature and return air temperature were monitored. The air set temperature was 1 ºC. However, as an input for the simulations, the measured delivery air temperature was i mposed at the inlet. T he fruit core temperatures in a carton are monitored in 6 pallets at 3 different heights. The shipment was monitored for about 21 days. The results of th ese validation experiments are discussed in the accompanying paper. Preprint 46 3 Results & discussion Validation of the porous -medium model We verify the accuracy of the physics -based porous media simulations. Therefore, we compare the porous -media simulation for cooling a stack of palletized cartons with a fully -resolved CFD simulation of this geometry (Thijs Defraeye et al., 2015b) . The average fruit pulp te mperature of each of the eight layers of fruit in a carton is shown in Figure 20 for both porous -medium and discrete object -resolved simulations. We included a comparison for simulations with a single -sphere correlation and a packed -bed correlation. The predicted airflow rate through the stack of cartons , for the imposed pressure difference over the stack, was below 1% of the discrete object -resolved simulations. When using the single -sphere correlation, a satisfactory agreement in cooling time and shape of the cooling curves is obtained . Still, the porous medium simulations cool down a bit faster than the CFD simulations . When using the packed bed correlation, we get a too-high cooling behavior using the porous medium approach, mainly of the bottom cartons of fruit. The reason s for the discrepancies between discrete and porous -medium simulations are the following: • The same CHTC -U correlation is assumed for all fruit in the entire pallet for the porous -medium approach . In contrast, the local relation between the CHTC and the airspeed will differ from carton to carton in the discrete model. • The porous -medium approach assumes a well- ventilated bed. In reality, fruit is packed. This packaging has little ventilation holes, by which several zones in the packaging are not well ventilated, and some of the fruit will cool less fast and also with different kinetics t han in a homogeneous packed bed. This could explain why the porous medium model predicts higher heat transfer rates. • The packed bed correlation we use d was derived from an entire packed bed based on the approach -flow temperature. However, in our implement ed porous -medium modeling approach, we use this correlation to calculate the local CHTC at each location in the packed bed based on the local airspeed . We also use the local air temperature to calculate the local heat flux from the fruit. A s such, using the packed -bed correlation for the CHTC of the full packed bed is likely less re presentative of the local CHTC and temperature we are modeling. The reason is that this correlation relates the heat flow from the fruit to th e approach flow air temperature and not to the local air temperature. All these differences are responsible for the different cooling times and kinetics (shape of the cooling curves). In this study, we use d the single -sphere correlation. This correlat ion accounts for the local temperature cl ose to the product . Since this single -sphere correlation gives an acceptable agreement with the discrete object -resolved simulations , we can conclude that it is acceptable to use it as an approximation for closely packed spheres in ventilated packaging. Figure 20. The v olume -averaged temperature of all fruit in a specific carton for each layer of cartons as a function of time using a CHTC correlation for (a) a single sphere and (b) a packed bed. Dashed lines are data from the dis crete object -resolved approach , and full lines are the porous -medium approach. (a) Single -sphere correlation for CHTC (b) Packed -bed correlation for CHTC Layer 1Layer 2Layer 3Layer 4Layer 5Layer 6Layer 7Layer 8 Airflow Discrete object resolved model-202468101214161820 0 20 40 60 80 100 120Temperature ( °C) Time (hours)Layer 1 Layer 2 Layer 3 Layer 4 Layer 5 Layer 6 Layer 7 Layer 8 -202468101214161820 0 20 40 60 80 100 120Temperature ( °C) Time (hours)Layer 1 Layer 2 Layer 3 Layer 4 Layer 5 Layer 6 Layer 7 Layer 8 Preprint 47 Exploratory findings from simulations This section gives an overview of our findings while setting up and troubleshooting the simulations. We report qualitative findings and do not provide additional data. Th is is because these simulations were exploratory , and we did not explicitly investigate each parameter in detail. However, the aim is to create awareness for modelers of specific hurdles we faced. Physics of airflow • 2D vs 3D . 3D simulations seem to have a roughly similar cooling behavior as our 2D simulations . • Laminar vs. turbulent flow . Turbulent flow significantly increases the pressure losses over the cargo compared to laminar flow. • Coand a effect . A larger inlet relaxes t he Coanda effect for turbulent flow and makes the flow remains less attached to the bottom wall. • Inlet geometry. The shape of the inlet did not significantly affect the cargo hold's flow fie ld significantly. • Baffle plate. The presence of a baffle plate did not significantly affect the overall airflow field. Locally , in the vicinity of the baffle plate , the airflow field was different. Simulation stability • T-bar floor . Modeling the T -bar floor as a thin porous medium (screen or perforated plate) gave rise t o instabilities in the simulation. The reason is likely that there is a significant change in airflow direction, namely from horizontal flow in the cavity be low the T -bar floor to vertical into the pallets of fruit. It is challenging for perforated -plate models to capture this change in the flow direction. As a result, e ven screens with very small resistances led to convergence problems. Using a larger height of the T -bar floor improved numerical stability. • Grid . It was more challenging to have proper convergence using a structured grid built out of quadrilateral elements than a triangular grid. However, quadrilateral boundary -layer elements were used on the walls as boundary -layer elements. • Porous -medium parameters . To increase numerical stability, it helped to ramp up the magnitude of the porous medium parameters, for example , by starting from 10% of the pressure resistance, in a few steps, to its full resistance. • Airspeed . Ramping up the airspeed in several steps helped to avoid numerical stability is sues. The p hysics of the recirculation zone predicted in the pallets for turbulent flow (Defraeye et al., 2022) . • T-bar and i nlet height. Increasing the height of the T-bar floo r and inlet to about twice its normal height (63.5 mm) removed the recirculation zone while keeping the airflow rate through the slot. The slot Reynolds number , however , stays the same as the airspeed decrease s with the same factor. • Slot Reynolds number. Reducing the airspeed , so the slot Reynolds number , reduces the size of the recirculation zone. • Void plug . The presence of the void plug did not impact the occurrence of a recirculation zone. Hence blocking the rear cavity to avoid airflow bypass near the door end has little impact. • Lamin ar flow . Creeping airflow (Stokes flow) instead of turbulent airflow does not lead to a recirculation zone in the cargo (Defraeye et al., 2022) . No Coanda effect (so jet attachme nt) is observed. • T-bar floor . Modeling the T -bar floor as a thin porous medium (screen or perforated plate) did not remove the recirculation zone and introduced challenges in convergence. Modeling the T -bar floor with discrete interior walls did not remove the recirculation zone but was more numerically stable. Others also modeled the T- bar floor discretely and found a recirculation zone (Getahun et al., 2017a; Tiamiyu, 2020) . Not modeling the T -bar floor reduced the size of the recirculation zone since the air jet is less restricted to remain inside the T- bar cavity and can move up more easily. The T -bar floor seems to enhance the presence of a recirculation zone. Modeling the T -bar floor as an isotropic porous medium of a few centimet ers thick removed the recirculation zone. However, this representa tion is not entirely correct as the airflow resistance in the horizontal direction is rather small in reality . • The g ap between the pallet and the wall at the refrigeration -unit end . The presence of th is gap did not affect the recirculation zone. Preprint 48 • The t urbulence intensity at the inlet did not affect the recirculation zone. An o utlook for future physics- based model development The current physics -based model can be extended to increase insights and improve accuracy. Turbulent airflow calculations are an essential future step to better capture the airflow, particularly in the channels upstream and downstream of the pallets. However, turbulent airflow was found in simulations to predict a recirculation zone in the first pallets at the refrigeration unit e nd. This could lead to a zone that cools down slower in these pallets (Defraeye et al., 2022) . The presence of this recirculation zone was not found in experiments explicitly yet. The current mismatch between simulated and experimental flow and thermal fields must be clarified first. We also ha ve a slight inconsistency when using the porous medium approach to calculate flow in the pallets of fruit is only valid under certain conditions (Defraeye et al., 2022) . The representative elementary volume (REV, computational cell) should be larger than the pore size or fruit diameter when usi ng the porous medium approach. This is not the case for large fruit in pallets in refrigerated container s. The alternative is explicitly modeling the thousands of fruit and hundreds of ventilated packages with vent holes . This scenario is computationally v ery expensive. One of our previous works required 40 million cells (W. Wu et al., 2019) to simulate all orange fruit in a single pallet. Various other physics -based model extensions include: (1) modeling the container in 3D, which is likely to affect the cooling at the sides of the pallets at the container walls predominantly ; (2) exploring the impact of buoyancy on the airflow and heat transfer in the cargo, which is only relevant at the start of cooling, since then the temperature differences between the fruit s urfaces and the air are the largest; (3) including the fan and other components, such as the evaporator heat exchanger, in the return air duct in the refrigerated container (Senguttuvan et al., 2021; Tiamiyu, 2020) ; (4) composing a library of porous media properties of different fruit packed in ventilated packaging, since the packaging type and loading density will affect the pressure resistance of the pa llets, and thereby also the airflow field and heat transfer, to some extent. We can extend the fruit quality predictions to include more relevant quality attributes . Examples are chilling injury incidence, pest mortality, mass loss of the fruit, occurrence of condensation , and the associated risk of microbiological decay , and for other fruit firmness and total soluble solids and acidity, for example . Note that our models do not aim to achieve very high absolute accuracy for fruit quality modeling . The model s are calibrated, but inherently we have variability in the models as the fruit properties vary with harvest year and location due to differences in breeding and growing conditions (Jedermann et al., 2017) . Another source of variability is that even fruit from the same orchard that are packed in the same container can exhibit differences in initial quality due to biological variability. The resulting induced i naccuracies lead to physics -based models not being readily adopted in commercial operations. Now, only a relative qualitative comparison of different shipments is possible. However, when stakeholders appreciate the value of this relative ranking, the results can still be very valuable for decision -making. The decision remains still at the discretion of the stakeholder, but they have more indicators to base their decision on, compared to, for example, a time -temperature curve only. Acknowledgments This work was supported by the Swiss National Science Foundation SNSF [200021_169372] and [200020_200629]. The funder was not involved in the study design, collection, analysis, interpretation of data, the writing of this article , or the decision to submit it for pu blication. This manuscript has been released as a pre- print . We acknowledge Joaquin Gajardo and Jörg Schemminger for creating or supporting the creation of some figures. We w ant to thank Richard Lawton from Cambridge Refrigeration Technology for providing us with detailed information on the specifications of refrigerated containers. Author contributions T.D. conceptualized the study and acquired funding; T.D. did the project admi nistration; T.D. performed the investigation, developed the methodology, and executed the simulations with key input of J.G. and C.V. ; T.D. developed the paper concept and simulations by critical discussions with C.V., J.G., L.L., T.B. a nd P.C.; T.D. wrote Preprint 49 the original draft of the paper; T.D. did the visualization; C.V., J.G., L.L., T.B. and P.C. performed critical review and editing. Declaration of g enerative AI and AI -assisted technologies in the writing process During the preparation of this work , the a uthor s used Grammarly (full paper) and ChatGPT (only section 2.3.7) in order to improve the spelling, grammar, and style of the text. No additional original content was generated using these AI -assisted technologies . After usin g this tool/service, the author s reviewed and edited the content as needed and take full responsibility for the content of the publication. Preprint 50