The scaling of kinetic and transport behaviors in the solution-phase chemistry of a plasma–liquid interface

The interactions between atmospheric-pressure, low-temperature, non-equilibrium plasmas and liquids have garnered attention due to promising applications including chemical synthesis,^1,2 nanoparticle and nanomaterial synthesis,³ wastewater treatment,⁴ fertilizer production (plasma agriculture),⁵ and medical therapeutics (plasma medicine).⁶ The mechanisms behind these applications often involve extremely reactive, short-lived species such as solvated electrons (⁠$eaq−$⁠),⁷ hydroxyl radicals (•OH),⁸ and hydrogen atoms (H•)⁹ introduced to or generated at the interface by the plasma. Further, longer-lived reactive nitrogen and oxygen species,¹⁰ such as hydrogen peroxide (H₂O₂),⁹ often arise from intermediate reactions of the short-lived species and are also important. These species react with chemical reactants in the liquid, referred to here as scavengers (S), leading to the eventual targeted outcome. This has led to an extensive study of how these plasma-generated species are formed¹¹ and what role they play in each application. Attempts to characterize systems have resulted in experiments and simulations that aim to explain their average lifetimes,¹² concentrations,^13,14 and temporal and spatial distributions,^14–16 all in an effort to more effectively design plasma–liquid systems. Simulations, in particular, have been useful to support the interpretation of experimental results. Some of these include modeling the concentration and penetration depth of $eaq−$⁠, and comparing them to spectroscopic measurements,¹⁵ studying the formation of plasmas in bubbles, including the formation of chemical species in the bubble and in the liquid,^17,18 and studying the transport of plasma species into droplets.¹⁹ In pulsed plasma systems, simulations have also clarified how species are generated between each pulse, and how long-lived species may accumulate in the system.²⁰ 

However, while characterization efforts are important for describing these systems, there is a need for analytical models that provide insight into general relationships between system variables and system outputs. Such generalizations, especially, when non-dimensional, reveal how variables are intrinsically connected and provide a framework for tailoring and optimizing a system. For example, it is often desirable to determine conditions that favor certain chemical reactions, i.e., selectivity, to increase the yields of a process as well as its efficiency. This is true for both synthesis and degradation of chemicals, in which increased conversion has been achieved by varying the system current,¹ or by designing multiphase reaction environments and igniting a plasma across bubbles or foam in the liquid.²¹ In the field of wastewater treatment, in particular, this issue has been extensively studied with a special emphasis on the power consumed in the process. Several reactor designs have been studied, which include bubbling,^21–23 flowing liquid through the gas-phase electrode,²⁴ increasing the area of electrodes and the area of plasma–liquid contact,²² or a combination of these,^22–24 all for the purpose of increasing the energy efficiency of contaminant degradation. Most of these improvements have relied on iterative improvements and heuristic design principles based upon physical experiments where issues like transport limitations are recognized but not quantified. In cases like these, it is desirable to understand, from a mathematical perspective, how the system variables affect the conversion efficiency and how these can be manipulated to optimize a plasma–liquid system, regardless of the application.

In our previous work, we studied the effect that the scavenger concentration and the current density had on the conversion and faradaic efficiency in a plasma electrolytic cell.^25,26 One takeaway was that, to increase the conversion efficiency, high reactant concentrations in the liquid phase should be used whenever possible. However, this is not always possible since the scavenger concentration may affect the purity or nature of the product, as in the case of nanoparticle synthesis.²⁷ Alternatively, the concentration might be predetermined, as in the case of contaminants in wastewater treatment. Additionally, those results were based on relatively long-time exposures of the plasma to the liquid (1–10 min) and assumed a uniform boundary flux over a defined area as a simplification of much more complex plasma profiles.^28,29 In this work, we take this analysis further by approaching the problem with finer resolution, both spatially and temporally, to understand how system variables affect the evolution of local reaction rates and conversion yields.

