Complete analytic solutions for convection-diffusion-reaction-source equations without using an inverse Laplace transform

Transient mass-transfer phenomena occurring in natural and engineered systems consist of convection, diffusion, and reaction processes. The coupled phenomena can be described by using the unsteady convection-diffusion-reaction (CDR) equation, which is classified in mathematics as a linear, parabolic partial-differential equation. The availability of analytic solutions is limited to simple cases, e.g., unsteady diffusion and steady convective diffusion. The CDR equation has been considered analytically intractable, depending on the initial and boundary conditions. If spatial adsorption and desorption of matter are super-positioned in the CDR equation as sink and source functions, respectively, then the governing equation becomes an unsteady convection-diffusion-reaction-source (CDRS) equation, of which general solutions are unknown. In this study, a general 1D analytic solution of the CDRS equation is obtained by using a one-sided Laplace transform, by assuming constant diffusivity, velocity, and reactivity. This paper also provides a general formalism to derive 1D analytic solutions for Dirichlet/Dirichlet and Dirichlet/Neumann boundary conditions. Derivations of the analytic solutions are found to be straightforward if a combination of the source function and the initial concentration provide a non-zero singularity pole of inverse Laplace transform.

In models of transport phenomena, the mass-transfer rate at a given location is ascribed to diffusion, convection (or advection), reaction, and sources/sinks, which are ubiquitous in a plethora of natural and engineered processes. Diffusion can be described mathematically by using the transition probability describing locally hopping molecules 1 , random fluctuating forces satisfying the dissipation-fluctuation theorem 2,3 , and random walk phenomena caused by irreversible increases in entropy 4,5 . Convective transport originates from the motion of solute-carrying fluid, as obtained from the continuity equation. A chemical reaction is depicted as a continuous transformation of solutes toward a chemical equilibrium in a bulk phase. Finally, source and sink functions are ascribed to transporting masses that are either created or annihilated at specific locations in the environment. Obtaining rigorous solutions of the coupled phenomena is of great importance in various scientific and engineering disciplines. In past decades, analytic approaches to the coupled transport phenomena has been limited to coupled advection -diffusion without reactions: 1D unsteady in an open channel with spatially varying velocity and diffusivity 6 , 3D steady in a planetary layer (semi-analytic) 7 , and 1D unsteady with variable coefficients in semi -infinite media 8 .
A full mathematical expression of the above four mass transfer mechanisms is known as the convection-diffusion-reaction-source (CDRS) equation, wherein a sink can be expressed as a negative source. Specifically, a 1D-unsteady CDRS equation may be written as where C x t ( , ) is the concentration at time t and position x, and D 0 and V 0 are the constant diffusion coefficient and convective velocity, respectively, K 0 is a first-order reaction constant, and S x ( ) is a source function. Because Eq. (1) is a parabolic partial differential equation, if Dirichlet boundary conditions (BCs) are assumed, a specific solution www.nature.com/scientificreports www.nature.com/scientificreports/ depends on an initial condition (IC) expressed as = = C x t C x ( , 0) ( ) I and two BCs: one at the inlet, = C x t C t ( , ) () 0 0 , and the other at the outlet, = C x t C t ( , ) ( ) 1 1 , where < x x 0 1 is given for convenience. For example, Dirichlet-type conditions of finite ≠ C 0 0 and = C 0 1 indicate a non-zero concentration at the inlet and a perfect sink of the transferring mass at the outlet, respectively. In contrast, if a semi-infinite domain is considered such as → ∞ x 1 , the Dirichlet BC of = C 0 1 is often switched to Neumann-type BC, i.e., C x , which is typically known as the zero-flux or exit BC. In the applied mathematics literature, a number of unsteady diffusion problems have been solved analytically using Green's function and Laplace transform (LT) techniques 9 with convection phenomena often discarded, especially, in long-term diffusion phenomena. For mathematical simplicity, Eq. (1) can be rewritten as by using dimensionless quantities defined as where ∞ C is a reference (often bulk or initial) concentration, L is the length scale of the spatial domain, Pe is the Peclet number, and the factor 2 of λ is for mathematical convenience. In groundwater research, 1D sub-surface domains are either bounded ξ < < (0 1) or semi-infinite ξ < < +∞ (0 ). In this dimensionless analysis, the inverse time scale τ −1 and the dimensionless reactivity κ quadratically increase with respect to L, i.e., τ κ ∝ ∝ − L L and 2 2 so that κτ is independent of the domain length scale, L. In the limit of → ∞ L , an effect of diverging λ in Eq. (2) is often nullified by the induced zero-flux condition at the outlet. In this work, a general analytic solution of the CDRS Eq. (1) is obtained by using a one-sided LT without conducting the inverse LT explicitly. The steady-state solution, denoted φ ξ ( ) ss , is obtained by using the residue theorem in complex analysis at the primary singular pole. The unsteady part of the full solution is obtained as a product of a spatial and a transient function, using the initial concentration, defined as µ ξ φ ξ τ ≡ = ( ) ( , 0) I and the source function σ ξ ( ), both in dimensionless forms. If µ ξ ( ) I and σ ξ ( ) have forms of hyperbolic or sinusoidal functions, then the general analytic solution is obtained simply by identifying the secondary singularity pole of the LT  φ ξ τ [ ( , )] in the complex domain. Several representative examples available in the literature were reproduced to confirm the analytical rigor and numerical accuracy of the current method. The current solution method, however, still contains a discontinuity due to the discrepancy between the IC and BC(s) at boundaries, which can be technically resolved by using a Fourier series.

