On Process Intensiﬁcation through Membrane Storage Reactors

: In this work, a dynamic, one-dimensional, ﬁrst principle-based model of a novel membrane storage reactor (MSR) process is developed and simulated. The resulting governing equations are rendered dimensionless and are shown to feature two dimensionless groups that can be used to affect process performance. The novel process is shown to intensify production of a desired species through the creation of two physically distinct domains separated by a semipermeable boundary, and dynamic operation. A number of metrics are then introduced and applied to a case study on Steam Methane Reforming, for which a parametric study is carried out which establishes the superior performance of the MSR when compared to a reactor operating at steady state (SSR).


Introduction
Over the course of the last forty years, process intensification (PI) has continued to develop as an area of active chemical engineering research, incorporating numerous considerations, including process safety [1] and process systems engineering [2]. PI encompasses any process design strategy that leads to a smaller, cleaner, safer, and/or more energy-efficient technology. Additionally, PI also includes system designs which reduce the number of devices employed [3]. To this end, chemical reactor design incorporating separation technologies continues to be an active area of PI research, with prominent examples being membrane reactors (MR) carrying out steam methane reforming (SMR)based hydrogen production at lower temperatures [4], dividing-wall columns combining multiple distillation channels into a single unit [5], and compact catalytic plate reactors for Fischer-Tropsch synthesis [6]. Our team has designed and analyzed intensified reactor configurations for hydrogen production via the Water Gas Shift (WGS) reaction. At the labscale, we have investigated the efficacy of carbon molecular sieve (CMS)-based membrane reactors for use in combined reactive separation networks [7]. We have also developed and experimentally studied a novel hybrid membrane reactor-adsorptive reactor (MR-AR) configuration. This technology, which is currently being field-tested, produces high-purity hydrogen with simultaneous carbon dioxide capture for application to Integrated Gas Combined Cycle (IGCC)-based power generation [8].
As advances in computational software have continued, there has also been a substantial increase in the number of tools developed for identifying new PI methodologies at the theoretical level. Advanced mathematical formulations, such as our Infinite DimEnsionAl State-Space (IDEAS) framework [9] and multi-objective optimization techniques [10], have helped introduce systematic approaches for developing and identifying PI pathways for various chemical systems. These advances have naturally led to the development of software for the generation of sustainable design alternatives to be used for PI purposes [11], and for using thermodynamic analysis to assess the viability of proposed technologies [12].
In addition to our experimental efforts noted above, we have carried out extensive theoretical research in developing PI models. These efforts include the design and sim-