The results here are derived from one-dimensional (1D) simulations of the non-linear reaction and diffusion equations for plasma-injected solvated electrons (⁠$eaq−$⁠) and solution scavenging reactants (S), where the primary competing reactions are second-order recombination of $eaq−$ and reduction reactions between $eaq−$ and scavengers. Our findings show that initially, these systems achieve 100% faradaic (reduction) efficiency and high reaction rates due to the large local initial concentration of scavenging reactants at the interface. However, as time progresses, the scavenger concentration depletes and the efficiency and rates become dependent on the diffusion of scavengers from the bulk to the interface. Theoretical analysis, confirmed by simulations, predicts a characteristic time for this transition that scales proportionally with the scavenger diffusivity (D_s) and the square of the scavenger bulk concentration (S_B²) and is inversely proportional to the electron flux squared (1/J²). The scaling proves to be remarkably accurate when the ratio of the scavenger reaction rate to that of second-order recombination of $eaq−$ is high (high S_B, low flux, or a combination of both). The time to depletion, however, decreases as the reaction rates become comparable in magnitude, making this an effective upper time limit for local depletion. Additionally, because this analysis is solely based on the reaction–diffusion of the species in the liquid phase, its corollary should apply to •OH as much as to $eaq−$ since their diffusivity and reactivity are similar.^30,31 Finally, this analysis is expanded to a three-dimensional (3D) plasma–liquid interface and the implications it may have in applications or reactor design are analyzed and discussed.

The simulations studied in this work revolve around the chemical reactions driven by $eaq−$ in the solution phase of a plasma–liquid system. Specifically, a direct-current, negatively biased plasma with a liquid water anode electrode is used as a model system (as illustrated in Fig. 1). The primary reaction of interest can be described by the general form,

where S represents any chemical scavenger being reduced by $eaq−$⁠. Although hundreds of chemical reactions have this form,³² here we take as reference the dissociative attachment to chloroacetate by solvated electrons,

which we studied experimentally in our previous work.²⁵ We do not include any competing •OH reaction, as those results showed that it did not affect the faradaic efficiency of the chloroacetate reaction. The main competing reaction, included here, is the second-order recombination of $eaq−$⁠,

which corresponds to hydrogen evolution in the electrolytic system.

One-dimensional (1D) behavior of the scavenging and $eaq−$ species into the solution (z-direction) can be described by two coupled partial differential equations, the reaction–diffusion equation for the chemical scavengers,

and the reaction–diffusion equation for the solvated electrons,

where [S] and [$eaq−$] are the scavenger and solvated electron concentrations, respectively, and a function of (z,t). Diffusion of each species is described by their respective diffusivities (D_s or D_e), and k_s and k_e are the reaction rate constants for scavenging (1) and second-order recombination (3) reactions.

Here, we solve Eqs. (4) and (5) in a simulated electrochemical cell with a solution depth of z_max = 1 cm to replicate the real length of the electrochemical cells used in our prior experimental work.²⁵ The boundary conditions for the scavengers [Eq. (4)] consist of two zero flux boundary conditions at $z=0$ and $z=1cm$⁠, representing the plasma–liquid interface and bottom of the reactor cell, respectively. The initial condition for the scavengers was set as a uniform value throughout the domain (S_B), which was varied between 1 and 100 mM for different simulated experiments. The boundary conditions for $eaq−$ [Eq. (5)] consist of a zero flux boundary condition at z = 1 cm at the bottom of the reaction cell, and a constant value flux applied at the plasma–liquid interface (z = 0), corresponding to the molar flux of $eaq−$ entering the solution from the plasma. Here, the molar flux of electrons is referred to as $J/F$⁠, where J is the plasma current density and F is Faraday's constant. For these simulated experiments, $J/F$ was varied between 10⁻⁶ and 10⁻¹ mol m⁻² s⁻¹ (10¹⁸ and 10²³ m⁻² s⁻¹). Initially, the $eaq−$ concentration is uniformly zero throughout the domain. The diffusivity for $eaq−$ is set to D_e = 4.80 × 10⁻⁹ m² s⁻¹,^30,33 and the diffusivity of the scavenger is set to D_s = 1.12 × 10⁻⁹ m² s⁻¹,³⁴ with sensitivity analyses going from 1.12 × 10⁻⁸ to 1.12 × 10⁻¹⁰ m² s⁻¹ as a way to explore different experimental conditions. Finally, the reaction rate for the scavenging reaction, the dissociative attachment to chloroacetate (2), was set to k_s = 1.48 × 10⁹ M⁻¹ s⁻¹, as measured in our previous work,²⁵ and was similarly varied between 1.48 × 10⁸ and 1.48 × 10¹⁰ M⁻¹ s⁻¹. The reaction rate for second-order recombination (3) was set as k_e = 1.1 × 10¹⁰ M⁻¹ s⁻¹.³² 

