Microscale Modelling of Flow, Heat and Mass Transport During Vacuum Cooling of Porous Foods: Effective Property Computation

A microscale modelling framework to compute effective properties related to flow, heat and mass transfer during vacuum cooling (VC) of porous foods was developed. A heterogeneous computational domain reconstructed from steamed bread (SB) was used for modelling, while the cellular water transport in SB investigated using nuclear magnetic resonance was applied for model validation. The computed porosity (63.47 ± 1.05%), effective permeability (1.91 ± 0.39 × 10−11 m2), effective thermal conductivity, (0.33 ± 0.08 W m−1 K−1), and effective diffusivity (5.56 ×  ± 0.24 10−8 m2 s−1) were in the same range as those measured from the experiment/literature. Also, the analysis revealed that microstructural variability significantly affected the estimated effective properties. The microscale model yielded results similar to the lumped formulation but provided details not visible in the latter. Therefore, the developed model provides a framework for multiscale modelling, which could lead to a better understanding of the underlying moisture loss mechanisms during VC.


Introduction
Bread is a cereal-based hygroscopic food product with cellular structures and sponge-like systems having interconnected gas cells/pores heterogeneously distributed in the solid matrix. These heterogeneous microstructural features strongly affect bread quality attributes such as texture, sensory, shelf life, and rheological properties (Luo et al. 2018;Li et al. 2020;Tian et al. 2020a). Moreover, the cellularly diffusive transport of liquid water in the solid phase and vapour transport in the gaseous phase depend on the morphological details of the sponge-like cellular structure (Aguilera 2005;Li et al. 2018;Alabi et al. 2020). Hence, evaluation of the coupled transport mechanism in both phases is important to understand and interpret the connection between the moisture transport and structural morphology of bread microstructure (Esveld et al. 2012;Zhan et al. 2018 ;Mahato et al. 2019). On the other hand, food processing also affects the microstructural stability, the relevant transport properties and the quality of the processed food (Aguilera 2005;Ajani et al. 2022a;Zhou et al. 2017); therefore, deconstructing the process-structure-function relationship is equally important. (Tian et al. 2020b(Tian et al. , c, d, 2021 Among these food processing methods, cooling is an important method used in removing heat from food products with vacuum cooling (VC) being a viable approach used for the rapid cooling of porous food products with large surface area and high moisture content Sun and Zheng 2006;Zhu et al. 2018Zhu et al. , 2020. The principle of VC entails the rapid evaporation and removal of moisture from within and on the surface of porous food products under a vacuum. Once the prevailing pressure falls below the saturated vapour pressure during VC, water evaporates by absorbing latent heat from the porous food products, thereby resulting in rapid cooling of the food products Sun and Zheng 2006;Zhu et al. 2018Zhu et al. , 2020. Among the porous food products, bread is a hygroscopic and porous staple food widely consumed around the world, in particular, steamed bread is a common Chinese staple with high porosity and initial moisture content (Kou et al. 2019;Gao et al. 2020), making it an ideal product for VC. However, despite the merits of VC in prolonging the shelf life and keeping the freshness of porous foods products Sun and Zheng 2006), high mass loss during VC is the major drawback of the technology (Wang and Sun 2006), which is a subject of both experimental and numerical studies Zhan et al. 2019aZhan et al. , 2019b.
Unlike experimental studies, numerical modelling is often used to better understand transport phenomena during food processing, and macroscopic and microscopic approaches are the two primary ways of modelling these food processes (Ajani et al. 2020). At the macroscale, several numerical models have been developed to investigate the VC process. These include analysing the VC systems (Wang and Sun 2002a), investigating the effects of VC on temperature distribution and moisture loss in meat (Wang and Sun 2002b, Sun and Wang 2006, evaluating the effects of VC on temperature distribution, moisture loss and bacteria infiltration in leaves (Song et al. 2016;Ranjbaran and Datta 2019;Huang et al. 2021). However, these macroscale models do not take into account the complexity of the cellular structures. Besides, since these lumped formulations are limited in describing the behaviours of moisture variation in the food materials and other state variables, it is important to investigate other possible transfer mechanisms that may occur during VC from a microscopic perspective. On the other hand, the effective properties obtained from the microscale numerical model using actual heterogeneous microstructure can be directly employed for lumped formulation at the macroscale. Moreover, studies have revealed that effective properties estimated from the analytical approach are less accurate since the models are based on either simplified assumptions or on experiments that might be ridden with errors (Zhou and Yang 2019;Esua et al. 2022). In addition, the variations in materials used and/or processing conditions are often excluded from the resulting effective properties obtained using the analytical approach (Zhou and Yang 2019;Cheng et al. 2021;Petrosino et al. 2022). Although microscale models have been developed for studying drying (Rahman et al. 2018), water sorption (Esveld et al. 2012) and dehydration during storage (Workneh et al. 2013), such models are not yet developed for the VC process.
Hence, the current study aimed at combining microscale flow, heat and mass transfer models to compute effective properties using a realistic two-dimensional (2D) cellular steamed bread (SB) structure during VC. Although the connectivity of a 3D microstructure is greater than its 2D counterpart, the latter is computationally economical for complex microscale models and easier to investigate. Moreover, the symmetry of the SB experimental samples also validates the use of a 2D model. In developing the 2D model, the liquid water sorption behaviours and vapour transport in the solid phase (solid matrix) and vapour phase (pores) were formulated, and the energy exchange in the cellular structure of SB was represented using a representative energy equation. In addition, the microscale model was validated by comparing effective properties derived from the model with those obtained experimentally using nuclear magnetic resonance (NMR) and from the literature. It is hoped that results from the current study should provide a better understanding of cellular transport mechanisms, and an accurate estimation of effective properties needed for multiscale modelling, which can be useful for guiding further industrial applications of VC technology.

Model Description and Assumptions
Two phases, namely the liquid water-filled solid phase (solid matrix) and the vapour phase (pores), were distinguished in the microscale model ( Fig. 1a) (Wang et al. 2011;Joardder et al. 2015;Ajani et al. 2022a). Also, the porous medium was assumed to be unsaturated (Voogt et al. 2011;Workneh et al. 2013). The VC process consisted of a depressurization/ pressure reduction stage denoted as stage 1 and a pressure holding or constant vacuum stage designated as stage 2 as shown in Fig. 1(b) (Ranjbaran and Datta 2019;Huang et al. 2021;Ajani et al. 2022b).
Unlike plant/fruit microstructures with turgor pressure (Workneh et al. 2013), porous bread structures are predominantly open pores with no turgor pressure; hence, the pressure gradient-induced flow field was applied only in the pores (Esveld et al. 2012;Wang et al. 2011). For mass transport, the water inside the gluten-starch solid matrix was assumed to be driven by diffusion. On the other hand, mass transport in the pores was driven by hydrodynamic pressure gradients and diffusion. A local equilibrium between the wet surfaces of the solid matrix (solid matrix wall) and vapour in the pores was assumed (Fig. 1b). Since NMR studies have shown that about 3% of water in the pores makes up the entire water in the SB microstructure and even some pores have been found to be completely filled with the gaseous phase (Rahman et al. 2018;Ajani et al. 2022b), as such, liquid water in the pores was ignored and the assumption of no accumulation between the solid matrix/ pore interface was made (Voogt et al. 2011;Workneh et al. 2013;Rahman et al. 2018). For energy conservation, a local equilibrium assumption was considered. Besides, the combination of local equilibrium between the two phases for mass and energy conservation enabled the connection between both gradients and the simplification of the microscale Fig. 1 a Coupled transport model description, b pressure history in the vacuum chamber during VC stages, and c flow chart of microscale to effective property model, d mesh of first microstructure, and e mesh of second microstructure models to effective models (Buckley 2015;Randrianalisoa et al. 2015). The information needed from the microscale model was achieved by imposing the measurable boundary conditions/variable drop obtained from the macroscale model developed using similar VC conditions (Ajani et al. 2022b), and computing the corresponding variable flux (Voogt et al. 2011;Workneh et al. 2013;Wagner et al. 2021). Unlike previous models on gas exchange, storage and dehydration of food products (Ho et al. 2009;Voogt et al. 2011;Workneh et al. 2013), a transient study was employed in the current work due to the time-varying changes in the variables over time during VC. Additionally, previous food models have been limited to estimating the effective properties related to the mass transport such as effective diffusivity (Voogt et al. 2011;Workneh et al. 2013). As such, the porosity associated with the porous food microstructure, the effective permeability related to the flow field, the effective thermal conductivity linked to the heat transport and the moisture diffusivity analogous to the mass transport during the VC process were also estimated from the current microscale model formulation. It is worth noting that except for the effects of microstructural variabilities on the effective properties, specific parametric study to evaluate geometrical properties such as pore connectivity was not considered in the current work.

Flow in the Vapour Phase (Pores)
A momentum conservation equation was used to model the hydrodynamic pressure gradient induced in the bread pores. The Stokes equation is appropriate for modelling the momentum in bread pores in which viscous flow is dominant (Wagner et al. 2021). Hence, the Stokes equation was employed to describe the flow in the pores and subsequently homogenized to compute the effective permeability (Sect. 2.5.2). The time-dependent Stokes equation was defined as: where , , and P are the density (kg m −3 ), dynamic viscosity (Pa s), velocity (m s −1 ) and pressure (Pa) of the fluid, respectively, and ∇ 2 is the Laplacian operator. The system of Eqs.
(1-2) was supplemented with time-varying pressure drop by applying the boundary conditions presented in Table 2 at the top and bottom with inflow and outflow pressure, respectively. Symmetry/slip boundary condition (n·u = 0) was applied on the other two sides (Petrosino et al. 2022).

Energy Conservation in Solid Matrix
The following is the energy equation for the solid phase: where C ps , T s and k s represent the heat capacity (J/kg K), temperature (K) and thermal conductivity of the solid matrix (W m −1 K −1 ), respectively. The term on the right-hand side of Eq. (3) is the conductive heat flow in the solid matrix.

Energy Conservation in the Pores
The energy equation for the pores is given by: where p is the effective density in the pores (kg m −3 ), C pp is the effective specific heat capacity in the pores (J/kg K), is the gas velocity (m/s), k p is the conductivity in the pores (W m −1 K −1 ), and T p is the temperature (K) of the pores. On the right-hand side of Eq. (4), the first and second terms are the conductive and convective heat flow, respectively.
The heat flux (W m −2 ) was related to the thermal conductivities in each phase and the temperature gradient, and due to the complex microstructure of porous food materials, the numerical approach was employed since it was straightforward (Randrianalisoa et al. 2015). Hence, to compute the effective thermal conductivity, a single representative energy equation was solved over the entire representative sample while including relevant domain properties of the respective phases and equilibrium assumed at the interphase. The temperature drop was applied at the lower and upper boundaries while an adiabatic-type boundary condition was applied at the other two edges (Randrianalisoa et al. 2015; Zhou and Yang 2019).

Liquid Water Transport in Solid Matrix
The transport of bound water in the solid matrix was modelled using diffusion law (Bird et al. 2007;Workneh et al. 2013). The liquid water diffusive flux through the solid matrix was defined using Fick's first law (Bird et al. 2007): where J ws is the water flux in the solid matrix (kg m −2 s −2 ), D s is the water diffusion coefficients in the solid phase (m 2 s −1 ), and ws is the solid phase density on wet basis (kg m −3 ). Also, ws was further expressed as: where ds , X s and a w are the solid matrix density (kg m −3 ) on dry basis, moisture content on dry basis and water activity, respectively. Also, X s a w is derivative of the moisture content on dry basis as a function of the water activity corresponding to the sorption isotherm. Using Eqs. (5) and (6), the unsteady diffusion in the solid matrix resulted in:

Gaseous Transport in Vapour Phase
The vapour diffusive flux through the vapour phase was defined as (Bird et al. 2007): where J v is the diffusive vapour flux in the solid matrix (kg m −2 s −2 ), v is the density of the vapour (kg m −3 ), and D v is the diffusion coefficient of water vapour in the air (m 2 s −1 ). In addition, v was further expressed based on liquid water activity as (Voogt et al. 2011): Here, vsat and a w represent the density of the saturated vapour (kg m −3 ) and water activity, respectively. Given the diffusive vapour flux (Eqs. (8) and (9)), and in addition to the convective vapour flux, the transient model of the gaseous transport within the pores was defined as (Bird et al. 2007): where is the velocity in the vapour phase (m s −1 ) and the second term in Eq. (10) denotes the vapour convective flux.
For the mass transport, the time-varying boundary conditions were imposed over the bottom and top boundaries of the domain using Dirichlet boundary conditions, no flux boundary conditions were applied at the lateral sides, while local equilibrium boundary condition ( v = ws ) was applied at the interphase between the vapour phase and the solid matrix (Voogt et al. 2011;Prawiranto et al. 2018).

Initial Conditions
The initial condition for the bread temperature was 77 °C and the initial velocity in the pores was 0 m/s. Also, the initial moisture contents in the solid matrix and pores were estimated from the NMR experiment, and the values including other relevant parameters are presented in Table 1. Specifically, the effective values of parameters such as the specific 1 3 heat capacity and thermal conductivity were estimated using mass fractions of SB compositions (Huang and Miskelly 2016) and relevant relationships (Sahin and Sumnu 2006).

Solid Matrix Water Activity
The Ferro Fontan model was used for the sorption isotherm in the solid matrix during VC and expressed as (Besbes et al. 2013): where a w is the water activity of the liquid water bound to the solid matrix and the adjustable isotherm parameters were obtained from Besbes et al. (2013).

Equilibrium Relationship
The relationship between the water vapour in the vapour phase and liquid water in the solid phase liquid phase was expressed as: where a w is the water activity, M w is the molar mass of water (kg mol −1 ), R is the gas constant (J mol −1 K −1 ), T is the temperature (K) and P sat v (T) is the saturation vapour pressure (Pa) at a reference temperature defined as (Sun and Hu 2003):

Macroscopic Property Estimation of Bread Microstructure
Generally, obtaining dynamic macroscopic properties/information from detailed structures without using empirical formulations could be achieved from the microscale model either by imposing the measurable variable flux (e.g. velocity) over the microstructure and computing the corresponding variable gradient (e.g. pressure drop) or setting the variable gradient and calculating the relevant flux change (Petrosino et al. 2022). Hence, information for computing the effective properties from the microscale model was calculated for the effective permeability (κ) (m 2 ), effective thermal conductivity (W m −2 ), and effective diffusivity (m s −1 ), by applying the relevant variable drop/boundary conditions and measuring the corresponding flux.

Porosity
The porosity of the microstructures was calculated as the fraction of the surface area of the pore space to the total surface area of the entire microstructure: where θ is the porosity and A p and A T are the area of the pore space (m 2 ) and the total area of the microstructure (m 2 ), respectively.

Effective Permeability
Porous media's permeability is its resistance to fluid flow when pressure gradients are applied over the porous media (Chaunier et al. 2008). Hence, permeability is important in studying pressure-driven mass transport in porous media such as during vacuum cooling of porous foods (Song et al. 2016;Ranjbaran and Datta 2019;Ajani et al. 2022b). Effective permeability (κ) in m 2 is an effective parameter related to the intrinsic property of the microstructures used to quantify and characterize fluid flow through porous media at the macroscale (Wagner et al. 2021). Hitherto, relevant information on the permeability of different porous foods is inadequate (Chaunier et al. 2008). Moreover, as input parameters for macroscale modelling of porous foods during vacuum cooling, estimation of permeability has been limited to complex empirical formulations or experiments (Wang and Sun 2002b;Wang and Sun 2006;Song et al. 2016;Ranjbaran and Datta 2019;Huang et al. 2021;Ajani et al. 2022b). Hence, a microscale numerical technique was employed in the current study to compute the effective permeability.
Since Darcy's law is valid at the macroscale (Darcy scale), the bulk flow through the entire porous region was defined using Darcy's law. The permeability of the bread porous region for each microstructure was estimated using volume-averaged quantities as follows: where κ, Q, , L, A, and ∆P are the effective permeability (intrinsic) in the flow direction (m 2 ), mass flow rate (Kg s −1 ), dynamic viscosity (Pa s), thickness of the microstructure (m), density (kg m −3 ), cross-sectional flow area (m 2 ) and pressure drop (Pa), respectively.

Effective Thermal Conductivity
The numerical approach allows for providing specific results based on different parameters and numerical methods such as the finite element method (FEM) are more appropriate for irregular geometries having spatially variable conditions and properties because the discretization process is not limited to a regular node positioning (Zhou and Yang 2019). Hence, FEM was used to model the temperature field of the SB microstructures which included the respective thermal conductivity of the solid matrix and the pores. The temperature drop was applied to obtain the heat flux Q T and the effective thermal conductivity of the microstructure was computed thus: where k ef ,,n , Q T , L and ΔT are the numerical effective thermal conductivity (W·m −1 ·K −1 ), heat flux (W·m −2 ), distance along the y-axis (m), and the temperature difference (K), respectively.

Effective Diffusivity
Hitherto, effective diffusivities employed in most macroscale modelling of food processes involve the use of simple empirical relationship derived from curve fitting of experimental data. However, perfect geometry is often considered in the curve fitting techniques which could lead to a discrepancy between experimental data and predictive model (Esveld et al. 2012;Workneh et al 2013;Zhou and Yang 2019;Zhang et al. 2018). Conversely, the macroscale properties estimation using the microscale numerical approach considers the material's heterogeneity in its formulation.
To evaluate the equivalent macroscopic gradient over the SB microstructure, the magnitude of the parallel diffusive fluxes can be compared. The overall flux of water over cellular solid food material such as SB is dominated by the flux in the vapour phase because the water migration distance in the vapour phase (pores) is large compared with the cell wall thickness of the solid matrix (Voogt et al. 2011;Esveld et al. 2012). Additionally, vapour diffusivity in the vapour phase is several orders of magnitude larger than the liquid water diffusivity in the solid phase (Voogt et al. 2011;Esveld et al. 2012). Hence, the macroscale moisture transport in the bread matrix was assumed to be dominated by the vapour phase (Bird et al. 2007;Voogt et al. 2011;Esveld et al. 2012) and can be represented by a single effective diffusion coefficient using Fick's law. According to local mass balance, the volume-averaged diffusive flux is equal to the water locally accumulated in the condensed and gaseous phase (Bird et al. 2007): where θ is the porosity, a w is the volume-averaged water activity, X s is the moisture content on dry basis, vsat and ds represent the density of the saturated vapour (kg m −3 ) and solid matrix density (kg m −3 ) on dry basis, respectively, while J is the volume-averaged flux (kg m −2 s −1 ) and t is the time (s).
Assuming a local equilibrium of the solid and vapour phase and neglecting the small mass accumulation in the vapour phase, the effective diffusivity Eq. (17) can be expressed as (Bird et al. 2007;Voogt et al. 2011;Esveld et al. 2012): D eff is the effective diffusivity (m 2 s −1 ), is the porosity, X s is the moisture content on dry basis, ds represents solid matrix density (kg m −3 ) on dry basis, and J is the volume average flux (kg m −2 s −1 ) in the vapour phase. Since the transient study was considered in the current work, the volume average fluxes were computed at various temporal instants during the vacuum cooling process. The time-varying volume average flux J (kg m −2 s −1 ) was calculated by integrating the water flux over the surface of the top boundary. The effective moisture diffusivity ( D eff ) in Eq. (18) only contains terms related to the moisture diffusivity in the vapour phase since the mass transport in this phase was the main contributor to the macroscopic transport and the equation is valid based on the assumption of equilibrium between the liquid water and vapour (Voogt et al. 2011;Esveld et al. 2012).

Model Implementation
In the current work, a representative element volume, which was reconstructed in our previous study (Ajani et al 2022a), was used as the computational domain for the microscale modelling. Details about the geometrical development are also contained therein (Ajani et al. 2022a). Three different geometries with different microstructural distributions were employed to investigate the microstructural variability (Ajani et al. 2022a). The SB microstructures were used for the transport modelling and imported from MATLAB (Version 2019b, MathWorks Inc., Natick, USA) into COMSOL Multiphysics 5.6 (COMSOL AB, Stockholm, Sweden) for the numerical simulation. Furthermore, the input parameters and initial and boundary conditions were applied to the model. Moreover, for the governing equations, the momentum was solved using the Creeping Flow module, while the heat and mass transport was solved using the Heat Transfer and Transport of Diluted Species modules, respectively (Fig. 1c). Unlike previous models that arbitrarily assumed relevant variable potentials for the spatial variation in microstructure, and accurately apply the macroscale VC conditions on the microscale model for realistic property estimation, the boundary conditions applied across the microstructures were derived from the macroscale model under similar VC conditions (Ajani et al. 2022b). In that way, an accurate and robust evaluation of the spatiotemporal variation can be considered to calculate the effective properties. The relevant boundary conditions at corresponding temporal instants were applied across the geometry to solve the transient model at each time step. The variable drops/boundary conditions are detailed in Table 2. The numerical approach allows for providing specific results based on different parameters and numerical methods such as the finite element method (FEM), are more appropriate for irregular geometries having spatially variable conditions and properties because the discretization process is not limited to a regular node positioning (Zhou and Yang 2019). Hence, FEM was used to model the temperature, vapour pressure and water distribution of the SB microstructures and calculate the corresponding effective properties. The effective properties including porosity, effective permeability, effective thermal  (18), respectively. For each microstructure, grid dependency tests were conducted using five physics-controlled mesh sizes namely coarse, normal, fine, extra fine, and extremely fine. The sufficiency of the mesh sizes was investigated using the temperature and it was observed that percentage errors between the temperature variation of the three largest elements were less than 1%; hence, the mesh sizes tagged "fine", with mesh sizes of 33,412 and 58,302 for M1 and M2, respectively, were used for the simulation to reduce computational space and time. Furthermore, a fully coupled Parallel Direct Sparse Solver (PARDISO) with an automatic nonlinear solver was used for the simulation. Besides, for the whole system of equations, the absolute and relative values were set to 0.001, and the formulated model was simulated on a 64 GB RAM and 3.6 GHz Core i9 PC (DESKTOP-V0OSV1R, Lenovo Group Ltd., Beijing, China).

Sample Preparation
Fresh SB samples were obtained from a local bakery shop at South China University of Technology (Panyu District, Guangzhou, China). The steamed bread was the Guangdong type and its formulation consists of wheat flour (100 g) (Feng Shun Flour, Feng Shun Flour Co. Ltd., Foshan, Guangdong, China) having a composition of protein (8.6 g), fat (0 g), carbohydrate (77.3 g), and sodium (0.012 g), and with sugar (5 g), baking powder (1 g), dry yeast (2 g) and water (80 mL). Moreover, the dough was proofed for about 1 h and a one-stage fermentation process was used for the steamed bread production. Furthermore, to apply similar conditions as in industrial practice, the temperature of freshly steamed bread samples was maintained at 85 °C before vacuum cooling. Due to the small size of the vacuum cooling system for the in situ investigation, the steamed bread with a rounded shape on the top, and a square base of about 2 cm 2 was used.

Experimental Setup
The experimental setup consisted of a small vacuum chamber, a vacuum pumping system (VRI-2, Zhejiang Value Mechanical & Electrical Products Co., Ltd., Zhejiang, China), and a data acquisition system. The data acquisition module (DAM-3232, Beijing Art Technology Development Co., Ltd., Beijing, China) was connected to a computer and a pressure transmitter (PT124B-212, Shanghai Zhaohui Pressure Apparatus Co., Ltd., Shanghai, China) with 100-10 5 Pa pressure range, 4-20 mA output current mode and 0.05% accuracy. A detailed schematic diagram of the experimental setup can be found in our previous study (Ajani et al. 2022c).

Experimental Procedure
During the experiments, an exponential decay function was used to model the pressure reduction in the vacuum chamber from the atmospheric pressure until the final value of approximately 600 Pa and is expressed as  where P vac is the instantaneous pressure (Pa) in the vacuum chamber, P i is the initial pressure (Pa), t is the vacuum cooling time (s), and (s −1 ) denotes the coefficient that characterizes the pressure reduction rate and the pressure profile over time in the chamber is depicted in Fig. 1b. The SB sample was placed in the miniature vacuum chamber, which was then inserted into the NMR probe for the dynamic investigation. Furthermore, Type-T thermocouples (Omega Engineering, Inc. CT, USA) were employed to study temperature profiles during the VC process. Moreover, the sample was weighed before and after the VC process, and the mean values of at least three sample replicates were recorded and used for analysis.

NMR Measurement
The NMR measurements were performed with a low-field nuclear magnetic resonance (LF-NMR) analyser (NM42-040H-I, Suzhou Niumag Analytical Instrument Corp., Shanghai, China) equipped with a 1 T permanent magnet at 32 °C. The Carr-Purcell-Meiboom-Gill (CPMG) sequence was used to assess the transverse relaxation times (T 2 ) with the following critical parameters: echo time (TE) of 0.25, waiting time (TW) of 1000 ms, frequency (SW) of 200 kHz, echo number (NECH) of 5000, and 2D slices with 2 mm thickness each. Multi-exponential fitting analysis of the relaxation data was performed using MultiExp Inv Analysis software (Version 3.0, Suzhou Niumag Analytical Instrument Corp., Shanghai, China).

Moisture Content Measurements by Gravimetric Method
The bread samples' dry weight was determined using an oven method at 105 °C with a halogen moisture meter (XY-100 MW-T, Shanghai Jiezhun Instrument Co., Ltd., Shanghai, China). The total MC of SB on a dry basis during the VC process was determined using a gravimetric approach with the following expression: where W 0 , W i and W d are the initial, instantaneous and dried weights of SB, respectively.

Moisture Content Measurements by NMR-Based Method
On a dry basis, the MC of free, loosely, and tightly bound water in the SB was estimated using its mass (g) and amplitude of the T 2 spectrum, assuming that only water was present and there was a linear relationship between the amplitude and water mass (Enninful et al. 2017). To estimate water samples, the bulk water parameters were first assessed and the water amplitude index was determined as follows: The NMR water mass (WM) was the water mass of the sample at each VC time and was calculated as (Enninful et al. 2017): (21) Amplitude Index(AI) = NMR amplitude∕unit mass of liquid 1 3 The corresponding MC was quantified using the following equation: where MC iNMR , Area i , ΣArea , and OW are the free, bound or total NMR-based MC, the amplitude area of free, bound or total water, the total amplitude area, and the oven-dry weight of the sample, respectively.

Effective Diffusivity Measurements
Effective diffusivity is the overall diffusion coefficient which describes the effective diffusion of moisture transport through the entire porous structure of the cellular food. The effective moisture diffusivity of the SB was calculated using the time-varying moisture contents on a dry basis ( MC iNMR ) obtained from the vacuum cooling experiment. The analytical equation for the evolution of moisture contents with vacuum cooling time of an infinite slab was given by Besbes et al. (2013): where MR NMR is the moisture ratio which is related to the MC iNMR (g/g), n is the number of terms of the Fourier series, D efe is the experimental effective moisture diffusivity (m 2 /s), t is the time (s), and, L is the thickness (m). Equation (24) can be simplified as: D ef ,e was related to the slope obtained from a plot of ln MR with vacuum cooling time.

Statistical Analysis and Validation
The relevant figures were plotted, the relationship between LF-NMR parameters and total moisture content was established, and validation of the developed model with experimental results was investigated using Origin software (Version 9.0, OriginLab Corp., Massachusetts, USA). The deviation was used to evaluate the numerical model accuracy and the accuracy was assessed by comparing the numerical effective properties with those estimated from experiments. The deviation was calculated from the following equation: where D ev is the deviation, E p,n is the relevant effective properties from the numerical approach, while E p,a is the effective properties obtained analytically or experimentally.
(22) WM = Amplitude of fluid in sample/Amplitude Index

Distribution and Variation of Pressure-Driven Flow
Although three microstructures were used for analysis, surface plots of the first microstructure denoted as M1 and the second microstructure designated as M2 are presented and their respective geometry meshes used for modelling are shown in Fig. 1d and e. Also, the spatiotemporal variation of the pressure drops (Pa) and corresponding flux/velocity (m/s) across the porous region of the two microstructures are displayed in Fig. 2. Generally, the spatial distribution of the pressure developed across the microstructures for the different VC times revealed that the pressure in the bread core was considerably higher than that of the surface as shown in Fig. 2a and b. The pressure gradients mostly exist from one pore to another along the direction of the applied gradients and bread pores at the same position in the gradient have almost uniform pressure drops. Similar variations in the gradients of relevant variables have been observed during microscale modelling of the gas exchange of fruit during storage (Ho et al. 2009) and dehydration (Workneh et al. 2013). The sharp gradient developed across the microstructures was due to their permeabilities which led to high resistance to the vapour flow field (Huang et al. 2021). As a result of the hydrodynamic pressure variations during vacuum cooling, the pressure drop across the microstructures varies over time and space. Hence, during the initial phase of stage 1 (Fig. 1b), the pressure gradients across the microstructures were minimal. As VC progressed, the pressure drop became most prominent at about 100 s of VC time ( Fig. 2a and b), thus implying that the vapour pressure gradients from the bread surface to the core were largest at this temporal instant. The high gradient developed across the microstructure just before the beginning of stage 2, could be attributed to a rapid drop in the prevailing vapour pressure below the equilibrium vapour pressure (flashpoint), which increased the evaporation and triggered the removal of more water vapour from the porous structure Sun and Zheng 2006;Ranjbaran and Datta 2019). However, during stage 2, the vapour pressure gradients at 300 s and 600 s reduced significantly as shown in Fig. 2a and b. This reduction was due to the prevailing vapour pressure becoming closer to the equilibrium vapour pressure during the constant vacuum pressure stage (Fig. 1b), thereby reducing the evaporation. Similar findings based on the macroscale VC models were also reported Sun and Zheng 2006;Ranjbaran and Datta 2019). Additionally, it was observed that the pressure gradient across the two microstructures was directly proportional to the flux/ velocity of the fluid at various temporal instants. Specifically, the maximum velocity across both microstructures ( Fig. 2c and d) was developed at 100 s when the pressure gradient was most prominent ( Fig. 2a and b), and similar relationships were observed at 300 and 600 s. What's more, the velocity was higher where the flow path was narrow as shown in Fig. 2c and d. Comparatively, although the same gradient was applied on both microstructures, the pressure drop at the same temporal instant varies. For example, at 100 s, the pressure drops for M1 and M2 were 25.1 Pa and 23.2 Pa, respectively. The microstructural variabilities including differences in permeability and porosity, between the two microstructures could be responsible for the variations in the pressure drop observed (Curcio et al. 2018). Also, the average time-varying fluxes (m/s) developed at the outlet of M1 and M2 are presented in Fig. 2e and the fluxes had a linear relationship with the pressure gradients. Additionally, the results revealed that the fluxes obtained from both microstructures were slightly different due to the variations in their respective microstructural distributions (Ho et al. 2009;Workneh et al. 2013). Furthermore, at 100 s, the highest difference between the total flux of M1 (6.6 × 10 -3 m/s) and M2 (2.8 × 10 -3 m/s) was observed as shown in Fig. 2e. This variation in the total flux between M1 and M2 could be particularly attributed to the differences in their respective porosities as further elaborated in Sect. 4.5.1. Fig. 2 a microscale pressure gradients in the vapour phase at different temporal instants of first microstructure, b microscale pressure gradients in the vapour phase at different temporal instants of second microstructure, c microscale flow field at different temporal instants of the first microstructure, d microscale flow field at different temporal instants of the second microstructure, and e macroscopic total flux at the outlet

