Fluid model of plasma–liquid interaction: The effect of interfacial boundary conditions and Henry’s law constants

The interaction between cold atmospheric plasmas and aqueous solutions has attracted increasing attention because it plays a crucial role in many promising applications such as in biomedicine,^1–3 agriculture,^4,5 and water purification.^6,7 The aqueous reactive species, directly or indirectly produced by gas plasmas, are deemed to serve as the main cocktail for those applications,^8,9 so a detailed understanding of the plasma–liquid interaction is very important for not only gaining insight into the production mechanism of aqueous reactive species^10,11 but also finding the approach to optimize the production of aqueous reactive species for applications.¹² The physical and chemical processes occurring at the plasma–liquid interface play a dominant role in the production of aqueous reactive species. However, so far, the experimental study on the plasma–liquid interface has rarely been reported. This is because, on the one hand, the thickness of the interface is only in nano- and microscales as estimated in different literature studies,^13–16 and it is too thin to be measured by most experimental methods; on the other hand, several reactive species have their mass transfer and heterogeneous reactions strongly coupled at the interface,¹⁷ and it is too complex for experiments to clarify their relationships. Therefore, as pointed out by The 2012 Plasma Roadmap, understanding the physical and chemical processes at the plasma–liquid interface is one of the major challenges in the research field of plasma–liquid interaction,¹⁸ which strongly hinders the development of plasma devices for water-relevant applications.

Given that great difficulties lie in the experimental studies, numerical models have been developed as an alternative approach for the study of plasma–liquid interaction. However, many models are zero-dimensional, and hence, the plasma–liquid interface is neglected,¹⁹ leading to a large inaccuracy of the simulation results. In recent years, one-dimensional (1D) and two-dimensional (2D) fluid models have been developed in which the plasma–liquid interface was considered,^16,20–22 but since much is unknown of the nature of the interface, the interfacial boundary conditions and their relevant coefficients needed for the simulation have a large uncertainty. For example, the thermodynamic, diffusive, and kinetic boundary conditions are used to describe the mass transfer of neutral species in the models, among which the thermodynamic boundary condition assumes a consistent concentration equilibrium and a flux continuity between gas and liquid phases,^10,20,22,23 while the diffusive and kinetic boundary conditions set the diffusivity and kinetic flux of each species as a function of its gaseous and aqueous concentrations, respectively.^11,16,19 All three boundary conditions have reasonability in theory, but the equilibrium states obtained by the models are different, leading to an uncertainty in the simulation results. Moreover, the Henry’s law constant, as an index of solubility and a key parameter of the boundary conditions, has not been measured for all the neutral species. Some of the Henry’s law constants are estimated theoretically with an uncertainty of 1–2 orders of magnitude,²⁴ and even some others are roughly speculated from those of other species with similar molecular structures.²⁵ Different values of Henry’s law constants used in the model would also lead to an uncertainty in the simulation results.

How much uncertainty in simulation results could be present with respect to different interfacial boundary conditions and Henry’s law constants used in the model? If it is big, then the simulation should be incapable of describing the plasma–liquid interaction unless a proper boundary condition and accurate Henry’s law constants are determined beforehand. This would certainly be a bad news for the fluid model simulation of plasma–liquid interaction, which has just arisen in the last few years. If it is small, then the simulation results are acceptable regardless of the different interfacial boundary conditions and Henry’s law constants. This would suggest that the simulation of plasma–liquid interaction can break through the lack of basic data to some extent and the doubts about the validity of simulation would be reduced. The uncertainty is possible to be small; for example, Lindsay et al. found that the uncertainty in electron dissolution is small as a function of the electron loss coefficient, which is unknown at the plasma–liquid interface. They altered the electron loss coefficient from 10⁻⁴ to 1 in a fluid model, and interestingly, the electron density hardly changed in the plasma-activated water (PAW).²⁶ 

In this paper, three 1D fluid models of plasma–liquid interaction are developed to investigate the uncertainty in simulation results with respect to different interfacial boundary conditions and Henry’s law constants. The neutral species are focused instead of the electrons as done by Graves et al., including diverse kinds of reactive oxygen species (ROS) and reactive nitrogen species (RNS) that are strongly relevant to the applications.^27,28 The plasma is a plasma jet with a feeding gas of He + 0.5 vol. % air, and the gas doping from the ambient air is also considered. This gas mixture is commonly used in plasma biomedicine.^29–31 The liquid is deionized water that is in contact with the plasma plume, and the plasma jet and deionized water are fully coupled in the model for calculation. The three fluid models are the same except for their interfacial boundary conditions of neutral species, i.e., the thermodynamic, diffusive, and kinetic boundary conditions are used in models I, II, and III, respectively. Therefore, a comparative study among the models could quantify the uncertainty in simulation results with respect to the interfacial boundary conditions. Moreover, the Henry’s law constants are altered according to their uncertainty reported in the literature. Based on this, three sets of Henry’s law constants are used in model I for a comparative study, and the uncertainty in simulation results with respect to the Henry’s law constants could be quantified.

This paper is organized as follows: the description of the fluid models is given in Sec. II, in which the interfacial boundary conditions and Henry’s law constants are provided and discussed. The density distributions in gas and liquid phases are given in Sec. III with respect to the different interfacial boundary conditions and Henry’s law constants, and discussions are made on the uncertainty in the simulation results and the production mechanism of the aqueous reactive species. Finally, a brief conclusion is given in Sec. IV.