Since real plasma–liquid systems are three-dimensional, this 1D representation of the reaction and diffusion of chemical species in the liquid phase will differ from a physical system primarily by ignoring lateral or radial diffusion (i.e., orthogonal to $z$⁠). If the initial and boundary conditions were uniform throughout the boundary, this would have no effect; however, most plasma–liquid systems, including the one on which this study is based, have a non-uniform boundary flux, J(r), which might lead to behavior not captured in this model. In many cases, including our previous work, radial differences in molar electron flux occur over distances on the order of millimeters,²⁵ whereas the penetration depth of electrons is on the order of tens of nanometers.³⁵ Because of this, the simplicity offered by a 1D model was prioritized over simulating the effects of radial diffusion, and this is discussed further in Sec. III D.

The coupled partial differential equations given by Eqs. (4) and (5) were solved numerically by employing a modified version of MATLAB's native pdepe function. The pdepe function discretizes the equations in space, yielding a set of ordinary differential equations (ODEs). These resulting ODEs are integrated using MATLAB's native differential-algebraic equation solver ode15s, which integrates using a variable time step. However, the solution is always returned on the fixed position-time mesh specified by the user.

To account for the large range of scales of interest, a non-uniform spatial mesh was used. The spatial grid was divided into three sections with fixed mesh parameters, spanning the full 1 cm depth of the domain. The finest mesh (Δz = 1 nm) was set near the plasma–liquid interface (z = 0 to z = 1 μm) in order to capture in detail the steep concentration gradients formed there. Simulations always began with a very small time step (⁠$Δt=1μs$⁠) to capture the rapid depletion of the scavenger at small times. While ode15s adjusts the time step dynamically, the step sizes were kept small by requiring the solution at short intervals (⁠$Δt=1μs$⁠). We confirmed that the solution was mesh independent and also compared simulations to our simplified analytical model that linearizes the second-order terms, as reported in prior work.³⁵ Good agreement under conditions where the simplification is valid indicated the simulations were accurate. More details on both benchmarking studies can be found in Sec. S.1 of the supplementary material.

To understand what happens at a specific point at the interface, with a defined incoming $eaq−$ flux, the coupled Eqs. (4) and (5) were solved to determine the spatial and temporal variations of the concentrations $[eaq−(z,t)]$ and $[S(z,t)]$⁠. For elementary reactions, the reaction rate per unit volume can be calculated by the multiplication of each of the reactant concentrations, elevated to their respective stoichiometric coefficients, and integrated over that volume. Here, the instantaneous reaction rate per unit area (⁠$R′′(t)$⁠) for the scavenger reduction [Eq. (1)] in a 1D system is calculated as the integration of $[eaq−(z,t)]$ and $[S(z,t)]$ over the entire domain,

The yield of the scavenger reduction reaction is then defined as

where $R′′$ is divided by the molar electron flux (J/F). We use these two parameters as our primary outputs in this analysis.

Building upon our previous work where we studied the effect of system variables on the faradaic efficiency,²⁵ we wanted to understand the way reactions behave on a more localized level and from a transient perspective. Figure 2 shows how the reaction rate per unit area and the yield vary in time for different molar electron fluxes (J/F). Notably, the resulting reaction rates are approximately equal to the molar electron flux itself at early times but decay from this value as time progresses [Fig. 2(a)]. Similarly, the local yield, which is normalized by the molar electron flux [Eq. (7)], is close to 100% at these early times and then falls below 100% as the reaction rate decays [Fig. 2(b)]. This behavior arises as the scavenger concentration is depleted to near zero at the plasma–liquid interface z = 0 nm (see Sec. S.2 of the supplementary material for the time evolution of S concentration profiles). The depletion time, which we arbitrarily defined as the time at which the yield falls to 99% (t₉₉), is identified by green circles in Fig. 2. It is dependent on the molar electron flux, with depletion occurring faster for higher fluxes.

Nondimensionalization of differential equations is often used to extract characteristic quantities intrinsic to a system and to describe how they scale with other important variables. This leads to parameters that are useful in analyzing the system and in understanding their relative importance. Nondimensionalization may also be used to transform the equation into one which is independent of most system parameters, yielding a solution which works not only for a particular experiment but for a large range of operating conditions.

