Development of a validated 2 D model for flow , moisture and heat transport in a packed bed reactor using MRI experiment and a lab-scale reactor setup

• A submitted manuscript is the author's version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.


Introduction
Climate change and depletion of fossil fuel resources are among the main issues under investigation all over the world. In Europe, the building sector is the largest consumer of energy and accounts for approximately 35% the total primary energy consumption, of which 75% is used for space heating and domestic water heating [1]. Solar energy is one of the most promising sustainable energy sources for replacing fossil fuels. However, to reach high solar fractions, long-term storage of solar energy is necessary.
A storage option with a promising potential is thermochemical heat storage utilizing thermochemical materials (TCM), by which heat can be stored almost loss-free over a long time [2]. Heat is stored into an endothermal dissociation reaction, splitting the thermochemical material into two components (charging), and, at a later time, the energy can be retrieved from the reverse exothermal reaction between the two components (discharging) according to the reaction AðsÞ þ BðgÞ $ ABðsÞ þ heat. For low temperature thermochemical heat storage, adsorption of water vapor on sorption materials [3] and hydration of salt hydrates [4] are frequently studied. Heat generated by a solar collector during summer can be employed to desorb water from the material, and the energy stored in this way can be released during winter by introducing water vapor to the dehydrated material. In this study, Zeolite 13XBF [5] is used as sorbent. Zeolite is a good candidate to be used in scientific reactor studies because of its high stability [6]. The adsorption kinetics is frequently modeled by the Linear Driving Force (LDF) model [7] [8] [9], because it is simple, analytical, and physically consistent [10].
In the literature, both open and closed systems are investigated for long-term thermal storage of solar energy [11]. In an open system, both sorbate and energy are exchanged between the system and the environment, while in a closed system, only energy is exchanged between the system and the surrounding environment. In this work, an open system is considered, because the open system concept seems more promising for seasonal heat storage because of robustness and low cost [12]. An important part of the thermochemical heat storage system is the reactor, where the charging and discharging of the thermochemical material take place. In this work, a packed bed reactor design was chosen because of low need of auxiliary energy in comparison with other types of reactors such as fluidized bed or screw reactors [13]. A disadvantage of the packed bed reactor concept is the risk of a nonuniform flow leading to a lower energy storage density. This can be avoided by consideration of specific measures in the design of the reactor, which require an in-depth understanding of the physical processes inside the reactor which can be obtained by developing a validated numerical model.
For better understanding of the influence of transport phenomena on the performance of the thermochemical heat storage reactor, a detailed model, which considers radial terms as well as axial terms, is needed. It is concluded from previous work [14] that the reactor wall plays an important role in the radial effects occurring in the reactor, and accordingly in the thermal performance of the reactor. This wall effect was illustrated by a strong radial temperature profile in the reactor. This may be caused by a combination of three different effects: Heat loss to the ambient causes a lower temperature near the wall, which is directly correlated to the heat transfer coefficient at the inside wall of the reactor and the thermal mass of the reactor wall. In the region near the wall, a higher bed void fraction is observed in the bed, because of the interaction of the particles with the rigid wall, since the spatial distribution of the particles must conform to the shape of the wall. The higher bed void fraction near the wall causes a higher local superficial velocity (called wall channelling), hence the produced heat by the reaction is depleted faster in this region. In addition, the higher bed void fraction near the wall results in less material near the wall, thus less reaction heat and a lower temperature as a result. Since the state of charge is a function of temperature and partial vapor pressure, the lower temperature near the wall during dehydration (due to heat loss) leads to a lower initial state of charge near the wall for the subsequent hydration.
An extensive investigation on the water vapor adsorption into zeolite 13X in an open system has been done by Mette et al. [15], and a 2D numerical model has been developed. In their work, the higher velocity in the near-wall region is implemented by using the extended Brinkman model developed by Vortmeyer and Schuster [16]. The radial heat and mass dispersion coefficients are estimated based on the so-called K r -model developed by Winterberg and Tsotas [17]. In these models, an effective viscosity, introduced by Giese et al. [18], is employed for the wall region, which is exponentially dependent on the particle Reynolds number in the bed. However, the effect of reactor diameter (or the Reynolds number) is not considered in this correlation.
The purpose of this investigation is to understand the significance of the above mentioned effects on the thermal performance of a thermochemical heat storage packed bed reactor. First, a literature survey is done on the different models expressing the transport phenomena in the packed bed, and the best representative model is chosen. Then, the model is validated by experimental results measured in a lab-scale test setup by comparing the pressure drop over the bed, velocity profile below the bed and temperature profile inside the bed. In addition, the profile of adsorbed water concentration in zeolite is compared with experimental results measured by MRI. Next, an investigation on a thermochemical heat storage reactor is done with the validated numerical model. From this work, predictions of the dynamic thermal performance of an adsorption bed on reactor scale can be achieved, which can be used for further studies on the design and  optimization of a thermochemical heat storage system. Suggestions are given in order to optimize the charging and discharging process times, hence, to improve the performance of the reactor.