Solution for Dirichlet Boundary
where β λ = +z 2 . Equation (3) indicates that Φ has two singular poles at = z 0 and λ = − z 2 (equivalent to β = sinh 0). The residue of the first pole (at = z 0) is equal to the steady-state solution, denoted φ ss : , is the final solution satisfying the IC, but it does not provide the inlet BC of φ 0 = 1 for arbitrary τ. Therefore, a proper conditional representation of the full transient solution is which resolves the discrepancy of the φ ξ τ = = ( 0, 0) value from µ = 0 0 and φ 0 = 1. To avoid the above conditional representation, an alternative expression was obtained by using a Fourier series 11 www.nature.com/scientificreports www.nature.com/scientificreports/ as an improved analytic form from the previous work 10 . Figure 1 shows the time-evolution of φ along ξ from the initial τ = ( 0) to the steady states. The initial profile of µ = 0 0 is shown as a horizontal line on the ξ axis, having a discontinuity with the BC of φ τ = (0, ) 1. This concentration jump at the inlet (ξ = 0) prevails until the steady state is reached. As φ increases from 0 to the steady-state profile over the entire domain except the boundaries shown in Fig. 2, the dimensionless time scale required to reach the steady state is graphically observed between 0 τ = and 10 (of an order of λ O(2/ ) 2 ). The conditional expression of Eq. (6) has a simpler transient dependence as a product of the steady-state solution and exponential transient decay, but the Fourier series of Eq. (7) resolves the discontinuity of an L-shaped initial concentration profile at the inlet boundary by using an infinite series. This problem is extended for non-zero values of φ 1 , 0 µ , and κ, for which the analytic solution is obtained as (8) e e e e e e e e (0 1, ) sinh ( (1 )) s inh ( ) sinh [1 ] sinh ( (1 )) s inh ( ) (See the Supplementary Information for detailed derivations.) In Eq. (8), the first term on the right-hand-side represents the spatial variation of φ from the two boundary values of 0 φ and 1 φ , which becomes the steady-state solution as τ → ∞; and the second and third terms show the variation in φ ξ τ ( , ) originating from the initial concentration µ 0 . At τ = 0, the first and second terms disappear and the third term confirms the initial concentration µ 0 . Equation (8) converges analytically to the previous Carslaw case by setting φ = 1 Figure 3 shows a similar concentration profile to that of Fig. 1 Figure 4 shows a 3D plot of φ versus τ and φ, emphasizing the independence of the steady state on the initial conditions. After τ exceeds around 1.0, the transient solution start converging to that of the steady state. Figure 5 considers a more general case in which the exit boundary value is non-zero, e.g, . A finite φ 1 may have better practicality in heat transfer applications, and therefore it is investigated here to confirm the mathematical robustness of the analytic solution of Eq. (8). In Figs. 5 and 6, the concentration φ increases from its initial value µ = .
0 2 0 , meeting the BCs of finite φ 0 and φ 1 and reaching the steady state. The time-evolution trend of Fig. 5 is similar to those in Figs. 1 and 3, showing φ ss monotonically decreasing from the inlet to the outlet. Figure 6 also indicates that φ ξ τ ( , ) starts converging to the steady state when τ passes around 1.
The inverse Laplace transform (iLT) is often challenging to obtain both analytically and numerically, owing to the sophisticated nature of singularity identification 9,10 and numerical sensitivity depending on specific algorithms applied. The uniqueness of the current approach is that the IC is used not only for the LT of the dimensionless governing Eq. (36) but also for the derivation of the unsteady part of the analytic solution without directly performing an iLT.

