A simple model for desulphurisation of flue gas using reactive filters

Desulphurisation of flue gas is essential before it can be released safely into the atmosphere. One way of removing sulphur dioxide is to use a purification device incorporating a reactive filter, in which the flue gas stream passes in front of a porous-catalyst-filled structure which converts the gaseous sulphur dioxide into liquid sulphuric acid. In this paper, we build and solve a simple mathematical model to describe the operation of a paradigm reactive filter. Our model captures the transport of sulphur dioxide through the device via advection in the main “outer” flow and diffusion through the catalyst structure along with the production of sulphuric acid. This sulphuric acid gradually accumulates in the filter rendering it less efficient. We determine the clogging time for an individual channel (that is, the time at which the entrance to the channel becomes completely filled with liquid) and explore how the concentrations of sulphur dioxide and oxygen and the thickness of the sulphuric acid layer change as the key dimensionless parameters are varied, comparing numerical and asymptotic results where appropriate. We then turn our attention to the device scale and solve our model numerically to determine the overall lifetime of the device. We vary the key dimensionless parameters and explore how they affect the efficiency of the device. In the physically relevant parameter regime, we find an explicit solution to the outer flow problem which agrees well with numerical solutions and provides a formula for the lifetime of the device. Finally, we propose a formula for determining the catalyst reaction rate, given data on the concentration of sulphur dioxide exiting the device.


Introduction
In the drive to protect the environment, reducing the concentrations of harmful chemicals that are released into the atmosphere has become a priority for industry. One key example is the removal of sulphur dioxide (SO 2 ) from flue gas. SO 2 is formed in vast quantities in industrial processes and power plants as a by-product of processing raw materials, such as crude oil and various ores [1][2][3]. Furthermore, it is a highly toxic gas that can cause acid rain and is linked to respiratory illnesses [3,4]. In order to decrease the concentration of sulphur dioxide in flue gas to be released into the atmosphere, various techniques are used, including "gas scrubbing" (both wet and dry) [5][6][7] [8], packed-bed absorption [9], and reactive filtering using catalysts [10]. However, the majority of these methods require high input power and a specifically situated operation site. Furthermore, they produce a lot of waste, such as gypsum containing a variable amount of impurities in the case of gas scrubbing, that needs to be disposed [2].
Our interest is in the use of reactive filters, since these are less costly than the other options in terms of maintenance, are easy to install in a factory, and do not generate toxic waste. These filters also produce sulphuric acid which can be easily captured, stored, and used for other purposes; the sulphuric acid also acts as a natural cleanser for the device. The filters are often made from stackable modules, each of which consists of a series of porous sheets containing catalyst material. Flue gas flows in front of these sheets from one end of the device to the other and diffuses into the porous sheets. When sulphur dioxide, oxygen (O 2 ), and water vapour come into contact with the catalyst-covered surfaces of the pores, liquid sulphuric acid (H 2 SO 4 ) is produced, as indicated in Fig. 1. In reality, there are multiple intermediate reactions that take place on the catalyst surface before acid is produced (see, for example, [11,12]), but the overall reaction is given by where r is the effective reaction rate taking into account all intermediate steps.
As the chemical reaction proceeds, liquid sulphuric acid accumulates in the pores and inhibits the transport of sulphur dioxide, which results in a reduction in device efficiency over time. The evolution of the sulphuric acid layers near the surface of the sheets results in the formation of surface droplets. When these droplets become sufficiently large, they run down the surface under gravity and accumulate at the bottom of the device. Our aim is to understand the transport of sulphur dioxide and oxygen in the device in order to gain insight into its operation and to determine how the control parameters alter device efficiency and longevity.
Although various models exist for other purification methods, including gas scrubbing [5,6,13,14] and absorption [15,16], reactive filters are less well studied. In [12], the removal of sulphur dioxide using activated carbon fibres, and the production of sulphuric acid, is studied. They develop a power-law model to describe the steady-state concentration of sulphur dioxide at the outlet of the device and find that this concentration is proportional to the weight of the catalyst and to a specific power of the initial concentrations of sulphur dioxide, oxygen, and water vapour. Furthermore, a detailed kinetic model for the evolution of the concentration of sulphur dioxide in an activated carbon fibre reactor is presented in [17]. They observe that increasing the inlet concentration of sulphur dioxide increases the outlet concentration, whereas increasing the concentrations of either oxygen or water vapour decreases it. Their model is based on a particular functional form for the concentration of sulphur dioxide in the fibre pores and, in addition, the governing equations make use of an effective uptake of sulphur dioxide and are derived from a simple averaging over the pore domain. Neither model involves the evolution of the liquid sulphuric acid produced by the catalytic reaction. In [18], a more sophisticated model for the operation of a catalytic device is derived for the case where spherical catalyst pellets are embedded in a porous matrix. Homogenisation techniques are used to upscale a model for the physical and chemical processes occurring on the pellet scale to one that is relevant for describing device operation. The model is used to predict the device-scale behaviour up until the point when the porous matrix becomes saturated near the inlet. A summary of models describing the related area of break-through-time predictions for gas removal using catalysts can be found in [19].
In this paper, we will build a simple mathematical model for the operation of a reactive filter device. We will assume the simplest possible geometry for the filter, namely a set of finite-length, thin, parallel, channels where the walls have been covered with catalyst material. In Sect. 2, we will derive a mathematical model for the device, which will incorporate an "outer" gas flow coupled with diffusion of sulphur dioxide and oxygen into the reactive channels and the consequent build up of sulphuric acid. In Sect. 3, we will nondimensionalise the model and determine the key dimensionless parameters that control the process. In Sect. 4, we will look at the physically relevant thin-catalyst-channel limit of the problem, which will enable us to reduce the complexity of the model. In Sect. 5, we focus our attention on the catalyst channel at the top of the device, and we will solve (numerically) for the sulphur dioxide concentration and thickness of the sulphuric acid layer in the channels. We will explore how the system behaves at the entrance to the channel and use asymptotic reductions to provide explicit expressions for some of the key variables. We will also consider a reduced model in the case when the liquid layer is very thin. In Sect. 6, we will solve the whole-device model which incorporates the "outer" gas transport as well as the activity in the channels. Finally, in Sect. 7, we will draw our conclusions and discuss extensions.

