A model for the lifetime of a reactive filter

In many industrial applications, reactive filters, which remove toxic substances using chemical reactions, prove very efficient and versatile. One such example is the removal of sulphur dioxide from flue gas by converting it into liquid sulphuric acid. In this paper, we extend a model developed in Kiradjiev et al. (SIAM J Appl Math 82:591–619, 2021) for the operation of a reactive filter to consider the long-time behaviour and determine how the filter becomes clogged with liquid and becomes inefficient. We explore how the parameters that govern the operation of the filter affect its efficiency and discuss optimal filtration strategies. We also consider a distinguished asymptotic limit in terms of the key dimensionless parameters and obtain an explicit solution that compares well with the numerical solution to the full problem. Finally, we compare the results of our model with those from a simplified model for a reactive filter considered in Breward and Kiradjiev (J Eng Math, 129(14):1–28, 2021).

Three modules of the filter device and a schematic of the filtration process (from [8]) that are currently in place to process this sulphur dioxide and prevent its release into the atmosphere, such as gas scrubbing [5][6][7], are costly, time consuming, and produce a large amount of waste.
W. L. Gore and Associates, Inc. have designed the Gore Mercury Control System (GMCS) filter, which converts gaseous sulphur dioxide into liquid sulphuric acid, which may be drained off and used for other purposes. The filter consists of microscopic catalytic carbon pellets that are held together by a polymer-composite network of fibres to form flexible sheets [8]. These sheets are then folded in a concertina-like fashion; the gaseous sulphur dioxide is passed through the triangular gaps and diffuses into the sheets where it reacts on the surface of the catalytic pellets, with oxygen and water vapour from the surroundings, to form liquid sulphuric acid. These concertina sheets are embedded in stackable modules and a schematic of the filtering device and the filtration process is shown in Fig. 1.
The conversion of sulphur dioxide to sulphuric acid takes place via a series of intermediate reactions (see, for example, [9,10]). However, these can be summarised by a single effective chemical reaction: Here, k is the overall reaction rate, taking into account all intermediate steps, and ( ) denotes gaseous state (g) initially, before a liquid layer has formed around the catalytic pellet, and aqueous solution (aq) afterwards. The sulphuric acid that is produced in the GMCS device acts as a natural cleanser of the filter by removing contaminant particles such as dust, and can be used for other purposes once drained off [8]. However, one drawback of the GMCS is that, if the liquid sulphuric acid does not drain away fast enough, it builds up in the void space between the catalytic pellets and the fibre network. If left for too long, the sulphuric acid will completely clog the filter, resulting in a dramatic drop in the processing rate of sulphur dioxide. To understand how to mitigate these clogging effects, we aim to build a mathematical model that describes the device behaviour. The model should describe the diffusive transport of sulphur dioxide within the void space of the filter and the subsequent accumulation of sulphuric acid on the pellets as a result of the catalytic reaction. This model may then be used to predict how the removal efficiency of the filter evolves with time as a result of the sulphuric acid build-up, and ultimately provide guidance on the optimal operating strategy that maximises the device efficiency.
A range of mathematical models exist in the literature to describe conventional methods for sulphur dioxide removal, including gas scrubbing [5,6,11,12] and absorption [13,14]. However, reactive-pellet chemical filters that involve a phase change of the contaminant, such as the GMCS device, are less well studied. Mochida et al. [10] study the removal of sulphur dioxide using activated carbon fibres to produce sulphuric acid. They develop a power-law model to describe the steady-state concentration of sulphur dioxide at the outlet of the device. They find that this concentration is proportional to the weight of the catalyst and to a specific power of the starting concentrations of sulphur dioxide, oxygen, and water vapour. In addition, they assert that the rate-limiting step of the reaction is the dissociation of the sulphuric acid that is produced into an aqueous solution around the carbon. Similar findings appear in Gaur et al. [15], where they develop a detailed kinetic model and describe the evolution of the gas concentration. They observe that increasing the inlet sulphur dioxide concentration or decreasing the oxygen and water-vapour concentrations increases the outlet concentration. However, in neither of these models is the evolution of the liquid sulphuric acid explicitly modelled.
A simple model for the behaviour of a reactive filter, in which the catalyst lines a series of parallel channels, is presented in Breward and Kiradjiev [16]. Here, we derive the equations governing the evolution of sulphur dioxide and oxygen in the gas and in the layer of sulphuric acid that forms on the surfaces, as well as the thickness of the acid layer in the channels, and couple this with the flow of gas through the filter. Under typical operating conditions, we find that the oxygen concentration within the device is constant. We postulate that, once the entrance to a channel has been filled with liquid sulphuric acid, no further transport of sulphur dioxide can take place in that channel. In doing so, we find that a travelling front propagates through the device, behind which all the channel entrances are saturated. We use the model to make predictions on the lifetime of the device, defined as the time it takes for the front to reach the exit of the device, and how the lifetime depends on the operating parameters. A more complicated model for a reactive filter, in which the catalytic pellets are assumed spherical, is presented in Kiradjiev et al. [17]. In this paper, we develop a homogenised model for a reactive filter in which we capture the microscale chemistry in a pair of coupled macroscale equations that describe the evolution of the sulphur dioxide concentration and the gas volume fraction in the filter. The model takes into account the complex microstructure of the porous medium through systematically derived effective quantities, such as the effective diffusivity, which are obtained by solving a corresponding cell problem. Oxygen and water vapour are assumed to have constant concentration due to their abundance under normal operating conditions. The key deficiency of this model is that it is valid only until the liquid layers in two neighbouring cells of the domain microstructure first coalesce; after this point, the topology within the unit cell changes. Our aim in this paper is to extend the model derived in [17] beyond the point of first coalescence by using the idea, introduced in [16], that once a continuous layer of liquid forms somewhere along the filter channel, there is no more uptake of sulphur dioxide at that location, and thus a clogging front travels through the filter channel until the filter stops working. This will enable us to examine the behaviour of the filter until it completely clogs with liquid, and to compare the behaviour with that predicted by the simple model in [16].
The outline of this paper is as follows. In Sect. 2, we will present an extended version of the mathematical model from [17] and derive an effective diffusivity that holds until the point when there is no single continuous gaseous pathway in the unit cell around a catalytic pellet. In Sect. 3, we will present the dimensionless model and, based on the values of the key dimensionless parameters, arrive at a reduced model for the filter. In Sect. 4, we will present numerical solutions to the model and vary the parameters that govern the performance of the filter. We will present an asymptotic limit of the model in Sect. 6 and obtain a Volterra integral equation for the position of the clogging front. In Sect. 7, we will explicitly compare the results from our model with those in [16]. We will draw conclusions in Sect. 8.

