In this section, we first provide numerical solutions to the ADEs and review existing approximations. We then develop new approximations for the full concentration profile.

A. Numerical solutions and approximations

Assuming abrupt connection to a reservoir at concentration CL0, the boundary conditions in the linear case are CL(x,0)=0; CL(1,t′)=CL0 and ∂CL/∂x|x=0=0 and the solution requires integration over 0≤x≤1, 0≤t′tmax′. Boundary conditions are similar in the radial case, replacing x with r. Numerical solutions are obtained in Matlab using the function “pdepe,” which can integrate parabolic and elliptic partial differential equations. The blue lines in Figs. 2 and 3 show concentration profiles on a logarithmic scale at different normalized times Δt′=t′−tF′ for 1D and 2D concentration, respectively, assuming that Pe=500. Here, tF′ is the filling time, discussed below.Near the reservoir (x,r≫0) all curves tend to a steady profile. In this region, approximate solutions CLx and CLr for the two cases can be obtained by neglecting time variation and diffusion, as CLxCL0≈1x;CLrCL0≈1r2.(9)These solutions are shown in green in Figs. 2 and 3. They match the numerical results well in the ramp regions but poorly elsewhere. The solution for the 1D geometry is known as the hyperbolic ramp.21–2421. J. Leng, B. Lonetti, P. Tabeling et al., Phys. Rev. Lett. 96, 084503 (2006). https://doi.org/10.1103/PhysRevLett.96.08450324. M. Schindler and A. Ajdari, Eur. Phys. J. E 28, 27 (2009). https://doi.org/10.1140/epje/i2008-10419-y In each case, concentration rises gradually with distance from the reservoir, but unfortunately tends to infinity at the origin. The rise is clearly faster for the radial case.Near the origin (x,r≪1), the profiles are Gaussian, with a peak value that increases linearly with time. Here, the following approximate solutions are valid, CLxCL0≈(t′−t′Fx)√(2Peπ)e−x2Pe/2,CLrCL0≈(t′−t′Fr)Pe4e−r2Pe/4.(10)Here, tFx′ and tFr′ are normalized filling times. Once again, the solution for the 1D geometry is already known.2121. J. Leng, B. Lonetti, P. Tabeling et al., Phys. Rev. Lett. 96, 084503 (2006). https://doi.org/10.1103/PhysRevLett.96.084503 These solutions are shown in red in Figs. 2 and 3. They match the numerical profiles near the concentration peak but are clearly inaccurate elsewhere. The peak is clearly wider in the radial case but, despite this, concentration factors are much larger.Filling times may be estimated by matching the peak concentrations CLMx=max(CLx) and CLMr=max(CLr) from Eq. (10) to numerical results. The full and dotted lines in Fig. 4 show numerical and analytic variations for 1D (LH) and 2D (RH) concentration and two values of Pe. In each case, the numerical variations gradually tend to the linear approximations, and the longer filling times but larger concentration rates achieved in the radial case should be noted.The points in Fig. 5 show the variations with ln⁡(Pe) of normalized filling time obtained by this matching, together with straight-line fits that allow the following analytic estimates: tFx′≈ln⁡(Pe);tFr′≈ln⁡(Pe/2).(11)Again, the expression for tFx′ is known.2525. J.-B. Salmon and J. Leng, J. Appl. Phys. 107, 084905 (2010). https://doi.org/10.1063/1.3354084 The estimates are accurate unless Pe is small. For Pe>4 (a trivial value), the filling time is longer for the radial case.

B. New analytic approximations

The approximations above are not satisfactory since they do not allow complete concentration profiles being drawn. We now develop new solutions that avoid this problem. To do so, we assume at the outset that the profile after filling can be written as the sum of time-varying and static parts. In the linear case, we, therefore, write for t′>tFx′, CLx(x,t′)≈(t′−t′Fx)C1(x)+C2(x).(12)Substitution into the ADE then yields the ordinary differential equation, C1≈(t′−t′Fx)+1Ped2C2dx2+xdC2dx+C2.(13)Since the time-dependent terms must vanish separately, C1 must satisfy 1Ped2C1dx2+xdC1dx+C1=0.(14)Direct substitution shows that C1(x)CL0=Aexp(−Pex22).(15)To match the filling rate for unit strip width, solute conservation requires that ∫01C1dx=CL0⋅1 so that A=2Pe/π. This result implies that the time-varying solution corresponds exactly to the upper of Eq. (10). The function C2 must then satisfy 1Ped2C2dx2+xdC2dx+C2=CL02Peπexp(−x2Pe2).(16)Surprisingly, this inhomogeneous equation can be integrated directly. Assuming the boundary conditions C2=∂C2∂x=0 on x=0, we obtain C2(x)CL0=Peexp(−x2Pe2)∫0xerf(yPe2)exp(y2Pe2)dy.(17)

It is simple to show that this result approximates a hyperbola as x→1, but the overall variation is now a realistic bounded function that tends to zero as x→0.

