Comparison of Finite Difference and Finite Volume Simulations for a Sc-Drying Mass Transport Model

Different numerical solutions of a previously developed mass transport model for supercritical drying of aerogel particles in a packed bed [Part 1: Selmer et al. 2018, Part 2: Selmer et al. 2019] are compared. Two finite difference discretizations and a finite volume method were used. The finite volume method showed a higher overall accuracy, in the form of lower overall Euclidean norm (l2) and maximum norm (l∞) errors, as well as lower mole balance errors compared to the finite difference methods. Additionally, the finite volume method was more efficient when the condition numbers of the linear systems to be solved were considered. In case of fine grids, the computation time of the finite difference methods was slightly faster but for 16 or fewer nodes the finite volume method was superior. Overall, the finite volume method is preferable for the numerical solution of the described drying model for aerogel particles in a packed bed.


Introduction
Supercritical drying of wet gels, such that they become porous aerogels, is the crucial step in aerogel production. The drying process uses supercritical CO 2 as extraction medium and is; thus, conducted under moderate temperature but elevated pressure and; therefore, may be cost-intensive. A comprehensive overview of experimental studies and mass transport models investigating the kinetics of supercritical drying of gels is given byŞahin et al. [1]. Within the last two years, a several model of the supercritical drying for spherical gel particles were developed [2][3][4][5][6][7]. All these works show that the gel drying process for particles is significantly faster than for monoliths of comparable size [2][3][4][5][6][7]. This observation highlights the fact that aerogel particles are highly promising for future applications in industrial products. Selmer et al. developed their physical mass transport model for the supercritical drying of gel particles in a packed bed and analyzed mass transfer steps that limit the overall drying kinetic to optimize the process [2,3]. Hatami et al. developed a similar model and improved the optimizing procedure to allow a fast drying process with low CO 2 consumption [6].
All reported mass transport models result in a set of coupled partial differential equations if the mass transport in the outer bulk fluid is considered additionally to the mass transport within the gel particles [3][4][5][6].Şahin et al. and Hatami et al. discretized the spatial dimensions via the finite difference method and solved the resulting coupled set of differential-algebraic equations using the built-in ordinary differential equation solver in Matlab [4][5][6]. Selmer et al. discretized the spatial dimension of the gel particles via the finite difference method [8,9] and the spatial dimension of the bulk fluid in Gels 2020, 6, 45 2 of 26 the autoclave via the finite volume method [3,10]. The finite volume method can preserve discretely important physical quantities such as mass [11].
During the supercritical drying process, the physical properties of supercritical CO 2 and its mixture with ethanol or other solvents strongly depend on the composition of the mixture which varies during the drying process from pure ethanol (or another solvent) in the gel to almost pure CO 2 [3]. Therefore, especially in mass transport models describing the supercritical drying process of small particles, a high spatial resolution of the domain with gel particles is required for accurate results. This, in turn, leads to high computation times.
For this reason, we set the aim of this work so as to investigate the behavior of numerical solutions of gel particles drying models according to their accuracy and efficiency, to increase the confidence in the simulation results reported so far.
In the present paper, two possible numerical discretization techniques of the resulting partial differential equations of Selmer et al.'s model [2,3] are implemented, namely the finite difference method [8,9] and the finite volume method [10]. Both discretization techniques are compared numerically by studying the quality of the approximations via Euclidean norm (l 2 ) and maximum norm (l ∞ ) errors, as well as the efficiency via the condition numbers of the linear systems to be solved.
The supercritical drying of spherical, ethanol-filled silica gels packed in a cylindrical bed serves as the example system. Two initial ethanol concentrations in the bulk fluid, which represent the pressurization using CO 2 with or without addition of excess solvent to the wet gels, are evaluated. Both pressurization techniques are applicable and implemented in supercritical drying processes.