Mathematical model
We begin by stating the homogenised device-scale model, derived in [17], describing the operation of a flue-gas filter. In this model, flue gas flows along a filter channel of length L and thickness 2d at a constant speed U and diffuses into a porous filter sheet of thickness H , width W , adjacent to the channel, where sulphur dioxide is catalytically converted into liquid sulphuric acid. A schematic diagram of the filter is shown in Fig. 2, where we also introduce our Cartesian coordinate system, within whichX measures distance into the filter sheet,Ŷ measures distance across the sheet, andẐ measures distance along the sheet.
We track the concentration of sulphur dioxide in the filter channel and in the filter sheet, denoted byŜ andŝ, respectively, and the gas fraction, or porosity, in the filter sheet, denoted by V, where, for clarity, we use hats to denote dimensional variables. As explained in [18], we assume a uniform flow in the filter channel and negligible diffusion in theẐ -direction, since the associated Reynolds and Péclet numbers are large. We assume further that there are no variations in theŶ -direction, and we work in the distinguished limit where the diffusive transport of sulphur dioxide on the macroscale balances the production of sulphur dioxide due to the chemical reaction. In this limit, the concentration of sulphur dioxide does not vary across the sulphuric acid layer. Our two-dimensional model reads where (2) is a cross-sectionally averaged equation for the sulphur dioxide concentration in the filter channels in which advection along the filter channel balances diffusion into the filter sheet, while (3) is derived by homogenising the transport equation for the sulphur dioxide concentration around the microstructure in the porous filter sheets, assuming that the surface reaction (1) follows the law of mass action, and that H 2 O and O 2 are present in abundance and so their concentration is constant. Equation (4) is the macroscale equation describing the change of porosity due to the production of liquid sulphuric acid on the surface of the pellets. Here,D(V) is the effective diffusivity, which incorporates the effect of porosity in (2) and (3) and we will present later, V m is the molar volume of sulphuric acid, K is the effective reaction rate, which depends on the microscale properties and is given by K = 4π k R 2 /β 2 s l 3 , where k is the rate of the chemical reaction, R is the radius of the pellets, β s is Henry's law constant for sulphur dioxide, and l is the inter-pellet distance, Θ is the Heaviside unit step function, and V crit is a critical gas fraction, below which there is no continuous path for the gas to flow through. We discuss the specific form of V crit later in this section. We will extend the model presented in [17] by following the approach in [16] and holding the gas fraction constant as soon as it reaches the critical value V crit at a givenẐ location on the surface of the filter sheet. At this point we stop solving (2) at this particularẐ location and setŜ = S 0 , the inlet concentration of sulphur dioxide. We note that we anticipate that this will introduce a propagating front located atẐ =Ĝ(t), behind which the sulphur dioxide and the gas fraction at the surface of the channel are constant.
The associated boundary and initial conditions (see [17] for details) arê where (5) represents continuity of sulphur dioxide concentration across the surface of the filter sheet, noting that s andŜ represent intrinsic concentrations, i.e. per unit volume of available space, (6) represents no flux of sulphur dioxide through the back of the filter sheet, the concentration at the inlet of the filter is given by (7), and, through (8), we assume that the sulphur dioxide concentration at the outlet is the same in the channel and in the sheet, which is due to the gas diffusing laterally as it leaves the channel and spreads over the end of the filter sheet. We also assume that, initially, there is no sulphur dioxide nor liquid sulphuric acid present in the device, as given by (9), and we introduce an initial gas fraction, V 0 , which is related to the microscale parameters by which follows from computing the volume of gas in a cubic cell of length l with a catalytic pellet in the centre, as shown in Fig. 2. We summarise the physical parameters and their typical values in Table 1. Using these values gives V 0 ≈ 0.98. We note that this model has been derived by formally homogenising a microscale model where the liquid sulphuric acid layer grows uniformly around spherical catalytic pellets that are assumed to be arranged in a periodic array. This assumption is an idealisation of the structure of the porous sheets, but we expect no qualitative difference in the results, and only minor quantitative discrepancy, since the periodic arrangement represents the average structure of the filter sheets. More specifically, our microscale model involves an advection-diffusion equation for the transport of sulphur dioxide in the gas and through the sulphuric acid layer coupled with an equation describing the growth of the layer due to the reaction on the surface of the pellets. The overall dynamics of the system are governed by four  Note that k is an effective rate, and its value is determined indirectly from experimental data, as described in [18] dimensionless parameters: the dimensionless reaction rate, the ratio of the sulphur dioxide diffusivities (in the liquid and gas), the solubility of sulphur dioxide in the liquid, and the ratio of the microscale length to macroscale length (there are two further dimensionless parameters: the product of the sulphur dioxide concentration and the molar volume of sulphuric acid, and the size of the pellets relative to the spacing between them). The upscaled model (2)-(4) is derived in the physically relevant parameter limit in which the reaction rate and the ratio of the diffusivity ratio to solubility are both small and of similar magnitude. In this limit, we find that the leading-order sulphur dioxide concentrations in the gas and in the liquid are the same and independent of the local microscale spatial variables. As part of the homogenisation, we derive two effective quantities, namely, the effective diffusivity and the effective sink term in (3), the latter of which is proportional to the square of the sulphur dioxide concentration. These two quantities capture the essence of the microscale processes and will be used in the macroscale model. We note that the form ofD presented in [17] and found by solving the corresponding cell problem is only valid until two neighbouring liquid layers coalesce, which first happens near the channel inlet. At this point, V = 1 − π/6, which corresponds to the volume of the space between a unit cube and an inscribed sphere of unit diameter (see Fig. 3a). After the first clogging occurs, there will still be continuous gaseous pathways from the channel into the sheet, unlike in a two-dimensional scenario, where coalescence of two neighbouring layers will block diffusion of sulphur dioxide through the gas into deeper areas of the filter. Thus, we calculate the effective diffusivity of sulphur dioxide beyond the point when V = 1 − π/6 by numerically solving a cell problem on the microscale geometry represented by a unit cube ω (with boundary ∂ω) around a spherical catalytic pellet, as shown in Fig. 2. The gas is located in the region ω g , and the gas-liquid interface is located at Γ i . Denoting the microscale concentration of sulphur dioxide in the gas (normalised by the inlet concentration S 0 ) by s g , we find that this is related to the normalised macroscale concentration s =ŝ/S 0 by , the cell function, satisfies the dimensionless cell problem given by (see details in [17]) where x = (x, y, z) are coordinates in the cell problem, related to (x,ŷ,ẑ) in Fig. 2 by (x,ŷ,ẑ) = l(x, y, z). This is subject to where e j are the unit vectors in the x, y, z directions, respectively, e r is the unit normal vector to the gas-liquid interface (see [17] for details), and Γ i is defined by where Σ r = {x 2 + y 2 + z 2 = r 2 |x, y, z ∈ R}. Here, unlike in [17], we assume that, after the liquid spherical layer touches the sides of the unit cubic cell, the evolution of the gas-liquid interface Γ i follows that of a capped sphere until the continuous gas phase is broken up into eight isolated gas pockets, as shown in Fig. 3b. We can readily calculate the gas fraction at this point to be With this assumption we neglect any local surface-tension effects that would minimise the total gas-liquid surface area around neighbouring pellets. The effective diffusivity tensor is given bŷ where D s is the molecular diffusivity of sulphur dioxide in the gas. We note that, due to isotropy, there is no preferred direction in calculating the effective diffusivity, and thus it is proportional to the identity matrix, i.e.D =DI. We solve (11)-(13) using COMSOL Multiphysics for each potential location of Γ i and then use the solution to calculatê D.
In Fig. 4, we show how the effective diffusivity, normalised by its maximum valueD = D s attained when V = 1, varies as a function of the gas fraction. This agrees with the results derived in [29] and [30], for example. We denote the point of first clogging when V = 1 − π/6 by a dashed line and observe that this coincides with a qualitative change in the behaviour of the diffusivity, marked by a slight change in the concavity of the profile. This can be explained by the change in topology of the domain, which limits diffusion of the gas around the sphere. We also note that the effective diffusivity becomes zero when V = V crit , as previously discussed.
We will now non-dimensionalise the model and identify the key dimensionless parameters.