Models for fluid flow in packed bed
Packed bed reactors belong to the most widely applied reactor types, originating from their effectiveness in terms of performance as well as low capital and operating costs. Because of their popularity, the modeling of flow through porous media attracts the interest of engineers and researchers alike due to the complexity of the modeling involved. However, the modeling can be considerably simplified if one is to consider a homogeneous porous medium where the possible void fraction does not vary significantly and a uniform flow distribution within the bed can be assumed. In this case, flow through packed beds can be modeled by analogy with flow in pipes when the bed void fraction is uniform and low (below 0:65 [19]). The pressure drop through packed beds is the result of frictional losses characterized by the linear dependence upon the flow velocity (Darcy term) and inertia characterized by the quadratic dependence upon the flow velocity (Forchheimer term). Adding these two contributions results in the well-known Ergun equation [20]: where p is the pressure, z the axial coordinate, the bed void fraction, g the viscosity, d p the particle diameter, q the density and u the velocity. The two first terms in the right-hand side, respectively known as Darcy and Forchheimer term, describes the pressure loss caused by the particles. Ergun determined the constants a and b to be 150 and 1:75 for the viscous term (often referred to as Blake-Kozeny-Carman constant) and the inertial term (often referred to as Burke-Plummer constant), respectively. The use of the simple semi-empirical Ergun equation is now generally accepted as a satisfactory prediction of pressure drop in packed beds. However, this is true for infinitely extended packed beds composed of particles that do not differ much in shape from that of spheres [31]. Corrections considering the effect of wall, bed void fraction and particle shape can be applied by modifying the constants in the Ergun equation. Some of the studies on the modified Ergun equation (in a reactor with diamter of d i filled with particles with diameter of d p ) are presented in Table 1. In the pressure drop relations, the average bed void fraction () is used. Some of the relations for average bed void fraction are presented in Table 2.
Existence of an annular wall zone, where the average bed void fraction is greater than in the core of the bed, forms flow channelling in this region. The influence of the wall becomes more significant as the reactor-to-particle diameter ratio d i =d p decreases. In order to take the near-wall flow channelling into account, a radial distribution of the bed void fraction ðrÞ can be used. A typical profile for packed beds consisting of particles with small deviations from spherical shape can be approximated by the following exponential expression after Giese [32]: where c is the void fraction at the core of the bed (far enough from the wall region), r the radial coordinate, R the radius of the bed.The wall region can provide conflicting effects on pressure loss because, on the one hand, the increase in void fraction leads to a reduction in bed resistance to the flow, whereas, on the other, wall friction itself leads to an increase in resistance. Cohen and Metzner [33] highlighted the twofold nature of the wall effect for the first time and have given a qualitative explanation by developing a model which distinguishes the fluid flow in the near wall region and the core region. This is further developed by Nield [34] by including an intermediate layer between the core and wall region. The model is further generalized by DiFelice and Gibilaro [35] for different reactor-to-particle diameter ratios. This type of model can circumvents problems of uncertainties with regard to conditions in the wall region.
The viscous friction at the wall, which increases the pressure drop in the bed, can not be neglected for small d i =d p ratio, due to the fact that the friction surface of the wall increases relative to the total bed surface corresponding to particles. The Darcy-Brinkman-Forchheimer equation, usually used to describe the flow profile within a porous medium bounded by rigid walls, is an extension of Darcy's law by a viscosity term in order to include the viscous forces near the wall; however its validity is restricted to low flow rates. The extended Darcy-Brinkman-Forchheimer equation incorporates the Ergun pressure loss relation as well, to allow for higher flow rates. Vortmeyer and Schuster [16] have quantified the wall effect with incorporating the Brinkman term in the Ergun equation: where g eff is the effective viscosity. The third term in the right-hand side accounts for the pressure loss resulting from viscous friction in near-wall region, known as the Brinkman term. Giese et al. [32]   Dixon [26] ¼ 0:4 þ 0:05ðd i =dpÞ À1 þ 0:412ðd i =dpÞ À2 2 Ribeiro et al. [27] ¼ 0:373 þ 0:917 exp À0:824ðd i =dpÞ À Á developed an exponential relation between g eff and Re p in packed beds of spherical glass particles, by comparing the calculated velocity profile and the measured one using Laser-Doppler-Velocimetry. Bey and Eigenberger [30] developed a nonlinear relation between g eff and Re p and Re in packed beds of spheres, rings and cylinders, by comparing the calculated velocity profile and the measured one by anemometry. The functions for particles with negligible deviation from spherical shape are presented in Table 3.
The dynamics of the adsorption process in a packed bed can be reliably predicted by means of a two dimensional model. However, this presupposes the implementation of both the thermal and the flow effect in the model [36]. In this work, the mass and heat transfer transport phenomena are modeled in axial and radial directions, by including the extended Darcy-Brinkman-Forchheimer equation (3) to model the fluid flow in the reactor. Best representing constants (2) are obtained by comparing the calculated pressure drop to the experimental results in Section 5.1. The most representative viscosity correlation from the different effective viscosity correlations (Table 3) is chosen, by comparing their effects on the velocity profile in Section 5.2.