Mass Transport Model of Supercritical Drying
The aim of supercritical drying is to extract a solute (in our case ethanol) from a wet gel using supercritical CO 2 (component 1). The mass transport model presented here, and discussed in Selmer et al. 2018 from the physical point of view [2], results in a system of coupled partial differential equations. The partial differential equations are written here in terms of concentration of ethanol (component 2), c 2 . In the drying process, ethanol is transported by diffusion from the inside of the porous gel particles (Equation (1)) to the bulk fluid ( Figure 1). The flow of the bulk fluid through the packed bed is modeled by convection (Equation (2)) ( Figure 1). The packed bed leads to non-ideal flow behavior, which is included through a time dependent dispersion factor D L (Equation (2)). Thus, two domains are considered: The spherical gel particles (subscript g), described by radial mass transport along the radial axis r; and the bulk fluid (subscript f ), described by axial mass transport along the axial axis z. Both domains are connected due to the mass transfer of ethanol across the boundary. This connection appears as a boundary condition for the mass transport within the particle (Equation (8)) and as a source term source 2, f for the mass transport of the bulk fluid (Equation (2)). It is assumed in our model that the gel particles are monodisperse and evenly distributed along the cylindrical autoclave. Since the mass transport within the gel particles depends on the mass transport of the bulk fluid, which in its turn varies along the axial packed bed height z, the mass transport within the gel particles should also depend on their position z in the packed bed. This fact is explicitly stated in (Equation (1)). The drying process is a dynamic process and thus depends on the time denoted with t. The source term source 2, f describes the ethanol being extracted from the gel particles (at axial position z in the packed bed) into the bulk fluid. The term is calculated as the temporal change of the amount of ethanol (in moles) within all gel particles divided by the accessible volume of the bulk fluid V f = A ac ·ψ·L (Equation (4)). A ac stands for the cross-sectional area of the autoclave, ψ for the porosity of the packed bed and L for the length of the packed bed. The temporal change of the ethanol amount within the porous gel sphere is negative during the course of the drying since ethanol concentration with time decreases, but the source term itself is positive. 3.09 • 10 and 7.18 • 10 • for the viscosity of the bulk fluid and between 2.34 and 1.66 • 10 for the resulting Schmidt number.
Parameters that are independent from the mixture composition are written in bold font throughout the text and as follows: The surface of a single spherical gel particle , the number of particles in the packed bed , the cross-sectional area of the autoclave , the length of the packed bed , the porosity of the packed bed and the mass flow rate of the ethanol-CO2 mixture . radial coordinate of spherical gel particle, r = 0-gel particle center, r = R-gel particle radius, ttime). V-volume flow, ρ f -fluid density (mixture of CO 2 and ethanol), r-radial coordinate of spherical gel particle, r = 0-gel particle center, r = R-gel particle radius, t-time).
Physical properties of supercritical CO 2 and its mixture with ethanol or other solvents strongly depend on the system parameters such as temperature, pressure and composition of the mixture [2]. The temperature and pressure are treated as constant during the drying process. The composition of binary mixture x 2 is directly related to the mixture concentration c and ethanol concentration c 2 (Equation (5)) and varies during the drying from pure ethanol (x 2 = 1) to pure CO 2 (x 2 = 0). Physical properties, such as the effective diffusion coefficient D g , the mixture concentration c g and the fluid density ρ f , as well as the parameters that are functions thereof, such as the interstitial fluid velocity u (Equation (3)) and the mass transfer coefficient β (Equation (8)), are highly dependent on the radial, axial and temporal distribution of the mixture compositions in the gel body x 2,g (r, z, t) and the bulk fluid x 2, f (z, t). For example, at a pressure of 10 MPa and a temperature of 321 K, the parameters being directly dependent from the mixture composition can vary during supercritical drying between 1.41·10 −9 and 7.68·10 −9 m 2 s for the effective diffusion coefficient D g (assumption: Gel porosity ε g = 0.93, gel tortuosity τ g = 3.48), 9.67·10 3 and 1.77·10 4 mol m 3 for the mixture concentration c g and the fluid density ρ f , 5.54·10 −9 and 2.87·10 −8 m 2 s for the diffusion coefficient in the bulk fluid, 3.09·10 −5 and 7.18·10 −4 Pa·s for the viscosity of the bulk fluid and between 2.34 and 1.66·10 2 for the resulting Schmidt number.
Parameters that are independent from the mixture composition are written in bold font throughout the text and as follows: The surface of a single spherical gel particle A p , the number of particles in the packed bed N p , the cross-sectional area of the autoclave A ac , the length of the packed bed L, the porosity of the packed bed ψ and the mass flow rate of the ethanol-CO 2 mixture . m.
Initial and boundary conditions are presented in Equations (6)-(8) for the gel particle domain and in Equations (10)- (12) for the bulk fluid domain. For each gel sphere (at a certain packed bed height z), the boundary conditions were given by the Neumann condition at the center (Equation (7)) and by the molar flux across the boundary layer at the surface (Equation (8)).