Dimensionless model
We non-dimensionalise (2)-(9) usinĝ where we have chosen the timescale associated with liquid growth, by balancing the terms in (4), assuming that this is the relevant timescale for the long-term operation of the filter. The dimensionless equations become subject to where we have the following five dimensionless parameters: Here, β represents the ratio of the diffusive flux of sulphur dioxide entering the filter sheet to the advective flux of sulphur dioxide passing through the filter channel and indicates how easy it is for sulphur dioxide to enter the sheet, δ is the aspect ratio of the filter sheet, κ measures the relative strength of the reaction to diffusion across the thickness of the filter sheet, and τ is the ratio of the advective timescale to the timescale associated with liquid growth.
Before we proceed, based on the parameter values in (27), since δ 1 and τ 1, from now onwards we will neglect the first term in (18) and the second term in (19). This means that S will evolve quasi-statically and diffusion in the Z direction is negligible, 1 which is due to the slenderness of the device (δ 1). We note that taking the quasi-static limit means that we lose the start-up problem where the concentration rapidly rises from zero as the sulphur dioxide fills up the channel for the first time. For completeness, the reduced dimensionless model is subject to In Sect. 7, we will compare this model to a similar one derived in [16]. The key difference between these two models is the presence of a non-linear effective diffusivity in (28)-(29), which systematically captures the effect of the microstructure on the transport of sulphur dioxide and liquid sulphuric acid within the filter sheet.
For future convenience, we also define the clogging time, t clog , to be the time at which V first equals V crit (which happens at Z = 0) and the complete clogging time, T clog , to be the time at which V equals V crit at Z = 1 (or, equally, when G = 1).
Before we solve the problem numerically, we calculate the clogging time explicitly, noting that it depends only on the behaviour of the filter at the inlet X = Z = 0, where s = S = 1. Evaluating (30) at this point, we obtain which we solve with (35) to find Thus, t clog is found by equating (37) to V crit , and therefore Redimensionalising (38) and substituting for V crit , the dimensional clogging timet clog is given bŷ We note that the clogging time depends on the reaction rate, the properties of the filter sheet, such as the radius of the pellets and inter-pellet distance, and the inlet concentration of sulphur dioxide. In particular, for less contaminated flue gas, where S 0 = 5 × 10 −3 mol m −3 ,t clog ≈ 63 days. We also note that the clogging time does not depend on the gas speed nor the dimensions of the filter, as previously noted. 5 Graphs of the numerical solution to (28)-(35) for a the spatial profile of S along the channel for t = 0, 1, 2, 3, 4, 5 (complete clogging time), b the spatial profile of V at X = 0 along the channel for the same times, c the temporal profile of S at the outlet Z = 1, d the spatial profile of V at Z = 1 across the filter sheet for the same times, and e the temporal profile of the front G(t), where the dashed lines correspond to t clog and T clog , respectively. In all plots, β = 6.7, κ = 0.88, V 0 = 0.98