The numerical model is developed based on the schematic diagram of the plasma–liquid interaction shown in Fig. 1(a). It contains a hollow high-voltage electrode, a Petri dish of deionized water, and a ground electrode under the Petri dish. The gas gap between the high-voltage electrode and water is 5 mm, the depth of water is 2 mm, and the thickness of the Petri dish is 0.5 mm. The model is one dimensional (1D) with similar physical dimensions as in the experimental setup, shown in Fig. 1(b). A sinusoidal voltage with a peak-to-peak value of 2 kV and a frequency of 13.56 MHz is applied on the high-voltage electrode, while the bottom of the dielectric layer is grounded. A similar radio-frequency plasma jet was reported previously.³² The working gas is helium, and the gas flow rate is 1 SLM (standard liters per minute). In addition, 0.5 vol. % air impurity is implemented in order to enhance the chemical reactivity of the plasma.^29–31 The influence of water evaporation on the background gas composition is also considered, but the reduction in water volume due to evaporation is neglected.

In this work, a fixed gas–liquid phase interface is considered and the phase change is distinguished numerically by computational domains and meshing. Overall, it solves the mass conservation equations [Eq. (1)] with the drift-diffusion approximation for both gaseous and aqueous species, the electron energy conservation equation [Eq. (2)] in the gas phase, and Poisson’s equation [Eq. (3)] in both phases as follows:

where the subscript e represents the electrons, k represents the working gas (He, N₂, O₂, and H₂O), and i and j represent the ith species and jth reactions, respectively. The symbols n, $Γ$⁠, S, m, μ, and D are the density, flux, gain/loss rate, mass, drift coefficient, and diffusion coefficient, respectively. E is the electric field, ΔE is the electron energy loss of inelastic collision, R is the reaction rate, and R_el is the momentum transfer collisional rate between the electrons and the working gas. ε is the mean electron energy, T is temperature, k_B is the Boltzmann constant, and q_i is the charge of species i. For more details on these governing equations, refer to our previous publications.^33–35 

The gas composition would significantly influence the discharge characteristics and the produced reactive species. Thus, an incompressible Navier–Stokes equation [Eq. (4)] and the convection–diffusion equation [Eq. (5)] are first solved in order to obtain the distributions of the working gas composition,

where u represents the velocity vector, ρ represents the overall mass density, p represents static pressure, μ represents the dynamic viscosity, k represents the working gas (He, N₂, O₂, and H₂O), and D represents the corresponding diffusion coefficient.

The 2D cylindrically symmetric distributions of the working gas consisting of H₂O, O₂, N₂, and He are shown in Figs. 2(a)–2(d), and their corresponding distributions along the axial direction are extracted and shown in Figs. 2(e) and 2(f), respectively. They are set to be the working gas composition in the following simulations and are regarded to be unchanged. These distributions included the influence of the water solution on air humidity by means of setting the saturated water vapor concentration at 300 K on the interface as a constant via the Antoine equation.^23,36

The chemical reactions have been carefully chosen in this model. Noting that the chemical processes are very sensitive to the admixture of reactive gas from either feeding gas or ambient air, it would be rather computationally unaffordable if the complete chemistry set, usually hundreds of species and thousands of reactions, is incorporated. Therefore, the key pathways with altered composition ratios, representing different locations in the gas phase, are calculated and extracted via a global model as reported previously.³⁷ The gaseous chemistry is taken to include all the key species and reactions at the considered working gas composition along the axial direction as obtained by Eqs. (4) and (5). This simplified approach for chemistry has already been applied for He + H₂, He + O₂, and He + air models in our previous works, and for details, refer to Ref. 37. As for the liquid phase chemistry, most of them are the same as in Ref. 20 and reactions of electrons and some short-lived species are appended according to Refs. 11, 12, and 19 in the meanwhile. Overall, 51 species and 223 reactions in the gas phase and 47 species and 148 reactions in the liquid phase are incorporated in the following models. The species in these models are listed in Table I. The gaseous and aqueous reactions calculated in the models are listed in Tables IV and V in the  Appendix, respectively, and their references are given in the tables. It is worth noting that, currently, the quantitative measurements for reactive species are limited, and it is very hard to achieve their temporally resolved distributions, especially in a very miniature targeted area as the interfacial layer. However, the characteristics of gaseous helium plasmas with air mixtures have been qualitatively validated to some extent in our previously reported work.³³ 

For the fluid models of plasma–liquid interaction, the gas and liquid phases are usually separated into different domains artificially with a fixed phase interface, and the mass transfer processes, such as solvation and absorption, are described by reasonable boundary conditions. Three kinds of commonly used interfacial boundary conditions for neutral reactive species are listed in Table II,^{10,11,16,20–23} where n_, Γ, H, D, and vth are the density, inward flux, Henry’s law constant, diffusion coefficient, and thermal velocity, respectively. The subscript interface represents the location of gas–liquid interface and gas and aq represent gas and liquid phases, respectively.