∀r, ∀z,
For the top of the packed bed, the Danckwerts' boundary condition [12] was chosen so that no loss of ethanol across the upper boundary occurs due to dispersion of ethanol (Equation (11)). It was useful to define the flux of ethanol f 2 as follows: In order to simplify the boundary condition at the top of the packed bed: The flux of ethanol across the upper boundary is set to zero (Equation (11)) according to the Danckwerts' boundary condition [12].
At the bottom of the packed bed, a Neumann boundary condition was used (Equation (12)).

Numerical Solution of the Mass Transport Model
The coupled partial differential equations were solved numerically. Equation (1), which describes the mass transport within the porous gel particles, was always solved using the finite difference method (FDM), whereas Equation (2), which describes the mass transport of the bulk fluid in the packed bed, was solved using either the FDM (Section 3.1) or the finite volume method (FVM) (Section 3.2), (Supplementary Materials). Equation (1) was discretized using an explicit scheme, since it is easy to implement. Without its dependencies on r, z , t, Equation (1) is written as follows: ∀z ∂c 2,g ∂t = 1 r 2 ∂ ∂r r 2 c g D g ∂x 2,g ∂r (13) which gives, after applying the product rule, ∀z ∂c 2,g ∂t = 2c g D g r ∂x 2,g ∂r + D g ∂c g ∂r Equation (14) was discretized using a forward difference for the time derivative and central differences for the spatial derivatives where the super index k refers to time, s to the axial index in the packed bed/autoclave domain (cf. Figure 2), and n to the radial index of the gel particle domain. N stands for the number of nodes used to discretize the gel particle domain. The temporal as well as both spatial discretizations were based on equidistant nodes resulting in constant time ∆t = t k − t k−1 and constant radial steps ∆r = r n − r n−1 . We can rearrange Equation (15) to n = 1, . . . , N ∀s, ∀k c k+1 2,g,n,s = c k 2,g,n,s + ∆t ∆r 2 1 − 1 n x k 2,g,n−1,s − 2x k 2,g,n,s + 1 + 1 n x k 2,g,n+1,s D k g,n,s c k g,n,s + 1 4 x k 2,g,n−1,s − x k 2,g,n+1,s D k g,n−1,s − D k g,n+1,s c k g,n,s + c k g,n−1,s − c k g,n+1,s D k g,n,s .
The initial condition of the gel particle domain was set to ∀n, ∀s, k = 1 c 0 2,g,n,s = c 2, g,start .
The boundary conditions of the gel particle domain were discretized using central differences for the spatial derivatives. Two artificial nodes (n = 0, n = N + 1) were placed beyond the boundaries n = 1 and n = N and calculated as follows: n = 0, ∀s, ∀k c k 2,g,0,s = c k 2,g,2,s , x k 2,g,0,s = x k 2,g,2,s , c k g,0,s = c k g,2,s , D k g,0,s = D k g,2,s The artificial node n = 0 was used to solve Equation (16) for c k+1 2,g,1,s and the artificial node n = N + 1 for c k+1 2,g,N,s . The domain which describes the mass transport in the bulk fluid (Equation (2)) was solved with both numerical methods (Sections 3.1 and 3.2). Without its dependencies on r, z, t, Equation (2) was written as follows: Using this definition of the flux (Equation (9)) allows Equation (20) to be rewritten as The source term source 2, f (Equation (4)), which describes the ethanol being extracted from porous gel particles to the bulk fluid in the packed bed with time, was discretized for both numerical methods in the same way (Equation (22)): It was calculated as the difference between the accumulated ethanol within a single porous gel particle (placed at axial position s ) at time t = k and at time t = k + 1, divided by the time step ∆t = t k+1 − t k and multiplied by the particles per bulk fluid volume (22)). The accumulated ethanol within a single porous particle was calculated by the summation of the ethanol concentration c k 2,g,n,s in each spherical shell element n, times the volume of each spherical shell element V g,n times the gel porosity ε g (Equation (22)): c k 2,g,n,s ·V g,n ·ε g − N n=1 c k+1 2,g,n,s ·V g,n ·ε g .
Gels 2020, 6, x FOR PEER REVIEW 6 of 31 ethanol within a single porous gel particle (placed at axial position ) at time = and at time = + 1, divided by the time step Δ = − and multiplied by the particles per bulk fluid volume = • • (Equation (22)). The accumulated ethanol within a single porous particle was calculated by the summation of the ethanol concentration , , , in each spherical shell element , times the volume of each spherical shell element , times the gel porosity (Equation (22)):

