Mathematical analysis of a Sips-based model for column adsorption

We investigate the dynamics of a contaminated fluid flowing through an adsorption column. We derive a one dimensional advection–diffusion equation coupled to a sink term that accounts for the contaminant adsorption. The adsorption rate is modelled by the Sips equation, where the order of the exponents is obtained through analysing the chemical reaction (as opposed to data fitting). By applying a travelling wave substitution the governing equations are reduced to two coupled ordinary differential equations. After non-dimensionalising and imposing typical values for the operating parameters we are able to identify negligible terms and so reduce the system to one where explicit solutions can be obtained. Distinct solution forms are provided for a range of Sips exponents. The approximate solutions are verified against numerics and experimental data. It is shown that if the breakthrough data is plotted in the form suggested by the analytical solutions then the result is a straight line when the correct Sips exponents are employed. This provides a clear indicator of the correct adsorption model. The straight line form permits a trivial optimisation procedure to determine the adsorption rate. This is the first time that analytical solutions have been obtained for sink terms beyond the standard Langmuir and Linear Driving Force models. ©

Graphical abstract and Research highlights will be displayed in online search result lists, the online contents list and the online article, but will not appear in the article PDF file or print unless it is mentioned in the journal specific style requirement. They are displayed in the proof pdf for review purpose only.

a b s t r a c t
We investigate the dynamics of a contaminated fluid flowing through an adsorption column. We derive a one dimensional advection-diffusion equation coupled to a sink term that accounts for the contaminant adsorption. The adsorption rate is modelled by the Sips equation, where the order of the exponents is obtained through analysing the chemical reaction (as opposed to data fitting). By applying a travelling wave substitution the governing equations are reduced to two coupled ordinary differential equations. After non-dimensionalising and imposing typical values for the operating parameters we are able to identify negligible terms and so reduce the system to one where explicit solutions can be obtained. Distinct solution forms are provided for a range of Sips exponents. The approximate solutions are verified against numerics and experimental data. It is shown that if the breakthrough data is plotted in the form suggested by the analytical solutions then the result is a straight line when the correct Sips exponents are employed. This provides a clear indicator of the correct adsorption model. The straight line form permits a trivial optimisation procedure to determine the adsorption rate. This is the first time that analytical solutions have been obtained for sink terms beyond the standard Langmuir and Linear Driving Force models.