Reactor setup
The system studied is a packed bed reactor filled with spherical beads (average diameter of 1:6 mm and 3:0 mm) of zeolite 13XBFK (CWK Chemiewerk Bad Köstritz GmbH), in which the adsorption of water takes place. Schematic view of the setup is shown in Fig. 1. The inflow of the reactor is prepared in a controlled way, in respect to the humidity and flow rate, by means of the GFC (Gas Flow Controller) and CEM (Controlled Evaporator Mixer). For hydration, air flow passes through the CEM where the air is mixed with the water flow coming from the water vessel regulated by a LFC (Liquid Flow Controller). In a practical real system with application in the built environment, a borehole system can be employed as the humidifier [37]. For dehydration, almost dry air flow passes through a heater, simulating the heat source (e.g. solar collector), heating up the air to 180 C to dehydrate the material.
The reactor is a cylindrical tube, with an inner shell made of Teflon, because of its low thermal conductivity, and an outer shell made of stainless steel (dimensions are presented in Table 5). In addition, a layer of glass wool insulation is used around the reactor in order to reduce the heat loss. The air enters the reactor from the top side and leaves the reactor at the bottom. The zeolite is placed on top of a filter at the bottom of the reactor. The positions of the 17 thermocouples attached to the reactor are shown in Fig. 2. Temperature is measured with thermocouples placed at the inlet (T in ) and outlet (T out ) of the reactor. Temperature is measured at five different heights numbered from 1 to 5 (h ¼ 90; 70; 50; 30; 10 mm), and three radii in the middle of the bed named M (r ¼ 0 mm), at the inside surface of the reactor body named B (r ¼ 35 mm) and at the outside surface of the stainless steel reactor wall named W (r ¼ 52:5 mm). Temperature and humidity of the outflow of the reactor is monitored by a humidity-temperature sensor. Since large errors in the absolute humidity can occur on low relative humidity measurements at high temperatures, the sensor is positioned after a chiller that cools down the outflow, by which the relative humidity rises and, therefore, the absolute humidity can be determined more accurately.
A manometer (Digitron P200UL) with an accuracy of AE1 ½Pa is used to measure the pressure drop for a range of flow rates (0-150 l/min) over the bed. The glass spherical particles of 1:5 mm, and two sizes of zeolite particles (with average diameter of 1:6 mm and 3:0 mm) with small variation from spherical shape are used. A hot wire anemometer (TSI 8386) is used to measure the velocity 15 mm below the bed, filled with zeolite particles with two different sizes, at different radial positions.

MRI setup
In order to verify the moisture transport model, another similar setup is used. The moisture content in the reactor is measured as function of time and position using an MRI scanner. The reactor has the same dimensions as the reactor used in the reactor setup. The influence of the flow rate and particle diameter (1:6 and 3:0 mm) on the adsorption is studied. The inflow air introduced to the system is conditioned in temperature and humidity, which is subsequently verified with the sensor before the packed bed.
The applicability of MRI measurement to follow the water density during the sorption/desorption process is previously tested and verified [38]. Magnetic Resonance Imaging (MRI) technique is used for non-destructively and quantitatively measuring proton densities in a sample. For constructing an image a Fast Field Echo (FFE) pulse sequence is used [39]. The intensity-time signal is detected by a coil in the MRI scanner and the signal is subsequently Fourier transformed to frequency, which is associates with spatial locations, and an image can be constructed [40].
The signal of water attenuates in an MRI experiment is correlated to the hydrogen density, the echo time and the transversal relaxation time [39]. Hence, it is possible to link the measured signal intensity of the MRI to the water density as the hydrogen density is directly related to the water density. As we used a medical scanner without Faraday coil inside, the coil detunes by variations in water content inside the coil. To counter for this effect, a reference sample is placed next to the sample to minimize the detuning and if necessary correct for signal variations as a result of detuning effects. In addition, signal variations as a result of temperature will be expected to be maximal 10% according to Curie's law [41], based on a DT of 30 K. As a temperature correction is not performed these MRI measurements are semi-quantitative and the signal can slightly vary as a result of temperature variation.