Mathematical model for the purification device
We focus on the simple two-dimensional device shown in the schematic in Fig. 2, where the flue gas flows along the outer region at the left-hand side of the device and diffuses into the catalyst-lined channels on the right. We neglect the transport of water vapour (assuming it is abundant) and also the hygroscopic effect of sulphuric acid (see [20]), and focus instead on the transport of sulphur dioxide and oxygen.

Model for the outer region
We let S(z, t) and C(z, t) denote the concentrations (measured in mol m −3 ) of SO 2 and O 2 in the gas in the outer region, respectively, where z denotes distance from the bottom of the device and t is the time. We assume that the device has length H and that the gas is transported in the outer flow by convection with constant speed U . Appealing to conservation of mass for SO 2 and O 2 , and assuming a priori that this outer flow is advection dominated, we write where Q s (z, t) and Q c (z, t) represent the flux (per unit depth, measured in mol m −2 s −1 ) of SO 2 and O 2 entering the porous structure, respectively, and l is the width of the outer region. We suppose that SO 2 and O 2 are supplied with a constant concentration at the inlet located at z = 0, and so we impose the boundary conditions: (3) Finally, we assume that there is no gas in the outer region initially, so we write As mentioned earlier, we assume for simplicity that the filter comprises numerous small parallel catalyst channels aligned perpendicularly to the direction of the outer gas flow. We focus on a single channel with the aim of determining the fluxes Q s and Q c of SO 2 and O 2 drawn out of the gas stream in the outer region. The geometry under consideration is shown schematically in Fig. 3. We assume that the surfaces of the catalyst channel are located at y = 0 and y = 2W and that the centres of the channels are a constant distance L apart. We use the coordinate x to denote distance into the filter, with the surface open to the gas flow located at x = 0 and the back of the filter located at x = L. We denote the concentration of SO 2 and O 2 inside the flue gas in the channels by s(x, y, z, t) and c(x, y, z, t), respectively and, assuming that the gas transport in the channel is diffusion dominated, we write where D g s and D g c are the diffusivities of SO 2 and O 2 in the flue gas, respectively. We suppose that the concentrations of SO 2 and O 2 are prescribed (by the outer region model) at the front of the filter, and that there is no flux of either species through the back of the filter. Thus, we write At the surface of the channel, the reaction turns the gaseous sulphur dioxide, oxygen, and water vapour into liquid sulphuric acid. Assuming that the water vapour is present in excess and, for simplicity 1 , that the reaction follows the law of mass action, the velocity v * of the sulphuric acid produced at the surface is where the concentrations of SO 2 and O 2 in the liquid layer that forms above the surface are denoted by S(x, y, z, t) and C(x, y, z, t), respectively, λ (measured in m 3 /mol) is the molar volume of sulphuric acid, r (measured in m 7 /mol 2 s) is the reaction constant, and the factor of two reflects that two molecules of H 2 SO 4 are produced in the reaction. The concentrations of SO 2 and O 2 dissolved in the liquid H 2 SO 4 layer evolve according to where D l s and D l c are the diffusivities of SO 2 and O 2 in the liquid, respectively, and u and v are the liquid velocities along the channel and across the channel, respectively. Appealing again to the law of mass action, the fluxes of dissolved SO 2 and O 2 removed from the liquid layer at the surface of the channels are At the moving interface, h(x, z, t), between the gas and the liquid H 2 SO 4 , we impose continuity of flux of SO 2 and O 2 , namely where u = (u, v), n is the normal to the interface, ∂/∂n is shorthand for n · ∇, and V n is the velocity of the interface. We also apply Henry's law, which states and comes from assuming that the chemical potentials for each species are in equilibrium at the interface. Here, β s and β c are partition coefficients for sulphur dioxide and oxygen, respectively. We further assume that no H 2 SO 4 molecules pass through the gas-liquid interface, and so the interface satisfies the kinematic condition u · n = V n , which can be written as: We note that (11) simplifies considerably under this assumption and reduces to conservation of diffusive flux across the interface. The liquid flow in the sulphuric acid layer is influenced by the (non-uniform) production of H 2 SO 4 at the channel surface and the surface tension of the air-liquid interface. To fully capture the motion, we would need to solve a model such as that presented in [21]. In this paper, however, our focus is on species transport, and so we use the simplest possible model for the liquid motion, in which we assume that there is no flow in the x direction, i.e., we set u = 0. With this assumption, conservation of mass requires that v is independent of y, and thus, using (8) and (13), we find that Equation (14) holds true until h = W , which we anticipate first occurs at the entrance to the catalyst channel. Subsequently, the channel entrance becomes blocked with liquid and the SO 2 extraction pathway switches to be solely through the liquid layer. We will discuss this further in Sect. 2.3. We assume symmetry at the centre of the channel, and so we write Finally, we assume that there is no gas in the channels and no liquid layer present at the start, so that we write