All the settings are the same in the three models except for their interfacial boundary conditions. Specifically, there are three kinds of commonly used interfacial boundary conditions, which are adopted in models I, II, and III, respectively. (1) Thermodynamic condition (model I) (Table II):^10,20,22,23 at the mesh point of the interface, the aqueous concentration is always in equilibrium with that of the gas phase; meanwhile, the flux continuity between them is applied to ensure the mass conservation. (2) Diffusive condition (model II) (Table II):^11,19 the diffusion coefficient at the interface varies with the density relation between the gas and liquid phases in the form of (Hn_gas − n_aq)/(Hn_gas). When the solvation balance is achieved at the interface, the diffusion coefficient decreases to zero and the mass transfer from the gas to liquid phase consequently terminates. (3) Kinetic condition (model III) (Table II):¹⁶ the inward flux of species into the solution is scaled to the kinetic flux of the species in the gas phase, 0.25nvth, with the same factor of (Hn_gas − n_aq)/(Hn_gas) as in model II. All three interfacial conditions calculate the overall relationship between the gaseous and aqueous concentration for neutral species, and the mass transfer would eventually end up at the solvation equilibrium at the interface. However, the rates required to achieve the equilibrium state given by these three interfacial conditions are different. It is evident that in model I, the phase equilibrium is assumed consistently, but in models II and III, the equilibrium is gradually obtained and the mass transfer of reactive species would be cut off once the equilibrium is achieved at the interface. However, the fluxes across the two phases are different and the equilibrium is supposed to be achieved at different moments. All three interfacial boundary conditions are reasonable, but there might be distinct different relationships between the gas and liquid phase concentration at the interface. Considering that for a system of plasma and liquid, multiple factors dynamically contribute to the mass transfer of reactive species across two phases¹⁷ and in the scope of the fluid model, the dissolution and reactions are carried out simultaneously. Therefore, it is quite uncertain whether the competition between mass transfer and reactions, or any other phenomenon, would shift the spatial–temporal character of reactive species.

For charged species, the descriptions of mass transfer across the interface are the same in three conditions. Briefly, the positive ions are supposed to dissolve immediately when they impact the interface because their potential energy is larger than any activation energy barrier required to cross over.³⁸ That of negative ions is rather less clear, and they are schematically regarded to freely transport from gas to liquid because they are far more polarizable than either neutrals or cations. For electrons, the measurement of accommodation at the interface is inaccessible, but numerical work finds out that the reflection or absorption would not much affect their distribution in the liquid phase.²⁶ Thus, electrons are set to pass through as soon as they reach the interface.

Henry’s law characterizes the equilibrium between the concentration of solute in an ideal dilute solution and its corresponding gaseous partial pressure, where the larger the dimensionless Henry’s law constant is, the stronger the solubility of a given gas would be. For some species such as H₂O₂ and O₃, their Henry’s law constants could be measured by experiments. Others such as OH and HO₂, which would transform into conjugate ions or other species in solution, could be theoretically estimated using their Gibbs free energy, reduction potential, or solvation energy, usually with an uncertainty of ∼10 times.²⁴ There are also species such as atomic O, N, and H that possess high chemical reactivity but are not very common in atmospheric chemistry nor chemical industries. However, they might be precursors in plasma chemistry and initially produced through collisions between electrons and surrounding air during the discharge and then transformed into other reactive species. Their lifetimes are relatively short, and the chemical processes are complex. Hence, it is hard to determine their Henry’s law constants by experiments or potential theory. However, they could be quite prevalent in plasma. For instance, O is found to be the most abundant ROS in He + air³³ and He + O₂ plasma⁴⁹ and N is found to be the key species that connects the chemistry of ROS and RNS.³³ Likewise, H is also an abundant neutral species in helium plasma with water evaporation.⁵⁰ They are supposed to be important in the gas phase during the treatment of plasma to liquid, but their Henry’s law constants for modeling are not adequate up to our knowledge. That of atomic H is obtained by interpolation and translation according to the linear relationship between the atomic radius and the solvation energy of noble gases, which have the same geometry of a single atom.^25,39 Evaluation of that of O and N is even rare. In some works, they are set to be the same as that of O₂ or N₂ because they are composed of identical elements.²⁰ Others set the constants to be 1 owing to their high reactivity in water and consequently the inaccessibility of saturation.¹⁹ It is clear that Henry’s law constants are essential parameters in the modeling of plasma–liquid interaction, but their values would be different by orders of magnitude in different numerical works, leading to different simulation results. Evaluating the uncertainty in simulation results caused by the use of different Henry’s law constants is of importance. For this reason, a comparative study is conducted and three sets of Henry’s law constants are applied in model I, which represent the minimum, medium, and maximum values according to their corresponding uncertainties, respectively. The aqueous distributions of reactive species are compared, and the effect of Henry’s law constants on the results is concluded. The Henry’s law constants are listed in Table III. It is worth mentioning that due to the lack of reference data of that of O and N, they are estimated with the same method as that for H.^25,39

Concerning the short period of radio frequency (RF) discharge, an acceleration approach is proposed to attain the desired discharge time, which is similar to those given in Refs. 51–53 in principle, and the calculation is addressed in the following manner as illustrated in Fig. 3. It is designed based on the character that radio frequency discharge usually works on the continuous discharge mode and the concentration of short-lived species would not change much during one RF cycle. After obtaining the working gas distributions, the plasma dynamics is calculated with the mass-conservation equations for all species as well as electron energy conservation and Poisson’s equation. When the phase-averaged density of any gaseous reactive species changes no more than 1% between two adjacent voltage cycles, the simulation is regarded to be in the electric steady state. Then, the calculation of short-lived species, whose density shows an apparent periodic pattern, is turned-off and only that of long-lived species⁵¹ as well as Poisson’s equation in liquid is calculated. The distributions of gaseous short-lived species and electron energy are set as their phase-averaged spatial distributions because the time step would be several times larger than the discharge cycle at this stage. When the concentration change in long-lived species increases by more than 50%, the model would again calculate the full plasma dynamics for five RF cycles. Since the lifetime of short-lived species is usually no more than one RF cycle, it would be enough to achieve a new steady state. The circulation would continue until it reaches the desired time point. This approach enables a large number of RF discharge cycles but to some extent neglects the electric character within one cycle. However, the plasma dynamics is adjusted at intervals to ensure the calculation results in a reasonable error range.