Transport phenomena equations
The equations which describe the dynamics of the system are formulated under the assumptions that: (a) the adsorbent beads have identical properties; (b) the bed properties are uniform; (c) the gas phase behaves as an ideal gas; (d) the gas and solid phases are in thermal equilibrium. The extended Darcy-Brinkman-Forchheimer equation (3) together with the effective viscosity (from Table 3) are used to model the fluid flow in the reactor and to obtain radial profile of the axial velocity. This velocity profile is then used in the mass and heat balances presented in (4) and (5), respectively. The governing mass and heat balance equations inside the bed are presented along the vertical coordinate (z), the radial coordinate (r) and the time (t). These two PDEs contain accumulation, convection, radial and axial dispersion, and source terms. ðrÞ @c @t þ uðrÞ @c @z À 1 r @ @r rD r ðrÞ @c @r where D z and D r are the axial and radial mass dispersion coefficient, respectively, q p is particle density and q water loading in the solid phase.
Heat balance for temperature of the gas and solid phases (T): qC P @T @t þ q g C P;g uðrÞ @T @z À 1 r @ @r rK r ðrÞ @T @r where K z and K r are the axial and radial heat dispersion coefficient, respectively, and DH is the reaction enthalpy. The heat capacity used in the accumulation term of the heat balance is the overall volumetric heat capacity for the gas and solid phases and can be estimated by: The overall volumetric heat capacity for the gas and solid phases consists of the heat capacity of air in the bed voids between particles (first term), the heat capacity of air in the pores inside particles (second term), the heat capacity of solid (third term), and the heat capacity of adsorbed water (forth term). It is dominated by the third and forth terms because of the higher density of the solid phase compared to the gas phase.
The wall temperature (T w ) in each layer is calculated by another PDE for conductive heat transport in a solid wall: This equation holds for all the layers in the wall (Fig. 1), namely Teflon (tef), stainless steel (ss) and insulation (ins). Their characteristics, i.e. the density (q w ), heat capacity (C P;w ) and conductivity (k w ), differ for each material. Cooling at the external surface of the outer layer (insulation) is calculated with _ Q loss ¼ h o ðT w À T amb Þ. The external heat transfer coefficient, h o , can be calculated from a Nusselt number relation at the outside wall of the reactor for natural cooling [42].

Mass dispersion and thermal dispersion coefficients
In (4) and (5), the third and fourth terms represent the radial and axial dispersion, respectively. In reaction engineering, the topic of dispersion (longitudinal and lateral) is treated in detail, and it is generally observed that the data for liquid and gases do not overlap, even in the appropriate dimensionless representation. In addition, dispersion coefficients are a function of different parameters such as length and diameter of bed, particle size, and density and viscosity of fluid. Therefore, finding good representative correlations for dispersion coefficients is important. For very low fluid velocities, dispersion is the direct result of molecular diffusion. As the velocity is increased, the contribution of convective dispersion becomes dominant over that of molecular diffusion. It is commonly assumed that the diffusive and convective components of heat and mass dispersion are additive, in both axial and radial direction [43]. In this section, correlations for mass and heat dispersion coefficients in both axial and radial directions are presented.

Axial dispersion
The axial mass dispersion coefficient D z and axial thermal dispersion coefficient K z can be calculated by [44]:  where Pe m and Pe h are the Péclet numbers for mass and heat transfer, respectively: Pe h ¼ The effective diffusion coefficient of the packed bed can mod- 45]; here d is the diffusivity. The effective thermal conductivity of the packed bed without fluid flow, k eff , can be estimated by the model of Zehner and Schlünder [46], which is satisfactory over a broad range of solid-to-fluid thermal conductivities and solid fractions, as follows: with B the shape factor which for spherical particles can be approximated as B ¼ 1:25 1ÀðrÞ ðrÞ 10=9 and b ¼ k p =k g the particle to gas heat conductivity ratio.
In (8) and (9), the value K z ¼ 2 is used [47] and is estimated theoretically from the equivalence of the dispersion model to a CSTRcascade by Aris and Amundson [48].

Radial dispersion
The radial mass dispersion and thermal dispersion coefficients can be presented based on the velocity in the core region of the bed (u c ) as follow: with f m ðR À rÞ ¼ RÀr K 2;m dp nm ; 0 < R À r 6 K 2;m d p 1; The functions mentioned in (13) and (14) have been proposed by Cheng and Vortmeyer [49] for fully developed forced convective heat transfer in a packed bed bounded between parallel plates and was experimentally approved by Vortmeyer and Haidegger [50] for both heat and mass transfer in a wall-cooled packed bed reactor in the presence of an exothermic chemical reaction. Winterberg et al. [17] have re-evaluated the functions and determined the parameters (presented in Table 4) by comparison of the model to the experimental data available from literature. They interpreted the quantities for the heat transfer case as follows (the parameters for the mass transfer case are analogous): K 1;h determines the slope of rise in the effective thermal conductivity with flow velocity; K 2;h sets the damping parameter after the coefficient begins to decline towards the wall; the exponent n h determines the curvature of the damping function.