Numerical results
We solve (28)-(35) numerically in MATLAB using an explicit Euler scheme and Newton iterations with 800 grid points in both directions. 2 In Fig. 5a, we show a spatial profile of the sulphur dioxide concentration, S, along the channel for different times using the parameters from (27).
We see that, as time increases, a propagating front forms, behind which S = 1. In Fig. 5b, we show the corresponding plot for the gas fraction V at the surface of the channel. We see that the gas fraction decreases to its minimum value, V crit , first at the inlet of the filter, and then until the sheet is clogged along the whole channel. In Fig. 5c, we show a temporal profile of the sulphur dioxide concentration at the outlet of the device. Before t = t clog , the outlet concentration changes very little. After t = t clog , it starts to increase more rapidly until it reaches S = 1 at t = T clog . This is often referred to as the "breakthrough curve" for S and is commonly used to assess the performance of a filter (see, for example, [31]). In Fig. 5d, we show a spatial profile of the gas fraction at the outlet of the channel, V| Z =1 , as a function of depth into the filter, since this shows us the extent to which the available filter has been used. We see that, while the porosity decreases rapidly near the channel surface, it changes very little deeper into the sheets, suggesting that these portions of the filter are not fully used. In Fig. 5e, we show the temporal profile of the position of the front G(t) behind which S = 1 and V = V crit . As predicted in (38), t clog ≈ 0.95. For t > t clog , the filter progressively clogs along the channel, decreasing its efficiency until the filter is completely clogged at T clog = 5, approximately five times longer than the clogging time. For the operating parameters we have used, this corresponds toT clog ≈ 13 days, whereas, for less contaminated flue gas, where S 0 = 5 × 10 −3 mol m −3 , t clog ≈ 325 days. We note that the propagation speed of the front is almost linear but slightly convex (G (t) > 0). This is because, as the filter clogs, the part that is still operational is decreasing in size, which speeds up the clogging in a positive-feedback mechanism. We further elaborate on this, and provide mathematical insight, in Sect. 6, where we obtain a relevant explicit asymptotic solution.
Following the approach in [16], we now vary the parameters of the model by doubling or halving their values to determine their effect on the behaviour of the device. We recall from (17) that the timescale of the problem is given by which changes, in particular, with k, l, and S 0 . Thus, in order to explore the behaviour of the solution as these (dimensional) parameters vary, we choose to plot all the results on a fixed time axis where t 0 = 2.65 days, which comes from evaluating (40) using the parameters from Table 1 (thus, t =t in Fig. 5). This way, we are able to compare all solutions on a common time axis that does not change regardless of which parameters we vary. We note that, after rescaling time with (40), (30) becomes is the ratio of t 0 to the liquid growth timescale for the particular set of parameters under consideration; α = 1 in the "standard" parameter case. In all figures that follow, we plot the outlet concentration S(1,t) and the position of the front G(t) as functions of time to measure the performance of the device. In particular, our measures for efficiency will be the value of the outlet concentration and the complete clogging time of the filter. To aid ourselves in the comparisons, we will also assume that we are given a pre-set outlet concentration threshold, say, 10 −3 mol m −3 , that should not be exceeded due to external control requirements, which is often the case when filtering flue gas.
In Fig. 6, we vary the dimensionless parameter β. This can be thought of as changing the half-width, d, of the filter channels, the speed, U , of the gas, or the length, L, of the filter channels, since these dimensional parameters appear only in the definition of β. As we increase β, the filter efficiency improves, since the outlet concentration decreases, and the complete clogging time increases, extending the lifetime of the filter. This reflects the fact that decreasing the width of the filter channels or the speed of the gas reduces the flux of sulphur dioxide through the filter, while increasing the length of the filter channel increases the residence time, and hence the interaction between the sulphur dioxide and the filter material. However, we note that, even though decreasing the channel width seems beneficial, in practice, the shear on the channel walls increases and results in the need for a higher pressure drop in order to maintain the same speed. Similarly, decreasing the speed too much will slow down the industrial process, while increasing the length of the filter channel will require more material. Thus, normally, there will be pre-set minimum channel width and gas speed, and a maximum length of the filter channel. Thus, fixing two of these, we (b) (a) Fig. 6 Graphs of a the temporal profile of S at the outlet Z = 1, and b the temporal profile of the front G varying β = 3.35, 6.7, 13.4. Here, κ = 0.88, V 0 = 0.98 Fig. 7 Graph of the temporal profile of the front G varying L by changing β =3.35 (solid green), 6.7 (dashed red), 13.4 (dotted black). The two asterisks denote the point of clogging for each of the two first two curves, respectively. Here, κ = 0.88, V 0 = 0.98 appeal to our convention of a given pre-set outlet concentration threshold, and we use our results to calculate the maximum gas speed, the maximum channel width, or the minimum channel length (depending on which ones we fix) in order not to exceed this threshold in the initial stage of the device operation. For example, if we fix the channel width and length, we find that the gas speed should not exceed 0.2 m s −1 in order to keep the outlet concentration below 10% of the inlet concentration, i.e. 10 −3 mol m −3 , during the initial stage of operation before the clogging front is formed. We see that the clogging time does not change when we vary β (see (39)). This is because it depends only on how fast the reaction proceeds at Z = 0, where S = 1. We also note that, if we attribute variations in β to variations in L, we need to appropriately scale the solution for G in Fig. 6b to reflect the change in the length of the channel, whereas Fig. 6a will remain the same. We show the corresponding plots of the position of the clogging front when varying L in Fig. 7. The resulting curves overlap each other, which is expected, since varying the length of the filter channel does not affect the dynamics of the front but only the time it takes to reach the outlet of the channel.
Since κ depends linearly on k and S 0 , whereas α depends linearly on k and quadratically on S 0 , we next explore varying κ and α in two different ways in order to mimic altering k and S 0 . In Fig. 8, we present solutions in which we have doubled or halved κ and α from their standard values in order to mimic varying the reaction rate k. We see that decreasing κ increases the outlet concentration of sulphur dioxide for early time and decreases the outlet concentration for later times. We observe similar behaviour in [16], which arises due to the change in the rate of clogging with time: initially, decreasing the reaction rate results in less sulphur dioxide extracted at the outlet of the filter, i.e. an increase in the outlet concentration. However, the lower rate of removal reduces the rate of liquid production, and hence clogging, which eventually decreases the outlet concentration. We also note that decreasing κ always increases the complete clogging time, since the reaction rate is decreased. The clogging time also changes in this case (cf. (39)), since we are altering the reaction rate. We also see that, as we decrease the reaction rate, the convexity of G(t) increases, which agrees with the observation in [16].  1.76, 4). Here, β = 6.7, V 0 = 0.98 In Fig. 9, we present solutions in which we have doubled or halved κ while multiplying α by 0.25 or 4, in order to explore the effect of changing the inlet sulphur dioxide concentration, S 0 , on the performance of the device. We note that we have scaled the solutions for S in Fig. 9a so that we can more clearly see the effect of altering S 0 . Decreasing the sulphur dioxide concentration decreases the outlet concentration and increases the lifetime of the filter as expected. We note that halving the sulphur dioxide concentration quadruples the clogging time, in agreement with (39), whereas the complete clogging time is almost tripled. This suggests that the lifetime of the filter is sensitive to changes in the inlet concentration of sulphur dioxide, as also found in [16]. In addition, the speed of propagation of the front decreases, as the concentration of sulphur dioxide is decreased, as expected.
In Fig. 10, we simultaneously vary κ, V 0 , and α in order to explore changing the inter-pellet distance l and thus the effect of varying how densely the catalytic pellets are packed within the filter sheet. We see that increasing l increases the outlet concentration for early times but decreases the outlet concentration for later times, which is reminiscent of the behaviour we observed in Fig. 8 when we decrease κ. This is because increasing l decreases the number of catalytic pellets in the filter sheet, and thus reduces the available area for reaction. This has the effect of decreasing the overall reaction rate, which can also be achieved by decreasing κ. We also note that increasing l increases the clogging time and the complete clogging time, since the pellets now are further apart. In particular, doubling the inter-pellet distance increases the clogging time almost eight times (cf. (39)) and triples the complete clogging time.
In Fig. 11, we simultaneously vary β and κ to show how varying the sheet thickness H alters the behaviour of the filter. Decreasing H decreases the outlet sulphur dioxide concentration and increases the complete clogging time, but the changes become smaller as H is continuously decreased. This suggests that there is a saturation type of behaviour, where, beyond a certain point, decreasing H does not have an effect on the system behaviour, as seen in [16].
In order to explore this further, in Fig. 12, we show a plot of the complete clogging time as a function of κ, keeping κβ 2 constant so that this corresponds to varying H only. We see that beyond κ ≈ 14, i.e. H ≈ 4 mm, there is a negligible change (less than 0.1%) in the clogging time, and, therefore, there is no benefit from making the filter sheet thicker than this value. We note that we observe the same behaviour in [16], with a prediction there of H = 2.2 mm.