Flux drawn out of the outer gas flow
We relate the fluxes, Q s and Q c , of SO 2 and O 2 drawn out of the outer gas flow to the activity inside the filter by writing where χ 0 = 2W/L is the void fraction at the surface of the filter. Once the entrances to the channels become blocked (i.e. h(0, z, t) = W ), the dominant SO 2 and O 2 transport route switches to being through the liquid. Given the differences in the diffusivities and the solubilities given later in Table 1, we anticipate that it is much more difficult for SO 2 and O 2 to diffuse through the liquid layer spanning the channel entrance than to diffuse into a channel that is not blocked with liquid. Thus, once h = W at x = 0 for some z, we stop solving for s, c, S, C, and h at that z location and instead simply set the fluxes Q s and Q c to zero. This is represented in (17) by

Nondimensionalisation
We are interested in the behaviour of the device on the timescale over which the liquid H 2 SO 4 layer has grown sufficiently that it is of the same order of magnitude as the thickness of the catalyst channels, found by balancing the terms in (14) when h ∼ W . 2 We scale distance along the device with the device length, distance along the channel with the channel length, and distance across the channel with the half the channel width, which we also use as the scaling for the thickness of the sulphuric acid layer. We scale the gaseous concentrations of sulphur dioxide and oxygen with their initial concentrations, and the concentrations in the sulphuric acid layer with the concentrations that are in chemical equilibrium with those in the gas. We thus use the following scalings: The dimensionless model reads, dropping the primes, δR with boundary conditions and initial conditions We have introduced nine fundamental dimensionless parameters: γ is the ratio of the width of the outer channel to the length of the catalyst channels, is the aspect ratio of the catalyst channels, Pe is the Péclet number associated with transport in the outer region, R is the dimensionless reaction constant, α is the ratio of the inlet concentrations of O 2 and SO 2 , δ is the fraction of the molar volume of sulphuric acid realised by the local concentration of SO 2 that is related to the input concentration, and the D i are diffusivity ratios. These are given by We present the sizes of the physical parameters for a typical operating situation in Table 1, along with the dimensions and operating parameters of a typical device. We note that we do not know the size of the reaction constant r . However, we infer r by noting that it takes a few days for droplets of H 2 SO 4 to appear on the front of the filter and use our scaling for t and the other parameters to calculate an estimate. The corresponding dimensionless parameters are shown in Table 2.
We see that the model depends on various combinations of the fundamental dimensionless groupings, namely The first set is taken from data (references shown). The second set is typical operating parameters. The value of r is inferred We immediately see that, in normal operation, we may drop the time derivatives everywhere except in the equation for the evolution of the sulphuric acid layer (25), and we do this henceforth. Further, the parameters associated with flux transfer across the gas-liquid interface in (31)-(32) are all small, which motivates us to simplify the problem further.

Reduction of the model in the catalyst channel
Exploiting the fact that the square of the aspect ratio of the catalyst channel, 2 1, we reduce the complexity of the model by looking for asymptotic solutions for the dependent variables of the form X ∼ X 0 + 2 X 1 + · · · . We work in the richest asymptotic limit of the remaining parameters, namely along with D 1 , χ 0 , are all O(1). Here, B 1 is the ratio of the diffusive flux of sulphur dioxide entering a catalyst channel to the convective flux of sulphur dioxide passing the channel entrance and indicates how easy it is for sulphur dioxide to enter the channel, B 2 is the ratio of the diffusive flux of sulphur dioxide through the liquid sulphuric acid to the diffusive flux of sulphur dioxide through the gas, multiplied by the inverse of the square of the aspect ratio and indicates how easy it is for sulphur dioxide to pass across the gas-liquid interface, B 3 is the ratio of the diffusive fluxes of sulphur dioxide to oxygen reaching the catalyst surface and indicates how oxygen is depleted at the catalyst surface, and B 4 is the ratio of the diffusive flux of oxygen through the liquid sulphuric acid to the diffusive flux of oxygen through the gas, multiplied by the inverse of the fourth power of the aspect ratio and indicates how easy it is for oxygen to pass across the gas-liquid interface. We will present their sizes, based on the values in Table 2, in (95) in Sect. 6. Given that we are unsure about the precise size of r , we proceed assuming further that the dimensionless reaction rate R = O(1).

Solution in the liquid layer
In this limit, diffusion dominates across the liquid layer, and the relevant leading-order (in 2 ) versions of (23)- (25) and boundary conditions (30)-(33) are Writing P = B 3 S − C, we solve the resulting (trivial) system for P to find that P is independent of y, and thus where s| h := s(x, y = h, z, t) (and equivalently for c| h ). We solve (39) and (42) to find that where b = S| y=0 and substitute (43) in to (41) to find that b satisfies Thus, h satisfies We use the solution for S and C to replace the flux conservation boundary conditions (31) and (32) on y = h with which we can then use to determine the gas flow in the catalyst channels without having to further consider the transport problem in the liquid layer. We note that, in this limit, the right-hand side of (47) is O( 4 ), and thus we do not expect much O 2 to be removed from the gas by this process.