All equations described above were solved by using the mathematics module of a time-dependent finite-element partial differential equation solver, COMSOL Multiphysics®. The total number of mesh elements is 1650 with the minimal value of 0.1 nm at the gas–liquid interface and 10 nm at other boundary nodes.

A comparative study is conducted in this work to unravel the dependence of the profiles of neutral species on the interfacial boundary conditions and the uncertainty in Henry’s law constants. For interfacial conditions, three cases are calculated, and each one has a unique interfacial boundary condition for neutral species, i.e., the thermodynamic, diffusion, or kinetic boundary conditions. For the uncertainties in Henry’s law constants also, three cases are calculated, and each one employs one set of Henry’s law constants representing the maximum, medium, or minimum value within their uncertainties. No other parameters are changed in the calculations. The gaseous distributions of reactive species are given. Then, the distributions of dominant neutral species nearby the gas–liquid interface are calculated by models I, II, and III (given in Table II) with the Henry’s law constants of set 2 (given in Table III), and those calculated by Henry’s law constant sets 1, 2, and 3 in model I are illustrated, and the corresponding effects are demonstrated. Finally, the aqueous chemical profiles are presented.

The gaseous distributions are shown in Fig. 4. The spatial–temporal distribution of electron density and electron temperature are plotted in Figs. 4(a) and 4(b), respectively, where the vertical axis represents the gas gap from the interface to the nozzle and the horizontal axis represents the time within one RF cycle. Figures 4(c) and 4(d) show the phase-averaged distribution of ROS and RNS, respectively, where the vertical axis represents the gas gap and the horizontal axis represents the average density.

Since the boundary conditions for charged species are the same in models I, II, and III, and their distributions are found to be very similar, only the results of model I (Table II) are presented. From Fig. 4(a), the spatial–temporal distribution of electron density reveals a similar periodically oscillation pattern as it is in gas phase RF discharge,³⁵ whereas it is not symmetrical between the upper and lower sides of the gas gap. This is because the water evaporation downstream traps electrons there. In Fig. 4(b), the maximum electron temperature appears at the upper side of the nozzle. According to the bipolar characteristic of RF discharge, there is also a rising of electron energy near the interface but much lower than that near the nozzle side. When helium flows out of the tube, the surrounding ambient air would mix into the working gas, especially when it comes near to the interface. The inelastic collisions are enhanced and especially more electrons would collide with increasing water evaporation to produce atomic H species,⁵⁴ which consumes much energy.³⁷ 

After a 10 ms discharge, the phase-averaged density of these species changes no more than 0.02% between two adjacent voltage cycles; therefore, they are regarded to be in steady states [Figs. 4(c) and 4(d)]. As for the gas compositions, this result shows the same tendency as that in the He + air parallel electrode discharge in a global model³⁷ and fluid model³³ where O and O₂(a) are the dominant ROS and NO, NO₂, and HNO₂ are the dominant RNS. As for the spatial distributions, the peak value of H-contained species appears near the gas–liquid interface,⁵⁵ while that of H-free species appears near the nozzle outlet. This difference is reasonable considering the spatial character of water evaporation and electron energy.

Regarding the distributions of reactive species calculated by three different interfacial boundary conditions (results of models I, II, and III, respectively), it is found that the major difference appears nearby the interface and the magnitudes of the differences vary with their Henry’s law constants distinctly. For this reason, the spatial distributions of six reactive species, i.e., O₃, N₂O₅, OH, HO₂, H₂O₂, and HNO₃, 10 µm above the interface and 3 µm beneath, are shown in Fig. 5. Their dimensionless Henry’s law constants are 0.25, 51, 942, 1.7 × 10⁵, 2.1 × 10⁶, and 5.1 × 10⁶, respectively, and they cover the typical range of solubility of common species in plasmas. The parameters for all species could be referred to the Henry’s law constants of set 2, given in Table III. The gaseous distributions of these six species are shown in Figs. 5(a)–5(f) and those of the liquid phase are shown in Figs. 5(g)–5(l), where the vertical axis represents density and the horizontal axis represents the distance to the gas–liquid interface. It is worth noting that there might be relatively large calculation error when the spatial scale is close to the collision mean free path, which is around tens of nanometers. Generally, the differences in distribution induced by different interfacial boundary conditions is seen in species with relatively large Henry’s law constants, such as HO₂, H₂O₂, and HNO₃, whose dimensionless constants are ∼10⁵–10⁶, and the discrepancy might cross several orders of magnitude above the interface. N₂O₅ is the only exception, which would be discussed later. Regarding the spatial range, the distinction only appears within micrometers above the interface and the aqueous distributions are hardly affected. The discrepancy in concentration induced by different interfacial boundary conditions is evaluated using Eq. (6), where the subscripts i and j represent the types of species and the model number, respectively. The relative deviation (RD for abbreviation) of concentrations to the minimal value in three results is calculated for each species,