Temperature Distribution and Variation
The temperature variation developed across the pores and solid matrices for M1 and M2 are presented in Fig. 3. The spatial variations of the temperature across the two microstructures revealed that the internal temperature in the bread core was higher than that at the surface and the temperature gradients along the same spatial planes were similar to each other as displayed in the surface plots in Fig. 3 (a) and (b). The gradient developed across the microstructures could be attributed to the internal heat generated at the bread core due to the evaporation induced during the VC process (Huang et al. 2021). On the other hand, the temporal variations revealed that temperature gradients for M1 and M2 were sharpest during the 100 s of VC implying that the internal heat generated in the bread core at this temporal instant was the highest (Sun and Wang 2006;Huang et al. 2021;Ajani et al. 2022b). As the VC progressed into stage 2, the temperature gradients significantly reduced across both microstructures ( Fig. 3a and b) due to the subsequent reduction in the amount of internal heat generated. Overall, it can be inferred the temperature distribution/gradient developed across the microstructure is directly proportional to the amount of evaporation and the ultimate moisture loss during the VC process. For example, the highest temperature gradients for M1 and M2 realized at 100 s implied that the moisture loss would be highest during 100 s of VC. The spatiotemporal temperature gradients observed in the current microscale model are similar to those of the macroscale counterpart except for the length scale of the former (Song et al. 2016;Ranjbaran and Datta 2019;Huang et al. 2021;Ajani et al. 2022b).
The volume-averaged heat flux developed at the outlet (surface) for both microstructures is presented in Fig. 3c, and the results showed a direct relationship with the temperature gradients at various temporal instants. Comparatively, the heat flux of the respective microstructures differed due to their microstructural variabilities and the highest heat flux difference was realized at 100 s, for both samples. Overall, it can be observed from the pressure drop and temperature gradients, that there is a similar spatiotemporal pattern that could be attributed to the relationship between the saturation vapour pressure and the temperature of the product (Eqs. 13).