Solution in the gas in the catalyst channel
In the region of the catalyst channel filled by gas, diffusion across the channel again dominates. Remembering the limit we are working in (given by (38)), we can easily see that , t), and thus the gas concentrations do not vary across the channel. In order to close the model, we need to go to O( 2 ) in the equations and boundary conditions. The relevant parts are with while (46) and (47) become and we remember that the lateral conditions are We first consider the O 2 problem. Integrating (48) across the layer, applying (49) and (51), integrating along the layer and applying (52) and (53), we find that i.e. oxygen is present in such a high concentration that the reaction does not reduce the leading-order amount of O 2 present in the catalyst channels. Turning our attention to the SO 2 problem, performing similar integration across the layer we find that, dropping the zeros, s satisfies which, in summary, must be solved alongside The behaviour of the channel model is governed by the three parameters B 2 , B 3 , and R and is influenced by the concentrations of SO 2 and O 2 at the inlet to the channel. Of these parameters, once the device is manufactured, B 2 is fixed by the chemistry and the geometry, while B 3 and R alter depending on the amount of SO 2 supplied.

Solution to the catalyst-channel model at the entrance to the device
We focus on the channel at the entrance to the device, where S = 1 and C = 1 (we will see in Sect. 6 that, in fact, C = 1 everywhere), since here the behaviour of the catalyst channel decouples from the behaviour of the outer flow.
Our aim is to explore the model behaviour and to determine the clogging time for the catalyst-channel entrance.
We solve (55)-(59) numerically and, in Fig. 4, we show how s and h vary with x, at various times, assuming all the parameters are set to 1. We see that SO 2 adsorption along the surface of the channel causes the SO 2 concentration to decrease into the filter, as expected, but that the behaviour at a particular x value is non-monotonic as t increases, with the concentration increasing everywhere in (0, 1] until t = 0.345, after which it starts to decrease. We also see that the thickness of the sulphuric acid layer is largest at the inlet and decreases monotonically into the filter. However, in contrast to the SO 2 concentration, the thickness of the layer increases monotonically with time until the channel clogs at the inlet at t = t clog = 1.544. In order to understand why the concentration changes non-monotonically with time, we consider the case where R 1 and B 3 > 1. We scale t = RT and write the system of equations as: with Numerical solutions in this limit suggest that We focus on the situation where b ∼ O(1) everywhere in the domain, and thus s − b = 1/B 3 . We are then able to solve the leading-orderin-1/R system, which reads Thus, differentiating s given by (65) and using ∂h/∂ T given by (64), we find that which is positive for all we conclude that (in this limit) s initially rises, then falls, as seen in Fig. 4 (albeit for O(1) parameters) and that this is caused by the interplay between the thinning gaseous layer (the (1 − h) factor in (64)) and the weakening of the effect of the reaction as the liquid layer grows (the (1/ h) factor in (64)).

5.1
Behaviour at x = 0 The main quantities of importance from the channel model for the operation of the device are the flux q of sulphur dioxide removed from the outer flow by the catalyst channels, defined by and the evolution of the thickness of the liquid layer at the channel entrance, h(0, 0, t). In Fig. 5, we show how these quantities vary with time. We see that the thickness of the liquid layer increases monotonically, and that the flux decreases monotonically, as t increases, with a very rapid switch off near the clogging time.
We now explore in more detail the behaviour of the system at x = 0 in order to make predictions for the clogging time t clog . At x = 0, Eqs. (44) and (45) are where we use the notation b 0 (t) and h 0 (t) here to indicate values of b and h at x = 0. We are interested in solving (68) and (69) until the entrance to the channel clogs; the time at which this occurs is given by t clog such that h 0 t clog = 1. Before proceeding, we immediately note that B 2 cannot influence the clogging behaviour since it Equation (70) may be solved to get an implicit formula for b 0 , which then givesh 0 through (68). However, it is just as informative to solve (69) and (70) numerically. In either case, we use the initial conditionh 0 (0) = 0 along with b 0 (0) = 1, since this is the solution to (68) when t = 0. In Fig. 6, we show how the evolution of b 0 andh 0 depends on B 3 (we recall that B 3 depends on the ratio S * /C * , and thus altering B 3 corresponds to varying the ratio of the concentrations of sulphur dioxide to oxygen in the gas flow). We see in Fig. 6 (left) that b 0 decreases as τ increases and, for B 3 > 1, the solution approaches the steady state b 0 = 1 − 1/B 3 while, for B 3 ≤ 1, the solution tends to zero as τ → ∞. The approach to zero is relatively slow; for example, solving the B 3 = 1 case, we find that We also show the value of b 0 for whichh 0 = R for R = 1 (red dots) and R = 5 (blue dots). We see that increasing R moves the dots to the right in Fig. 6 (left), with consequence that, as R increases, b 0 becomes increasingly close to its steady-state value. In Fig. 6 (right), we also see that increasing B 3 slows the rate of increase ofh 0 , and that the clogging time, for a particular value of R, increases as B 3 increases.
In order to see the trends in the clogging time as we vary the parameters, we first recall that the time scaling in our nondimensionalisation includes the parameters r and S * . In order to assess the influence of varying r and S * , we remove this hidden dependence by considering t clog /RB 3 . We show how t clog /RB 3 varies with B 3 for R = 1, 10, and 100 in Fig. 7 (left). We see that increasing B 3 decreases the scaled clogging time, and that, for large B 3 , this behaviour is independent of R. To find this behaviour, we make the anzatz for  respectively, which we solve to find that Therefore, we conclude that which we see in Fig. 7 (left) is in excellent agreement with the numerical results. Conversely, when B 3 1, we find that the leading-order solution to (68) and (70), written in original variables, is Thus, at t = t clog , we find that Again, we see excellent agreement between the asymptotic predictions and the numerical results, for all three values of R shown in Fig. 7 (left). In Fig. 7 (right), we explore how t clog /RB 3 varies as we alter R, for various values of B 3 . We see that the scaled clogging time decreases as R increases, as expected, since an increase in the reaction rate generates more liquid which means that the inlet clogs more quickly. We show the small-R asymptotic solution (which we will establish in the appendix) as dashed lines and we see excellent agreement. Moreover, we note that the small-R limit holds over the whole domain for large values of B 3 , as shown in Fig. 7 (right) where we see that the asymptotic line lies on top of the B 3 = 100 line. We also show the large-B 3 asymptotic solutions as dotted line and also see excellent agreement there too. In the inset to Fig. 7 (right), 4 we show the behaviour of t clog /R versus R and include the curve for B 3 = 0. We see that this follows the B 3 = 1 curve reasonably closely, and that t clog /R → 0 as R → ∞ for both B 3 = 0 and 1. We note that, when B 3 1, (77) holds over the whole domain and provides an exact formula We note that the behaviour of t clog given by (77) and (78) are identical to the numerical results and are not shown separately in Fig. 7.