Introduction
Achieving a toxic free environment is essential to improve both population health and to fight against climate change [1][2][3]. Nowadays, fixed bed adsorption is one of the most widely used systems to remove contaminants from liquids, such as wastewater [4][5][6], or gases, such as greenhouse gases and specifically CO 2 from air [7], so it is of the utmost importance to improve and optimise these processes.
Adsorption is a surface phenomenon where a species attaches to the surface of a bulk adsorbent material by a physical process, such as the weak Van Der Waals force or electrostatic attraction, or by chemical bonding, see [8] for example. Frequently adsorption occurs as a result of both physical and chemical processes. The chemical reactions taking place between the free contaminant and the adsorbent matrix can be complex but they are often known and well understood, so the global orders of the reaction can be inferred in advance.
There exists a wealth of experimental, computational and mathematical studies in the literature dealing with the many forms of contaminants and adsorbents. However, the complexity of the underlying adsorption reactions can lead to the use of phenomenological laws which result in inconsistencies between * Corresponding author.
E-mail address: abel.valverde@upc.edu (A. Valverde). models [8,9]. Many studies involve a number of coupled differential equations which can only be treated numerically. In this paper our aim is to provide a simpler mathematical model which can accurately and consistently predict results from a range of adsorption experiments and so may be used to understand the effect of operating parameter regimes and, ultimately, lead to improvements in the process.
It is standard for mathematical models of column adsorption to be based on a one dimensional (in space) system of partial differential equations which describe the evolution of the (cross-sectionally averaged) concentration of contaminant in the fluid, c(x, t), and adsorbed contaminant, q(x, t). The variation of concentration has previously been derived by applying a mass balance, see [9][10][11] for example, coupling the advection of contaminant to a mass sink which represents the rate of attachment to the adsorbent. The form of this mass sink is key to the success of the model. If the sink term is represented by A classic adsorption model is attributed to Langmuir [16]. This defines Q (c, q) = k ad c(q m − q) − k de q . (2) It is based on the assumption that the adsorption mechanism is caused by physical interactions between a monolayer of adsorbate and the available sites on the surface of the adsorbent. The rate of adsorption is then proportional to the contaminant concentration, c, and the available sites, q m − q (where q m represents the maximum available sites). The constant of proportionality is denoted k ad . Desorption is proportional to the amount of material already adsorbed, with a constant k de .
In column adsorption tests there are two key experimental measurements which provide information about the process. The first is the isotherm, which is obtained experimentally by allowing the adsorption process to continue for a sufficiently long time so that the outlet concentration matches the inlet value. The second measurement is the outlet concentration, providing what is termed the breakthrough curve, c = c(L, t), where L denotes the position of the end of the column. For the isotherm, the column may be weighed at the end of the experiment (and compared with the original weight) to determine the equilibrium adsorbed mass q e , or equivalently integrating the measured outlet concentrations can provide this value. Mathematically the isotherm corresponds to the steady-state of (1), such that Q (c e , q e ) = 0, where subscript e denotes the equilibrium value. In column sorption c e is the inlet value, i.e. equilibrium occurs when the concentration matches the inlet value throughout the column. Experimental measurements of q e for different inlet concentrations provide the isotherm. Adsorption isotherms are discussed in detail in [8,[17][18][19]. Setting Q = 0 in (2) determines the Langmuir isotherm q e = q m K L c e where K L = k ad /k de . Comparison of the isotherm equation against experimental values permits the evaluation of q m , K L .
The well-known Linear Driving Force (LDF) model uses the Langmuir model to define This has a constant isotherm, q e , independent of the concentration. Clearly this is unrealistic in regions of low or zero concentration where adsorption is predicted even without the presence of contaminant, see the discussion in [9]. The Freundlich isotherm is an empirical relation between the adsorbed fraction and the concentration in the fluid at equilibrium [20][21][22]. It is expressed as where k F is a constant (with respect to c e ) and l an exponent. The downfall of Freundlich is the lack of a slowing down mechanism.
Increasing c e always increases q e . In an attempt to rectify this, Sips [23] proposed which corresponds to a sink of the form The exponents m and n are determined by the partial orders of the chemical reaction (and hence do not form part of any parameter fitting process). The method to identify m, n from the chemical reaction is discussed in Section 2.1. In the case m = n = 1 the Langmuir isotherm is reproduced, while with K L c m e ≪ 1 the Freundlich isotherm is retrieved. The LDF model corresponds to a constant q e , this appears in the limit K L c m e ≫ 1. Provided isotherm data has been employed to determine q m , K L there should remain only a single unknown in the system of equations, the adsorption rate, k ad , which may be determined by fitting the model to breakthrough data. If a single parameter permits close agreement with a large number of data points it may be viewed as a good indicator that the model and underlying assumptions are correct. If k ad remains relatively constant with respect to changes in c in then the sink model, which defines the k's as constants (with respect to c, q) of proportionality may be correct. If k ad varies with c in then the initial assumption has failed and a new sink must be sought. One of the main issues discussed in [9] is the variation of k ad with inlet concentration predicted by numerous investigations: a clear indication that an incorrect model has been applied. However, good agreement with data is often achieved using incorrect models by setting known values as unknowns to provide extra fitting parameters. The review of Patel [4, Table 3] summarises how model ''constants'' vary with inlet concentration, flow rate, temperature and column length. One conclusion of [9] is that it is crucial to compare against a number of data sets to verify that k ad does indeed remain approximately constant as operating conditions vary. The Sips model coupled to the advection-diffusion equation investigated in the following sections has previously been verified against multiple data sets for the adsorption of Cr(III) and toluene [9]. Here we add another system to this list, CO 2 adsorption on a polyethyleneimine which is a primary amine rich polymer. However, our goal here is not to prove the applicability of the Sips model, that has been achieved in [9], rather it is to determine analytical solutions for physically relevant exponents and to exploit these in identifying the correct adsorption model and unknown parameters.
The paper is organised as follows. In Section 2 we derive the three dimensional governing equations for the adsorption process and then perform the cross-sectional averaging which results in a one-dimensional advection-diffusion equation (1). Non-dimensionalising the system permits us to identify the dominant (and negligible) terms. In Section 3 we obtain travelling wave solutions to the reduced system for a range of physically realistic exponents for the Sips model. The approximate formulae are verified through comparison with the numerical solution of the full system. Then, in Section 4, we compare the analytical results against the experimental data of Sujan et al. [24] which deals with the adsorption of dry CO 2 . Finally, in Section 5 we present the conclusions of this work.

Derivation of the governing equations
A typical form for an industrial adsorption tank is presented in the left image of Fig. 1. Over most of the domain they are cylindrical and filled with some porous, adsorbent material capable of capturing contaminants. The adsorbent consists of randomly packed pellets or rough spheres of slightly varying size leading to an unknown, heterogeneous distribution. Consequently the packing and therefore fluid flow past the adsorbent is complex and so, in order to obtain a tractable mathematical model, several simplifications and assumptions must be imposed. In the following derivation we will make the six, generally standard, assumptions: (A1) The column is regarded as a circular cylinder (A2) The adsorbent distribution is approximated using an ensemble average, that is one over a number of cross-sections, so that we may define an average cross-section throughout the column defined by a constant void fraction. (A3) The fluid velocity field is parallel to the axis of the container and taken as constant.     Area denoted by |B i | and with boundary ∂B i .

n:
Normal vector to ∂A.  Diagonal diffusion matrix, with diagonal D x , D y , D z . u: Velocity field of the fluid. c: Average density of the contaminant over each cross-section. c in : Average density of the contaminant at the inlet. ϵ: Void fraction.
(A4) The walls of the reactor are impermeable, so the fluid enters the container at the inlet and exits only at the outlet. (A5) The concentration of contaminant is sufficiently low such that its removal does not affect the flow. (A6) The adsorption process is isothermal.
To derive the mathematical model we define the x axis as coinciding with the axis of the container. The void region between the solid is defined as A while B i , i = 1, . . . , n represents the set of areas occupied by the solid. Following assumptions (A1), (A2) each cross-section is equivalent, that is, independent of x and so one may think of the average solid region as a set of small cylinders oriented along the x axis, as depicted in the right of Fig. 1. The notation employed in the derivation is shown in Table 1.
We denote by ρ c the density of the contaminant carried by the fluid and D = (D x , D y , D z ) is the diagonal diffusion matrix. Since all sections are equal we shall assume that D x , D y , D z are constant along the column.
A mass balance over the fluid region A gives Due to the random nature of the packing we average over the column cross-section and define the average value of the contaminant density by c(x, t) where The left integral in (10) can be split into two terms, where D (y,z) denotes the 2 × 2 diagonal matrix with entries D y , D z . We note that, since A is independent of x (assumption (A2)), The Gauss-Green theorem [25] permits the second integral on the right hand side of (11) to be written as Since the column is impermeable, (assumption (A4)), the summation involves fluxes purely into the adsorbent material and does not include a term at the column boundary. Since the velocity field u is constant and parallel to the x-axis, assumptions (A3) . (14) The rate at which mass leaves the fluid region A, represented by the summation term must balance the mass gained by the solid occupying B. If we denote by m ad (x, t) and ρ c ad (x, t) the averaged mass and concentration density that has been adsorbed at a given section x at time t, then where |A col | is the cross-sectional area of the column. The rate of mass loss between fluid and adsorbent is then Theoretical models typically deal with a dimensionless contaminant density where ρ b is the density of the bulk adsorbent material (defined as the ratio of the mass of adsorbent to the column volume).
Replacing ρ c ad by q in the sink term of Eq. (14) and noting |A| = ϵ|A col | results in Eq. (16) is a standard mass balance: the amount of contaminant changes over time due to advection, diffusion and adsorption as it is converted to q. The key question then is, how is the contaminant adsorbed?

Sips adsorption model
As discussed in the introduction, several sink models are common in the literature: a number of which are reviewed in [9]. For the present analysis we will employ the Sips model given by Eq. (7). The exponents m ≥ 0, n ≥ 1 may be defined as integers and, rather than appearing as additional model fitting parameters, are calculated through the order of the chemical reaction. The Sips sink term may be related to adsorption kinetics of the form where C represents the contaminant in the fluid, A the available sites for reaction in the adsorbent surface (q m − q [mol/kg]) and m and n are the stoichiometric coefficients of the reaction. Note that in order for the Sips sink term to correctly describe the kinetics, the partial orders of the reaction must be equal to the stoichiometric coefficients. This might not always be true depending on the mechanism that leads to the global reaction (17). Moreover, the products should also result into an order n of reaction.
An example of this situation is the reaction between dry CO 2 and an amine group present on the surface of adsorbent polymers such as polyethylenimine (PEI). This kind of adsorbent has important applications in the removal of CO 2 from air and we will investigate an experimental study of this situation in Section 4. Primary and secondary amine reactions for this problem are discussed in [26][27][28][29][30]. In the case of a primary amine where R represents organic chains in the PEI impregnated adsorbent. Identifying the contaminant c as CO 2 , the available sites q m − q as RNH 2 we see that one gas molecule reacts with two solids to produce two solids (identified as q). Consequently we write m = 1, n = 2 and the sink term (7) becomes However, as we will demonstrate later, the kinetics of the reaction may also be inferred from the experimental isotherm or breakthrough data. This may be preferable when the details of the reaction are unclear or unknown.

Initial and boundary conditions
Initially the column is free from contaminant and no adsorption has taken place, indicating At the inlet we write a mass balance in terms of the original variables where ρ in is the concentration of the contaminant entering the filter, which we assume is constant, and D B stands for Brownian diffusion. To integrate over the column inlet we note that since contaminant only occupies the void region where c in = ρ in (since the region before the column is all void). The velocity inside the column is related to that outside by u in = ϵQ in /|A| = uϵ where Q in is the fluid flux just upstream of the inlet. Integrating Eq. (21) over the cross-section at For typical column adsorption operating conditions dispersion is orders of magnitudes greater than Brownian diffusion, see for instance [31], so we shall neglect D B . Since D x is constant it may be taken outside the integral and then we obtain the well-known Danckwerts condition [32] At the outlet, a similar argument suggests where c b is the breakthrough concentration. However c b (t) is the main quantity to be predicted (since the intention is to apply the model beyond the experimental study). Danckwerts therefore proposed a simpler, approximate condition Pearson [33] states that this is the correct condition although arrived at intuitively rather than rigorously. The indeterminacy arises due to the imposition of discontinuous dispersion coefficients. Consider the inlet region, far upstream the flow is constant and Brownian diffusion dominates, within the column the flow is diverted around numerous particles causing dispersion to dominate. In the immediate vicinity of the inlet the upstream flow must adjust slightly as the fluid moves to enter the void regions and this results in augmented diffusion. Consequently, rather than a discontinuity in diffusion/dispersion, there will be a narrow boundary layer around the inlet (and outlet) where the conditions adjust. Pearson reconciles the problem by imposing a piecewise continuous diffusion/dispersion coefficient over the boundary layer. Assuming continuity of contaminant concentration at the boundaries, in the limit that the boundary layer thickness tends to zero, conditions (24), (26) are retrieved. These boundary conditions will be employed in the numerical scheme of Section 3.3.
Eq. (16) is an advection-diffusion equation: the diffusion term indicates an infinite speed of propagation and hence c, theoretically, extends to infinity. However, due to the presence of the mass sink, the concentration could become negative (in which case the physical solution is cut off at the point where c = 0) or asymptote to zero. The issue of whether the mass sink results in a finite (when c is cut off) or infinite domain is discussed in detail in [9]. In Section 3 we will examine a travelling wave approximation to the system. A travelling wave is a specific form of wave moving with constant speed and retaining its shape for all time. This, together with the fact that Eq. (16) is an advection-diffusion equation, prevents the application of boundary conditions (24), (26) which apply at the (finite) fixed points x = 0, L. Instead, when we examine the travelling wave form we impose along with vanishing derivatives at both plus and minus infinity. These are consistent with Eqs. (24), (26) when extended far upstream and downstream. The conditions as x → −∞ indicate that sufficiently far from the wave front the concentration approaches the inlet value while q tends to the corresponding equilibrium value. As x → +∞ both quantities tend to some, as yet unspecified, constant value.

Non dimensional form
Setting c = c inĉ , q = q mq , x = Lx and t = τt, the non-dimensional model reads where the adsorption length and time-scales and the parameters Da, Pe -1 and δ are given by The parameter Da is the Damköhler number, noting that L/u is the advection time-scale it may be interpreted as the ratio of advection to reaction time-scales. Pe -1 is the inverse Peclet number, representing the relative importance of diffusion to advection. δ is an indicator of the relative importance of desorption to adsorption. In an efficient adsorbing system δ is expected to be small. However, we note that desorbing systems are also of practical interest, for example when regenerating an adsorbent or in erosive/leaching systems [14,15], in which case δ would be large.
The initial conditions in non dimensional form read, For the numerical scheme we will impose the boundary conditions These conditions will be used in the numerical resolution of Eqs. (28). For the travelling wave we employ along with vanishing derivatives at both extremes.

Travelling wave solution
After an initial transient the physical process of column adsorption follows a fixed pattern. Contaminant enters at the inlet and travels a certain distance until a region where adsorbent material is available to capture the contaminant. As adsorbent is exhausted the process slowly moves forward. On physical grounds this suggests a travelling wave form. This has been verified in a number of papers through comparison with numerical solutions for adsorption with a Linear Driving Force model [10,11], desorption with an LDF based on the contaminant concentration in the fluid [14] and adsorption with a Langmuir sink [9]. For the Sips model it will be verified in Section 3.3.
We now introduce the similarity variable η =x − s(t), where s(t) helps track the position of the wave. In models with a strong sink, such that the concentration reaches zero over a finite distance, s(t) is typically defined as the point whereĉ(s(t),t) = 0.
For weaker sinks the concentration may never reach zero, hence we apply a more general definition such thatĉ(s(t),t) =ĉ f wherê c f ∈ (0, 1). The speed of the travelling wave is given by v = ∂ˆt s and so s(t) = vt + s 0 (since the travelling wave does not hold at t = 0 we cannot assume that s 0 = 0). We now define wheret f ,x f are such thatĉ(x f ,t f ) =ĉ f . This choice merely defines the origin of the η-axis, that is η(x f ,t f ) = 0.
Changing variables in Eq. (28) and writingĉ(x,t) = F (η) and q(x,t) = G(η) leads to As stated earlier, in contrast to thex,t system, the travelling wave holds over an infinite domain. The inlet condition is then replaced by the first of conditions (32) which indicates The amount adsorbed tends to the equilibrium value (in nondimensional form) hence Applying these far-field limits to the second equation in (34) determineŝ which is simply the non-dimensional form of (6). The final two conditions of (32) indicate proportional to some power of c: as c → 0 so does the strength of the sink [9]. The LDF, which is independent of c, is an example of a strong sink but we do not examine that here.
Integrating the first equation in (34) and applying conditions (38) results in a first order equation Applying conditions (35), (36) the velocity of the travelling wave is obtained This demonstrates that the larger the value of F ∞ the faster the wave moves.
In practice adsorbent materials are chosen to actually adsorb a contaminant. Materials that leave a positive concentration over an infinite length are of little practical use, so from this point on we set F ∞ = G ∞ = 0 and we therefore choseĉ f = 1/2. This results in the position η = 0 defining the mid-point of the travelling wave F (0) = 1/2. However, we note that the more general form may be relevant to studies of desorption or regeneration of adsorbent materials. The mass balance and velocity may now be expressed as Eqs. (40) may be used to remove G from Eq. (34) to obtain a mass balance in terms of F alone,

Approximate solutions for the travelling wave
In most applications advection dominates over diffusion and therefore the inverse Péclet number Pe -1 ≪ 1, see for example [9,11,15,[34][35][36] and also Section 4 later in this paper. So we proceed by writing F = F 0 + Pe -1 F 1 + O(Pe −2 ). After noting that δq n e = (1 −q e ) n (see Eq. (37)) the leading order of Eq. (41) may be written We define the polynomial in F where a = 1/q e > 1, so that the solutions F 0 of (42) with initial condition F 0 (0) = 1/2 satisfy Eq. (44) can be implicitly solved for a number of realistic m, n combinations. The solutions are found after obtaining explicit expressions for the integral on the right hand side of (44), which depends on the number and location of the real and complex roots of the polynomial f (u). For example, if m ≤ n or In Appendix we provide proofs and results pertaining to the derivation of physically realistic solutions of Eq. (44). In Table 2 we provide the forms of Ψ (x) for different values of n and m.
Then we take the expressions from Table 2, vary F 0 ∈ (0, 1) and generate the corresponding η values. The three curves with m = 1 closely coincide for η > η 1/2 while for η < η 1/2 they diverge to a greater extent as n varies. The curve with m = 2 shows a much slower decay to zero as η → ∞. From Eq. (40) we have G 0 =q e F 0 , so the leading order solution for the adsorbed mass takes exactly the same form as that shown in Fig. 2 but slightly compressed on the vertical scale due to the factorq e < 1.

Breakthrough curves
In column adsorption tests there are two key experimental measurements which provide information about the process. The first is the breakthrough curve, which is the concentration at the outlet. The second is the isotherm, which measures the equilibrium adsorbed mass for a given inlet concentration. Breakthrough may be easily measured experimentally and can be used to verify (or not) model predictions. To convert the solutions of Table 2 to the physical domain and hence determine the breakthrough curve we start with Eq. (46). To convert back to thex,t system we note η =x − vt − s 0 where s 0 is unknown. Since the travelling wave does not hold at t = 0 the initial conditions cannot be applied to determine s 0 , consequently it is typical to specify the experimental timet 1/2 when the outlet concentrationĉ = 1/2. The outlet corresponds tox =L and we have defined η = 0 as  Table 2 Analytical forms of Ψ (x) defined in (45) obtained by using the formulae provided in Lemma 2 for m̸ =n and in Lemma 3 when m = n, see Appendix.
n, m values whereĉ b (t) =ĉ(L,t). This expression provides an exact solution, in implicit form, for the variation oft withĉ b , whereĉ b ∈ [0, 1], for a range of m, n combinations.
Using the results in Table 2, where a = 1/q e , one can then obtain the different forms of breakthrough curves: • n = m = 1: • n = 2, m = 1: • n = 2, m = 2: • n = 3, m = 1: In this case, as stated in Lemma 2 and shown in Table 2 the expression differs depending on whether 1/q e is less or greater than 4. However, according to expression (37) 1/q e = 1 + δ 1/n , so in most applications where desorption is small this number will be close to one or at least 1/q e < 4.
Hence all of the results for F 0 may be applied to G 0 but adjusted by the factorq e .

Verification of the travelling wave against the full numerical solution
The numerical solution of the system (28) was obtained using the method of lines (MOL). For this the spatial coordinate was discretised using the finite differences method (FDM) with a central second-order scheme. This produced a set of coupled ODEs which were solved using the MATLAB routine ode15s. The ode15s function is a variable-step, variable-order solver based on the numerical differentiation formulae of orders 1 to 5 and is recommended for stiff problems. The stability and convergence of the temporal integration was controlled by the MATLAB function, while the stability of the spatial discretisation (convective stability) was controlled through the cell Reynolds number R c = ∆x/ Pe -1 , where ∆x is the spatial step size [37].
The central second-order spatial discretisation between two positions x i at time t j reads The boundary conditions (31) at x = 0, L may be represented by The initial conditions areĉ i,j =q i,j = 0.
In Fig. 3 we present numerical and analytical predictions for the breakthrough curve for three combinations of m, n. Two curves are shown in each plot, the dotted line representing the numerical solution of the full system, the solid line representing the analytical expression using the results of Section 3.1 which neglect Pe -1 , Da. The numerical solution of the travelling wave ODE including Pe -1 , Da, Eq. (41), is not shown because it is indistinguishable from the full numerical solution. All three figures employ the same conditions and correspond to the experimental operating conditions described in Section 4, such thatq e = 0.85, Da = 10 −3 and Pe -1 = 0.2. Since L, τ depend on m, n, see (29), their values change between the curves despite constant operating conditions. We have taken Da, Pe -1 slightly higher than the values presented near the end of Section 4 so that at least some difference between the curves may be observed. Even with the high values it is clear that the travelling wave, which neglects Da, Pe -1 , provides excellent agreement with the numerical solution of the full system. These results demonstrate that the approximate travelling wave may be used to accurately represent the evolution of the concentration and so henceforth we will only employ the approximate travelling wave solution.

Validation against experimental data
As discussed in the introduction there are two key pieces of experimental information: the isotherm and the breakthrough curves. So now we will investigate the correspondence of the isotherm models against experimental data and then compare travelling wave solutions against experimental breakthrough curves. From here on we will work in dimensional form, so that the results follow the convention of the experimental studies.
It is almost obligatory to present breakthrough curves in studies of adsorption and much less common to present the full isotherm data. In this section of the paper we will validate our model against the breakthrough data of Sujan et al. [24]. This paper deals with the adsorption of dry CO 2 by fibres composed of cellulose acetate (CA) and silica (SiO 2 ) which have been enriched with a low content of polyethyleneimine (PEI) (0.7 g PEI/g SiO 2 ), i.e. a primary amine rich polymer. The adsorbent will be referred to as PEI-CA-SiO 2 from now on. It is assumed that attachment is through chemisorption where the CO 2 molecules react with the amine functional groups [24]. In keeping with the majority of papers in this field they omit the isotherm data, which we instead take from work of the same group reported in Kwon et al. [38]. The experimental set-up is slightly different in the two papers. The adsorbent in [24] is PEI-CA-SiO 2 , while in [38] it is PEI-SiO 2 (i.e. without cellulose acetate). Since the reaction is between the PEI and CO 2 and the same dose of PEI is used in both experiments (0.7 g PEI/g SiO 2 ) at a similar temperature (30 • C) it is reasonable to assume that the isotherms will be similar. In terms of the model this means that the ratio q e /q m should be approximately the same between experiments.
The operating conditions for the breakthrough curve shown in [24, Fig. 3a] are detailed in Table 3. The process is carried out in a circular cylinder with a constant average column void fraction under isothermal conditions at 35 • C, all in line with assumptions (A1, A2, A6). The inlet concentration (0.0154 mol/m 3 or 395 ppm) is sufficiently low that mass loss effects on the flow are negligible and the assumption (A5) of a constant fluid velocity therefore holds. The reaction between CO 2 and PEI was discussed in Section 2.1 where we obtained m, n = 1, 2 however in the following we will examine other combinations to demonstrate the different possible behaviours.

Isotherm analysis
The isotherm represents the equilibrium adsorbed mass for a given inlet concentration. The mathematical formulae for the where K L = k ad /k de is the Langmuir constant. Eq. (56) corresponds to Eq. (37) in dimensional form. Note that the shape of the curve is determined by the ratio m/n rather than a certain m or n. Thus, the study of the isotherm will only provide the relation between these partial orders, their actual value must be determined by some other method, such as the chemical reaction or the breakthrough curve.
To demonstrate the different forms for different values of m/n in Fig. 4 we present the fitting of the nonlinear isotherm expression (56) to the data points of [38] (shown as circles). For the four combinations m/n = 1, 1/2, 1/3, 2/3 we calculate the best fit to 1/q m , 1/(q m K 1/n L ) using Matlab's fit function. For large c e the best fit is provided by m/n = 1/2: in order to estimate the maximum q m it is essential to study the high concentration behaviour. The value of q m ranges between 1.65 mol/kg for m/n = 1, 1.73 for m/n = 1/2 and m/n = 1/3, to 1.93 for m/n = 2/3. This indicates that the variation of q m is moderate, with a maximum change of around 14%, and so does not strongly depend on the factor m/n. Thus, a very approximate value for q e /q m can be obtained regardless of the chosen sink term. Motivated by the fit with m = 1 and n = 2 we choose q m = 1.73 mol/kg and then reading off from Fig. 4 we obtain q e = 1.455 at 0.0154 mol/m 3 to Table 4 Results from the fittings in Fig. 5. give q e /q m ≈0.84. We assume that this value of q e /q m is typical for the PEI-CO 2 adsorption process. After substituting for q e /q m , m, n = 1, 2 and c e = 0.0154 the isotherm determines The high value of K L = k ad /k de is indicative of the dominance of adsorption over desorption. Taking q e /q m = 0.84, k de = k ad /1816 we now match our model to the breakthrough data of Sujan et al. [24] which now has only a single unknown, namely k ad .

Fitting results
We begin by writing the implicit expression for the breakthrough curve, Eq. (47), in dimensional form and the function Ψ is provided in Table 2 for different values of m, n. The only unknown is τ , which depends on k ad . To calculate τ we consider the set of experimental data at the breakthrough curve given by where (t i , c bi ) are the experimental breakthrough data points. Eq. (58) represents a straight line with gradient τ . If we plot the set of experimental data points in the form defined by S then, with the correct model, the data points should also follow a straight line and the gradient determines τ . The results of the linear fitting t − t 1/2 = τ S are presented in Fig. 5.
The solid line in Fig. 5 represents the analytical form (58) while the circles represent the data points S. Fig. 5(b) is perhaps the strongest vindication of the choice m, n = 1, 2: only when plotted in this form does the data follow a straight line. This suggests that if the chemical reaction is not known, or fully understood, then by simply plotting the data according to the various solution forms we may infer the ratio m/n. Table 4 shows the results of the fitting of Fig. 5. Although apparent from the figure the accuracy of the m = 1, n = 2 model is confirmed by the lowest value of the Sum of Squares Error (SSE) and highest R 2 . In [9] it is stated that a high R 2 is not a sufficient test of a model's goodness of fit: the value R 2 = 0.93 for the case m = n = 1 may suggest it as a possible contender for the isotherm however the SSE is significantly higher than for m = 1, n = 2 and from Fig. 5 it is clearly not the appropriate isotherm. Taking c in =0.0154 mol/m 3 , q m = q e /0.84 = 0.61/0.84 = 0.726 (see Table 3) the obtained slope τ = 1/(c in q m k ad ) = 605.5 corresponds to k ad = 0.145kg·m 3 /(mol 2 s). Using the value calculated in the previous section, K L = k ad /k de = 1816, results in k de = 8 × The length-scale L = 9.43 × 10 −4 m indicates that the reaction occurs over a region of width of the order 1 mm. The low value of the Damköhler number demonstrates that the reaction timescale is much greater than the flow time-scale and consequently the neglect of the time derivative term in the advection-diffusion equation is justified. The inverse Péclet number has been calculated using a dispersion coefficient D x =10 −5 which is based on the experimental findings on dispersion in fluids flowing through packed beds reported by Levenspiel [31]. The neglect of terms of order Pe -1 suggests possible errors of the order 11%. However, since the boundaries are fixed (hence the solution is exact at either end of the domain) the true error is lower. This is standard with asymptotic reductions, fixing the boundaries results in smaller errors than predicted by the size of neglected terms alone. As demonstrated in Fig. 3 even with Pe -1 = 0.2 the error is negligible. The small value of δ suggests a further reduction of the system may be possible, however, neglecting δ equates to neglecting desorption which plays a small but important role in the process and in quantifying the isotherm. Fig. 5 represents the breakthrough curve in a rather obscure form. In the literature breakthrough is invariably plotted as c(L, t)/c in , hence in Fig. 6 we present the results in this more standard form using the parameter values calculated above.
As expected the solution with m, n = 1, 2 provides excellent agreement with the data however the curve for m = n = 1 also shows good agreement for t ≤ t 1/2 . It is then possible that a curve of this form may be accepted as a suitable fit whereas the parameterisation used in Fig. 5 makes it quite clear that it is not correct. There are many published works where the agreement is good before or after (depending on the fitting criteria) the half-way point. The use of the incorrect isotherm may be one explanation for this.
Another way to place the present models in a more standard context is to write down the dimensional form of breakthrough.
For the case m = n = 1, which corresponds to the standard Langmuir isotherm, writing Eq. (48) in dimensional form we obtain This may be rearranged to give c(L, t) = where a = q m /q e > 1. This equation has not been reported previously in the literature. The results for m = 1, n = 3 and m = 2, n = 3 may be derived from Eqs. (51), (52).

Conclusions
In the literature it is standard to compare a mathematical model for the breakthrough curve with data and, from this comparison, determine unknown system parameters. The result can often lead to a mismatch with certain regions of the data and, much more importantly, the fitted parameters break the model assumptions, specifically that they remain constant with respect to the inlet concentration.
In this paper we take as our starting point the model developed in [9], which uses the Sips adsorption model, and generalise it to obtain analytical solutions for a range of exponents of the Sips model. The Sips isotherm has the advantage of being more general than the Langmuir, Linear Driving Force and Freundlich isotherms. It reduces to each of these models in certain limits and can be linked to chemical reactions if they occur.
The main conclusions of our work are: 1. The governing equations may be reduced to a pair of ordinary differential equations. Neglecting dispersion leads to a first order system, which accurately represents the full system, but permits analytical solutions. For certain combinations of the Sips exponents, m and n, the first order system is found to have heteroclinic orbits that connect the initial value of the concentration to a state where all contaminant is removed: these correspond to a good adsorbent material. Exact solutions in implicit form may be found for physically meaningful combinations of m, n. 2. Plotting the breakthrough data in the form provided by the relations between t and Ψ (expression (58)) should result in a straight line when the correct Sips exponents are employed. This is a clear indicator of the correct adsorption model. The more standard approach, plotting the outlet concentration against time, is much less clear and can easily result in the use of incorrect sink parameters. The straight line form permits a trivial optimisation procedure to determine the adsorption rate. 3. It is important to use as much data as possible when examining an adsorption process. Knowledge of the isotherm and chemical reaction can be used to evaluate the Sips exponents, the ratio of adsorption to desorption and the maximum amount that may be adsorbed.
Published papers tend to focus on breakthrough curves, whereas invaluable information contained in the isotherm or chemical reaction may be wasted. By combining all available information a much more reliable, generally applicable result can be obtained. The combination of information can also act to check consistency of predicted values, for example if the chemistry suggests values for m, n, does the corresponding plot of the data points for t against Ψ result in a straight line?
In summary, in this paper we have developed a novel set of solutions which correspond to a variety of Sips adsorption models. By combining information from the chemistry, isotherm and breakthrough we found remarkable agreement with experimental data for CO 2 capture at low levels. In practice not all information may be available, in fact we could have obtained all unknowns from the breakthrough data but the added information provided a valuable consistency check. Clearly now the aim is to apply the model to a wide range of processes, with the aim of better understanding and improving adsorption techniques. The ultimate goal is, of course, to make the technology sufficiently efficient so that it may be used extensively in the fight against climate change and environmental pollution.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The experimental data used in this paper has been reported by Sujan et al. [24] and Kwon et al. [38].

Lemma 1.
Let m, n be positive integers such that 1 ≤ m < n, and let f (x), x ∈ R, be the polynomial of degree n + m defined by where a ∈ R, a > 1. Then: (iii) if n − m is odd, f has no real negative roots; (iv) if n − m is even, there exists a value a * (depending on n and m) such that for a ∈ (1, a * ) f has no real negative roots, for a = a * , x 1 = (m−n)a * m is the unique negative zero of f with multiplicity 2, and for a > a * f has exactly two real negative roots, x 1 , x 2 ; (v) all the complex roots of f have multiplicity 1.
Therefore, the first statement is proved. In order to prove the second and third statement, we need to count the number of real roots of the polynomial p defined in (67). We write p(x) = p 1 (x) + p 2 (x), where p 1 (x) = (a − 1) n x r and p 2 (x) = −(a − x) n , are both real analytic functions on any interval I = [c, d] ⊂ R.
We denote by N I (p) the number of roots of p on a given interval I, counted according to multiplicity. Then, using Theorem 2 of Voorhoeve and Van Der Poorten [39] we have that N I (f ) ≤ 1 + N I (W (p 1 , p 2 )), where W (p 1 , p 2 ) denotes the Wronskian determinant associated to p 1 and p 2 . A straightforward calculation gives W (p 1 , p 2 ) = (a − 1) n x r−1 (a − x) n−1 (mx + ra).
Notice that the roots of W (p 1 , p 2 ) are x = 0, x = a and x = −ra/m < 0.
For a given ϵ > 0 small enough and M > 0 big enough, we define the following three intervals: I 1 = [−M, −ϵ], I 2 = [ϵ, a−ϵ] and I 3 = [a + ϵ, M]. We will count the number of roots of p in each one of the intervals I j , j = 1, 2, 3.
Notice that p(a) > 0 for any value of n. Therefore: -If n is even, lim x→+∞ p(x) = −∞, so that for M big enough there is one change of sign and N I 3 (p) = 1.
-If n is odd, using the fact that p(x) > 0 for x > a we obtain N I 3 (p) = 0.
As the number of roots is valid for any ϵ > 0 and M > 0, the second statement is proved.
With respect to the claim of the third item, if r = n − m is odd, p(x) < 0 for x < 0, so there are no other real roots and the statement is proved.
With respect to the claim of the fourth item, we consider r = n − m even and the interval I 1 . Due to the fact that N I 1 (W (p 1 , p 2 )) = 1, we have that N I 1 (p) ≤ 2.
In fact, we will prove that, depending on the parameter a, p can have none, one (double) or two (simple) negative roots. We write so that we look for negative solutions of g(x) = (a − 1) n .
It is straightforward that Using the fact that g(−ra/m) is an increasing function with respect to a, and comparing with (a −1) n it is clear that there exists a unique value a * > 1 such that g ( −ra * m ) = (a * − 1) n , and that for a ∈ (1, a * ), g(x) > (a − 1) n , and for a > a * , g(x) = (a − 1) n has two negative solutions. This ends the proof of the claim of the fourth item.
In order to prove the final statement, let z 0 ∈ C\R be a zero of the polynomial p and N z 0 (p) its multiplicity. Then, using Theorem 1 of Voorhoeve and Van Der Poorten [39], N z 0 (p) ≤ 1 + N z 0 (W (p 1 , p 2 )).
(iii) if n is even and m − n is odd f has exactly three no null real roots, x 1 , x 2 ∈ (0, a), x 3 > a. If a < m/(m − n), x 1 < a(m − n)/m < x 2 = 1; if a > m/(m − n), x 1 = 1 < a(m − n)/m < x 2 ; if a = m/(m − n), x 1 = x 2 = 1. (iv) if n is odd and m − n is even f has exactly three no null real roots, x 1 < 0, x 2 , x 3 ∈ (0, a). If a < m/(m − n), (v) if n and m − n are odd, f has exactly two no null real roots x 1 , x 2 ∈ (0, a). If a < m/(m −n), Furthermore, all the complex roots are simple.
Proof. The first claim follows from writing where r = m − n. Furthermore, p(1) = 0 and p ′ (x) = 0 if and only if x = ra/m. The claims on the number and multiplicity of real roots of f follow easily from the behaviour of x r (a − x) n . Finally the last statement can be proved using similar arguments as in Lemma 1. □