Solving the Diffusion-Advection Equation Using the Finite Difference Method
Finite difference methods are a flexible tool to discretize partial differential equations and easy to implement. Figure 2 shows the equidistant axial discretization of the bulk fluid/packed bed domain for the FDM and the FVM (Section 3.2). For FDM, the first discretization node was placed at = 0 and the last one at = . Therefore, both boundary volume elements/cells (grey boxes) were half of the length of the inner volume elements: Additional to the number of nodes S in the packed bed domain, two more artificial nodes ( and ) outside the domain were needed to discretize the boundary conditions via the FDMs. In the following, we will distinguish between two versions of possible discretizations using the FDM to solve numerically the diffusion-advection equation. The difference is given by the discretization of the advective term in Equation (20).
x 2 x s x s-1 x s+1 x 2 x 3/2 Δz s x 3/2 Figure 2. Axial discretizations of the packed bed for FDM and FVM. S-number of nodes, x s -s th node (blue), x 0 , x S+1 -artificial node (green), ∆x-distance between neighboring nodes, x s+ 1 2 -cell interface (in the middle of two neighboring nodes), ∆z s -length of volume element (grey boxes, boundary boxes in dark grey).

Solving the Diffusion-Advection Equation Using the Finite Difference Method
Finite difference methods are a flexible tool to discretize partial differential equations and easy to implement. Gels 2020, 6, 45 7 of 26 Figure 2 shows the equidistant axial discretization of the bulk fluid/packed bed domain for the FDM and the FVM (Section 3.2). For FDM, the first discretization node was placed at z = 0 and the last one at z = L. Therefore, both boundary volume elements/cells (grey boxes) were half of the length of the inner volume elements: Additional to the number of nodes S in the packed bed domain, two more artificial nodes (x 0 and x S+1 ) outside the domain were needed to discretize the boundary conditions via the FDMs.
In the following, we will distinguish between two versions of possible discretizations using the FDM to solve numerically the diffusion-advection equation. The difference is given by the discretization of the advective term in Equation (20).

Discretization of Advection Term: Version A
In version A, the advective term − ∂(c 2, f ·u) ∂z was discretized. The discretization of Equation (20) using an explicit forward difference in time, a first order backward difference and a second order central difference in space results into