Solution for Mixed Dirichlet-Neumann Boundary Conditions.
The approximate semi-unbound system is investigated by considering a Neumann-type BC at the outlet (instead of φ = 0 1 ) such as is set for simplicity, for which the LT is nullified. A set of the inlet-Dirichlet and outlet-Neumann BCs gives a matrix form to obtain coefficients B 1 and B 2 of the complementary part: Then, the iLT is simply written as   . The denominator of Eq. (12) increases with κ (through α) having a fixed λ, thus confirming that the steady -state concentration φ ξ ( ) ss reduces for high reactivity κ. The time scale to reach the steady state is of an order of α − O( ) 2 not only for this zero outlet flux case but also for other previous cases. A full solution with non-zero flux J 1 and constant µ 0 is calculated as follows sinh sinh cosh (18) whose detailed derivation can be found in the Supplementary Information. Figure 7 shows the unsteady profiles of φ ξ τ ( , ) having the inlet -Dirichlet φ = ss 0 of which physical meanings are as follows. The concentration φ is maintained at the level of the inlet concentration φ 0 if the reaction is absent (of Eq. (20)) or the convection is predominant over diffusion and reaction (of Eq. (21)). The reaction process significantly decreases the outlet concentration due to the zero-flux in the absence of convection (of Eq. (22)). The influence of the exit BC on the evolution of φ to the steady-state is shown in Fig. 8. As the outlet concentration is not forcefully bounded by constant φ 1 , the overall profile φ rapidly increases from zero within a time interval, τ < ≤ 0 1 , instantaneously satisfying the zero outlet-flux.
Solution for a full CDRS equation with spatially varying source function. As discussed above, the analytic solution of the CDRS equation is subject to finding the steady-state solution corresponding to ξ Λ = z ( , 0) and calculating the unsteady part by using the IC. This argument seems to be, however, restricted to the governing CDRS equation with constant initial concentration without any source or sink functions. In this section, a sophisticated case is discussed for  having ˜= − a a 1 2 and γ = − z 2 3 4 . Note that the zero BCs relate the particular and the complementary solutions because the combined form should satisfy required conditions. Following the one-sided LT described in the Method section, the particular steady-state solution is obtained as φ with the two coefficients determined as 3 4 . The iLT of Φ C will produce the complementary solution φ C,ss in real space. The iLT of the full solution ξ Possible negative values of φ do not indicate that the solution φ ss has unacceptable physical meaning, but instead indicate that it is ascribed to the zero initial concentration. Although the source function σ does not explicitly depend on φ, finite κ (regardless of its sign) does not allow translational invariance of φ; i.e., φ and φ ε + , where ε is a constant, are not governed by the identical transport equation. The singularity of Eq. (24) occurs at 2= − − − ± , which is found by setting the denominator of Φ to zero. Then, the coefficient of the ξ b sin( ) term becomes~γ weighted to the steady-state solution, thus yielding the final solution of , a new axis is defined as χ πξ = and used for Figs. 9-11. The initial concentration µ = 0 0 is superceded on the horizontal axis of Fig. 9. The numerical solution for φ seems to reach zero at χ π = n /5, which was expected by the functional form of σ of Eq. (23), but the analytic solution of Eq. (30) clearly indicates that a phase shifts as much as θ 0 , which is calculated as 0.1365 in radians. At the first and second extrema at χ π = /10 and π 3 /10, respectively, the numerical solution deviates from the analytic solution, especially for τ = 1. Although the theoretical and simulational results show small discrepancies, it is worth noting that the differences are dominant where σ ξ ( ) decreases exponentially with sinusoidal fluctuations. Figure 10 shows variation of φ ξ τ ( , ) with respect to τ at three different locations. Both solutions show the rapid convergence of φ toward its state-state value after 0 2 τ > . . The rapidly fluctuating phenomena in the initial state ( 01) τ < . were not captured by the numerical solution, mostly because the time interval was not short enough. Zhong et al. indicated that slight mismatches between FEM and SGM results were due to unequal time-step sizes, however, specific grid size or time interval were not reported. In the presence of the source function such as σ of Eq. (23) decreasing stiffly from ξ = 0 for τ > 0, constant grid size may not be the best practice; instead, a better approach can be a variable grid size being inversely proportional to the concentration gradient in Scientific RepoRtS | (2020) 10:8040 | https://doi.org/10.1038/s41598-020-63982-w www.nature.com/scientificreports www.nature.com/scientificreports/ magnitude. Figure 11 clearly shows the initial fluctuation of φ with respect to τ along the χ direction, which was not visually observed in the numerical work of Zhong et al. 12 in Fig. 9. The spatial variation of φ is primarily due to the variation pattern of the source function σ and the transient behavior is ascribed to Eq. (35) of diminishing sinusoidal variation. Figure11 also implies the necessity of the variable grid-steps (in ξ direction) based on the concentration gradient.
A more complex form of a source function can be represented using multiple terms or a Fourier series. As the general solution for a source function of Eq. (23) is obtained in this study, it is straightforward to extend the current method for cases of functionally sophisticated source terms as long as the complexity of the functional forms provide no more than one non-zero singularity pole. In other words, because a linear combination of σ ξ ( ) and µ ξ ( ) I is included as an integrand to obtain I ξ p ( , ) and Φ ξ p ( , ), mathematical procedures to obtain the final analytic solution φ of the CDRS equations are straightforward if σ and µ I are expressed with sinusoidal functions and/or infinite Fourier series.