Moisture distribution and variation
Since the macroscopic transport is dominated by the vapour phase, the predicted cellular water distributions developed during the VC for M1 and M2 in the vapour phase are presented in Fig. 4a and b, respectively (Esveld et al. 2012). From the figures, it can be observed that pores at the same position tend to have similar and uniform water density. Also, the spatial variations mainly exist from one pore to another along the direction of the applied gradient which could be attributed to the hydrodynamic pressure gradients developed across the microstructures (Ho et al. 2009;Ajani et al. 2022b). Additionally, the vapour density at different temporal instants differs as shown in Fig. 4a and b and this could be attributed to the differences in the gradients developed during the VC process . In the same vein, some variations in the vapour density distribution were observed between M1 and M2 and the differences in the microstructures might be responsible (Ranjbaran and Datta 2019;Huang et al. 2021;Ajani et al. 2022b). Also, the distribution of the water vapour density was similar to that of vapour pressure as discussed in Sect. 4.1, thus signifying the relationship between them. In addition, the moisture variation and removal increased with an increase in the temperature gradients at different 1 3 temporal instants. Specifically, the moisture removal was more prominent at 110 s which corresponds to the same temporal instant with the highest temperature gradients as discussed in Sect. 4.2. Similarly, moisture variation/removal was directly proportional to the vapour pressure gradients developed at a corresponding temporal instant during the VC process.