Varying the parameters
We now vary the parameters to determine their influence on the flux, q, of sulphur dioxide drawn into the catalyst channel. In Fig. 8, we show the effect of varying B 2 (this is analagous to altering the aspect ratio of the channels), B 3 (this is analagous to altering the ratio of the concentration of sulphur dioxide to oxygen), and R. We see that increasing all these parameters increases the flux of sulphur dioxide drawn into the catalyst channel, although we note that increasing R decreases the actual clogging time, as can be seen in the inset to Fig. 8 (right). In Fig. 9, we show how the total amount of SO 2 extracted by the channel before it clogs , given by − t clog 0 q dt, varies with RB 2 . We see that the total amount of SO 2 removed increases with RB 2 and that the behaviour scales like RB 2 when RB 2 1 and like √ RB 2 when RB 2 1. We note that it is possible to obtain an asymptotic solution to (55)-(59) in the case when R 1; we present this solution and compare it with a numerical solution for R small, in the Appendix.

The thin-liquid-layer (small-time) approximation
For small time, we expect the thickness of the layer of liquid sulphuric acid to be small compared with the width of the catalyst channels. When h 1, h = εH say, Eq. (57) reduces to s = b, and thus (55) and (56) reduce to where T = t/ε. We note that RB 2 , and thus, the solution, only depends on the solubilities, and not the diffusivities, of SO 2 and O 2 in the H 2 SO 4 layer. We integrate (79) to find that s satisfies where s e = s(1) and we have picked the negative gradient on physical grounds. In the thin-liquid-layer limit, the flux q ≈ −∂s/∂ x| x=0 . Unfortunately, the flux given by (81), evaluated at x = 0, depends on s e , and thus, we have to calculate s e . In principle, we could integrate (81) to give an implicit formula for s e , but it is just as instructive to solve (79) with boundary conditions (58) numerically. We show the resulting concentration profiles from a numerical simulation in Fig. 10 (left) for various values of RB 2 . We see that increasing RB 2 decreases the concentration of SO 2 in the channel and increases the flux of SO 2 drawn out from the gas. In order to quantify the affect of altering RB 2 , we rescale (79) and (58) using x = ξ/ √ RB 2 , and write the model as the first-order system Remembering that we expect p to be negative, we show in Fig. 11 (left) the phase plane for (82), for various values of p(0). The trajectories split into two sets. For p(0) > −2/ √ 3, the solutions reach zero slope at some finite value of s and at some finite value of ξ and represent the set of realisable solutions. For p(0) < −2/ √ 3, the solutions reach s = 0 at some finite slope; we reject these as unphysical since they do not reach p = 0. We show the dividing trajectory as a red-dashed line. This trajectory has s, p → 0 as ξ → ∞ and, thus, corresponds to p(0) = −2/ √ 3; this is in line with an appropriately scaled (81) with s e = 0. We solve (82)-(84) for a given value of p(0) and calculate ξ e , backing out the relationship between ξ e and − p(0), which we plot in Fig. 11 (right). We see that, as ξ e increases, the flux increases (agreeing with (121) for small ξ e ) and then tends to 2/ √ 3 for large ξ e , in keeping with the results presented in Fig. 10 (right) which we will describe on the next page.
Returning to Fig. 10 (left), we also see that the value of s e decreases as the value of RB 2 increases. We can, thus, find an approximation for the flux of SO 2 entering the channel by assuming that s 3 e 1. Neglecting s 3 e in (81), we arrive at the simple formula We show how q varies with RB 2 in Fig. 10 (right); our numerical results are presented as dots and the asymptotic results from (85) as a dotted line. We see excellent agreement between our simple formula and the numerically calculated flux, even for RB 2 = O(1). We also note that q ∼ √ RB 2 for large RB 2 , as previously seen in Fig. 9 albeit for the total SO 2 removed. The solution to (85) with boundary condition (58) (left) is This solution does not satisfy the boundary condition at x = 1. In fact, there will be a region near x = 1 where the terms involving s e that we have previously neglected are important. However, since the flux at x = 0 is accurately captured, we do not worry about this further. We next show how the (very thin) layer of liquid evolves. In Fig. 12 (Left), we show the spatial variation of the rate of increase of the thickness for various values of RB 2 . We see that increasing RB 2 localises the rate of film growth, and hence the liquid, near the entrance to the channel, as we might expect. Finally, we investigate how the flux into the catalyst channels behaves in the thin-layer case as we move along the device, i.e. as we alter S, remembering that we expect it to be directly proportional to S 3/2 in cases where we can neglect s e . We show the resulting behaviour as we alter S in Fig. 12, for RB 2 = 1 and RB 2 = 10. We see that the flux decreases as the inlet concentration decreases. We further see excellent agreement between the approximate solution given by (85) and the numerical solution when RB 2 = 10, and reasonable agreement even when RB 2 = 1. We note that (85) provides a very simple non-linear relationship between the flux of SO 2 removed from the gas stream and concentration of SO 2 in the gas stream that would have been difficult to predict had we not considered the catalyst-channel-scale model. We will use this relationship in Sect. 6.2, where we propose a simple method for determining r .