Discretization of Advection Term: Version B
In version B, the advective term −u· . Thus, the discretization of Equation (20) using an explicit forward difference in time, a first order backward difference and a second order central difference in space is given by The corresponding discretized initial and boundary conditions for versions A and B are Equations (24) and (25) can be rearranged to a linear system (Equations (29)-(33)) using the discretized boundary conditions (Equations (27) and (28)) and corresponding cell sizes ∆z s (Equation (23)).
∀s, ∀k c k+1 Gels 2020, 6, 45 Version A: Where a FDM,s is nonzero only for the indices s − 1, s, s + 1 and given by Version A:

Solving the Diffusion-Advection Equation Using the Finite Volume Method
Under certain conditions, finite volume methods conserve mass discretely and thus are suitable here. In the FVM, the ethanol flux f 2 (Equation (9)) is evaluated at the cell interfaces x s+ 1 2 and the ethanol concentration c 2, f at the nodes x s (Figure 2, Equation (35)). To discretize the boundary conditions appropriately (cf. Equations (42)-(44)), the first discretization node was placed at z = x 1 and the last node at z = L . Thus, the bottom boundary volume element was halved ( Figure 2): The grid definition of the FVM leads to a slightly different grid compared to the FDMs (including different ∆z s ) using the same number of nodes S. This is important for the subsequent comparison of both methods.
To solve Equation (21), we adapted the time-invariant complete flux scheme version developed by Farrell and Linke [11], which is based on Voronoï meshes [13]. They considered a numerical scheme for stationary problems, which we extended to time-dependent problems via integration: Equation (21) was integrated on a discrete cell ( Figure 2) at time k yielding Equation (35). Here f k 2,s+ 1 2 describes the molar flux which crosses the interface between two neighboring cells. It can be divided into a homogeneous part and an inhomogeneous part (Equation (36)), which were calculated using Equation (37) to Equation (41) according to Farrell and Linke [11]. Their idea is based on complete flux schemes [14][15][16]. Thus, we used the supercritical drying model of gel particles in a packed bed as application example for the complete flux scheme developed by Farrell and Linke [11].
The boundary conditions of the FVM (Equations (42)-(44)) were chosen according to Equations (10)- (12). The difference to the FDM (Equations (26)-(28)) is that the boundary conditions are written in terms of ethanol fluxes f k 2,s . Equation (43) implies that no ethanol crosses the upper boundary which is analogous to the Danckwerts' boundary condition [12]. At the bottom and at the end of the packed bed the ethanol flux is only influenced by the transported ethanol in the bulk fluid (pure convection without any dispersion) (Equation (44)), since the spatial derivative of the ethanol concentration in the bulk fluid c 2, f stays zero at z = L (Equation (12)).