Experimental Moisture Variation
The dynamic MC during VC was evaluated using the T 2 inversion curves presented in Fig. 4c. Initially, four peaks of T 21 , T 22 , T 23 and T 24 with corresponding signal amplitudes of A 21 , A 22 , A 23 and A 24 were noticed in the inversion curve as shown in Fig. 4c. The different relaxation times represent the respective molecular environments of the water molecules. In addition, T 21 is the tightly bound water (TBW) attached to the protein structures, Fig. 4 a Water vapour density variation at different temporal instants of the first microstructure, and b water vapour density variation at different temporal instants of the second microstructure, c T 2 inversion curves at different times during the VC process, and d experimental moisture content variation on dry basis T 22 is the loosely bound water (LBW) attached to the polysaccharides, while T 23 and T 24 are the free water components found in the pores (Kou et al. 2019). The T 21 and T 22 components were less mobile and harder to be removed during the VC process than the T 23 and T 24 counterparts. This could be attributed to the binding of T 21 and T 22 components to the protein structures and polysaccharides in the bread, respectively (Zang et al. 2017;Kou et al. 2019). Furthermore, the higher relaxation times of the T 23 and T 24 components made them free and were easily removed as the VC progressed (Fig. 4c).
To calculate the MC from the experimental data, the integral of the NMR signal peaks (Fig. 4c) was correlated with the amount of water in the SB samples measured before and after the VC processes (Ajani et al. 2022a, b, c). What's more, the MC calculated using the NMR approach was in good agreement with that calculated using the gravimetric-based approach with an R 2 of 98.78%. Besides, the experimental MCs during the VC process are presented in Fig. 4d and the MC variation was employed in calculating the effective diffusivity.