Calculated using Eq. (6), species with larger Henry’s law constants tend to have larger RD in the gas phase and the RD of species with low solubility is negligible. However, the discrepancy would quickly diminish with the increasing distance away from the interface. Specifically, at the range of 1 µm above the interface, HO_2gas is the most affected species with RD(HO_2gas) = 98%, followed by HNO_3gas and H₂O_2gas, whose RDs are 95% and 92%, respectively. However, when the evaluated gap expanded to 20 µm above the interface, the largest RD is only ∼6.7% for all gaseous species, which indicates that the difference in species concentrations caused by different interfacial boundary conditions quickly vanishes within several micrometers. For aqueous species, their differences of concentration between the three models are negligible. Even within 1 nm depth under the interface, the RDs are generally no larger than 2% regardless of their solubility.

The effect of Henry’s law constants could be explained with the two-film theory in principle, which is widely applied to evaluate the mass transfer between two phases. Specifically, the mass transfer between two phases is jointly determined by both the gaseous and aqueous thin films clinging to the interface. The mass transfer rate could be described using the total gaseous uptake coefficient, K_G, according to the following equation:

where H represents the Henry’s law constant and k_L and k_G are the mass transfer coefficients in the liquid and gas phase, respectively, which are proportional to their diffusion coefficients, D, in the corresponding phase.

As observed in Fig. 5, the discrepancy between three interfacial boundary conditions exists for species with a relatively large Henry’s law constant. Taking H₂O₂ as an example, its Henry’s law constant is ∼10⁶ and its k_L is thousands of times smaller than k_G owing to the proportionality to the diffusion coefficient in each phase. Hence, the first item on the right hand of Eq. (7) is far larger than the second item, i.e., 1/Hk_L ≪ 1/k_G, and it is supposed to work in the gas-film control mode. In other words, the capacity of the solution is large enough and the mass transfer of H₂O₂ is rate limited by the supplement from the bulk gas phase. Regarding the production efficiency of reactive species by plasma, it may be not effective enough to sustain an instant solvation equilibrium at the interface, not to mention the whole liquid bulk. This results in different profiles of species between three models with diverse boundary conditions. For model I, the throughout solvation equilibrium is artificially installed at the interface, i.e., c_aq = Hc_gas. However, this is not the case for models II and III, where the solvation is progressively developed and it is not required to reach the equilibrium state throughout the discharge. When the Henry’s law constant is too large to sustain the proportion between gaseous and aqueous concentration as described in model I (or there is a fast process to reduce the aqueous species), a large negative concentration gradient is numerically required. This explains the remarkable concentration drop of H₂O₂ above the interface in model I. This pattern also exists for HO₂, HNO₂, HNO₃, etc., whose Henry’s law constants (no less than 10⁴) are large enough to work in the gas-film control mode according to Eq. (7). Moreover, no matter which interfacial boundary condition is used, the aqueous profiles are identical because the production of species in gas bulk is the same in three models. As for species with a relatively small Henry’s law constant, the solvation equilibrium is readily achieved in all three models, and thus, both gaseous and aqueous profiles are the same regardless of the interfacial boundary conditions.

It is worth noting that multiple reactive species present in PAW (plasma-activated water) and the liquid phase chemistry also largely contribute to the aqueous concentration. Fast destruction reactions could further enlarge the demand for mass transfer from the gas phase and enhance the concentration gradient above the interface in model I. In other words, the reactions might dominate in the reaction–diffusion competition⁵⁶ and Henry’s law constant should not be the only index used to evaluate the possible difference between different interfacial boundary conditions. This is the case of N₂O₅ as shown in Figs. 5(b) and 5(h), whose transportation across the interface is a typical reactive uptake process.⁵⁷ It would quickly hydrolyze or react with water to form HNO_3aq/NO_3aq⁻ once it enters the solution. This procedure makes the solvation saturation very hard to achieve even though its Henry’s law constant is only 51. It results in a concentration discrepancy with several times of magnitudes between the three models.

To estimate the spatial range of gaseous distribution of species under the influence of interfacial boundary conditions, EL, the following equation is proposed:

where D_g is the gaseous diffusion coefficient, c_g is the concentration in gaseous bulk, and R_p is the production rate of the given species, respectively. The maximum range of the gas gap being affected is estimated to be tens of micrometers above the interface, which is consistent with our observation in Fig. 5.

Spatially, the range of deviation caused by different interfacial boundary conditions, i.e., tens of micrometers, is much smaller than the realistic scale of a plasma jet interacting with a flat interface. However, one should still be careful if the volume of the interfacial gas layer is comparable to the bulk volume of gas phase in the studied situation. Interaction through microbubbles is a typical example. Some studies found out that the production of PAW with microbubbles (∼hundreds of micrometers in diameter) could enhance the efficiency of activation owing to the large surface-to-volume ratio, and the smaller size of bubbles is beneficial to the accumulation of reactive species in the liquid phase.⁵⁸ In this circumstance, as described by the thermodynamic interfacial boundary condition in model I, the dramatic decrease in highly soluble species near the interface may result in an underestimation of its gaseous density and consequently mislead the gaseous chemistry.