Adsorption
In the mass balance for the gas phase in (4) and heat balance for the gas and solid phases (with thermal equilibrium) in (5), the source terms are related to the kinetics of the adsorption reaction. For many adsorption systems, the diffusion-controlled kinetics may be satisfactory represented by the Linear Driving Force (LDF) approximation first time introduced by Glueckauf [51]: where q and q eq are the adsorbed and equilibrium concentrations in the solid phase in moles of water per kilograms of dry zeolite, and k LDF is the mass transfer coefficient between the fluid phase and the solid phase. In previous work [14], an extensive discussion on the Linear Driving Force (LDF) approximation is presented. The Langmuir-Freundlich equation has been used to estimate the equilibrium loading of water in the adsorbed phase on zeolite (q eq ): where q max is the maximum amount of adsorbed water at high water vapor pressures, p water is the partial pressure of water vapor, and b and n are temperature-dependent parameters, which can be described as [52]: where T 0 is the reference temperature and is assumed to be 0 C. Parameters b 0 and n 0 are the adsorption affinity constant and exponent constant at reference temperature, respectively. DE is the activation energy for desorption. By using the van't Hoff equation, a concentration dependent form of the isosteric heat of adsorption can be obtained as follows [52]:

Parameters
The complete model as presented, is developed and solved in COMSOL Multiphysics, which is a Finite Element Method (FEM) based program. The PDEs are discretized in space by Lagrange linear shape functions. Both the bed and wall domains are meshed with an element size of 1 mm (chosen based on a mesh convergence study). The solver is a fully coupled time-dependant solver. A Backward Differentiation Formula (BDF) scheme (maximum order of 5 and a minimum order of 1) with an initial time step of 0:001 s is used as the time stepping method. The relative, absolute and event tolerances are set to 1EÀ2. The other parameters used in the model (i.e. properties of the material, characteristics of the reactor, and operational conditions) are presented in Table 5.

Pressure drop
The pressure drop over the bed of particles in the reactor setup is measured for glass particles of 1:5 mm diameter and zeolite particles of 1:6 and 3:0 mm average diameter, for a range of flow rate 0-150 l/min with increments of 5% (superficial velocity 0-0.65 m/s). The measured pressure drop per unit of reactor length are compared to the one calculated using the relations discussed in Section 2.
In Fig. 3(a), the pressure drop measured over the bed filled with glass particles is compared with the ones calculated by using different relations for friction factors a and b in the modified Ergun equation. It can be seen that using cor #5 from Table 1 (Cheng [25]) gives the best fit, for the whole range of flow rate. The estimated bed void fraction is more accurate for the bed filled with glass particles, than for the one filled with zeolite particles, since the range of size distribution is wider for zeolite. The small and big zeolite particles are in a range of 1:4-1:8 mm and 2:5-3:5 mm in diameter, while the glass particles are almost uniform in size. In addition, the lower sphericity of zeolite particles has an effect on both the bed void fraction and the friction factors (a and b) in the modified Ergun equation. In Fig. 3(b), the effect of the average bed void fraction on the calculated pressure drop over the bed is shown. For both small and big zeolite particles, cor # 2 from Table 2 (Ribeiro et al. [27]) gives a better estimation than cor # 1 [26]. For large zeolite particles, the deviation is larger and a better fit is achieved by ¼ 35:8%.

Velocity profile
Measured radial velocity profiles below the bed filled with small and big zeolite particles for flow rate of 30; 60 and 120 l/min are shown in Fig. 4. Maldistribution of the flow near the wall can be clearly seen for all the cases. Additionally, wall channelling becomes less pronounced for higher Reynolds numbers. Since the inertial term in the bed overrides the viscous term near the wall and results into more dispersion of momentum, thus the velocity profile becomes smoother. This effect is the same for both sizes of the zeolite particles. However, the peak in the velocity profiles for larger particle size is measured to be at a larger distance from the wall. Furthermore, the velocity peak near the wall is wider and of a smaller magnitude for the larger particle size. That is caused by the fact that the void space near the wall is larger for a pack of larger particles compare to one of small particles.
Different relations for effective viscosity are presented in Table 3. The different models for predicting effective viscosity result in a large spread in simulated velocity profile in packed bed. It confirms that the models should not be used outside the range over which they are measured. Therefore, the ones proposed by Einstein [28], Brinkman [29] and Givler and Altobelli [53] are not suitable here. In addition, it can be concluded from the velocity profile measurements that the shape of radial profile depends on the Reynolds number (or superficial velocity) and the particle size (Fig. 4). Therefore, the relation proposed by Bey and Eigenberger [30] is chosen, which gives a more general estimation of effective viscosity in pack bed compared to Giese et al. [18].
In Fig. 5, the effects of the Reynolds number and particle size on simulated velocity profile are shown over a small region near the wall. Unfortunately, the velocity distribution in the bed itself could not be measured. It should be noted that the modeled velocity peak in the bed occurs in a very narrow region near the wall, while the velocity profiles measured below the bed show wider peaks. That is because the velocity profile changes rapidly from the shape in the bed (simulated velocity profile in the bed) to the shape related to flow in an empty pipe. Vortmeyer and Schuster [16] and Bey and Eigenberger [30] claimed that it is possible to calculate this effect by solving numerically the twodimensional Navier-Stokes equation for the developing flow inside the empty tube with the artificial flow profiles of the packed bed as inlet condition. Despite the fact that the measured velocity profiles are in the transition state, the effects of the Reynolds number and particle size are still clear. These effects are similar between measurements below the bed and simulated velocity profiles, when the relation proposed by Bey and Eigenberger [30] is used for the effective viscosity.