Solution Procedure of a CDRS Equation.
To solve the dimensionless governing Eq. (2), the IC is given , indicating a specified inlet concentration scaled by ∞ C and the perfect sink of mass at the outlet, respectively. Second, the outlet condition can be of the Neumann-type, by using a concentration  ) . Then, Eq. (2) is rewritten as and sup indicates the supremum. This Tauberian theorem for the LT indirectly supports the current derivation to obtain the steady state solution as both limits of → ∞ T and → ∞ R are analytically implemented. In the original work of Talbot 18 , the Bromwich integral 26 representation was used to transform the integration path into a parabolic parameterized curve that give fast convergence. An improved algorithm -at least twice faster than Talbot's work -was proposed by Trefethen et al. 27 with fewer parameters, which are more difficult to compute. A universal observation in the iLT to solve the parabolic partial differential equation such as CDRS equation is as follows: first, numerical algorithms for the iLT are available in the literature, but results are subject to specific parameters to be determined; second, as the iLT deals with theoretically infinite imaginary domain, the numerical accuracy cannot be improved by simply increasing the number of divisions of a real-space integration; and third, a numerical result of the iLT is hard to validate unless it is compared with numerical solutions independently obtained using different methods.
The current method is uniquely applicable to solving the CDRS equation -especially when the non-trivial steady state exists -of which mathematical generality needs to be further studied rigorously. Applications of the current method to conventional examples successfully reproduced analytical and numerical results found in literature, especially for cases that the initial concentration and a source function are represented as complex Fourier series or a simpler combination of sinusoidal and hyperbolic functions, of which combination provides a non-zero singularity pole during the iLT.