Results and Discussion
First, we numerically analyze the accuracy of the presented discrete schemes by studying the convergence behavior and the closure of the mole balance. Second, we discuss the computational efficiency of all three methods by examining the condition number as well as the duration of the computation.
The following drying conditions were utilized for all calculations: System pressure P = 10 MPa, system temperature T = 321 K, gel sphere radius R = 3.175 mm, gel porosity ε g = 0.93, gel tortuosity τ g = 3.48, ethanol molar fraction within the gel sphere at the start of the drying x start 2,g = 1, ethanol molar fraction within the gel sphere at the end of the drying (depends on system temperature and defines the drying end) From Figure 3 it is obvious that increasing the total number of nodes S increases the quality of the solutions for all three methods. For S = 2 and x 2, f ,start = 0.05, the FVM overshoots slightly the reference solutions in the first drying minutes due to a sharp increase of the molar fraction x 2, f ,S at this position (bottom of the packed) in the beginning of the drying process. The FDM version B prolongs the drying process compared to the reference solutions.   In the previous paper [3], 20 nodes were taken for the same calculation using the presented FVM and x 2, f ,start = 0.05 as initial condition. From Figure 3 it is obvious that these 20 nodes (marked as In a second step, the solutions of the numerical methods were analyzed not only at the final node s = S (bottom of the packed bed) for varying drying times (as presented above), but also for all other nodes ( s = 1, . . . , S ) in order to have a spatial error distribution ( Figure 5).     In a third step, the overall errors calculated for all spatial nodes of the packed bed height ( = 1, … , ) and for all drying time steps ( = 1, … , ) are presented in Figure 6 for different discretizations/grids ( = 2; 4; 8; 16; 32 ) and both investigated initial conditions. It is obvious that the overall errors are higher for the initial condition , , = 0.95 compared to the calculations using , , = 0.05 . For both initial conditions and both errors, and , the FVM shows the best convergence behavior, meaning the error decreases the fastest for an increasing total number of grid points (or smaller grid sizes ∆ ). To be able to compare the molar fractions on coarse meshes with the fine reference solution x REF 2,f , the molar fractions on the finer reference mesh were interpolated to the coarse grid. The deviations from the reference solution were calculated for each time step k = 1, . . . , K and for each spatial step s = 1, . . . , S. Thus, the following mole fraction vector was defined to calculate the l ∞ and l 2 errors: x 2,f := x 2,f,1 , . . . , x 2 Figure 5 shows the resulting l ∞ and l 2 errors calculated at each single node s = 1, . . . , 16 for all time steps k (as example for a total amount of S = 16 nodes). The l ∞ error depicts the maximum value of the absolute error, the l 2 error corresponds to the l 2 norm of the error.
As observed already in Figure 4, the highest maximum absolute error (l ∞ error) at the last node s = S = 16 for the case S = 16 nodes and x 2, f ,start = 0.05 was given by the FVM due to the approximated boundary condition at z = L. Nevertheless, the FVM shows small l ∞ errors for all other nodes (similar to FDM version A) for both initial conditions ( Figure 5). Overall, the smallest l 2 errors for all nodes (for all time steps and both initial conditions) were identified for the FVM. Version A of the FDMs shows generally smaller l ∞ and l 2 errors compared to Version B (as already seen for the single node s = S at the bottom of the packed bed (Figure 4)). For both errors, l ∞ and l 2 , and all three numerical methods, the errors are generally higher for the initial condition x 2, f ,start = 0.95 compared to the case using x 2, f ,start = 0.05 due to an almost stepwise change of the molar fraction at the top of the packed bed (s = 1) in the beginning of the drying process. The molar fraction x 2, f ,1 changes here at the beginning of the drying process from nearly pure ethanol (component 2) to the autoclave entering pure CO 2 (component 1). The errors are; therefore, higher at node s = 1 compared to the other nodes.
In a third step, the overall errors calculated for all spatial nodes of the packed bed height (s = 1, . . . , S) and for all drying time steps (k = 1, . . . , K) are presented in Figure 6 for different discretizations/grids (S = 2; 4; 8; 16; 32) and both investigated initial conditions. Next to the above analyzed molar fractions, the calculated drying times as well as the dimensionless number 1 (both parameters were defined in the previously published paper [3]) are influenced by the total number of grid points in the region of packed bed (Figure 7). The calculated drying time is the time from start to end of the drying process, whereas the drying process ends when the ethanol molar fraction is less than x , ( = 321 ) = 0.0109 within all gel particles  It is obvious that the overall errors are higher for the initial condition x 2, f ,start = 0.95 compared to the calculations using x 2, f ,start = 0.05. For both initial conditions and both errors, l ∞ and l 2 , the FVM shows the best convergence behavior, meaning the error decreases the fastest for an increasing total number of grid points S (or smaller grid sizes ∆x).
Next to the above analyzed molar fractions, the calculated drying times as well as the dimensionless number K1 mean (both parameters were defined in the previously published paper [3]) are influenced by the total number of grid points S in the region of packed bed (Figure 7). The calculated drying time is the time from start to end of the drying process, whereas the drying process ends when the ethanol molar fraction is less than x end 2,g (T = 321 K) = 0.0109 within all gel particles of the packed bed [3]. x end 2,g depends slightly on the system temperature (see previously published paper [3]). The dimensionless number K1 mean is a relative measure that shows which mass transport step is the limiting one for the overall drying kinetic, being either diffusion in the gel spheres or mass transport in the bulk fluid. It can be used to find a fast-drying process at low CO 2 consumption [3].
From Figure 7a it is obvious that a minimum can be reached for the calculated drying time by increasing the number of nodes S. The black crosses in Figure 7 mark the settings of the FVM used in the previously published paper [3]. The calculated drying time from the previously published paper [3] is close to the found minimal value.
An increase of the total number of nodes S leads to an increasing K1 mean number (Figure 7b). The result of the previously published paper [3] is again close to the K1 mean using the FVM and S = 32. Nevertheless, it is unclear if a maximum can be reached by increasing the total number of grid points S further.

Mole balances
During the drying process no chemical conversion occurs and thus the balance in terms of mole fraction of ethanol and carbon dioxide were considered. The closure of the mole balance of ethanol is shown in Figure 8 in the form of the relative error (Equation (52)) As expected, the FVM shows the best closure of the mole balance among all discretizations due to its requirement by definition to close the mass balance for each volume element (Figure 8a,b). The

Mole Balances
During the drying process no chemical conversion occurs and thus the balance in terms of mole fraction of ethanol and carbon dioxide were considered. The closure of the mole balance of ethanol is shown in Figure 8 in the form of the relative error (Equation (52)) As expected, the FVM shows the best closure of the mole balance among all discretizations due to its requirement by definition to close the mass balance for each volume element (Figure 8a,b). The relative error of the FDM version A shows a better closure of the mass balance compared to FDM version B (Figure 8a). A maximum (corresponding to an absolute error minimum) was observed for the FDM version A (Figure 8b). The black crosses in Figure 8  Comparing the closure of the mole balance of the FVM (for fixed Δ and total number of nodes ) depending on the dimensionless number 1 (Figure 9) shows that the relative error is increasing almost linearly with increasing dimensionless number 1 . An increasing dimensionless number 1 stands for an increased convective term in the packed bed [3]. For a constant total number of nodes , an increase of the time step Δ leads to an almost linear increasing relative error (not shown in the diagrams). Comparing the closure of the mole balance of the FVM (for fixed ∆t and total number of nodes S) depending on the dimensionless number K1 mean (Figure 9) shows that the relative error is increasing almost linearly with increasing dimensionless number K1 mean . An increasing dimensionless number K1 mean stands for an increased convective term in the packed bed [3].
For a constant total number of nodes S, an increase of the time step ∆t leads to an almost linear increasing relative error (not shown in the diagrams). Comparing the closure of the mole balance of the FVM (for fixed Δ and total number of nodes ) depending on the dimensionless number 1 (Figure 9) shows that the relative error is increasing almost linearly with increasing dimensionless number 1 . An increasing dimensionless number 1 stands for an increased convective term in the packed bed [3]. For a constant total number of nodes , an increase of the time step Δ leads to an almost linear increasing relative error (not shown in the diagrams).  (49)). The condition number of the FVM is the lowest, whereas the condition numbers of the FDMs become closer to each other for increasing number of nodes. The condition number of the FDM version B is up to more than four times higher than the other condition numbers on coarse grids. Additionally, it fluctuates in the beginning of the drying process. The condition number of the FVM used in Selmer et al. 2019 [3] with S = 20 nodes is between the condition numbers at S = 16 and S = 32 of the FVM. The condition number of FDM version B is more sensitive to the initial condition x 2, f ,start than the condition numbers of the other two methods.
Next to the above presented drying time dependencies of the condition numbers, the corresponding maxima of the condition numbers are shown as a function of the reciprocal of number of nodes S in Figure 11. The previously discussed results can be easily rediscovered here: The FVM has the lowest maximal condition numbers, the maximum condition numbers of all investigated numerical methods increase with increasing number of nodes or decreasing grid sizes and the condition number of the FDM version B is most sensitive with respect to the initial condition. Next to the above presented drying time dependencies of the condition numbers, the corresponding maxima of the condition numbers are shown as a function of the reciprocal of number of nodes in Figure 11. The previously discussed results can be easily rediscovered here: The FVM has the lowest maximal condition numbers, the maximum condition numbers of all investigated numerical methods increase with increasing number of nodes or decreasing grid sizes and the condition number of the FDM version B is most sensitive with respect to the initial condition.