Moisture content
The MRI technique is used to investigate the moisture content q, in the zeolite bed. The MRI signal intensity in a cross section at a fixed height of the bed (L z ¼ 90 mm) is shown in Fig. 6(a) (at t½s ¼ 1225; 2042; 2453; 3270; 3755; 4572 from top figure to bottom, respectively) and over a longitudinal section at the center of the bed is shown in Fig. 11(b) (at t½s ¼ 410; 1636; 2864; 4166; 5394; 6622 from top figure to bottom, respectively). The moisture content is indicated ranging from blue, low moisture content, to red, high moisture content. The earlier occurrence of the reaction in the wall region compared to the core region is clearly visible in the longitudinal section view. As can be seen also in the cross section view, the higher moisture content happens first in an annular region near the wall, and later in the core region.
The MRI experiment is performed for different flow rates and particle sizes. For each case, the normalized integration of the measured signal intensity over a cross sectional slice located at a fixed height (L z ) during time is calculated. The normalized signal intensity directly represents the total moisture content in the specific slice of the bed. In Fig. 7, the total adsorbed moisture is shown vs. time (t À t 0 ), where t 0 is calculated by t 0 ¼ ðL z =uÞ ðq b Dq=c in Þ. The slope of the curves is an indication of the reaction front sharpness. For example, as an extreme case, if the curve was a step function it means that the reaction front is flat. It is observed that decreasing the flow rate promotes the wall channelling, and for a larger particle size the intensity increases slower compared to a smaller particle size.
The location of adsorbed moisture front can be extracted from the MRI experimental results. From each of the longitudinal section images presented in Fig. 6(b), the location of the moisture front is extracted. Fig. 8 shows the comparison between the extracted front from the MRI measurement and the one calculated by the numerical model. An excellent fit is found by choosing K 1;m ¼ 2; K 2;m ¼ 0:44 and n m ¼ 2 for parameters used in (13) and (15). However, in order to generalize this result measurements on more different conditions are needed. The amount of adsorbed water at different heights of the bed over time can also be derived from the cross section images presented in Fig. 6(a). Fig. 9 shows the comparison between the experimental and numerical concentration of adsorbed water in the solid phase (q). The concentration at each height initially increases slowly due to the peak near the wall and later increases sharply due to the flat core region.

Hydration
The temperature measured and simulated at the location of M thermocouples at different heights of the bed (1-5 as explained in Section 3.1) during hydration of the bed is shown in Fig. 11. In addition, the concentration of water at the outlet of the reactor is shown in Fig. 12. The water introduced at the inlet of the reactor is adsorbed during hydration and the temperature increases as a result of released adsorption energy. When the reaction front passes the location of a certain thermocouple and the temperature drops to the temperature of the cold humid inflow air, the material at that location is completely hydrated, hence no further reaction occurs at that specific place in the reactor and no energy is released. The temperature at the outlet stays high till all material has reacted, after which the temperature at the outlet drops and the water concentration increases to the inlet concentration.
If the reaction front is flatter in the radial direction, it will lead to a sharper transition in the temperature and concentration at the outlet of the reactor. In Fig. 10, the concentration of water in the gas and solid phases, and the temperature of the bed at a specific     bottom, the water concentration of the outflow air start to increase. Conversely, when the temperature (T) front passes the bottom, the temperature of the outflow air start to decrease. The sharpness of this transition depends on the flatness of the fronts. Effects of different considered parameters on the transition are studied in Section 6.1.

Dehydration
The bed is dehydrated in the reactor setup. The temperature program used for dehydration consists of a ramp from ambient temperature to 180 C with a heating rate of 0.67°C/min, and a isotherm at 180°C. The measured and simulated water concentration at the reactor outlet during 20,000 s dehydration are shown in Fig. 13. As can be seen, the comparison between experimental and numerical curves is far from perfect. The experimental result shows a bimodal shape, while the numerical result shows only one peak. The first peak can be caused by escape of loosely bonded water molecules to zeolite at lower temperatures, and the second peak by release of strongly bonded water molecules to zeolite at higher temperatures. However, this should be investigated by measurements in small scale samples. In spite of different profiles from experiment and simulation, the total amount of desorbed water from the material over time is almost the same. After 20,000 s of dehydration, the total amount of desorbed water calculated from experimental and numerical results is 87.8 g and 86.5 g, respectively. However, the amount of desorbed water and the State of Charge (SoC) in the material are dependent on the duration of the charging process.
In Fig. 14, the SoC in the material as a function of time is shown. The SoC after 16,000 s is still high in the wall region at the bottom of the reactor, while it is getting drier after 20,000 s. Finally, after 24,000 s, the average water loading in the reactor is around 1.09 mol/kg. It is almost the highest achievable SoC with this charging temperature (around 180°C). As can be seen, the improvement in the SoC from 20,000 s to 24,000 s after the start of charging process is minor. Different duration of the charging process leads to a different final SoC for the charging process and, subsequently, a different initial SoC for the succeed discharging process.   Therefore, the final state of charge of the zeolite bed as a function of the radial and axial coordinate obtained after a dehydration simulation is used as an initial condition of the subsequent hydration simulation, and contrariwise. It concludes the importance of the durations of the charging and discharging processes. The optimum durations of the charging and discharging processes are studied in Section 6.2.