Optimal relationship between the gas speed and the inlet concentration of sulphur dioxide
Finally, having observed a large variation in the results when altering the inlet sulphur dioxide concentration in Fig. 9, it would be convenient to have a simple prediction for how the gas speed can be chosen in response to changing the inlet concentration, so that the complete clogging time is the same. In order to do this, we make use of the previous observation that filter sheets thicker than 4 mm exhibit the same behaviour. This means that, if we are given a sheet of this thickness, or larger, for mathematical convenience we may assume that the sheet is infinite in thickness without altering the resulting behaviour of the filter. We then scale X = βX to eliminate β from (28). In order to be able to compare results on the same timescale, however, we again choose the timescale given by t 0 in Sect. 4 and defined by (17). Similarly, since we want to vary the inlet concentration, we need to rescale s appropriately, as in Fig. 9. Denoting the new inlet concentration by S , we write (s, S) = (S 0 /S )(s,S). The governing equations become with corresponding boundary and initial conditions (31)-(35). We see that the solution depends on the single parameter Thus, assuming the filter sheet is thick enough, if we increase the inlet concentration by a factor of n, we need to decrease the gas speed by a factor of √ n so that B remains the same, and thus we maintain the same efficiency and lifetime of the device.

The limit β 1 and κ 1
To gain further insight into the behaviour of the system, we examine the limit of thin filter channels (or slow gas flow) and thin filter sheets (or small reaction rate), i.e. when β 1 and κ 1. Specifically, we consider the distinguished limit β → ∞ and κ → 0 with βκ = λ = O (1). In this case, after expanding S, s, and V in powers of κ, in the form we obtain an explicit solution for the leading-order contributions of S and V until t = t clog , following the approach in [16]. Since κ 1, the leading-order version of (29) (dropping the superscripts) is We integrate (48) and apply (32) to find We note that this also satisfies the leading-order version of (28), where β 1. The first-order versions of (28) and (29) are Multiplying (51) through by λ, integrating the result with respect to X from 0 to 1, and using (31), (32), and (50), we obtain a single ordinary differential equation for S, namely which we solve alongside (30), together with The solution reads We note that, using (56b), we recover the expression for t clog , given by (38) by setting V = V crit and Z = 0. For t ≥ t clog , there is a moving front at Z = G(t), behind which we have Ahead of the front, we solve (52) and (53), together with the conditions where (59) follows from (56b). We find that the solution for t ≥ t clog and Z ≥ G(t) is In order to find G(t), we evaluate the expression for V in (60) at Z = G(t), rearrange and use (57) and (38) to obtain the following non-linear Volterra integral equation of the second kind for G(t): The advantage of (61) is that we have an implicit expression defining the front, which we can manipulate to obtain further insight about the front profile, and in particular its convexity. Differentiating (61) with respect to t (using Leibniz' rule), we obtain Fig. 13 Graphs of the numerical solution (solid black) to (28)-(35) and the asymptotic solution (dashed red), given by (56a,b), (57), (60), and (61), for a the spatial profile of S along the channel for t = 0, 1.5, 2, 2.45 (complete clogging time), b the spatial profile of V at X = 0 along the channel for t = 0, 0.5, 1, 1.5, 2, 2.45, c the temporal profile of S at the outlet Z = 1, and d the temporal profile of the front G(t). In all plots, β = 10 3 , λ = 1, V 0 = 0.98. We note that there are fewer times depicted on Fig. 13a than on Fig. 13b, since the solution for S is static until t = t clog , as evident from (56a) and Fig. 13c where we have used the definition of t clog (38). Since G(t clog ) = 0, evaluating (62) at t = t clog , we find that Differentiating (62) with respect to t again and rearranging the resulting identity, we obtain Evaluating (64) at t = t clog , using (63), and rearranging, we obtain This shows that, at least initially, the profile of the front is convex, which supports our observations in Sect. 4. In Fig. 13, we show a comparison between the numerical solution to (28)-(35) for S, V, and G and the asymptotic solution, given by (56a,b), (57), (60), and (61).
We see that there is excellent agreement in all plots, but we note that we had to use a very large value for our asymptotic parameter β, namely, β = 10 3 . We find that the clogging and complete clogging times for this regime are t clog = 0.95 and T clog = 2.45, respectively. In Fig. 13d, we see that the speed of the front Z = G(t) increases slightly with time. The profile nonetheless is close to linear, which is reminiscent of the findings in [16], where linear profiles are observed for large reaction rates and slight deviations appear when the reaction rate becomes small.
We also use the solution to (61) to explore how the complete clogging time, obtained by finding when G = 1, depends on the remaining parameters of the system in the limit β → ∞ and κ → 0 with βκ = λ = O(1). In Fig. 14a, we show how the complete clogging time depends on λ on a log-log scale. This can be thought of as changing the reaction rate, k, which is the most uncertain parameter in the model. For an easier physical interpretation, we rescale dimensionless time t with λ (i.e. t = λT , or in terms of the original dimensional time, t = (H L/dU S 0 V m )T ) so that, when varying this parameter through k, the time axis will remain unscaled, since k does not appear in the definition of this new timescale. We see that, as λ → 0, T clog → ∞, as expected, and T clog monotonically decreases with λ. We also note that the decrease becomes very gradual after about λ = 2, which suggests that the reaction rate does not dramatically affect the complete clogging time beyond that point. In Fig. 14b, we show how the complete clogging time depends on the inter-pellet distance l by varying the initial gas fraction V 0 and λ simultaneously (keeping (1 − V 0 )/λ constant). We note that V 0 ≥ 1 − π/6 ≈ 0.47, which corresponds to the case when neighbouring catalytic pellets touch each other. Again, since we have rescaled time with λ, varying l does not have an effect on the timescale, since it does not appear in the definition of the new timescale. We see that, as the inter-pellet distance increases (i.e. the initial gas fraction increases), the complete clogging time increases, as expected, since the pellets are farther apart. We observe an approximately slow linear increase in the clogging time until about V 0 = 0.95, and a rapid increase afterwards.
In both Fig. 14a and b, we see approximately linear profiles (in the log-log scale) for λ 1 and 1 − V 0 1, respectively, which is indicative of power-law dependence of the complete clogging time on these parameters. In order to explore this further, we take the limit λ → 0 in (61) (remembering to rescale t with λ so that we do not alter the timescale when changing k and l, as previously mentioned), expand both sides to first order in λ, and differentiate the result to obtain a differential equation for G(t), namely, Applying G = 0 at t = t clog /λ (after rescaling time), we obtain Since, by definition, G = 1 at t = T clog , we find that In Fig. 14a, we show the asymptotic result (68a,b) with and without the exponential correction, and we see excellent agreement with the numerical solution.
Returning to the limit 1 − V 0 1, we recall that, when varying the inter-pellet distance l, we do so by keeping (1 − V 0 )/λ = 0.02 constant (using the parameters from Fig. 13). Thus, the limit of small 1 − V 0 corresponds to small λ, which we have already considered. Therefore, to obtain the corresponding asymptotic plots for small 1 − V 0 in Fig. 14b, we plot (68a,b) in terms of 1 − V 0 rather than λ, noting that there is an additional V 0 inside the definition of t clog (38). In particular, Again, we see excellent agreement between the asymptotic results and the numerical solution.

Comparison with the simple model from [16]
We now explicitly compare the results of our model with our simple model presented in [16]. This model applies to a rectangular geometry, where, instead of catalytic pellets, the filter sheet consists of parallel catalyst channels, as shown in Fig. 1. Thus, as the reaction proceeds, a liquid layer of sulphuric acid forms along the walls of each of the catalyst channels. When the liquid completely blocks an entrance to a particular channel, the filter is assumed to stop working (i.e. there is no more uptake of sulphur dioxide) at that location.
For ease of comparison, we present the equations from [16] after relabelling the variables and parameters in the equivalent form to (28)-(35). In particular, the dimensionless parameters β and κ are identical to their counterparts in [16], where we note that the grouping 4π R 2 /l 3 , which appears in the definition of κ, represents the ratio of catalyst surface area to available volume, and the corresponding quantity in the simplified model from [16] is equal to 1/W , where W is half of the height of the catalyst channels. To aid the comparison, we introduce the gas thickness, h g instead of the liquid thickness, h, and we have rescaled time to remove any numerical pre-factors in (72). Finally, in order to compare the two models, we need to take a particular distinguished limit of the model from [16] that corresponds to our model (28)-(35). This is a simplified sublimit of the model presented in [16] in which the concentration of sulphur dioxide does not change across the sulphuric acid layers; it is derived by taking R → 0, B 2 → ∞ with RB 2 = O(1), and B 3 → 0 in the original nomenclature. We obtain subject to where β and κ are as defined in (27), h g denotes the thickness of the gas phase above the liquid layer of sulphuric acid and can be identified with V from our model, and χ 0 is the initial gas fraction in the filter, which, in the rectangular geometry shown in Fig. 15, corresponds to the ratio of the height of a catalyst channel to the distance between the bottom of the catalyst channel and the top of the adjacent catalyst sheet layer. For the purposes of our comparison, we identify χ 0 with V 0 . We again note that the main difference between our model and the simplified version of the model from [16] is that the "effective diffusivity" in the simplified model is linear with the gas fraction, h g . We also observe that, in one of the limits considered in [16], where sulphur dioxide concentration is allowed to vary across the microscopic liquid layers in the catalyst channels, a linear profile of the clogging front with time is recovered in contrast to the convex shape we obtained here. In Fig. 16, we show plots of the temporal profiles of the outlet sulphur dioxide concentration and the position of the clogging front comparing our model and the simplified version of the model from [16]. We see that both models give approximately identical behaviour until the clogging time. The qualitative behaviour afterwards is the same for both models, but there is a small quantitative difference in the complete clogging time, where the model from [16] predicts a 10% increase in the complete clogging time for the operating parameters we have used. This can be explained by the fact that the model (70)-(77) uses an effective diffusivity that is linear with the gas fraction and thus larger in magnitude than the convex diffusivity depending on V for our model shown in Fig. 4. This, however, shows that a simplified model for the catalytic filter is appropriate to capture the qualitative behaviour of the device and be reasonably close to the quantitative predictions from the more rigorous homogenised model presented in this paper.

Conclusions
In this paper, we presented a model for a reactive filter that removes sulphur dioxide from flue gas through a chemical reaction producing liquid sulphuric acid. We used the homogenised model derived in [17] and solved it until the filter becomes completely clogged with liquid. To achieve this, we obtained the effective diffusivity of sulphur dioxide in the porous filter sheet as a function of the gas fraction, until there is no single continuous gaseous pathway from the channel into the sheet. We found that the gas fraction in the filter sheet at this point is approximately 3%. After that, a continuous layer of liquid forms in the filter sheet at the channel surface, and the transport of sulphur dioxide is greatly reduced, since its diffusivity in the liquid is much smaller than the diffusivity in the gas, and the (b) (a) Fig. 16 Graphs of a the temporal profile of S at the outlet Z = 1, and b the temporal profile of the front G obtained from solving our model (solid black) and the model from [16] (dashed black). Here, β = 6.7, κ = 0.88, V 0 = 0.98 filter becomes inefficient. We defined the clogging time of the filter to be the first time when clogging occurs (which happens in the filter sheet at the inlet) and introduced a propagating clogging front along the filter channel, behind which the sulphur dioxide concentration is given by its inlet value, and the gas fraction corresponds to complete clogging. We measured the efficiency of the device by considering the outlet concentration of sulphur dioxide and the lifetime of the filter, i.e. the complete clogging time (that is, the time at which the whole surface of the filter sheet is completely saturated with liquid sulphuric acid).
We then solved the full problem numerically and observed that, under normal operating conditions, the time until complete clogging of the filter was almost four times longer than the time to first clogging. Looking at the gas fraction across the thickness of the filter sheet, we found that it decreased very little near the back of the filter, which suggested that the filter clogged only near the channel surface. We also found that the clogging front speed increased a little with time, which agrees with the results in [16], where a simple model for a reactive filter is presented.
We varied key dimensionless and dimensional parameters and explored their effect on the device performance. We compared the results and the efficiency of the filter in different parameter regimes by looking at both the outlet concentration of sulphur dioxide and the complete clogging time. We found that, as expected, decreasing the channel width or the speed of the gas decreases the outlet concentration and increases the complete clogging time. However, in practice, this would mean increasing the pressure drop across the filter, or slowing the process altogether, which might be costly to a company. Thus, we assume a given threshold outlet concentration that cannot be exceeded due to external requirements, and then use our results, fixing the channel width, for example, to find the maximum flow speed of the gas in order to achieve this in the initial stage of operation (or, equally, fixing the speed, we can find the maximum channel width of the filter). For example, using the normal operating parameters, but varying the speed, we find that the gas can flow at a maximum speed of approximately 0.2 m s −1 in order not to exceed an outlet concentration of 10% of the inlet value. We also varied the length of the filter and observed that this increased the efficiency of the filter, as expected, by decreasing the outlet concentration and increasing the complete clogging time of the device. We found that decreasing the reaction rate increased the outlet concentration to begin with but then decreased it for later times, and increased the complete clogging time, as expected. We also observed that the smaller the reaction rate, the higher the convexity of the graph of the position of the clogging front with time was, indicating faster clogging closer to the time of complete clogging. We found that the behaviour of the filter is sensitive to the changes in the inlet sulphur dioxide concentration. In particular, halving the outlet concentration resulted in an almost tripled complete clogging time. This suggests that it would be beneficial to have filters tailored and optimised for specific concentration ranges, rather than using the same filter for varying inlet concentrations. We also found that increasing the inter-pellet distance increased the outlet concentration to begin with but then decreased it in the late stage of the device operation, and substantially increased the complete clogging time. In particular, doubling the inter-pellet distance increased the complete clogging time by a factor of eight and increased the outlet concentration initially by almost three times. This means that placing the catalytic pellets farther apart will increase the time of the filter operation but at a higher outlet concentration, which is maintained for a longer time compared to the case of small inter-pellet distance. We also varied the thickness of the filter sheet and found that increasing the thickness beyond 4 mm had no effect on the behaviour of the filter. We obtained similar results in [16] for a simpler filter design made of catalytic channels rather than pellets. We also explored the relationship between the gas speed and the inlet sulphur dioxide concentration while keeping the behaviour of the filter the same as a way of control that can be applied when using various inlet concentrations of sulphur dioxide. For a sufficiently thick filter sheet (such as 4 mm as above), this relationship can be found by considering a single dimensionless number in the governing equations, which revealed that increasing the inlet concentration by a given factor, n, may be balanced by decreasing the gas speed by √ n. We were able to obtain an asymptotic reduction of the model assuming that the filter channels were thin (or the gas speed was slow) and the filter sheets were thin (or the reaction rate was small). We derived an explicit solution for the sulphur dioxide concentration in the filter channel, which was static until the clogging time, and the gas fraction in the filter sheet. A travelling-front solution was found after the clogging time, in which the position of the clogging front was given by numerically solving a Volterra integral equation. We also explored the effect of changing the reaction rate and the inter-pellet distance on the complete clogging time and found that increasing the reaction rate or decreasing the inter-pellet distance decreased the complete clogging time, where the decrease became very slow as the parameters were continuously changed.
We also compared the results of our model to those of a simplified version of the model presented in [16], which represents the filter as an array of rectangular catalyst channels rather than pellets. In the relevant limit, the main difference between the two models was the functional form of the effective diffusivity, which was linear in the simplified model and a convex non-linear function in the homogenised model. Both models agreed on the initial behaviour of the filter, and the simplified model predicted a 10% increase in the complete clogging time of the filter.
There are multiple ways our model can be refined. When finding the effective diffusivity, we assumed that, once two neighbouring liquid layers touched, capped spheres of liquid formed around the pellets until the continuous gaseous pathway for sulphur dioxide was broken up. In practice, as soon as such layers touch, there will be an immediate change in the topology of the liquid phase that would result in a more complicated convex shape that minimises the surface energy. Thus, a more careful examination of this phenomenon will be useful in evaluating the precise effective diffusivity, even though we anticipate little quantitative or qualitative changes in the diffusivity profile due to these local effects. In addition, when introducing the clogging front, we ignored the diffusion of sulphur dioxide through the macroscopic liquid layer that forms along the surface of the filter sheet, since its rate is very slow compared to that in the gas. In a more detailed model, we could capture this effect by explicitly tracking the diffusion of sulphur dioxide through that layer. We showed a comparison between our results and the results from [16] by considering the model from [16] in the limit studied here. An alternative comparison would be to consider our homogenised model in the limit of [16], in which the sulphur dioxide concentration varies across the liquid layers. In our model, we assumed a one-step chemical reaction, but in reality there are multiple intermediate steps.
Incorporating these in the chemical kinetics will improve the model's predictive power and will allow for a more accurate measurement of the chemical rates, which will be known for each reaction. Another useful modification to our model will be adding the effect of liquid drainage from the filter, which affects the rate at which liquid is accumulated in the filter sheet. A starting point to achieve this would be to consider a filter composed of hydrophilic fibres and pellets, which will facilitate transport within the filter (see, for example, [32]). Finally, sulphuric acid is hygroscopic, i.e. it absorbs moisture form the surroundings. Modelling this effect (see, for example, [33]) is important, as absorption of water is another clogging mechanism that will affect the performance of the filter.
Our models provide useful insight into behaviour of reactive filters and into optimising their performance, enabling the reduction of sulphur dioxide released into the atmosphere. material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.