Porosity
The mean porosity was computed from the different microstructures and the results revealed that their variability affected the porosity. For example, the computed porosity for M1 and M2 was 62.89% and 65.01%, respectively. Additionally, the mean porosity from the numerical approach (63.47 ± 1.05%) was compared with that obtained experimentally by Gao et al. (2015) (64.81%) and the deviation was 0.02 as shown in Table 3. The comparison showed that the numerical results were found to be in the same range as the porosity of steamed bread measured experimentally by Gao et al. (2015). Conversely, the reported porosity for baked bread in the same study was 71.56% (Gao et al. 2015), and much higher porosities for five sandwich bread, ranging from 72.48 to 84.00%, were reported in another study by Wang et al. (2011). The variation in the porosities among the different bread types could be attributed to the respective formulation and processing methods (Gao et al. 2015). Overall, it can be inferred that the current numerical approach was accurate and the slight variations could be due to the differences in the microstructural variability of the experimental and numerical studies and the microstructural connectivity associated with using a 2D geometry (Workneh et al. 2013). Lastly, the accuracy of the computed porosities reflects the sufficiency of microstructures as representative samples.

Effective Permeability
The permeability computed from the cumulative microstructures is presented in Table 3 and the mean permeability of the microstructures was estimated as 1.91 ± 0.39 × 10 −11 m 2 . The variation among microstructures could be attributed to the differences in their porosities as discussed in Sect. 4.5.1. Similarly, reports from microscale model development by Curcio et al. (2018) revealed that differences in porosity affected the estimated gradient/ flux. Besides, the computed permeability from the cumulative microstructures was compared with the experimental data obtained from the literature. Experimental work conducted by Wang et al. (2011), revealed that the different bread permeability ranged from 1.35 × 10 −10 to 2.35 × 10 −13 m 2 . In another study, Chaunier et al. (2008)   permeability ranging from 3.61 × 10 −10 to 4 × 10 −12 m 2 . This showed that the computed values were within the permeability ranges reported for bread experimentally; thus, the accuracy of the proposed microscale framework in estimating the effective permeability was confirmed. Additionally, this also confirms the adequacy of the microstructural sizes employed as representative element volume (REV) for evaluating the effective permeability (Esveld et al. 2012). Also, the average permeability reported by Wang et al. (2011) was compared with the current study and presented in Table 3. The results showed a deviation of 0.05, which could be attributed to variations in experimental and numerical properties. Overall, the effective permeability, which hitherto has been limited to rigorous experimental investigation and empirical relationships, can be estimated using the microscale framework developed in the current study.