Computation time
The computation time as a function of the reciprocal number of nodes is shown in Figure 12. In principal, the computation time increases for decreasing mesh sizes (or increasing number of nodes ) for all three methods. The calculations using the FVM take less time than using the FDMs (except on the finest mesh). Additionally, the calculation time increases for the FDM version B for coarse meshes due to an overestimation of the calculated drying time (Figure 7a) and thus longer calculation durations using a fixed time step. Hence, a minimum in the calculation duration of FDM version B can be observed in Figure 12.
The black cross marks the settings of the previously published paper [3]. Its value is smaller than for the calculations of this paper due to a time step which is twice the time step of the calculations made here. Even though the time step is doubled for the calculations in Selmer et al. 2019 [3], the calculation quality is acceptable as shown in the previous sections.

Computation Time
The computation time as a function of the reciprocal number of nodes S is shown in Figure 12. In principal, the computation time increases for decreasing mesh sizes (or increasing number of nodes S) for all three methods. The calculations using the FVM take less time than using the FDMs (except on the finest mesh). Additionally, the calculation time increases for the FDM version B for coarse meshes due to an overestimation of the calculated drying time (Figure 7a) and thus longer calculation durations using a fixed time step. Hence, a minimum in the calculation duration of FDM version B can be observed in Figure 12.