In order to understand the nature of the depletion time t₉₉, nondimensional analysis was applied to the reaction–diffusion equation for the scavenger [Eq. (4)] using the following nondimensional variables:

where S_B is the bulk concentration of the scavenger in the liquid;

where t_c is an undefined characteristic time of the system; and

where z_c is an undefined characteristic length.

The last variable in Eq. (4) is the concentration of $eaq−$ and nondimensionalizing it is less obvious as it must account for two issues. First, the $eaq−$ concentration is not uniform across the plasma–liquid interface,¹⁵ and second, the $eaq−$ concentration is obviously dependent on system variables such as the electron flux. However, as shown in Fig. 2(b), before the interface is depleted of scavengers, the reaction yield is close to 100%, indicating that nearly all injected electrons react with the scavenger. Because scavenger reduction [Eq. (1)] dominates so strongly over second-order recombination [Eq. (3)] during this initial time period, we can treat the $eaq−$ concentration as a constant. This simplifies the scavenger reduction reaction rate to a (pseudo)first-order expression for [S]. Similarly, because the yield is so close to 100%, the reaction rate for the scavenger reduction can be equated to the entirety of the molar electron flux (J/F), as shown in Fig. 2(a), reacting over the characteristic length, z_c. With this assumption, $ks[eaq−]$ can be approximated as

Using these approximations, it is possible to define a characteristic time and characteristic length that render the reaction–diffusion equation independent of most system variables, including the scavenger concentration, scavenger diffusivity, and the molar electron flux. The resulting characteristic time is

and the resulting characteristic length is

A detailed derivation of these characteristic variables is shown in Sec. S.3 of the supplementary material.

In this system, the characteristic time is the time it takes to transition from a kinetics-limited system to a transport-limited system. Physically, it represents the time at which the scavenger concentration at the boundary is depleted and needs to be transported from the bulk to continue reacting with injected electrons (see Sec. S.2 of the supplementary material), and as such, physically represents the same thing as t₉₉. Similarly, at this time, the characteristic length represents the distance between the interface and where the scavenger approximates its bulk concentration (see Sec. S.3 of the supplementary material for more details).

We compared this analytical expression for the characteristic time t_c with the depletion time t₉₉ we extracted from simulations (e.g., green circles in Fig. 2). These are shown in Fig. 3 for a number of simulated conditions of S_B, J/F, and D_s as listed in Table S2 in Sec. S.4 of the supplementary material. All simulated measurements of t₉₉ collapse along a single line, suggesting that there is a linear relation between t₉₉ and the characteristic time, regardless of the scavenger bulk concentration, scavenger diffusivity, or molar electron flux. Note, however, that the slope of the line (illustrated with the solid line) is smaller than unity (illustrated with the dashed line). We term this the measured scale factor, and it is t₉₉/t_c = 0.66 for this data. The deviation from unity is due, in part, to the mathematical simplifications used to determine the analytical characteristic time, including Eq. (11), that are not present in the fully coupled, non-linear simulations that yield t₉₉. Differences can also arise because t₉₉ is arbitrarily defined as a drop in the yield from 100% to 99% and thus approximates but is not the exact time when the system transitions from kinetics-limited to transport-limited. Figure 3 does, however, confirm that t₉₉ is strongly associated with the transition from a kinetics-limited to diffusion-limited reaction at the interface, as well as with the depletion of the scavenger at the interface.

One of the assumptions made in obtaining the analytical characteristic time in Eq. (12) is that the scavenger reduction reaction dominates over second-order recombination at short times before scavengers are depleted from the interface. This raises the question of whether t_c has any validity when the rates of both reactions become comparable, either by having a much lower scavenger concentration, a smaller scavenger reaction rate constant, or a significantly higher molar electron flux.

To understand this, we tested how the ratio of t₉₉ to t_c behaved as a function of the ratio of the rates of both reactions, $rS=ks[eaq−][S]$ and $re=ke[eaq−]2$⁠. However, as mentioned earlier, a non-uniform solvated electron concentration profile greatly complicates any mathematical calculation. Here, we used an analytical approximation for $[eaq−]$ derived in our previous work,³⁵ 