Parametric study
In the previous section, the model is validated by the results from different experimental techniques performed in the reactor with the characteristics presented in Table 5. In this section, the validated model is used to study the effects of different parameters on the performance of an up-scaled reactor. A segment of a large scale reactor with a segment volume of about 62.5 L is defined in the model and a larger air flow rate (compared to the lab-scale case) is applied to get a high heating power in the order of 1 kW, corresponding to a full scale design. The inlet water concentration is chosen at 13 mbar water vapor pressure (being the vapor pressure at 10°C, which is a typical borehole temperature in the Dutch climate [37]), and the whole reactor is initially at ambient temperature. The characteristics of this large scale reference case are presented in Table 6. In this section, the dimensionless time is used, which is defined as the time divided by the residence time (h ¼ t=s), for dehydration and rehydration processes.

Effects
In this section, the effect of different parameter settings in the model on the reaction front in the reactor are evaluated. The water concentration and temperature at the outlet of the reactor during the discharging process for the different cases are presented in Fig. 15. Considered Cases are: case #0: original model (reference case) case #1: flat velocity profile case #2: flat initial SoC case #3: no heat loss through reactor wall case #4: combination of cases #1 and #2 By neglecting the radial dependencies of the velocity, initial SoC and bed void fraction (case #4), the transition in the water concentration and temperature at the outlet of the reactor is  sharper compared to the reference case (case #0). The cases #1; #2 and #3 show a intermediate transition compared to these two extreme cases. The cases #2 and #3 show a similar trend, which means that the main effect of decreasing the heat losses, is to obtain a flatter SoC during dehydration, and therefore a sharper transition compared to the reference case.

Duration of charging and discharging processes
The convective heat flux introduced to the reactor during the charging process (dehydration) and extracted from the reactor during the discharge process (rehydration) are calculated and shown in Fig. 16. The powers during the dehydration and rehydration are calculated by: The energy provided to the reactor during the charging process and energy extracted from the reactor during the discharge process can be calculated by integrating the powers over time: For the system presented in Table 6, the energies are calculated for the charging process till t ¼ 30; 000 s and the discharging process till t ¼ 35; 000 s. The calculated energies for the charging and discharging processes are 111.2 MJ and 45.6 MJ, respectively. The efficiency is calculated by (24), and is around 41% for this complete cycle of charging-discharging.
The effect of the charging and discharging duration on average power (power per flow rate P=Q ) during discharging and total efficiency of the reactor are investigated (Fig. 17).
For D=H ¼ 1 and s ¼ 1:56, increasing the duration of charging process (t deh ) from 20,000 s to 30,000 s will increase the power, but the change in power by increasing t deh form 30,000 s to 40,000 s is negligible. However, the efficiency will drop by increasing t deh . The efficiency is generally increasing by increasing the duration of charging process (t reh ). However, the change is negligible for t reh larger than 30,000 s (h reh P 19; 230). For the case t deh ¼ 30; 000 s, the maximum power is achieved by t reh being between 30,000 s and 35,000 s (19; 230 P h reh P 22; 440). The reaction front for this case at different rehydration times is shown in Fig. 18. The low total efficiency is caused mostly by the heat losses during the dehydration process (19.76 MJ) compared to the rehydration process (1.47 MJ), for t deh ¼ 30; 000 s and t reh ¼ 35; 000 s. The heat loss mechanisms are conduction through the walls (around 91% and 43% of the heat loss during the dehydration and rehydration processes, respectively) and sensible heat in the material (around 9% and 57% of the heat loss during the dehydration and rehydration processes, respectively).
For D=H ¼ 1 and s ¼ 0:31, both the efficiency and power per volumetric air flow rate are improved. This is caused by the fact that heat losses are less for this case with higher flow rate (s ¼ 0:31) compared to the case with lower flow rate (s ¼ 1:56), because the process is faster. The optimum for this case is also occurring for t reh being between 6000 s and 7000 s (19; 230 P h reh P 22; 440). However, the operating flow rates should be chosen based on the charging and discharging strategy developed according to the heat supply and demand in the house.
For D=H ¼ 2 and s ¼ 1:56, both the efficiency and power per volumetric air flow rate decline. The efficiency is the highest for t reh ¼ 40; 000 s (h reh ¼ 25; 640 s), while the power is at the lowest. The optimum power still occurs for t reh being between 30,000 s and 35,000 s (19; 230 P h reh P 22; 440). However, the reactor aspect ratio should be chosen based on the COP of the reactor, considering the pressure drop over the bed and the required fan power to blow air through the bed.
Considering both power and efficiency, the optimum suggested operating condition is to apply dimensionless duration of charging (h deh ) of around 19,230 and dimensionless duration of discharging (h reh ) between 19,230 and 22,440. Specifically, for a reactor with a bed volume of about 62.5 L and air flow rate of 0:04 m 3 =s (s ¼ 1:56 s), the optimum durations of charging and discharging processes are around 30,000 s and 35,000 s, respectively.