Effective Thermal Conductivity
The effective thermal conductivity estimated from the microscale model is presented in Table 3. Compared with the literature, the mean thermal conductivity calculated was 0.33 ± 0.08 W m −1 K −1 and ranged from 0.26 to 0.43 W m −1 K −1 . Also, Zuniga and Le-Bail (2009) reported an effective thermal conductivity of bread ranging from about 0.12 -0.42 W m −1 K −1 . In another study, Sumnu et al. (2007) reported a much lower effective thermal conductivity for baked bread (0.15 -0.30 W m −1 K − ) and the low values might be attributed to the differences in higher porosities and lower moisture contents compared to that of the current study. The differences between both approaches could be attributed to relevant assumptions made in both approaches. Also, challenges such as difficulties in measurements and measuring instruments during experiments might be responsible for the variations (Randrianalisoa et al. 2015;Zhou and Wang 2019).

Effective Diffusivity
To compare the effective diffusivity of the mass transport model with that of the experiment, the mean of the time-averaged effective diffusivity for the respective microstructures was calculated and presented in Table 3. The computed effective diffusivity of the microstructures revealed some variations in the values across the different samples which could also be attributed to the varying microstructural distributions such as connectivity, pores size, and distribution. Comparatively, the mean effective diffusivity of the model (5.56 × ± 0.24 10 −8 m 2 s −1 ) was different from the experiment's (8.51 ± 0.45 × 10 -8 m 2 s −1 ) and the deviation of 0.34 was reported between them. Moreover, the differences between the experimental and numerical values were likely due to the relevant assumptions made and differences in experimental and numerical conditions. Although the numerical approach assumed a local equilibrium including the domination of the vapour transport for its computation (Voogt et al. 2011;Esveld et al. 2012), the experimental estimation of diffusivity was also based on certain simplifications such as the assumption of perfect geometries and the exclusion of some time-varying conditions, hence, the discrepancies between them.

Conclusions
For the first time, a microscale model was developed using reconstructed microstructures of SB to predict the cellular moisture changes and temperature distribution during VC of SB and calculate the effective properties such as permeability, thermal conductivity and diffusivity. In addition, numerical results were validated using the effective properties calculated from NMR with reasonable accuracy. The diffusion model shows that morphological parameters have a more profound impact on water migration than material properties. Overall, the ellipse tessellation-based transport model incorporated detailed heterogeneity with greater morphological details of the sample than the network-based model, while using relatively lesser computational resources compared with using a complex direct voxel-based simulation. Results indicated that the current model could provide a framework for a multiscale model to better understand the process-structure-function relationship during the VC process. The main advantage of the microscale numerical modelling approach is the possibility of computing effective properties for macroscale models for VC processes without using the traditional approach of rigorous experimental investigation of individual samples to formulate relevant empirical correlations. Lastly, 3-D models developed with actual microstructures could be used to improve the current model accuracy with an increase in computational resources.