which corresponds to the $eaq−$ interfacial concentration in the presence of scavengers, where the scavenger reaction dominates and there is no scavenger depletion. This expression assumes that $ks[S]≫ke[eaq−]$⁠, but it should be possible to use it as a way to understand how the depletion time might deviate from the scaling presented in Sec. III B as the ratio between reduction and recombination decreases, as long as the inequality holds true. This approximation in Eq. (14) also proved to be very similar to that generated by simulations before the scavenger was depleted (see Sec. S.2 of the supplementary material). If the local scavenger concentration is additionally approximated as its bulk concentration $[S]=SB$⁠, the ratio $rS/re=ks[S]/ke[eaq−]$ can be approximated as

where every parameter is known (k_s, D_e, k_e) or is a system variable (S_B, J/F).

A plot of t₉₉/t_c as a function of $ks[S]/ke[eaq−]$ in Eq. (15) is shown in Fig. 4. As expected, at high values of $ks[S]/ke[eaq−]$⁠, the ratio of t₉₉/t_c approaches the measured value of t₉₉/t_c = 0.66 extracted from Fig. 3, although it does saturate at slightly higher values, which we attribute to uncertainty in numerically estimating t₉₉ for some of our longer simulations. However, at low values of $ks[S]/ke[eaq−]$⁠, the ratio of t₉₉/t_c decreases to the point where the characteristic time described in Sec. III B cannot be used to estimate the time of depletion. Interestingly, the decrease in (t₉₉/t_c) occurs in an orderly fashion, and regardless of the scavenger concentration, these ratios fall along the same curve. This behavior shows that the ratio between the two reactions, the scavenger reduction, and the second-order recombination, as estimated by Eq. (15), is the only thing required to know if the scaling from Sec. III B is valid.

From this analysis, several conclusions can be drawn. First, the characteristic time, t_c, can reliably be used to approximate the depletion time as long as the ratio of scavenger reduction to second-order recombination is greater than 5000, with reasonable order of magnitude estimates at ratios as low as 1000. Second, the characteristic time describes the maximum time it takes to deplete the interface, as t₉₉/t_c only decreases at lower ratios of $ks[S]/ke[eaq−]$⁠. And third, it shows that at very low ratios of $ks[S]/ke[eaq−]$⁠, there is a point where it would be impossible to have high yield (>99%) even at extremely short times. This indicates that electrons will be lost to second-order recombination near instantaneously from the start of reactions.

On a similar note, a recent publication by Zheng et al.¹⁶ correctly pointed out that Eq. (14), which is from Ref. 35, may not be valid at low scavenger concentrations. Using 1D simulations of 1 mM AgNO₃ as a model system, they show that the interfacial $eaq−$ concentration for molar electron fluxes between 0.01 and 0.4 mol m⁻² s⁻¹ is the same as without scavengers after 10 μs. It goes on to suggest that Ag⁺ and NO₃⁻ concentrations were not high enough to sufficiently scavenge $eaq−$ and that is why the interfacial $eaq−$ concentration followed the analytic model for a solution with no scavengers. Here, we posit that the question of whether the scavenger concentration is too low or cannot be described with the analytic model for scavengers can be answered with the scaling shown in Figs. 3 and 4. For example, even using their lowest reported molar flux, the ratio of the rates [Eq. (15)] would still be below 100, suggesting that the interface would in most cases be depleted from the start, i.e., at fractions of the characteristic time. Given that the characteristic time would fall between 0.04 and 70 μs, depending on the molar electron flux, it is no surprise that the interface is depleted after 10 μs, in which case the interfacial $eaq−$ concentration would be equal to that in neat water (i.e., without scavengers). However, at higher scavenger concentrations or a lower molar electron flux, the $eaq−$ concentration described by Eq. (14) should be valid before depletion, as shown in Fig. S2 of the supplementary material.

The analyses in Secs. III A–III C are based on the solution to the 1D reaction–diffusion equations, which reflect only a specific point at the plasma–liquid interface. However, this analysis can provide insight into how a physical 3D system behaves, and we do so here by analyzing the implications when the plasma takes a Gaussian profile at the liquid interface. We chose this profile for several reasons, including it being the shape of the plasma–liquid system on which we have conducted experiments [see Fig. 5(a)]^29,36 and because plasma species like •OH and atomic oxygen have been reported to have a Gaussian or a Gaussian-like profile.^13,37,38 Additionally, Gaussian profiles occur naturally in nature when there is diffusion associated with a localized source of a species,³⁹ which could occur when a species such as •OH is generated in the liquid phase.^8,40