Aspect ratio
The effect of the reactor aspect ratio on the efficiency and COP of the reactor is studied in this section. Volume of the reactor and air flow rate are 62.5 L and 0:04 m 3 =s, respectively. The optimum durations of charging and discharging process (around 30,000 s and 35,000 s, respectively), concluded from previous section, are used here. The COP of the reactor is calculated as the total energy extracted during rehydration divided by the total required fan energy to blow air through the bed during rehydration and dehydration (COP ¼ E reh =E fan ). The total required fan energy is calculated as E fan ¼ Q DPðt reh þ t deh Þ=g fan , which is the require energy to overcome the pressure drop over the bed DP (estimated by the fluid flow model), for an air flow rate of Q, during rehydration and hydration processes, with a fan with an efficiency of g fan ¼ 50%. Fig. 19 shows the efficiency and COP of the reactor for different aspect ratios. The maximum efficiency of around 41:5% is achieved for the aspect ratio of 1:5, however the change in efficiency for the other aspect ratios (ranging from 0:5 to 3:0) is negligible (ranging from 40:5 to 41:5%). While, as expected, the COP increases continuously by increasing the aspect ratio.

Conclusions
A literature survey is done in order to find the best model to evaluate pressure drop and velocity profile in a packed bed of porous media, considering the radial effects, specifically the effect of the rigid wall of the reactor. By comparing the measured pressure drop over the bed to the model results, good agreement is found for the Ergun equation with the modified constants presented by Cheng [25]. However, the particle diameter distribution must be taken into account, because a wider deviation from the average diameter results in a lower bed void fraction and higher pressure drop. The velocity is measured and wall channelling is observed. An increase in wall channelling was found for low flow rates and smaller particle diameters. The best representative model for this behavior is found to be the Brinkman extension to the Darcy-Forchheimer equation in combination with an effective viscosity expression by Bey and Eigenberger [30] and radial bed void fraction profile by Giese [32].
A 2D mass and heat transfer model is developed by assuming a homogeneous temperature for the gas and solid phases, and the radial profile of the axial fluid velocity in the bed from the aforementioned velocity model. In addition, the heat transfer in the solid wall of the reactor is modeled and coupled to the bed model at the inside wall boundary. The model is validated using experimental thermal analysis from the reactor setup, and the moisture content results from the MRI setup.
The MRI technique is proven to be an excellent addition to the experimental research into the reactor. The moisture content in a reactor of up to 3 l can be measured as a function of time with a spatial and time resolution of approximately 2 mm and 2 min (not simultaneously). The effects of flow rate and particle diameter on the shape of the moisture front is studied and compared to the numerical model for the benchmark case. From this analysis, the optimal value of the mass dispersion parameters presented by Winterberg et al. [17] are determined, and the dependence of the  Fig. 19. Efficiency and COP of the reactor for different reactor aspect ratios, for Q ¼ 0:04 m 3 =s and dp ¼ 3:0 mm.  dispersion parameter K 1;m on the Peclet number is validated. Further investigation is needed in order to generalize and quantify it for a range of conditions. Reasonable agreements between the temperatures and concentrations from the experiment and numerical model during hydration and dehydration are found. Measurements indicated that the influence of the dehydration cycle on the hydration cycle is significant. Therefore, the final state of charge of the zeolite bed as a function of the radial and axial coordinate obtained after dehydration simulation is used as an initial condition of the hydration simulation, and contrariwise. It concludes the importance of the durations of the charging and discharging processes.
Finally, the significance of each considered effect in the model on the performance of the reactor is assessed. An investigation is done on the effect of re/de-hydration duration on the average power and total efficiency of the reactor. The optimum suggested operating condition for a reactor with a bed volume of about 62.5 L and air flow rate of 0:04 m 3 =s (s ¼ 1:56 s) is to apply durations of charging and discharging processes of around 30,000 s and 35,000 s, respectively. The aspect ratio of 1:5 leads to the highest efficiency in the reactor, however the change in efficiency for the other aspect ratios (ranging from 0:5 to 3:0) is negligible. While, as expected, the COP increases continuously by increasing the aspect ratio. The insight gained in this study on the influence of flow rate and tube-over-particle diameter ratio on the velocity, mass and heat transfer can be used in further studies to design an optimal reactor for a large scale storage system.