Device operation
We now consider how the device operates by coupling the transport of gas through the device with the removal of sulphur dioxide by the filter. As we have seen in Sect. 4.2, the leading-order concentration of oxygen in the catalyst channel does not vary in the channel, and thus, there is no leading-order flux of oxygen into the channel. Thus, we see from (21) that the device does not extract sufficient oxygen to make a significant change to its concentration as the gas passes through, and hence C = 1 everywhere.
The model for sulphur dioxide removal can be summarised as follows: and we note that, in practice, we solve (88) and (89) until h(0, z, t) = 1, after which we simply hold h = 1 at that z location and do not solve (88) and (89) further. The relevant boundary and initial conditions are Using the data in Table 2, we calculate that the physically relevant parameter values are and we (arbitrarily) pick χ 0 = 0.5 in order to estimate the size of the parameters under normal operation. We recall that R = 150 (but that it depends on the unknown r ). We note that we can control U , the speed of the gas flow, and possibly H , the length of the device (by adding additional filter modules in series) once the system has been built, and that S * and C * can vary depending on the gas fed into the filter. Thus, the main controllable parameters in the system are B 1 , B 3 , and R, while B 2 remains constant and known once the device has been built. We solve (87)-(94) numerically using an explicit Euler scheme for (87), (89), (91), and (94) along with Newton iterations for the non-linear algebraic problem (90) and boundary value problem (88), (92), and (93), discretised using a central finite-difference scheme. We implemented the scheme in MATLAB using 800 grid points for x and z (beyond which altering the step size did not appreciably alter the solution), with the parameters given in (95) and show in Fig. 13 how S and h| x=0 vary with z, for various values of t. We see that S and h| x=0 decrease along the device, with both S and h| x=0 initially confined to a region close to to the upper surface, and that both increase monotonically with time. We see that, once h| x=0 = 1 at z = 0, a front, located at z = G(t) say, propagates down the device behind which S = 1 and h| x=0 = 1. In Fig. 14 (left), we show how the concentration of SO 2 remaining in the gas stream at the outlet of the device varies with time, 5 while in Fig. 14 (right), we show how the front z = G(t) propagates. We see that the concentration of SO 2 exiting the device (for release into the atmosphere) remains low (less than 10% of the input amount) until about t ≈ 100, which is the equivalent of 296 days, and then increases linearly until just before the time at which the device completely stops working, which we call the device lifetime, t L . Thus, we conclude that the device becomes less efficient as time progresses. For this particular set of parameters, we find that t L = 264.3, which is equivalent to 2.1 years. We also see that G = 0 up until t clog , since in this period, the liquid layer thickness is less than one everywhere in the device (as can be seen in Fig. 13 (right)). After t = t clog , we see that G moves linearly through the device.
We now look at the effect of varying the key parameters on the amount of SO 2 exiting the device. In each case, we compare solutions in which we use the parameters in (95) with ones in which we halve, or double, the relevant parameter(s). In order to ensure that we compare like-with-like, we plot the dependence against t/RB 3 , to account for the hidden r and S * in our nondimensionalisation for time. In Fig. 15, we vary B 1 ; B 1 is directly proportional to the length of the device H and inversely proportional to the gas speed U and the width of the outer channel l. In Fig. 15 (left), we see that increasing B 1 increases the device lifetime, as expected, since this corresponds to increasing the length of the device, or decreasing the gas speed or the width of the outer channel (these both reduce the flux of SO 2 passing through the device). We also note that there is a slight kink in the curve for B 1 = 4; this occurs at t = t clog and is associated with the switching off of growth at the inlet and the transition to travelling wave behaviour. The kink is, in fact, present in all the curves, but decreases in magnitude as B 1 increases. In Fig. 15 (right), we see that the time at which G starts to move is the same for each curve, as expected, since t clog does not depend on B 1 . After t = t clog , we see that G(t) is linear in all three cases, and that increasing B 1 decreases the slopė G. We next present, in Fig. 16, the effect of varying the inlet SO 2 concentration S * (which results in scaling B 3 and R simultaneously; both are doubled or halved) on the concentration of SO 2 leaving the device and the front position. We note that we have scaled the vertical axis in Fig. 16 (left) artificially so that a concentration of 1 corresponds with the normal value of S * . We see that increasing the inlet concentration increases the outlet concentration, as found in [17], and consequently reduces t L . This is to be expected, since increasing S * increases the amount of SO 2 passing through the device (per unit time) and, thus, increases the amount of SO 2 entering the catalyst channels. In Fig. 16 (right), we see that the front starts propagating earlier as S * is increased, in keeping with the fact that t clog varies with both R and B 3 . We note that B 3 1 and R 1 for all three parameter sets in Fig. 17 (and, indeed, for all the results, we are presenting), so t clog can be approximated by taking the large R limit of (77), yielding We see that the predictions from (96) are in good agreement with the numerical behaviour seen in Fig. 16. After t clog , we see that the front propagates linearly through the device for all three choices of parameters, and that the slope increases as S * increases, in keeping with our intuition that a higher concentration of SO 2 at the inlet results in faster film growth throughout the device.
Next we look at how varying the dimensionless reaction rate alters the behaviour of the device. Since we do not know r , we look at five different R values to assess inaccuracy in r of two orders of magnitude. We show how the outlet concentration varies as we change R in Fig. 17 (left). We see that increasing R results in less SO 2 being released into the atmosphere at early time (until t/RB 3 ≈ 7, and then increasing R increases the amount of SO 2 released, resulting in a decrease in t L . This behaviour can be explained as follows: at early time, increasing R results in more SO 2 extracted by the catalyst channels close to the inlet and, thus, lower levels of SO 2 at the outlet. However, a higher rate of SO 2 removal causes faster growth of the H 2 SO 4 layers and, thus, faster clogging of individual catalyst channels and a corresponding increase in SO 2 reaching the outlet. In Fig. 17 (right), we see a small change in the clogging time as R increases, in keeping with (96). We also see that, after t clog , the fronts propagate linearly through the device, except for R = 15, where the speed increases slightly as the propagation progresses. We also see that the speed of propagation increases slightly as R increases. The decrease in t L is about 14% for an increase in the reaction rate of two orders of magnitude, suggesting that device operation is not so sensitive to the reaction rate. Further, given that the slopes of G are approximately constant (for R > 75), the change in t L is in line with the change in t clog , which (when B 3 = 0) alters by 20% as R changes from 15 to 1500.
Finally, in Fig. 18, we show the effect of varying the length L of the catalyst channels (via altering both B 1 and B 2 appropriately; B 1 √ B 2 remains constant). We see that increasing L increases the lifetime of the device, but the rate of increase decreases as L increases. We explore this behaviour further by considering the case where RB 2 1. In this case, we scale x = 1/ √ RB 2 ξ , and the modified parts of (87)-(94) are and while the problems for b and h remain unchanged. In the scaled coordinates, we see that increasing RB 2 increases the device length, and we anticipate that, as we saw in Fig. 11 (right) for the single-channel-no-liquid case, the scaled flux p will become constant even for moderate values of RB 2 . Since B 1 √ B 2 is independent of L, we expect that the device operation will be completely independent of L for RB 2 1. To illustrate this behaviour, in Fig. 19, we show how t L /RB 3 varies with B 2 (R = 150 throughout). We see that, above B 2 ≈ 2, the device lifetime does not appreciably change. Thus, we calculate that there is no point making the filter thicker than L = 2.2 mm, for the operating parameters presented in Table 1 From (95) and the value of R in Table 2, we see that R 1, B 1 1, B 2 1, B 3 1 during normal device operation, and so we now explore that limit. We first note that B 3 1 and so neglect all B 3 terms in (88)-(90). We then solve (90) to find that Thus, scaling t = RT , we write (88) and (89) as We now assume that R 1, so that we may neglect the terms of O( √ R) in (100) and (101), and look for an asymptotic solution of the form X = X (0) + B 2 X (1) +· · · , for small B 2 . The solution to the leading-order version of (100) with boundary conditions (92) and (93) is s (0) = S (0) , and thus, s (0) does not vary with x. The leading-order version of (101) is and we see that h (0) also does not depend on x. Proceeding to the next order in (100), integrating, applying (93), and evaluating at x = 0, we find that We use the flux given by (103) and the relationship (102) to write (87) as a single partial differential equation for h (0) , namely where Θ[.] is the Heaviside function. We immediately see that h = 1 is a solution to (104). When h < 1, we integrate (104) with respect to T to find where we have used the fact that h (0) = 0 at T = 0. The solution to (105) is To determine f (T ), we first consider T ≤ T clog := t clog /R. In this case, we solve (102) at z = 0, where S (0) = 1, to find that in line with (77) expanded for R 1. For T clog < T < T L := T L /R, we have h (0) = 1 for z ≤ G(T ). We apply the boundary conditions A (0) = 1 and S (0) = 1 at z = G(T ) to find that Thus, the leading-order solution is with S (0) found by solving (102) to give we note that the solutions for h (0) and S (0) are identical for t > t clog . We use (109) to predict that the (leading-order) device lifetime, T L , is given by i.e.
We see from (112), in this limit and to leading order, that neither the dimensional device lifetime, t dim,L , nor the dimensional clogging time, t dim,clog , depends on the reaction rate. This might seem counter intuitive, but we explain it as follows. While Eqs. (10) and (14) each depend on the reaction rate, if we combine them we find that Thus, the scaling for the dimensional clogging time 6 is given by setting h d = W , which results in t clog ∼ β s W 2 /λS * D l s , and does not depend on r . We recall from (96), however that the next term in the asymptotic expansion for t clog /RB 3 (and thus, its dimensional equivalent) does depend on the reaction rate. We conclude that, with this particular choice of parameter sizes, the device is diffusion limited rather than reaction limited at leading order and weakly varies with the reaction rate. The second term in the expression for t dim,L given in (112) is the ratio of the moles of sulphuric acid needed to fill the whole void space (= 2W L(H/L)/λ) to the flux of sulphur dioxide supplied (=S * lU ) and, thus, represents the time it takes to supply the device with the amount of suphur dioxide required to produce the amount of sulphuric acid needed to completely fill the void space. We do not proceed further with our asymptotic expansion, since terms of are all likely to come in at the same time. Instead, in Fig. 20, we present a comparison between our asymptotic solution and a numerical solution where we take B 2 = 0.01 and R = 10,000. We see excellent agreement between the numerical results and the asymptotic predictions, except where S (0) is small.

Determining the unknown parameter r
Our model can be used to estimate the value for the (dimensional) reaction rate r by using the solution to the thin-liquid-layer problem given in (85) and applying it in the device-scale model at early time. As we have seen in Sect. 5.3, when the liquid layer is thin, the flux of sulphur dioxide removed by the catalyst channels is given by (85), evaluated at x = 0. Thus, we can replace (87) with which must be solved along with boundary condition (91). The solution is analagous to (86), and the redimensionalised version reads We can, thus, back out that which provides a simple formula for working out r given the physical parameters of the model and an early-time experiment to determine S(H ).
efficiency via the concentration of SO 2 released into the atmosphere), and then the efficiency decreases until the device completely stops working, at about 2.1 years. Varying the parameters, we saw that increasing the length of the device or reducing the flux of SO 2 entering the device resulted in increasing the device lifetime, while increasing the reaction rate decreased the device lifetime. We do not know the value of the reaction rate precisely, and we showed that a change in reaction rate of two orders of magnitude only altered the device lifetime by a factor of 14%. We found that the device was particularly sensitive to changes in the concentration of SO 2 that were passed through it; increasing the concentration significantly reduced the device lifetime. Further, we explored the behaviour of the device as we altered the depth of the catalyst channels. We found that the device lifetime plateaued above L = 2.2 mm (for the operating parameters presented in Table 1), indicating that there is not much point building a device with filter thicker than this value. We then explored the physically relevant parameter regime, solved the appropriate asymptotic limit of the model explicitly, and found excellent agreement between the asymptotic solutions and numerical simulations, although we had to make B 2 and R smaller and bigger, respectively, than their physically relevant sizes to show convergence. We were also able to derive an explicit formula for the device lifetime. Finally, we used the small-liquid-layer approximation in order to provide a formula for determining the catalyst reaction rate, given some data on the concentration of sulphur dioxide exiting the device.
In summary, we conclude that the device lifetime can be increased by making the device longer, decreasing the speed of the flow, and increasing the length of the catalyst channels (up to a certain length). Moreover, the device should be tailored to specific gas concentrations, to avoid the large changes seen as S * is varied. We also conclude that the device is not particularly sensitive to the uncertainty in the reaction rate.
There are numerous limitations with our model that need refining. For example, we have assumed that the catalyst channels stop taking in sulphur dioxide once the sulphuric acid layer spans the catalyst channel. However, in practice, sulphur dioxide will diffuse through this layer to replenish the trapped gas behind this blocked region, and the sulphuric acid layer will continue to grow, albeit even more slowly than previously. Extending our model to include the fully saturated region is our obvious next step. Another key step will be to compare our model to experimental data in order to fit the value of the reaction rate and to give predictive power. We have used a simplified one-step reaction in our model but, in reality, there is a multi-step process. Extending our model to include multiple reactions would be useful. Finally, sulphuric acid is hygroscopic and, thus, draws in water vapour which results in accelerated growth of the liquid layer. Since water vapour is also present in the gaseous phase, extending our model to incorporate transport and absorption of water vapour, using a model such as in [20], would also enhance our model's predictive power.
Our overarching aim is to improve the understanding of reactive filters, and this simple model provides a stepping stone in that direction.
Since s (0) = 1, the leading-order flux of sulphur dioxide into the channel is zero, and we find that the O(R) flux is constant. We proceed to O(R 2 ) and find that the flux satisfies and thus, which we see scales like RB 2 , as seen previously in Fig. 9 for the total SO 2 removed in the small RB 2 limit. We use (119) to find an expression for the clogging time, t clog , namely In Fig. 21, we compare our asymptotic solution to the numerical solutions for R = 0.001, B 2 = 1, and B 3 = 1. We see excellent agreement between the two solutions, even near the clogging time.