We consider a system that is axisymmetric about the center of the plasma (r = 0 in Fig. 1), and we seek solutions for $[eaq−(z,r,t)]$ and $[S(z,r,t)]$⁠. Under the assumption of a Gaussian current of electrons to the liquid, the molar electron flux (J/F) at the boundary (z = 0) can be described as a function of the radial distance (r),

where σ_p is the Gaussian radius of the plasma and $J0$ is the peak current density at r = 0. In 1D, the solutions of Eqs. (4) and (5) correspond to a line stretching from a single point at the plasma–liquid surface toward the bottom of the liquid (see Fig. 1). Here, we assume that radial diffusion of both the solvated electrons and scavengers is negligible compared to the dominant diffusion into the depth of the solution. Thus, rather than solving axisymmetric expressions of Eqs. (4) and (5) in cylindrical coordinates, we assume that $[eaq−(z,r,t)]$ and $[S(z,r,t)]$ can be effectively approximated by the superposition of 1D solutions for $[eaq−(z,t)]$ and $[S(z,t)]$ at every point r, where the boundary condition (J/F) varies as described by Eq. (16). We justify this assumption because radial concentration differences of many plasma–liquid systems, including that of our previous work, occur over a span of millimeters,^13,25,37,38 whereas similar concentration differences occur over a span of micrometers or even nanometers into the solution. (This can be observed in the [S] profiles in Sec. S.2 of the supplementary material.)

As we have shown that the yield is approximately 100% at very early times [Fig. 2(b)], the radial reaction rate profile $R′′(r,t)$ should take the Gaussian shape of the molar electron flux [Eq. (16)] at the beginning of plasma–liquid contact. As time progress, the profile should change shape as parts of the interface become depleted. Figure 5(b) shows the Gaussian molar electron flux (solid line) as well as the reaction rate for a bulk scavenger concentration of S_B = 100 mM at times of t = 0.001, 0.003, 0.01, and 0.1 s. It is clear that the reaction rate at the center of the plasma (r = 0) decreases rapidly in time. Interestingly, at the tails of the Gaussian profile, where the flux of electrons is lower, the interface depletes much more slowly such that in time the edges of the plasma yield the highest reaction rate and the center of the plasma the lowest reaction rate. Put another way, the curvature of the reaction rate profile changes from concave (t = 0.001 s) to convex (t = 0.003 s) due to depletion at the center.

As Fig. 5(b) shows, the normalized radius $rσp$ at which the yield is still 100% changes with time and by the very nature of the distribution, does so in a non-linear fashion. This position along the plasma profile, which we call the radius of depletion $rσpd$⁠, can be estimated from Eqs. (12) and (16). First, we consider Eq. (12) and the equivalent molar electron flux that occurs at the characteristic time when the interface is depleted of scavengers and the reaction transitions from kinetics-limited to diffusion-limited,

If we plug this molar electron flux into Eq. (16), we can rearrange and determine an expression for the normalized radius where depletion occurs as a function of time as

Figure 6 shows how $rσpd$evolves with time for different scavenger concentrations as predicted by Eq. (18) (lines) and as extracted from simulations (points). Initially, at very short times, the radius of repletion is near the centerline, where the molar electron flux is the highest and scavengers are depleted quickly. As time progresses, the radius of depletion moves outward and asymptotically plateaus, where the plateauing radius depends on the concentration of scavengers in the bulk (for higher concentrations, the plateauing radius is greater). Notably, the analytical model in Eq. (18) underpredicts the simulated data, and we attribute this to the simplifications inherent to the model for the characteristic time [Eq. (12)] discussed in Secs. III B and III C, where non-linear effects are neglected. An important takeaway from this behavior is that depletion occurs rapidly where the electron flux is the greatest, and the highest yield will often occur at points along the plasma–liquid interface distant from the greatest electron flux. Such insight is important for designing reaction systems to achieve high yield and potentially spatial uniformity. To that end, system design to prevent depletion from occuring at these short times may be needed.

The mathematical analyses of Secs. III A–III D, based on the reaction–diffusion equation, are not exclusive to $eaq−$ chemistry but rather applied to any plasma chemical species with extremely fast reaction rates in the liquid phase that can also be consumed by second-order recombination, such as •OH, which recombines to form hydrogen peroxide (H₂O₂). Here, we consider the implications of our analyses in the design of reaction systems for applications that rely on chemistry driven by solvated electrons and hydroxyl radicals.