Above all, the interfacial boundary conditions would only affect the gaseous distributions in an interfacial gas layer with several micrometer thickness for species with large Henry’s law constants (H > 10⁴) because their mass transfer is rate limited by their production rate there. Species with low Henry’s law constants are hardly affected because their mass transfer is rate limited by their solubility. No matter which interfacial boundary condition is used, there is little influence on the aqueous profiles. Since the affected interfacial gas layer is very thin, if the spatial scale of the gas phase is much larger than that of the interfacial gas layer, i.e., tens of micrometers, all three interfacial boundary conditions have little effect on the characteristics of gaseous and aqueous bulk. However, special attention is recommended if the spatial scale of the gas phase is comparable to the thickness of the interfacial gas layer, where the thermodynamic boundary condition (model I) may lead to an underestimation of gaseous density and may further affect the evaluation of the mass transfer process.

The Henry’s law constant represents the solubility of a given gas, which is a prerequisite to estimate its mass transfer from the gas to liquid phase and thus of great importance to the evaluation of PAW. To assess how much these uncertainties would affect the numerical results, three sets of Henry’s law constants are adopted in model I, respectively, noted as sets 1, 2, and 3, which stand for the lowest, middle, and largest value within their corresponding uncertainties.²⁴ The Henry’s law constants in three sets approximately cover the common values in the present numerical works. Details of the parameters are given in Table III. The distributions of aqueous species are compared in Fig. 6. Figures 6(a)–6(h) show the distributions of species with different Henry’s law constants in three sets to ascertain how much the uncertainty in Henry’s law constants would affect their aqueous distributions, while Figs. 6(i)–6(l) show those of species with constant parameters in three sets to ascertain if the uncertainty would cause ripple effects to them.

As can be seen in Figs. 6(a)–6(h), the change in Henry’s law constants only affects the distribution of species with relatively small Henry’s law constants such as atomic oxygen and nitrogen, H_aq and NO_3aq, whose solvation equilibrium at the gas–liquid interface is readily sustained. However, it hardly affects the profiles of HO_2aq/O_2aq⁻ or OH_aq. To evaluate the deviation of aqueous concentrations due to the change in Henry’s law constants, Eq. (6) is applied for each species. It is found that the largest RD appears on O and N(²D), which is 57.5 and 55.3, respectively, followed by 48.6 for NO₃ and 4.4 for N. RDs of other species are generally no more than 5%, apart from 70% of H. It should be noted that even though the average concentrations of O, N [or N(²D)], H, and NO₃ change for orders of magnitudes in three sets, their dissolution depths and maximum concentrations are limited, i.e., no more than micrometer level and hundreds of nM, which are far below the scope of detection in experiments. Apart from the above species, the constants of N₂O₅ also expanded for two orders in the three sets, but its aqueous profiles are almost unchanged, which would be discussed later.

Generally, the differences in aqueous distributions only exist in species with low Henry’s law constants (H < 1). Concerning the relationship between the gaseous and aqueous concentration at the gas–liquid interface, it is found that the impact of Henry’s law constants depends on the accessibility of the steady solvation equilibrium. If the transportation from the gas phase competes over the destruction by aqueous reactions or downward diffusion along the concentration gradient, the reduction of species at the interface could be sufficiently replenished by mass transfer from the gas phase to saturate at once. Therefore, the gas–liquid mass transfer is rate limited by the solubility of the given species. If the Henry’s law constants get relatively larger, as in the case of set 3, more species would transport into liquid owing to the stronger solubility and their concentrations would increase as well. This is the case for O, N(²D), H, N, and NO₃, whose Henry’s law constants are small enough (no more than 1) to work in the liquid-film control mode.

However, this is not the case for N₂O₅, OH, and HO₂. Note that many species present in PAW and the aqueous chemistry would much impact the transportation of species. The aqueous destruction reactions may be too rapid to sustain the solvation equilibrium at the interface, or it would take a much longer time to achieve the equilibrium than the calculated timescale. As mentioned above, the capacity of N₂O₅ in liquid is boosted because of its reactive uptake process⁵⁷ and the production rate of gaseous N₂O₅ is inadequate to maintain a solvation equilibrium. Similarly, aqueous OH is also a short-lived species and it is quickly transformed into H₂O₂, HO₂, or other species. As for HO₂, it works in the gas-film control mode owing to its strong solubility in all three sets. The increase in its Henry’s law constant in modeling cannot enhance its solvation due to the limited gaseous concentration. In a word, when the aqueous destruction process is strong enough to compete over interfacial mass transfer or when its Henry’s law constant is already large enough to work in the gas-film control mode, their mass transfer is rate limited by their gas phase production. For a given discharge condition, the bulk concentration is certain, and therefore, the aqueous profile is unchanged.