Mathematical Formulation
In this study, a composite one-dimensional model for the intensified MSR process is derived which captures and highlights several attractive MSR characteristics. Following non-dimensionalization, dimensionless groups governing the MSR behavior are identified. The MSR is considered to be a composite thermodynamic system (r) comprised of three domains that are spatially distributed and exclusive of each other (reactor voids (v), catalyst pellet (c), and storage pellet (s)), and a maximum of two phases (gas (g), and solid (so)). The volume fractions of the voids (v), catalyst (c), and storage (s) domains are denoted as ε v , ε c , ε s , respectively, while the volume fractions of the gas and solid phases in the storage domain are denoted as ε gs , ε sos , respectively. In the context of this study, the following are then considered to hold: The communication between the catalyst and reactor voids (gas phase) domains is quantified using an effectiveness factor approach, and the communication between the storage domain and reactor gas phase domains through a permselective layer on the storage domain's control surface is quantified using species transport equations obeying Sieverts' Law. The composite system (r) is considered to be isothermal, whereas the reactor voids domain (v) is also considered to be isobaric. The molar diffusion flux of gaseous species in the voids domain is considered to be negligible compared to the species convective molar flux. No reaction occurs within the storage domain. The reacting mixture is considered to be an ideal gas. The equation describing Mass Conservation (molar form) of species i in phase g within domain v on a system r volumetric basis is: , T + S i,gv,gs (t, z) i = 1, NC The Peclet number, Pe = v e f f L char D e f f , is typically used to determine whether a process is dominated by convective or diffusive mass transport processes, where v e f f is a gas velocity, L char is a characteristic length, and D e f f is the effective diffusion coefficient of the transferrable species. Although there is still considerable debate regarding the choice of these parameters, which depend on the structure and geometry of the system at hand [19], it is generally agreed that Pe 1 is sufficient to ignore diffusive effects. For example, the porosity of the system that is available for diffusion is generally smaller than the total porosity itself [20]. The parameter D e f f is evaluated as D e f f = ε v D AB using the binary diffusion coefficient, D AB , of the transferrable species, which in turn can be evaluated using the Chapman-Enskog solution to the Boltzmann Equation along with a Lennard-Jones potential for estimating the effects of intermolecular forces. Using the formulation established in our previous results [14], the values of D AB for the species present in the considered case study fall in the range of 0.02 ≤ D AB ≤ 0.15 cm 2 s . Values for the characteristic length can be obtained by considering L char to be equal either to the length or the radius of the reactor, which in turn yield Pe values that satisfy Pe ≥ 1 × 10 5 and Pe ≥ 100, respectively, and justify neglecting diffusion in the void domain." The equation describing Mass Conservation (molar form) of species i in phase g within domain s on a system r volumetric basis is: where the temporal and spatial dependence of all volume fractions, the effectiveness factor, and the catalyst pellet density are not shown, because they will next be considered constant.
Sieverts' Law has commonly been employed in quantifying species' molar fluxes through Pd membrane layers [18]. For such membranes, it has been shown, both theoretically [21] and experimentally [22], that the molar flow rate is proportional to a nth power differential of partial pressures, where 0.5 < n < 1. In this work we employ the same transport equation with n = 1, because it has been shown to hold true for high-temperature SiC membranes [23]. The resulting equation is then: where β i is the i th species molar permeance through the permselective layer, α s,v is the storage-void interfacial area per unit reactor volume, and P i,gv , P i,gs are the i th species partial pressures in the gas phase of the reactor void and storage domains. Considering an ideal gas mixture, the following relations hold: where R is the universal gas constant, and T is the reactor temperature that is considered to be constant over both time and space. Combining the above equations yields: Summing the gas phase species conservation equations in the void domain and utilizing the isobaric void domain assumption yields: Substitution of Equation (7) into Equation (6) then yields: Next, dimensionless variables are introduced, and dimensionless groups are identified governing MSR behavior: where P * is the operating pressure of the reactor, L * is the length of the reactor, v * is the inlet velocity to the reactor at time zero, and t * = L * /v * is the reactor's residence time at time zero, and r* is defined later in this paper. Given that z is the reactor's axial coordinate, L * is the length of the reactor, and z z L * , the reactor inlet and outlet can be designated as z = 0 , z = 1, respectively.
Given that the reactor voids domain (v) is considered isobaric, with operating pressure P * , then NC ∑ k=1 P k,gv = P * , which implies that NC ∑ k=1 P k,gv = 1. Then, the resulting dimensionless equations for the MSR model are: The above MSR dimensionless model suggests that two dimensionless numbers determine the MSR's dynamic behavior. The first dimensionless quantity, Θ = 1/Pe, provides a measure of how effectively the desired product species is being extracted from the reactor void domain into the gas storage domain, as compared to its axial convective transport in the reactor void domain. The Peclet number, P e , is a dimensionless group commonly used in the analysis of membrane reactors [24][25][26][27]. The Damkohler number, D a , captures the importance of the rate of convection compared to the reaction rate and can be thought as the ratio of the characteristic time for convective flow to the characteristic time for reaction. The use of these dimensionless numbers in assessing the performance of membrane reactor systems reduces the dimensionality of the design parameter space for the MSR process and is widely practiced in MR experimental and theoretical research [28,29].
As mentioned earlier, the proposed MSR process is operated in a dynamic (periodic) manner. A possible such dynamic operation is a "Partial Pressure Swing Operation" (PPSO) of the MSR, which keeps the reactor pressure and temperature constant, and involves three modes of operation. In the first Operating Mode (OM), with duration τ 1 (dimensionless duration τ 1 ), the MSR operates in a Loading-Reaction mode for phase (g) in domain (v) (where the reactants are being loaded (fed) into the MSR and the desired reactions are carried out in the presence of (c) and accounted for in (v)), and in a Storage mode for phase (g) in domain (s) (where the desired species is preferentially transported from (v) to (s), where they are stored). In the second OM, with duration τ 2 (dimensionless duration τ 2 ), the MSR operates in a Reactant-Flushing mode for phase (g) in domain (v) (where the reactants are being removed from (v)), and in Storage-Maintenance mode in phase (g) in domain (s) (where the desired species is maintained in storage within (s)). In the third and final OM, with duration τ 3 (dimensionless duration τ 3 ), the MSR operates in an Emptying mode for phase (g) in domain (s) (where the desired chemicals are emptied from storage within (s), and transported into (v)), and in Unloading-Production mode for phase (g) in domain (v) (where the desired species are removed from (v), to yield the desirable MSR products).
Comparing the performance of a periodic process with that of a traditionally steadystate process necessitates the introduction of several performance metrics. These metrics can be defined for each OM and for the process as a whole, namely for the duration of all three aforementioned OMs. To this end, the axial molar flowrate of the i th species F i,gv is expressed in terms of a reference molar flowrate F * and a dimensionless molar flowrate F i,gv , as follows: Then, the following dimensionless metrics are defined and explained below:

Limiting Reactant Conversion
An important metric for the assessment of MSR performance is the conversion of a limiting reactant R LI M . Typically, this reactant will be fed in the MSR during one OM, but may be removed from the MSR during multiple OMs. Thus, we define its conversion over all OM as follows:

Desired Product Ratio
Another important metric for the assessment of MSR performance is one that captures the extent to which the desired product is formed as the limiting reactant is converted. Thus, the following OM-dependent metric is introduced that is equal to the ratio of the i th species amount generated during OM k over the limiting reactant's R LI M amount fed throughout all OMs.
where τ 0 = 0, τ 0 = 0. An overall OM-independent metric is also introduced that is equal to the ratio of the i th species amount generated during all OMs over the limiting reactant's R LI M amount fed throughout all OMs.
This is often referred to as product yield. It then holds:

Desired Product Recovery Fraction
Finally, we define the fraction of product generated during a given OM over the total product generated over the duration of all OMs as the Product Recovery fraction, defined as: The system of non-linear first-order partial differential equations described by Equation (10) is used to simulate all three OMs of the MSR PPSO, and is solved using a direct finite element solver with implicit multistep free backward differentiation (BDF) with total time discretization consisting of 1 × 10 6 time steps for each simulation. A fractional initial step for backward Euler consistent initialization was used, and linear first-order Lagrange shape functions were used for space discretization and consistent stabilization based on full calculation of the residual was employed. Because of the changing boundary conditions, a minimum non-linear damping factor of 1 × 10 −9 was used during the Newton-Raphson iteration. The above solution strategy was coded and implemented using the COMSOL Multiphysics software platform. The MSR concept is next illustrated for the SMR process. We demonstrate, through parametric analysis, the impact of the previously identified dimensionless parameters D a , Θ on the performance of an SMR-MSR process operating under PPSO.

Steam Methane Reforming (SMR) Case Study
In this case study we propose intensifying SMR-based hydrogen production, using a MSR under PPSO whose design is pursued using the above developed dimensionless MSR model. It is shown through the use of quantifiable design metrics that the proposed SMR-MSR outperforms a conventional SMR-SSR. Approximately 95% of the hydrogen presently produced industrially in the United States is obtained via the SMR reaction [30]. Industrial reactors operate at or near equilibrium, and due in part to the high endothermicity of the SMR reaction, encompass a large portion of the operating cost of a refinery. Although membrane separation technology, using polymeric and inorganic membranes, has been investigated and shows promise for broad industrial applications [18], there is still a need for high pressure and temperature gas-phase applications.
According to Xu and Froment [31] SMR is carried out through three reversible reactions RX1, RX2, and RX3, with enthalpies of formation and kinetic rates as shown below: The kinetic and equilibrium constant values in the above expressions are listed in Tables 1 and 2.  Table 2. SMR reaction equilibrium constants.

Equilibrium Constant Units
The above kinetic rates can be brought in dimensionless form, as suggested in [17,32]. In the definition of the Damkohler number, the reference reaction generation rate r * is selected as the reaction rate for RX1 evaluated at the reactor entrance conditions during OM 1: The resulting dimensionless reaction rates become: The PPSO of the SMR-MSR is carried out in three distinct operating modes: OM 1 (Loading-Reaction/Storage), OM 2 (Decarbonization/Maintenance), and OM 3 (Unloading-Production/Emptying) each of which is characterized by different processes occurring in the (v) and (s) domains. These three OMs have a time duration designated as τ 1 ,τ 2 , and τ 3 , as described next and shown in the Figure 1 below.

OM 1: MSR Loading-Reaction/Storage Phase
At the beginning of this OM, phase (g) in domain (v) is largely composed of steam and minute amounts of hydrogen, phase (g) in domain (s) only contains hydrogen at a pressure slightly above the partial pressure of hydrogen in (v), and a mixture of methane and steam at a constant flowrate is fed into the MSR. During this OM, the hydrogen that is produced from the SMR reactions taking place begins to permeate into (s) as its partial pressure in (v) increases above the total pressure of (s). Under isobaric conditions in (v), the outlet flowrate varies, as described in Equation (10). The duration τ 1 of OM 1 is chosen so that a desirable level of limiting reactant conversion is attained. Reactant conversion, as defined by Equation (12), is defined over all three OMs. Typically, however, the limiting reactant is fed to the MSR only during OM 1, whereas any unreacted limiting reactant exits the MSR during OM 1 and OM 2. As discussed below, the duration of OM 2 is fixed and chosen to be equal to the reactor's OM 2 residence time. Thus, the duration τ 1 of OM 1 is the only degree of freedom that can be chosen so that the amount of unreacted limiting reactant that exits the MSR over τ 1 + τ 2 (i.e., over the duration of OM 1 and OM 2) is acceptably low, compared to the amount of limiting reactant fed to the MSR during OM 1, so that a desirable level of limiting reactant conversion can be attained.