The primary finding of our analyses is that the characteristic time where the reactions transition from kinetic- to transport-limited depends on the diffusivity and bulk concentration of scavengers and the molar electron flux at the interface. Furthermore, for non-uniform molar electron fluxes, this behavior can result in the position where the reaction rate and yield are highest changing over time. Based on this, there are a few possibilities to improve the yield of a transport-limited reaction: increasing the transport driving force by complementing diffusion with advection; modifying the reaction that causes depletion to occur; or changing the reaction system entirely, such as using a multiphase system.

While the first option seems to be a straightforward and has oft been studied in electrochemical systems,⁴¹ a simple analysis suggest that this obvious approach might be incapable of overcoming the fundamental limitations of a plasma–liquid system. Consider the case of a continually flowing channel just beneath the plasma–liquid interface as illustrated in Fig. 7(a); in principle, the continuous flow transports reacted scavengers away from the interface and refreshes scavengers to the interface. In order to overcome the rapid depletion that occurs, the flow speed would need to ensure that the solution traverses the diameter of the plasma (on the order of ∼1–10 mm) in a time on the order of t_c. The minimum required flow speed would be

where d is the diameter of the plasma. While this velocity depends strongly on the bulk scavenger concentration and the flux of the plasma species, simple substitution can be used to get an idea of how fast of a flow would be required. For example, for solvated electrons introduced through a DC plasma system with a liquid anode, the flux at the center of the plasma is about 0.01–0.1 mol m⁻² s⁻¹ based on reported currents,^1,25,42 which means that for an S_B of 1 mM, a minimum flow velocity of 10²–10⁴ m/s would be required to replenish the interface, which is unrealistically fast. Overall, the required liquid flows would be realistic only for reactions with a very high scavenger concentration, in which case depletion might not be an issue or a simpler technical solution could be implemented. This is supported by our earlier result that showed that stirring the solution, which creates a similar flow at the plasma–liquid surface, did not significantly increase the reduction efficiency of chloroacetate by $eaq−$⁠.²⁵ 

Alternatively, the liquid could flow in the direction of plasma impingement, as illustrated in Fig. 7(b), where advection normal to the interface would directly contribute to transporting scavengers to the depleted region. In fact, systems that incorporate advection normal to the interface have been studied,⁴³ primarily for analytical chemistry applications. It is not difficult to envision a similar system being built with a focus on optimizing an interfacial chemical reaction. In this case, the flow speed should be sufficient to replenish the depletion distance, as estimated by z_c in Eq. (13), over the characteristic depletion time t_c. That is,

or

A bigger problem with using advection to overcome depletion, however, may come in the form of a viscous boundary layer in the solution at its interface with the plasma. Given that any $eaq−$ or •OH concentration profile will form mostly in the first 100 nm of the liquid, the benefits of advection will be hindered if the viscous boundary layer in the liquid is significantly larger than the penetration depth of these species. To know if the boundary layer will be similar to that of the reaction zone near the interface, we can approximate the boundary layer thickness of a liquid flow toward a plasma interface as that of the spread of a liquid jet over a horizontal plane, δ, where $δ=Rνu1/2$⁠,^44,45 with R being the radius of the jet, ν being the kinematic viscosity of the liquid, and u being the velocity of the liquid jet. Using water as a solvent, a jet radius of 1 cm, and a flow velocity of 1 cm/s, the boundary layer thickness would be around 1 mm, which is many orders of magnitude larger than the penetration of plasma-generated species, indicating that the boundary layer would pose a problem when trying to use advection to increase the transport of the scavenger from the bulk of the solution to the interface. However, this analysis assumes that the boundary layer in the liquid phase of water flowing toward the plasma would behave as with a solid boundary. While this will not be the case, particularly, because of lower friction between water and the gaseous plasma, it points to a potential issue in the design of plasma–liquid systems trying to overcome transport limitations via advection.