Performing a similar procedure for the radial case but now assuming that CLr(r)≈(t′−t′Fr)C1(r)+C2(r) for t′≥tFr′ and requiring that ∫01C12πrdr=CL0⋅2πr⋅1/2 (where the factor of ½ arises from the radial velocity profile) leads to C1(r)CL0=Pe4exp(−rPe24).(18)This result implies that the time-varying solution corresponds exactly to the lower of Eq. (10). Substitution and integration then yields C2(r)CL0=(Pe4)exp(−r2Pe4).(19)Here, ei(x)=∫−∞xettdt is the exponential integral and γ≈0.57721 is the Euler–Mascheroni constant. Equation (19) tends to 1/r2 as r→1, but the overall variation has again been replaced by a bounded function.Figures 6 and 7 show these analytic solutions at different normalized times Δt′ for linear and radial concentration, respectively, again assuming that Pe=500. These results should be compared with Figs. 2 and 3.

The new solutions present physically realistic approximations for the entire concentration profile and (despite some inaccuracy for small Δt′) tend to the numerical variations as Δt′ increases.

Using the analytic solutions, the filling times may be estimated as tFx′=∫01⁡C2dx∫01⁡C1dx;tFr′=∫01⁡C22πrdr∫01⁡C12πrdr.(20)Numerical integration shows that these expressions agree well with Eq. (11) and Fig. 5.We have not been able to perform an analytic integration; however, we note that Eqs. (17) and (19) may be approximated as C2(x)CL0≈1−e−(xPe/2)3x;C2(r)CL0≈1−e−12(Per24)2r2.(21)These expressions highlight the limiting behavior near the source and stagnation point more clearly. Usefully, they may both be integrated analytically, to yield tFx′≈ln⁡(Pe)+;tFr′≈ln(Pe2)+.(22)In each case, the first term may be recognized as the corresponding filling time in Eq. (11), confirming the dependence on Pe. The second term is a small constant error, which has the numerical value −0.1541 in the 1D case and −0.5778 in the 2D case.Although the 1D geometry has a shorter filling time, a flatter ramp, and a sharper concentration peak, these results imply that the radial geometry concentrates more quickly after filling. In fact, the ratio of the peak values CLMr and CLMx for t′≫tF′ isThis ratio is unity for Pe≈10 and 10 for Pe≈1000. Radial-flow evaporation concentrators are, therefore, fundamentally more effective than linear concentrators, and this advantage is clear in the earlier Fig. 4. The only uncertainty is whether a suitably high Pe can be achieved experimentally before permeability effects limit capillary rise to an intermediate stagnation point.2929. R. Syms, Biomicrofluidics 11, 044116 (2017). https://doi.org/10.1063/1.4989627 We now confirm this is the case.

C. Model limitations

We now consider three limitations of the model. First, the analyte supply is likely to be finite. ADEs can always be solved numerically for such cases. For example, the blue lines in Fig. 8 show concentration profiles at different times for a 1D concentrator with Pe=1000, assuming injection of an analyte slug for a time Δt′=0.25 followed by pure solvent. Solute may be seen entering from the source until t′=Δt′; after this, it travels to the origin by solvent pumping. The concentration peak now tends to a steady state, which from previous analysis must clearly be CLx≈Δt′C1(x). The prediction of this expression is shown in green in Fig. 8; the agreement with the final numerical curve is exact. A similar result is obtained for radial flow, confirming retention of its geometric advantage for a finite supply.Second, the velocity spread inherent in porous media results in Taylor dispersion, often modeled by an effective axial diffusion coefficient Da=D(1+κVm), where D is the molecular coefficient, V is the velocity, and κ and m are constants. For the velocity profile in a 1D concentrator, this can be written as Da=D(1+αxm), where α=κ(XM/τe)m. The effect of a varying diffusion coefficient Da=Df(x) is then to modify the ADE to ∂CL∂t′=(1Pe)∂∂x(f∂CL∂x)+x∂CL∂x+CL.(24)The previous method allows solutions for C1 and C2 to be found for the infinite source case, for both m=1 and m=2. The expressions are complicated, so here we give only those for C1, m=2;C1=A(αx2+1)Pe2α.(26)Here, A is a constant chosen to conserve the quantity of solute in the peak and is most easily found by numerical integration. Figure 9 shows this solution for Pe=1000,m=1, and different α.

The effect of increasing α is to broaden the concentration peak and reduce the concentration rate. However, because the velocity is inherently low near the origin, large values are required to make a significant difference. A similar ADE may be constructed for radial concentrators; this can be solved numerically, but we have so far failed to find analytic solutions.

Third, non-uniform evaporation may give rise to a nonlinear normalized flow velocity profile g(x). In this case, the ADE modifies to ∂CL∂t′=1Pe∂2CL∂x2+∂(gCL)∂x.(27)The previous method again allows solutions for C1 and C2 to be found for the infinite source case. Assuming the polynomial variation g(x)=a1x+a2x2+a3x3⋯ the solution for C1 is the modified Gaussian concentration profile,

For evaporation profiles varying mainly near the source, the effect on the peak shape is also small.