Figures 6(i)–6(l) show the profiles of some species with the same Henry’s law constants in the three sets. Their distributions are similar in the three sets no matter whether their solubility is weak or strong. It is evident that the change in Henry’s law constants with generally strong solubility, such as in HO₂, does not impact its aqueous distribution because its gaseous supplement is limited and its mass transfer does not change much even though its Henry’s law constant is increased by 100 times. The variation in Henry’s law constant with weak solubility such as in H and NO₃ also impact little on other species even though their concentration may vary a lot between the three sets of Henry’s law constants. On the one hand, their mass transfer from the gas phase is little after all, and they are not pivotal species in the aqueous chemistry and contribute little to the production of other species, which would be discussed later. However, one should be careful to set the Henry’s law constant for species such as O, whose gas phase concentration is remarkably large and the aqueous chemical reactivity is supposedly high. O would quickly transform into OH in liquid, and the latter one is found to be the pivotal species in aqueous chemistry. For instance, concerning its fast destruction reaction, the Henry’s law constant of O is set to be 1 in Ref. 19 and its mass transfer would be much enhanced compared with this work, where its Henry’s law constant is set to be the same as that of O₂, owing to the identical element composition.^10,20 As a result, the production of OH increases significantly and then shifts the aqueous chemistry. This is because O is one of the most abundant species in the gas phase and it is supposed to be sufficient to alter the aqueous chemistry if its solubility is increased immoderately, even to work in the gas-film control mode. Hence, a detailed validation of its reasonable value by experiment is preferred in future works.

In a word, the uncertainty in Henry’s law constant would only affect the aqueous distributions for species with small Henry’s law constants (H < 1) because their mass transfer is rate limited by their solubility. Species with large Henry’s law constants are hardly affected because their mass transfer is mainly determined by their gaseous production. There is little effect on the whole picture of aqueous chemistry because the solubility of affected species is weak after all. However, special attention is necessary to the species with an estimated Henry’s law constant, which is dominant in the gas phase and highly reactive in the liquid phase.

Based on the investigation in Secs. III B and III C, the choice of interfacial boundary conditions and the uncertainty in Henry’s law constant have little effect on the spatial distributions of most reactive species in liquid in this circumstance. Therefore, a 10-s treatment of plasma to the deionized water is calculated with model I in Table II and with the Henry’s law constants as given in set 2 of Table III. The spatial profiles of main ROS and RNS at time = 0.1 s, 1 s, and 10 s are shown in Figs. 7(a)–7(c) and Figs. 7(d)–7(f), respectively, where the vertical axis represents the concentration and the horizontal axis represents the dissolution depth in the liquid phase. As shown in Fig. 7, H₂O_2aq is the most dominant ROS_aq. Its maximum local concentration and dissolution depth reach up to mM level and mm level after a 10-s treatment, respectively. In addition, even though the average concentrations of HO_2aq and O_2aq⁻ are much lower than that of H₂O_2aq, which are ∼72 nM and ∼2 nM, respectively, they could also dissolve deeply into ∼1 mm depth, with the maximum concentration of tens of μM. The distribution of OH_aq becomes steady quickly within 1 s. Its maximum concentration is around dozens of μM and the depth with concentration no less than 1 nM is of micrometer level. For RNS, HNO_3aq/NO_3aq⁻ and O₂NOOH_aq/O₂NOO_aq⁻ are the most dominant species, whose maximum concentration reaches up to mM level and dissolution depth reaches up to ∼1 mm. Compared with Figs. 7(e) and 7(f), the distribution of ONOOH_aq/OHOO⁻_aq becomes steady after 1 s, with the maximum concentration of ∼0.4 mM and the dissolution depth of hundreds of μm.

Interestingly, there are different patterns of the spatial distributions of reactive species. For hydrophilic species such as H₂O_2gas and HNO_3aq/NO_3aq⁻, their concentrations decrease monotonically with increasing dissolution depth, but this is not the case for HO_2aq and NO_2aq. Specifically, the concentration of NO_2aq would increase again when penetrating deeper into the liquid and the distribution of HO_2aq becomes step-like when the treatment time is longer than 1 s. These special profiles indicate that the localized aqueous chemistry competes over the mass transfer from the gas phase in their dissolution and the localized dominant chemical process may change with increasing treatment time or the dissolution depth. Thus, a spatially resolved long timescale modeling work is of importance.

To unravel the chemical process in liquid, the key chemical pathways of reactive species after a 10-s treatment are shown in Fig. 8, where the vertical axis roughly represents the depth of the solution. The key generation/loss pathways of one species are identified as the reactions whose rates dominantly contribute to the total generation/loss rate of that species in the calculation. Since it has been concluded from Fig. 7 that the key chemical processes may vary with depth, the key pathways are distinguished and calculated in two separated regions, namely, dissolution depth within 1 µm and deeper than 1 µm. The solid lines represent important pathways, while dotted lines represent unimportant pathways. Hollow arrows in the bars are for denoting that for these species, the diffusion dominates their dissolution toward the deep layer.

It can be concluded from Fig. 8 that taking the aqueous reactive species as a whole, the reactivity of liquid is introduced by gaseous OH_gas, H₂O_2gas, and HNO_2gas. The transportation of HO_2gas from the gas phase also contributes to its production in liquid, but it is less significant than the chemical reactions, which is in agreement with the findings from Fig. 7. HNO_3aq is the most hydrophilic species in the model, and the transportation from the gas phase dominates its aqueous production with no doubt. However, HNO_3aq and its conjugate ion NO_3aq⁻ are quite chemically stable in aqueous species that they contribute little to the aqueous chemistry but only to the acidity of the solution according to its low pKa.

For HO_2aq, apart from the mass transfer from the gas phase, it is primarily produced by the reaction between OH_aq and H₂O_2aq (R1). The reduction pathways of HO_2aq are different in the shallow and deep layers of the solution. Within the depth of no more than 1 µm, it is mainly eliminated by reaction with OH_aq to form O_2aq (R2). While transporting deeper into the solution and with the reduction of OH_aq, the combination of HO_2aq with NO_2aq to form O₂NOOH_aq (R3) ⁵⁹ gradually takes over control,