An alternative to trying to improve transport is to modify the way the reactions are induced by the plasma. For example, affecting the reaction at the interface in the form of pulsed voltages has been shown to increase efficiency with respect to a non-pulsed (DC) reactor.²³ While the chemical species generated are the same, pulsed operation allows for the species in the bulk to diffuse back to the interface in the time between each pulse, which is generally on the order of 0.01–0.1 s.^22,46 In principle, depletion could be avoided by allowing 2t_c between each pulse (the second half for allowing replenishing of the interface, $zc2Ds$⁠), corresponding to a pulse frequency faster than 1/2t_c. The actual time, however, might be lower depending on the actual system as the role of other species formed in between pulses could be important, including the potential products formed by the $eaq−$-scavenger reaction.

In a similar manner, decreasing the flux of the plasma species can increase the efficiency greatly as the time spent under depletion will be lower, particularly, if coupled with a pulsed system. In the case of solvated electron chemistry, this would imply reducing the current density, either by reducing the current or increasing the area of plasma contact. Ultimately, the characteristic time scales as the square of the current density [Eq. (12)], which indicates this would be a very effective lever in system design.

Finally, multiphase systems could also be employed to overcome transport limitations. For example, studies have shown that flowing bubbles or even foam through a reactor greatly increases the efficiency at which dyes are oxidized.^21,22 Similarly, spraying droplets greatly increases reaction efficiency when compared to a reference reactor.²³ While the simple explanation is that bubbles or droplets increase the surface area to volume ratio of plasma–liquid contact, thereby distributing the flux across a wider boundary closer to the bulk of the solution, it should not be discounted that bubbles rupture, often at very short times,⁴⁷ and both bubbles and droplets may continuously deform, changing the boundary with the plasma. Given that reactions by $eaq−$ or •OH are very fast, a high surface area to volume ratio would ensure that, even if the scavenger close to the interface depletes, a greater surface area would have been treated. It is also true that the smaller the droplet or bubble, the volume of the gas or liquid where the reaction occurs will be greater as a proportion of the total volume, leaving less of the solution untreated. Finally, it is possible that when bubbles or droplets deform, especially when they split into smaller bubbles or droplets, the time to depletion is reset, making the characteristic time from Eq. (12) an effective marker, but one that is continuously dialed back by the nature of that plasma–liquid system.

In this work, we analyzed a transient, 1D, interfacial reaction at the plasma–liquid interface between solvated electrons and their chemical scavengers. With the scaling of the scavenger reaction–diffusion equation, a characteristic time closely associated with the transition to a transport-limited system was derived, which is proportional to the scavenger diffusivity and to the square of its concentration and inversely proportional to the square of the molar electron flux. This characteristic time is shown to be very similar to the time it takes to deplete the interface, particularly when there is a high ratio between the scavenger reduction rate and electron second-order recombination but becomes much greater as this ratio decreases, making it an effective estimation of the maximum time it would take to deplete the interface.

Additionally, this analysis is extended to 3D to analyze how depletion would occur in a physical system where the plasma current takes the form of the Gaussian profile leading to a Gaussian molar electron flux across the plasma–liquid interface. The analysis shows that the scavenger at the center of the plasma (the point with the highest flux) depletes quickly, which leads to a reaction rate lower than the surface molar electron flux, while the reaction rate at the edge of the plasma continues to follow the original Gaussian profile. From this analysis, the radial point of transition from Gaussian to a depleted interface can be estimated as a function of time. Furthermore, given these results are based on scaling of kinetics and diffusion at the interface, their conclusions should be valid for other highly reactive plasma-generated species, including H• and •OH.

Finally, this analysis has important implications for reactor design. For the most part, we restate some of the conclusions from the plasma wastewater treatment literature while arguing from an interface depletion rather than energy efficiency perspective. Namely, pulsed voltage plasma systems can be effectively used to increase the yield, and consequently, the efficiency of the reaction, and multiphase systems, particularly, bubbles and foam in the liquid phase can increase efficiency by increasing the surface area to volume ratio and constantly modifying the boundary. Flowing the solution through the area of plasma–liquid contact is also analyzed, giving further motivation and quantification of its limitations, including high velocity flows and the persistence of a boundary layer at the interface.

See the supplementary material for additional information on the simulation characteristics and benchmark with other models, concentration profiles, analytical derivations, and data points presented in Figs. 3 and 4.

This work was supported by the U.S. Army Research Office under Award No. W911NF-17-1-0119. D.M.B. is supported by the U.S. Department of Energy Office of Science, Office of Basic Energy Sciences under Award No. DE-FC02-04ER1553.

The data that support the findings of this study are available from the corresponding author upon reasonable request.