Conclusions
In this work, we evaluated three different numerical methods, based on finite differences and finite volumes methods, to solve the advection diffusion equation arising in the context of the supercritical drying process of spherical gel particles in a packed bed.  The black cross marks the settings of the previously published paper [3]. Its value is smaller than for the calculations of this paper due to a time step which is twice the time step of the calculations made here. Even though the time step is doubled for the calculations in Selmer et al. 2019 [3], the calculation quality is acceptable as shown in the previous sections.

Conclusions
In this work, we evaluated three different numerical methods, based on finite differences and finite volumes methods, to solve the advection diffusion equation arising in the context of the supercritical drying process of spherical gel particles in a packed bed.
The FVM showed the best closure of the mole balance, best convergence behavior and lowest condition numbers. Therefore, this method is recommended to solve the coupled partial differential equations describing the supercritical drying kinetic in a packed bed developed in Selmer et al. 2018 [2]. A sufficient number of nodes is needed especially for high temporal changes of the molar fractions in the bulk fluid. The reason is that high spatial and/or temporal changes of the molar fractions result into high variations of the physical mixture properties and thus changes all particular mass transfer steps directly or indirectly. Especially the effective diffusion coefficient within the porous gel particle, the mass transfer coefficient and the interstitial fluid velocity within the packed bed are influenced. More detailed information is given in [2,3].
The FDM version A showed good accuracy and computational efficiency. It is easy to implement and should be used with a sufficient number of nodes.
The FDM version B showed the lowest accuracy and efficiency. Thus, it is not recommended to be used.

Acknowledgments:
We would like to thank Anna-Sophia Behnecke for developing the code for the computation of the diffusion within the porous gel spheres.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.  [−] Function V f m 3 Volume bulk fluid V g,n m 3 Volume element of spherical gel particle

Nomenclature
x 2, f mol EtOH mol EtOH +mol CO 2 Ethanol molar fraction in the bulk fluid x 2, f mol EtOH mol EtOH +mol CO 2 Vector ethanol molar fraction of the axial bulk fluid elements s and time indices k