OM 1: MSR Loading-Reaction/Storage Phase
At the beginning of this OM, phase ( ) g in domain ( ) v is largely composed of steam and minute amounts of hydrogen, phase ( ) g in domain ( ) s only contains hydrogen at a pressure slightly above the partial pressure of hydrogen in ( ) v , and a mixture of methane and steam at a constant flowrate is fed into the MSR. During this OM, the hydrogen that is produced from the SMR reactions taking place begins to permeate into ( ) s as its partial pressure in ( ) v increases above the total pressure of ( )

OM 2: MSR Decarbonization/Maintenance Phase
OM 2 begins at the final conditions of OM 1, and the reactor is then fed pure steam at a constant flowrate. If the reactor is operating under plug flow conditions with no diffusional effects, as is the case here, a sharp steam front is formed that propagates through

OM 2: MSR Decarbonization/Maintenance Phase
OM 2 begins at the final conditions of OM 1, and the reactor is then fed pure steam at a constant flowrate. If the reactor is operating under plug flow conditions with no diffusional effects, as is the case here, a sharp steam front is formed that propagates through the reactor and reaches the reactor outlet in time equal to the reactor's OM 2 residence time. In the reactor section ahead of the steam front, the SMR reaction continues to occur in the void domain, and hydrogen is transported from the void to the storage domain. Behind the steam front hydrogen is transported from the storage to the void domain. Thus, under plug flow conditions, the duration τ 2 of OM 2 is chosen to be equal to the reactor's OM 2 residence time, because at that time complete reactor decarbonization has been attained. (If dispersion was to be present in the reactor, the steam front would not be as sharp. In that case, the duration τ 2 of OM 2 may need to be chosen to be greater than the reactor's OM 2 residence time for the desired level of reactor decarbonization to take place).

Phase 3: MSR Unloading-Production/Emptying Phase
Similarly, OM 3 begins at the final conditions of OM 2, and the reactor is again fed pure steam at a constant flowrate. This action empties the contents of (s) back into (v), generating a mixture of readily separable, high-pressure hydrogen and steam as the SMR-MSR product. The duration τ 3 of OM 3 is chosen so that the hydrogen in the domains (v), (s) is below a desirable low level.
To compare the performance of the SMR-MSR under PPSO with that of a traditional SMR-SSR, we utilized the process performance metrics defined by Equations (12)-(16) for the SMR case study. To this end, identifying methane as our limiting reactant, Equation (12) becomes: Similarly, Equations (13)- (16) can be expressed for all products (carbon monoxide, carbon dioxide, and hydrogen), and are listed below for hydrogen.
Equation (24) captures the H 2 molar production ratio in OM 1 over the total amount of methane fed. Equation (25) describes the hydrogen yield over all three OMs, i.e., the total moles of hydrogen produced divided by the moles of methane fed to the reactor. Lastly, Equation (26) is used to calculate the hydrogen recovery ratio during the 3rd OM, which is of importance because it quantifies how much hydrogen is produced in readily purifiable form.
To determine the time duration τ 1 , τ 2 , and τ 3 for the three OMs, a criterion for terminating each OM must be selected. The reactor contents at the beginning of OM 1 consist largely of steam, while the reactor feed at that time switches to an appropriately chosen methane steam mixture, which in turn leads to a methane front propagating through the reactor. Thus, at each fixed reactor location, methane partial pressure initially increases, reaches a maximum value, and then decreases as the reactor feed is switched to steam at the beginning of OM 2. The maximum methane partial pressure attained becomes ever lower at downstream reactor length locations. With this behavior in mind, our selection τ 1 for the duration of OM 1 will be related to the conversion (utilization) of the raw material (methane).
The duration τ 2 of OM 2 is selected to be equal to the OM 2 residence time, ensuring that the dimensionless pressure of all carbon-containing gas phase components in the domain will be zero, as will be the level of carbon impurities in the hydrogen product generated during OM 3. Similarly, the duration of OM 3 is selected as the time τ 3 at which the dimensionless partial pressure of hydrogen in domain (v) at the reactor exit falls below a set value, e.g., P H 2 ,gv (τ 1 + τ 2 + τ 3 , 1) ≤ 0.01, which ensures that most of the hydrogen in both domains (v) and (s) has been removed and thus the reactor is largely filled with unreacting steam.
Simulation of the dimensionless Equation (10) of the MSR model requires that a number of model parameters first be specified. It is considered that the reactor pressure and temperature are P * = 25 bar and T = 900K, respectively, and the MSR feed at the beginning of OM 1 contains steam and methane at a molar ratio of 3, and minute amounts of all other species. Knowledge of these parameters enables the computation of the reference reaction generation rate r * using Equation (19), and expression of the dimensionless reaction rates R 1 ,R 2 ,R 3 in terms of the dimensionless species pressures, using Equations (20)- (22). Selection of the storage domain, catalyst domain, and gas phase in storage domain volume fractions as ε s = 0.5, ε c = 0.35, ε gs = 0.495, respectively, of the effectiveness factor as η = 0.5, and of the gas permeance ratios for all species (the first species being hydrogen) In this work, assessment of SMR-MSR performance was performed by carrying out sixteen combinations (trials) of the above five parameters, as listed in Table 3 below. The values for D a and Θ were chosen by selecting a range of reasonable values for the parameters shown in Equation (8). Fixing the operating pressure P * = 25(bar), temperature T * = 900(K), catalyst density ρ c = 2355.2 kg m 3 [14], and reference reaction rate , T = 13.606 mol kg·s , leaves residence time t * (s) as the only free parameter for calculation of D a . Because of the rapid rate at which the SMR reactions attain equilibrium, the range of considered residence times was 0.01 ≤ t * (s) ≤ 0.5, so that the concentration profiles along the reactor do not immediately reach equilibrium. The resulting D a range is 1 ≤ D a ≤ 6. Given the above range of residence times 0.01 ≤ t * (s) ≤ 0.5, and pressure P * = 25(bar), temperature T * = 900(K) values, evaluation of a range for Θ requires range information for the species membrane permeance β i The β 1 range is 1.05 · 10 −8 ≤ β 1 ≤ 1.21 · 10 −6 for hydrogen [23,33], and the α s,v range is selected so that the resulting Θ range is 1 ≤ Θ ≤ 50. Each trial was simulated in COMSOL for three cycles of operation in order to ensure that dynamic operation had reached its long-term periodic behavior. COMSOL implementation of the alternating feed compositions for each OM was carried out by approximating for each OM an OM-specific Heaviside function h k with an OM-specific Verhulst function, W k , as shown in Equation (27) below: Figure 2 below illustrates the dynamic behavior of the methane inlet partial pressure at the reactor inlet as a function of dimensionless time, considering τ 1 = τ 2 = τ 3 = 5, a = 1 and c = 8.5. In Figure 3 the time evolution of all species' (except, of course, for the hydrogen in the storage domain) dimensionless partial pressures at the reactor inlet is shown, over three operating cycles. It is assumed that the conditions at the inlet of the reactor at the beginning of OM 1 is the same as its condition at the end of OM 3. During OM 1, the reactor's feed consists mainly of steam and methane, at a 3:1 steam/methane ratio, and trace amounts of carbon monoxide and carbon dioxide. During OM 2 and OM 3 the reactor's feed consists only of steam. Similarly, Figure 4 shows the time evolution of species dimensionless partial pressures at the outlet of the reactor through the same three cycles of operation. Because the reactor's condition at the beginning of OM 1 is the same as its condition at the end of OM 3, at early operational times the reactor outlet consists mainly of steam. As the reaction proceeds, a mixture of methane, carbon dioxide, and carbon monoxide begins to exit the reactor. The hydrogen in the void domain continues to rise until it reaches its peak, which coincides with the end of OM 2. In Figure 3 the time evolution of all species' (except, of course, for the hydrogen in the storage domain) dimensionless partial pressures at the reactor inlet is shown, over three operating cycles. It is assumed that the conditions at the inlet of the reactor at the beginning of OM 1 is the same as its condition at the end of OM 3. During OM 1, the reactor's feed consists mainly of steam and methane, at a 3:1 steam/methane ratio, and trace amounts of carbon monoxide and carbon dioxide. During OM 2 and OM 3 the reactor's feed consists only of steam. Similarly, Figure 4 shows the time evolution of species dimensionless partial pressures at the outlet of the reactor through the same three cycles of operation. Because the reactor's condition at the beginning of OM 1 is the same as its condition at the end of OM 3, at early operational times the reactor outlet consists mainly of steam. As the reaction proceeds, a mixture of methane, carbon dioxide, and carbon monoxide begins to exit the reactor. The hydrogen in the void domain continues to rise until it reaches its peak, which coincides with the end of OM 2.   Figure 5 shows the time evolution of the dimensionless partial pressure of hydrogen in the void and storage (represented by R H2 and S H2) domains, during the last cycle of operation, at four different reactor locations (25%, 50%, 75%, and 100% of the total reactor length), and Figure 6 shows the dimensionless partial pressures of methane and carbon     Figure 5 shows the time evolution of the dimensionless partial pressure of hydrogen in the void and storage (represented by R H2 and S H2) domains, during the last cycle of operation, at four different reactor locations (25%, 50%, 75%, and 100% of the total reactor length), and Figure 6 shows the dimensionless partial pressures of methane and carbon  Figure 5 shows the time evolution of the dimensionless partial pressure of hydrogen in the void and storage (represented by R H 2 and S H 2 ) domains, during the last cycle of operation, at four different reactor locations (25%, 50%, 75%, and 100% of the total reactor length), and Figure 6 shows the dimensionless partial pressures of methane and carbon dioxide in the void domain, during the last cycle of operation, at the same four reactor locations. In Figure 5 it can be seen that at the 25% length location the hydrogen concentration in the void domain is initially practically zero; then, as methane reacts to generate hydrogen, it increases, first rapidly and then less so as it begins to transport into the storage domain; finally, the concentration decreases as steam is fed into the reactor. The same behavior is observed at larger lengths. locations. In Figure 5 it can be seen that at the 25% length location the hydrogen concentration in the void domain is initially practically zero; then, as methane reacts to generate hydrogen, it increases, first rapidly and then less so as it begins to transport into the storage domain; finally, the concentration decreases as steam is fed into the reactor. The same behavior is observed at larger lengths.  In Figure 6, at the 25% reactor length location during the last cycle of operation, the methane and carbon dioxide concentrations are shown to first increase, reach a maximum, and then decrease, in correspondence with the composition changes at the inlet of the reactor shown in Figure 3. The methane (carbon dioxide) maximum occurs at the latter (earlier) part of a plateau-like region in time. At points further down the reactor (50%, 75%, and 100% of the reactor length), the methane (carbon dioxide) pulse exhibits a progressively lower (higher) maximum, which is indicative of its role as a reactant (product) in the reaction that takes place. The magnitude of the methane pulse remains larger than that of the carbon dioxide pulse until it reaches the end of the reactor. In the Supplementary Materials section, Figures S1-S3 (Figures S4-S6) show the dimensionless pressure profiles for all species vs. time for three reactor locations (25%, 50%, and 75% of the reactor length), for the last cycle of the 7th (13th) trial D a = 2, Θ = 10 (D a = 4, Θ = 50), which are qualitatively similar to the behavior shown in in Figures 4 and 5 above.    Figures 7-9 show the dimensionless partial pressure profiles for all species along the length of the reactor at the end of OM 1, OM 2, and OM 3, respectively, for the 7th trial D a = 2, Θ = 10. In these simulations we selected the end of OM 1 to be when P CH 4 ,gv z=0.3 = 0.05 occurs for an operating time for the first OM of τ 1 = 0.566. From Figure 7 one observes that the contents of the reactor void domain up to point (z * = 0.3) contain a mixture of reforming products, whereas at larger reactor lengths there exists only water, validating that the reactor is operating in plug flow and that no species diffusion has occurred down the reactor. From the reactor inlet (z * = 0) to (z * = 0.3), the dimensionless partial pressure of the hydrogen in the storage shows a marked difference from the dimensionless partial pressure of hydrogen in the reactor void, which highlights the time-dependent nature of the process. As seen later in the analysis of the influence of the parameter Θ, reducing this difference during early operational times correlates to higher conversion and better overall reactor performance. Figures 8 and 9 show the corresponding profiles at the end of OM 2 with τ 2 =1.00, and at the end of OM 3 with τ 3 = 2.51, respectively. From Figure 8 one sees that at the end of OM 2 the reactor void contains a mixture of steam and hydrogen only, and that hydrogen has been successfully transferred to and stored in the storage domain. From Figure 9 one sees that, at the end of a cycle of operation, the reactor has returned to its initial state of operation, with the reactor void being filled with only steam and the contents of the storage medium having been emptied.     With the operating times for each OM set, it is now possible to calculate the metrics as described by Equations (23)- (26). These are shown for the different cases studied in Table 4. The SMR-MSR shows higher overall conversion than the traditional SMR-SSR for all cases studied. In addition, the SMR-MSR enables the in situ separation of the hydrogen product, which is a distinct advantage over the traditional SMR-SSR that requires an additional downstream separation step. For the design parameters employed in trial 7 (Da = 2, Θ = 10) for which the SMR-MSR behavior has been described in detail above, the SMR-MSR exhibits a 116% increase in methane conversion over the SMR-SSR from 0.2564 to 0.5553. Additionally, the hydrogen recovery is 83%, indicating that a significant portion of the total generated hydrogen exits the reactor during OM 3 in the form of ultra-high purity (UHP) hydrogen. For the case studies conducted here (see Table 4), a maximum increase in methane conversion of 235% was observed, indicating the significant benefit that can be gained by appropriately selecting design parameters. The substantial increases in conversion are accompanied by simultaneous significant gains in hydrogen yield. The hydrogen recovery ratios reported in Table 4 are comparatively higher than those by similar membrane reactors operating at steady state [34]. These high SMR-MSR metric values suggest that the SMR-MSR outperforms the SMR-SSR, and that the SMR-MSR operating at the selected temperature 900 T K = (which is lower than the temperatures traditionally used in SMR [35]), may be economically viable and realizable [36]. With the operating times for each OM set, it is now possible to calculate the metrics as described by Equations (23)- (26). These are shown for the different cases studied in Table 4. The SMR-MSR shows higher overall conversion than the traditional SMR-SSR for all cases studied. In addition, the SMR-MSR enables the in situ separation of the hydrogen product, which is a distinct advantage over the traditional SMR-SSR that requires an additional downstream separation step. For the design parameters employed in trial 7 (Da = 2, Θ = 10) for which the SMR-MSR behavior has been described in detail above, the SMR-MSR exhibits a 116% increase in methane conversion over the SMR-SSR from 0.2564 to 0.5553. Additionally, the hydrogen recovery is 83%, indicating that a significant portion of the total generated hydrogen exits the reactor during OM 3 in the form of ultra-high purity (UHP) hydrogen. For the case studies conducted here (see Table 4), a maximum increase in methane conversion of 235% was observed, indicating the significant benefit that can be gained by appropriately selecting design parameters. The substantial increases in conversion are accompanied by simultaneous significant gains in hydrogen yield. The hydrogen recovery ratios reported in Table 4 are comparatively higher than those by similar membrane reactors operating at steady state [34]. These high SMR-MSR metric values suggest that the SMR-MSR outperforms the SMR-SSR, and that the SMR-MSR operating at the selected temperature T = 900 K (which is lower than the temperatures traditionally used in SMR [35]), may be economically viable and realizable [36]. The effect of Θ (the inverse Peclet number) on reactor performance is also shown in Table 3 (operational times) and 4 (performance metrics), with an increase in all metrics shown for increasing Θ. The cases with the lowest Θ value (Θ = 1) studied are characterized by minimal increases in methane conversion and hydrogen yield over their traditional reactor counterparts, and low recovery ratios, suggesting a lower limit threshold of Θ > 1 for desirable SMR-SSR performance. For a constant value of Da (and thus a constant value of the SMR-MSR dimensionless residence time), as Θ increases, the operating time τ 1 in OM 1 can be increased, because reactor conversion then increases as Θ increases and it takes a longer time for the unreacted methane to reach the outlet of the reactor in substantial amounts. In contrast, operational times τ 3 in OM 3 decrease upon increase in Θ for each fixed Da. These beneficial changes in OM duration times, coupled with the increase in the hydrogen recovery ratio as Θ increases (for fixed value of Da), suggest that Θ should be increased as much as possible in order to achieve desirable reactor performance. Figures 10-13 show corresponding results for the 13th trial (Da = 4 and Θ = 50). Figure 10 shows the time evolution of the dimensionless partial pressure for the various species at the outlet of the reactor through three cycles of operation. Figure 11 shows the time evolution of the dimensionless hydrogen partial pressure in the void and storage domains during the last cycle of operation at four locations in the reactor (25%, 50%, 75%, and 100% of the total reactor length). Figure 12 shows the time evolution of dimensionless partial pressures for methane and carbon dioxide in the void domain during the last cycle of operation at the same four reactor locations. Figure 13 shows the dimensionless partial pressure profiles for all species along the length of the reactor during cycle 3 at the end of OM 1. Figures S4-S6       Comparing Figure 10 with Figure 4, we observe a large increase (decrease) in the amount of carbon dioxide (methane) exiting the reactor. Figures 5 and 11 show similar patterns of behavior for hydrogen in the reactor and storage void domains, with the hydrogen in the storage medium exhibiting a noticeably shorter response time in Figure 11. This is attributed to the higher Θ value, which allows the hydrogen to be more easily  Comparing Figure 10 with Figure 4, we observe a large increase (decrease) in the amount of carbon dioxide (methane) exiting the reactor. Figures 5 and 11 show similar patterns of behavior for hydrogen in the reactor and storage void domains, with the hydrogen in the storage medium exhibiting a noticeably shorter response time in Figure 11. This is attributed to the higher Θ value, which allows the hydrogen to be more easily Comparing Figure 10 with Figure 4, we observe a large increase (decrease) in the amount of carbon dioxide (methane) exiting the reactor. Figures 5 and 11 show similar patterns of behavior for hydrogen in the reactor and storage void domains, with the hydrogen in the storage medium exhibiting a noticeably shorter response time in Figure 11. This is attributed to the higher Θ value, which allows the hydrogen to be more easily transferred between domains. Comparing Figures 6 and 12, we see a dramatic change in the behavior for both methane and carbon dioxide. In Figure 6, the methane waveform retains its three distinct phases through all measured points, and is also always greater in magnitude until the exit of the reactor. This is not the case for Da = 4 and Θ = 50, where the magnitude of the carbon dioxide pulse becomes greater than that of the methane pulse at the 50% location. This behavior is also seen in Figures S4-S6, which show (for the case with Da = 4 and Θ = 50) the dimensionless partial pressures for all species during the last cycle of operation at locations equal to 25%, 50%, and 75% of reactor length, respectively. This is particularly highlighted in Figure 13, which shows the species profiles in the reactor at the end of OM 1. When compared with Figure 7, it can also be seen in Figure 13 that the dimensionless partial pressures of hydrogen in the storage and reactor void are much closer in value, indicating a much higher rate of transfer of separated hydrogen.
From a process design and operational perspective, there are several parameters that can be adjusted to alter Da and Θ. For a fixed inlet composition at the beginning of OM 1, the Da can be increased (decreased) by increasing (decreasing) the reactor's residence time, which in turn can be achieved by increasing (decreasing) the reactor's length and/or by decreasing (increasing) the reactor's residence time. The dimensionless parameter Θ can be increased (decreased) by increasing (decreasing) the residence time, storage-void domain interfacial area and hydrogen permeance, and/or by decreasing (increasing) the void domain volume fraction. Θ can be altered even while keeping the Da fixed, by altering the last three aforementioned design parameters. In particular, increasing the preferential hydrogen permeance through the storage medium's permselective layer will increase Θ and can be accomplished through appropriate selection of the layer's pore structure and material properties.

Conclusions
In this work, a novel, first principle-based, spatially dependent model capturing the membrane storage reactor behavior was presented and simulated. This novel intensified process, termed the MSR, is capable of overcoming operational limitations of traditional reactor design by combining multiple mass transport processes into a single unit, and increases the production rate through dynamic operation. The resulting desired products are delivered at high pressure and in readily separable form, which leads to additional energetic and operational cost saving. The assessment of MSR behavior was described by a one-dimensional dynamic model which, when cast in dimensionless form, revealed two dimensionless parameters, Da (the Damkohler number) and Θ = 1/Pe mem (the inverse Peclet number), that capture MSR behavior. Performance analysis of the MSR process was undertaken by introducing a number of metrics, and a case study on Steam Methane Reforming for the production of hydrogen was subsequently carried out. Numerous simulations of the developed spatially dependent dynamic model enabled parametric studies of the effect of the aforementioned Da and Θ = 1/Pe mem dimensionless groups, and established that maximizing both groups leads to improved MSR performance. A comparative analysis between SMR-MSR and SSR showed that the SMR-MSR process obtained a higher methane conversion, X CH 4 , and a greater yield of hydrogen, Ω H 2 . It was also shown that the SMR-MSR hydrogen recovery ratio, R H 2 ,3 , was comparable to or above those obtained by currently investigated membrane technologies.
There are several notable advantages of the proposed MSR process, the first being its ability to be modularized and to accommodate a variety of economic production scales. Although the current case study emphasized application to hydrogen production for use at the level of a refinery or a large plant, smaller scales of production are also quite feasible. The MSR process can conceivably be applied to hydrogen-fueling stations so as to facilitate the creation of sustainable decentralized hydrogen generation. Additionally, implementation of the MSR process into existing plants would help reduce many of the economic barriers faced by industrial operations when converting to new technology. Using SMR as an example, retrofitting of an existing plant would consist mainly of reactor repacking with catalyst and storage media. Because of the need for the MSR process to direct material flow, it is imagined that additional fluid lines (for both feed and effluent) need to be constructed. This additional cost, however, is not expected to be large, especially when considering that no additional exogenous chemical components are required. Additional component cost from MSR implementation can also be expected from the need to incorporate dynamic control equipment because, as noted earlier, most industrial reactors operate at steady state. All of this suggests that the MSR process has the potential to overcome many of the barriers that hinder many potential process intensification technologies and be implemented at an industrial scale.
The presented one-dimensional dynamic model was used to demonstrate the novel MSR concept. Because the MSR process is shown to offer significant performance advantages, future work will attempt to rigorously quantify energetic savings by incorporating energy balances into the MSR model. Additionally, development of a MSR model that can capture multi-scale effects and help quantify storage media structural parameters will be investigated.