Regarding OH_aq, three quarters of them are received from the gas phase¹⁹ and the rest is produced by the collisional decomposition of H₂O_2aq by solvated electrons (R4). Part of them is then combined back into H₂O_2aq (R5),¹⁹ and the rest is then transformed into HO_2aq with H₂O_2aq (R1) or converted back into O_2aq with HO_2aq (R2), which is also the key destruction pathway in the deeper region. As shown in Fig. 8, even though the dissolution depth of OH_aq with a concentration larger than nM is only ∼1 µm, it plays an important role in the consumption of ONOO_aq⁻ [(R6) and (R7)] and NO_2aq⁻ (R8) at micrometers under the interface, owing to the reversible reaction (R9). Through (R9), OH_aq and NO_2aq combine into ONOOH_aq in the shallow layer of liquid. As ONOOH_aq diffuses into the deeper layer, it again decomposes back into OH_aq and NO_2aq,¹⁹ which is also the dominant production pathway for NO_2aq inside the solution. This explains the step-like profiles of NO_2aq,

Owing to the strong solubility and chemical stability, most of H₂O_2aq are transported from the gas phase⁶⁰ and then diffuse downward. In addition, part of them are produced by the recombination of OH_aq via (R5) ⁶⁰ and most of them are consumed by the solvated electron to form OH_aq (R4) or by OH_aq to form HO_2aq (R1).

The rest of ROS and RNS are relatively less connected to the aqueous chemistry system according to Fig. 8, owing to their weak solubility or low gas phase concentration.

Atomic oxygen, O_aq, and ozone, O_3aq, are mostly produced by mass transfer from the gas phase. Most of O_aq are consumed by H₂O_aq and H₂O_2aq to form OH_aq^19,60 and HO_2aq (R10) and (R11). Half of O_3aq is oxidized by OH_aq into HO_2aq (R12), and the rest is transformed into HO_3aq by HO_2aq (R13). NO_3aq and N₂O_5aq are also mostly absorbed from the gas phase. NO_3aq is then converted into NO_3aq⁻ by HO_2aq (R14) or charge transfer reaction (R15), while N₂O_5aq quickly hydrolyzes into NO_3aq⁻ via (R16). For more details on the complete reaction pathways, refer to the  Appendix,^{50,55,61–139} 

Three one-dimensional fluid models are developed and compared in this paper for the interaction between a plasma jet and deionized water, considering the difference in interfacial boundary conditions and the uncertainty in Henry’s law constants as reported in the literature. The plasma jet is excited by a radio-frequency voltage, and the mixing of ambient air into the He + 0.5 vol. % air feeding gas is considered. For interfacial conditions, three cases (models) are calculated, and each one has a unique interfacial boundary condition for neutral species, i.e., the thermodynamic, diffusion, or kinetic interfacial boundary conditions. These interfacial boundary conditions are typical ones as reported in the recent literature. For the uncertainties in Henry’s law constants, three cases (models) are calculated as well, and each one employs a set of Henry’s law constants representing the maximum, medium, or minimum value within their uncertainties. Then, the dependence of density distributions of reactive species, both in the gas and liquid phases, on the interfacial boundary conditions and the Henry’s law constants is obtained to support the validity evaluation of the simulation method.

It is found that the density distributions of reactive species are nearly independent of the interfacial boundary condition in both the gas and liquid phases, except for the interfacial gas layer with a thickness of several tens of micrometers above water. In such a gas layer, the densities of high soluble reactive species (H > 10⁴) are different by several orders of magnitude among the interfacial boundary conditions. Although the densities in the gas layer are sensitive to the interfacial boundary condition, the particle fluxes on the water surface are not; hence, nearly nothing has changed in the plasma-activated water with the interfacial boundary conditions. When the Henry’s law constants are changed, it is found that only low soluble reactive species (H < 1) have their aqueous densities following the change. Such species mainly exist in the shallow layer of the water and their densities are very low, so generally the plasma-activated water has not changed. Therefore, it could be concluded that the simulation of plasma–liquid interaction is generally independent of the given interfacial conditions and the uncertainty in Henry’s law constants.

The density distributions and chemical pathways of reactive species in the deionized water are obtained by simulation. It is found that the original aqueous reactive species that dissolved from the gas phase are mainly H₂O_2gas, OH_gas, and HNO_2gas, and the dominant ROS and RNS obtained after the 10-s plasma treatment are H₂O_2aq and HNO_3aq/NO_3aq⁻, respectively. Regarding the liquid chemistry, it is interesting that the stepwise oxidation of HNO_2aq by OH_aq into NO_2aq and then into ONOOH_aq plays an important role in the shallow layer of the water, and ONOOH_aq could diffuse into a deeper region and decompose there, making the short-lived species NO_2aq exist at a depth of more than 1 mm underwater.

This work was supported by the National Science Foundation of China (Grant Nos. 51722705 and 51837008), the Fok Ying-Tong Education Foundation (Grant No. 161055), the State Key Laboratory of Electrical Insulation and Power Equipment (Grant No. EIPE19122), the fund of the State Key Laboratory of Laser Propulsion and Application (Grant No. SKLLPA-29), and the Fundamental Research Funds for the Central Universities (Grant No. xzy022019034).

The data that support the findings of this study are available within the article.

Tables IV and V are given in Appendix.