Modeling of reactive species interphase transport in plasma jet impinging on water

The interaction between low-temperature atmospheric pressure plasma and water is of primary relevance to an increasing number of applications, from water treatment to medicine. The interaction between an argon plasma jet and water is investigated using a three-dimensional (3D) time-dependent computational model encompassing turbulent gas flow and induced liquid motion, gas–water interface dynamics, multiphase species transport, and gas- and liquid-phase chemical reactions. A single-field approach based on the volume-of-fluid (VoF) method together with conditional volume averaging (CVA), is used to consistently describe the dynamics of the interface together with interfacial reactive mass transfer. Three CVA-based interface species transport models, based on arithmetic, harmonic, and unified mixture species diffusivities, are evaluated. Simulations of a plasma jet impinging on water at different gas flow rates are presented. The resulting deformation of the interface and the production and accumulation of hydrogen peroxide, reactive oxygen, and nitrogen species corroborate prior findings in the research literature showing that higher jet velocities and associated increased interface deformation led to the enhanced transport of reactive species across the plasma-water interface. The VoF-CVA approach appears promising for the modeling of general plasma-liquid multiphase systems.


Introduction
The interaction between plasma and liquid is gaining increased interest due to its principal relevance in a wide range of Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.applications, spanning from water treatment, agriculture, and plasma-activated water [1][2][3], to sterilization, disinfection [4][5][6][7], and medicine [8][9][10][11].Several applications exploit the interaction between an atmospheric pressure plasma jet (APPJ) and water.The need to advance such applications and understand the underlying physical-chemical processes due to the interaction of a plasma jet and liquid water has motivated diverse experimental and computational studies.
Experimental studies of the interactions between an APPJ and water have been mainly based on using helium, argon, and air as working gases.Uchida et al [12] experimentally investigated the characteristics of a helium APPJ operating with flow rates from 1.5 to 10 slpm impinging on deionized water.Their results showed that the radial expansion of the plasma afterglow after impinging on the water interface rapidly grows with increasing gas flow rate.This expansion produces a larger surface for interaction between plasmagenerated species and water, favoring plasma-water treatment processes.Brubaker et al [13] investigated the characteristics of a helium APPJ interacting with water for varying flow rates and applied voltages.They compared the cavity depth formed by the interface deformation due to the impingement of gas and plasma jets.Their results showed that the use of a plasma jet leads to significantly greater cavity depths and significantly greater evaporative cooling of the water compared to those using a gas jet.Also, their investigations showed that the acidity of treated water increases near the impingement point and propagates into the bulk liquid following the recirculation patterns induced within the liquid.The main contributor to acidification was the dissociation of plasma-generated reactive oxygen and nitrogen species.Park et al [14] demonstrated that a helium APPJ significantly increases the hydrodynamic stability of the interface cavity compared to that obtained with a pure helium jet.Their experiments were conducted for varying flow rates, as well as for different widths and heights (amplitudes) of the 10 kHz high-voltage pulses driving the electrical discharge.Their results revealed that higher widths and heights of the voltage pulses increase cavity depth and stability due to increased surface charge along the interface.Van Rens et al [15] conducted experiments on the impingement of a kINPen ® -generated argon APPJ onto water.They showed that, for a given flow rate, the depth of the cavity formed was more significant with a plasma jet compared to a gas jet (i.e.plasma on versus plasma off), corroborating the findings reported in [13].Liu et al [16] experimentally investigated an air APPJ impinging on water in penetrating mode (i.e.relatively large interface cavity depths).The investigators analyzed the effects of jet height and diameter, gas flow rate, and liquid sheet depth on the dry spot diameter under the penetration point and characterized the associated generation of surface waves.Winter et al [17] employed Fourier transformed infrared spectroscopy and laser-induced fluorescence spectroscopy for the investigation of the dependence of the concentration of gas-phase hydrogen peroxide (H 2 O 2 ) in an argon APPJ impinging on a liquid cell culture medium.Their results showed that the transport of gas phase-generated H 2 O 2 is the primary source of H 2 O 2 in the liquid compared to liquid-generated H 2 O 2 .
Computational modeling of plasma-liquid interactions complements experimental studies, primarily to provide insight into characteristics not practically accessible experimentally, but also to guide equipment and process design and optimization.Studies of the impingement of gas jets onto water provide the basis for investigating APPJ-water interactions.Hwang and Irons [18] performed computational studies of the impingement of a gas jet on a water bath.They obtained relations for cavity depth and width as function of the net momentum supplied by the gas flow.Nguyen and Evans [19] computationally investigated gas-liquid interface deformation under the impingement of a gas jet onto water.Their model included surface deformation captured using a volume-of-fluid (VoF) approach and the Young-Laplace equation to describe the pressure jump across the interface.Their model was implemented in Ansys Fluent ® and validated against analytical results.Ersson et al [20] performed computational modeling of a top-blown bath to analyze jet impingement depth at the water surface and the induced flow motion in the bulk water.Their numerical model was based on a VoF method to capture interface deformation.Their results presented good agreement with experimental measurements.Muños-Esparza et al [21] numerically investigated the impingement of an air jet onto water under different surface deformation modes using a two-dimensional (2D) transient VoF model implemented in Star-CCM+ ® .Their results were validated against experiments using particle-image-velocimetry and level detection and recording measurements, showing good agreement.
Computational investigations of APPJs interacting with liquid are scarcer than those dealing with gas jets.This is partly due to the relative novelty of the relevant applications and the increased complexity of describing plasma flows, especially within multiphase systems.Levko et al [22] presented the computational modeling of an atmospheric pressure air discharge over a liquid water cathode.Their model accounted for plasma generation in the air through electron impact, surface emission, and the induced motion of the liquid.Their results of the water electrode and gas breakdown dynamics corroborated experimental observations.Verlackt et al [23] performed computational modeling of an argon APPJ impinging on water.Their numerical simulations showed the relocation of reverse vortexes formed in the water towards the interface.The relocation effect was more pronounced for greater jet velocities, significantly affecting the accumulation of long-lived reactive species within the liquid.Lindsay et al [24] presented a 2D axisymmetric model of the interaction between an air APPJ and water.Their model accounted for conservation of mass, momentum, thermal energy, and chemical species within a static pre-defined jet-water interface.Their results emphasized the importance of convective transport induced in the liquid and of interface deformation on the distribution and the volume-averaged uptake of hydrophobic species.Heirman et al [25] computationally investigated plasma-treated liquids using a model that combined a 2D axisymmetric fluid dynamics model and a zero-dimensional (0D) chemical kinetics model.They showed that only long-lived RONS (i.e.H 2 O 2 , HNO 2 , HNO 3 , HO 2 , O 3 , and ONOOH) could accumulate in the bulk liquid before the plasma treatment.After treatment, only H 2 O 2 , HNO 2 , and HNO 3 concentrations remained relatively unchanged.Semenov et al [26] numerically investigated interphase transport phenomena in an APPJ impinging on water.Their model was validated with experimental findings of water cavity depth and hydrogen peroxide (H 2 O 2 ) accumulation within the liquid.The investigators demonstrated that the concentration of hydrophilic species in the liquid does not affect the transport in the gas phase under representative experimental conditions.Within investigations of plasma-liquid interactions, gasphase, and liquid-phase chemistry have been studied to a significantly greater extent than the processes happening at the plasma-water interface, particularly interphase transport, namely the transfer of reactive species across the gas and liquid phases.The present work addresses the modeling of species transport in plasma interacting with liquid-particularly, the coherent handling of the conservation of chemical species across phases within an evolving plasma-liquid interface.Most computational studies of plasma-on-liquid systems have relied on a pre-defined plasma-liquid interface (e.g.[22][23][24]).In such case, conventional (single-phase) species conservation equations can be solved in the plasma domain and in the liquid domain, coupled by a constitutive relation among the species fluxes between phases (such as Henry's law to describe the thermodynamic equilibrium of chemical species at the interface) imposed along the pre-defined interface.If the interface is not pre-defined, then a coherent approach to describe species transport considering an evolving interface has to be used.The present work addresses this challenge with a single-field approach based on the VoF method together with conditional volume averaging (CVA).The VoF-CVA approach has been shown to avoid known issues in multiphase transport models, such poor stability, accuracy and/or loss of conservativeness [27][28][29][30], and hence it is suitable for the simulation of more complex plasma-liquid systems, such as plasma-within-liquid or plasma interacting with bubbles or mist.Three interface species transport models based on CVA are evaluated, namely models based on arithmetic, harmonic, and unified mixture diffusivities [28][29][30].
This study focuses on the impingement of an argon APPJ from the kINPen ® plasma source onto water, as schematically depicted in figure 1, as a representative plasma-liquid system of practical relevance.Previous studies have shown that the uptake of reactive species from an APPJ into the liquid phase is a very effective process, even if the species' surface reactivity is minimal [31].In this context, the present work addresses the transport and formation of reactive chemical species due to an APPJ impinging on water.The computational model is developed using the OpenFOAM computational framework [27].The model encompasses gas jet dynamics, interface deformation, induced liquid motion, interphase transport of reactive species, and chemical kinetics in the gas and liquid phases.The jet flow rate is treated as the main control parameter, as it has a dominant effect on inducing liquid motion and, consequently, on the accumulation of reactive species within the liquid.This work aims to validate the singlefield VoF-CVA approach as a general and effective method for the modeling of plasma-liquid systems, as well as to provide insight into the effects of impinging APPJ velocity on surface deformation, the transport of reactive species across the interface, and their accumulation within the water.

Mathematical model
The mathematical description of the interaction between a plasma jet with liquid water requires the consistent coupling of a multiphase fluid flow model, a model to describe interface species transport, and a model of chemical kinetics in each phase.These models are described next, follow by their computational implementation.

Multiphase fluid flow model
To capture the phenomena of a dynamic two-phase interface, the model employs the VoF method.Two main interface reconstruction techniques are used to separate discretization cells into different domain phases in multiphase flow solversgeometric and algebraic.The former type, which includes the VoF approach, is based on explicitly reconstructing the interface's representation within each computational cell with some pre-defined functional form.Although geometric methods prevent the artificial proliferation of the interface due to numerical factors, their implementations tend to be complex and lead to the increased computational cost of simulations [32].In contrast to geometric methods, algebraic methods avoid the numerical reconstruction of the interface and instead rely on algebraic manipulations (e.g.compressive differencing schemes) to describe its effects on the flow.Although algebraic methods are more straightforward and computationally inexpensive, the interface obtained spans several cells instead of being defined at the sub-cell level as with geometric methods.
In the present model, the dynamics of a two-phase flow are described by the equations of conservation of momentum, mass, and of the VoF indicator function, in this case, the liquid fraction α (i.e.α = 1 implies liquid phase, α = 0 indicates gas phase, and 0 < α < 1 represents the region near the interface) [33].
Conservation of momentum is given by: where t represents time, ∂ t ≡ ∂/∂t the temporal derivative, ∇ ≡ ∂/∂x the gradient operator with x as the spatial position vector, u velocity, ρ density, p pressure, g acceleration due to gravity, τ the stress tensor, and the last term on the right-handside describes the surface tension force acting on the interface.In the last term in equation (1), the dynamic interface is denoted as S(t), σ is the surface tension coefficient between the two phases, κ is the local surface curvature, n the vector normal to the interface, and δ represents the Dirac delta function.Mass conservation, assuming argon and water as incompressible fluids, is expressed as: whereas the equation describing the evolution of the liquid fraction α is given by: (i.e. a material derivative of α equal zero) indicating that α gets transported with the flow.This study considers argon and liquid water as the primary fluids at standard temperature and pressure conditions (the concentration of reactive species from the plasma is relatively small).Hence, the gradient of the stress tensor is described as for an incompressible Newtonian fluid, namely: where µ is the dynamic viscosity, assumed constant for each phase (gas or liquid), and is related to the kinematic viscosity ν by ν = µ/ρ; and T indicates the transpose operator.
The last term on the right-hand side of equation (1) describes surface tension forces.These are described by the continuum surface force approach developed by Brackbill et al [34], namely: This approach leads to a description of the interface force that directly depends on the local variation of the liquid fraction α.
The model uses the auxiliary variable p rgh defined as: to incorporate the combined effect of hydrostatic pressure and gravity.The main purpose of using p rgh is to avoid severe pressure gradients due to hydrostatic effects.Substitution of equations ( 4), ( 5) and ( 7) into the momentum equation leads to the final form of momentum conservation used in the model, i.e.
Turbulence modeling is essential to describe APPJs operating under conditions leading to turbulence (i.e.large Reynolds numbers) and to correctly predict the deformation of the gas-liquid interface.The present model adopts Chien's low-Reynolds k-ε turbulence model [35].Chien's k-ε model is better suited for describing the impingement of circular jets than the standard k-ε model [36].The k-ε model is composed of two additional governing equations describing the evolution of the turbulent kinetic energy k: and the evolution of turbulent kinetic energy dissipation ε: with: and Chien's k-ε model uses the following coefficients' values: The turbulent kinematic viscosity ν t is determined from: Using as the wall distance y the distance to the plasmawater interface, as adopted in [37,38], the wall-bounded turbulent Reynolds number is given by: The modified wall distance y * , used in equation (11), is given by: Finally, to complete the definition of the model, the walldamping function f ν and the damping function f 2 are given by: and

Multiphase species transport
The conservation equation for a generic chemical species through a uniform medium is given by: where the subindex i denotes the ith chemical species, C i is the specie's molar concentration, D i is its total diffusion coefficient, and R i is its net chemical reaction rate.
If the gas-liquid interface is known, then an equation ( 16) can be solved in each phase (gas and liquid), coupled by a constitutive relation among the species fluxes between phases (such as Henry's law) imposed along the pre-defined interface.If the interface is not pre-defined, as it is the case in models that seek to capture the interface such as those adopting the VoF method, the description of species conservation within a multiphase system is not trivial.There are primarily two different approaches to model mass transfer in multiphase flows within VoF methods.These are based on arithmetic averaging and harmonic averaging of diffusivities across phases, as developed by Marschall et al [29] and Haroun et al [28], respectively.In the first model, Marschall et al introduced the concept of the continuous species transfer (CST) method, which is based on CVA to get a single-field formulation (i.e. a single-field to describe properties across both phases) for species transport, similarly as for the other fields involved in the VoF formulation (pressure, velocity, momentum, etc).The model by Marschall et al is largely analogous to the earlier model by Haroun et al [28], except for the manner of handling species fluxes across the interface.In the present work, both approaches are implemented, as well as their generalization in the unified model developed later by Deising et al [30], which is based on the CVA of diffusivities and accounts for the curvature effect on the gradients of concentrations.These models are described next.

Arithmetic mean mixture diffusion model
The multiphase species transport model developed by Marschall et al [29] is used to model species mass transfer across immiscible interfaces in free-surface flows.In this approach, the solubility is determined by Henry's law with a constant coefficient.The conditions that are satisfied at the interface are Henry's law for the ratio of liquid and gas species concentrations and the equality of mass fluxes at the gas and liquid sides of the interface, namely: In equation ( 17), H i is the Henry's law coefficient (solubility constant) for species i in the gas phase diffusing into the liquid, and the superscripts l and g denote liquid and gas, respectively.Equation ( 18) implies the conservation of mass of species C i across the interface (i.e. the interface does not store mass).
The single-field species transport equation is attained following the CVA method within the VoF approach, leading to the species interphase transport equation: In equation ( 19), the single-field expressions to describe species concentrations, diffusivity, viscosity, and density are respectively given by:

Harmonic mean mixture diffusion model
The model developed by Haroun et al [28] is based on using the harmonic mean of mixture diffusion coefficients to describe species transport across the interface.This model defines an interface species harmonic diffusivity as: This diffusion coefficient is derived from the definition of Henry's law solubility constant (equation ( 17)).
Similar to the arithmetic mean mixture diffusion model, the single-field species transport equation is attained following the CVA approach, leading to the single-field species conservation equation: It can be noted that, by contrasting equations ( 19) and ( 25), the model by Haroun et al does not have an extra term on the right-hand side involving the gradient of the diffusion coefficient (C i ∇D a i ).

Unified mixture diffusion model
More recently, Deising et al [30] presented a model based on CVA that extends the CST method and unifies the above two approaches for handling interface species transport.The CVA model leads to the following single-field species transport equation: A comparison of the above multiphase species transport models reveals that the unified mixture diffusion model has an extra term describing curvature effects on interface species transport (i.e. the fourth term on the right-hand-side of equation ( 26)).

Multiphase chemical kinetics
The net reaction term in the multiphase species transport equations (equations ( 19), ( 25) and ( 26)), namely R i , is given by [38]: where M w,i is the molecular weight of species i, N r is the total number of reactions, and R i,r is the Arrhenius molar creation rate of specie i in reaction r.The general form of the reactions system is: where N s is the number of chemical species in the system, ν ′ i,r and ν ′ ′ i,r are stoichiometric coefficients, and M i describes the number of molecules of species i.The chemical kinetics models for the gas and liquid phases are formulated with forward reactions only (instead of forward and backward reactions).Therefore, the molar rate of creation of specie i in reaction r is given by: where C j,r is the molar concentration of specie j in reaction r (in units of mol m −3 ), η ′ j,r and η ′ ′ j,r are rate exponents for reactant and product species j, respectively, in reaction r; and k f,r is the reaction rate constant.
The reaction rate constants have the Arrhenius form, modified to handle reactions near the interface following the VoF-CVA approach.Specifically, for reactions happening in the liquid phase, the reaction rate constants have the form: whereas for the gas phase reactions, they have the form: In the above expressions, f (α) is a smoothstep function that determines the extent of a reaction across each phase, A r is the pre-exponential factor, T is the absolute temperature, β r is the exponential temperature factor, and T a is the activation temperature.The smoothstep function used in the model implementation is:

Computational implementation
The computational implementation of the mathematical model is based on the finite volume method (FVM) discretization implemented in the interFoam solver within OpenFOAM.OpenFOAM [27] is an open-source framework based on the C++ language with serial-and parallel-computing capabilities for implementing transport models using FVM discretizations.The interFoam solver is designed to model two-phase incompressible isothermal flows based on the color function volume-of-fluid approach.Phases are described using an Eulerian formulation, while the interface is modeled using an indicator function (local volume fraction of one phase, which is taken as the liquid fraction α in the present work).The functionality of the interFoam solver is augmented by adding the low-Reynolds k-ε turbulence model and the interphase transport of species model described above.
The numerical discretization of the VoF model utilizes the multidimensional universal limiter with explicit solution (MULES) approach [39].MULES is a numerical interface compression method that limits the numerical representation of fluxes to preserve the sharpness of the interface.The technique adds a heuristic term to equation (3), i.e. a compression term that improves the numerical solution at the interface, namely: where u r is the relative velocity between phases defined by: with n f is the face-centered interface normal vector, C r is a compression factor, S f is the cell face surface area, and φ is the volume flux (i.e.rate of change in volume per unit area).The maximization operator (max) in equation ( 34) is applied over the whole domain, whereas the minimization operator (min) is performed locally at each cell face f.The third term on the left-hand side of equation ( 33) aims to minimize the numerical diffusion of the gas-liquid interface.More details of the model are presented in the work by Deshpande et al [33] and in the work by Rusche [40].
The temporal discretization of the model's transport equations is based on the first-order-accurate implicit Euler scheme.Spatial gradients are handled via second-order Gaussian integration with linear interpolation.Laplacians are discretized using the Gauss scheme with linear interpolation.Surface normal gradients are computed using the explicit nonorthogonal correction method.The linear upwind scheme is used to discretize the velocity advection, and the discretization of species concentrations uses the limited linear differencing scheme.
The following solvers are used to solve the linear systems of equations attained from the discretization of the model's equations: smooth linear using Gauss-Seidel smoother for the liquid fraction and momentum conservation equations; algebraic geometric multi-grid for the pressure equation (i.e.pressure correction equation derived from the equation of conservation of total mass); and the preconditioned bi-conjugate gradient for the species concentrations equations.
The computational implementation of the model has been validated with three benchmark cases, which assess the main aspects of the model, namely the description of interface dynamics, multiphase species transport, and multiphase chemical kinetics.The description of the benchmark cases and the attained results are presented in appendix.All the validation studies showed good agreement between the results obtained with the multiphase species transport model and the reference results.

Model set-up
The interphase species transport model, based on the VoF-CVA approach an implemented with a FVM discretization in OpenFOAM, is used to study the interaction between an APPJ and water following the work by Semenov et al [26].The setup consists of the impingement of the jet afterglow from a pulsed steamer discharge produced from a kINPen ® plasma source device placed perpendicularly on top of a cylindrical container partially filled with water.The computational model set-up, depicting the 3D discretization of the spatial domain and its boundaries is shown in figure 2. As shown in figure 2(a), the computational domain encompasses one quadrant (90 • ) of the complete cylindrical domain.The distance from the kIN-Pen's nozzle to the interface h is 10 mm.The water container has a diameter D of 3 cm and is filled to a height d of 2.5 cm leading to a total water volume of 17.6 ml.The domain boundaries in figure 2(b) include the symmetry planes given by x = 0 and z = 0.The cross-sectional view of the mesh in figure 2(b) shows the finer discretization used near the region of expected interface deformation (y < 0) and that the finite volumes near the center of that region are approximately equally-spaced.
The model described in section 2 allows simulations encompassing turbulent gas flow and induced liquid motion, air-water interface dynamics, multiphase species transport, and gas-and liquid-phase chemical reactions.The actual plasma jet from the pulsed streamer discharge within the kINPen ® is not explicitly modeled due to the high computational costs associated with capturing all relevant phenomena.Instead, the model describes the impingement of an argon gas jet seeded with reactive RONS depicting the afterglow from the kINPen device.Furthermore, the discharge environment is assumed to be composed mostly of argon (instead of air).The species included in the model for both gas and liquid phases are OH, H 2 O 2 , NO, NO 2 , N 2 O 4 , and HNO 2 .The set of species and chemical reactions used in the simulations are adopted from Lindsay et al [24] and are shown in table 1.The set of seven reactions describes the evolution of representative RONS within the gas and water domains.The model considers a uniform temperature T equal to 300 K throughout the gas and liquid phases, and hence values of the reaction rate constants used in the simulations are explicitly indicated in table 1.

Boundary conditions
The boundaries conditions for the model of a APPJ impinging on water are listed in table 2. To replicate the experimental conditions in [17], seven values of jet flow rate Q are considered in the simulations.The values of Q range from 1.1 to 1.7 slm with steps of 0.1 slm.Given the diameter of the nozzle equal to 1.5 mm, the average inlet velocity u inlet range from 7.2 to 11.2 m s −1 with steps of 0.66 m s −1 .The velocity profile imposed over the Inlet boundary is set as u inlet = [0 −u y , inlet 0] T , where u y , inlet (r) corresponds to a parabolic velocity profile with mean value u inlet and r = (x 2 + z 2 ) 1/2 is the radial coordinate.Pressure is set equal to atmospheric pressure p atm = 1 atm at the Outflow radial and the Outflow top boundaries.The liquid fraction α is set to 0 at the Inlet boundary, corresponding to a pure gas phase.
The k-ε turbulent flow model requires appropriate definitions of inlet quantities to describe the level of turbulence in the inlet gas jet.The following empirical relations obtained from [20,38] are used to determine inlet turbulence quantities, k inlet and ε inlet : where d device is the diameter of the jet and is equal to 3 mm.Since the Reynolds number at the inlet ranges from 1670 to Inlet concentrations for all chemical species C i,inlet in the model are obtained from Lindsay et al [24].The resulting chemical species concentrations are specified at the Inlet boundary using the values listed in table 3.For all variables at boundaries without specified values, zero-gradient (Neumann) boundary conditions are employed.
Transport properties for all the species in the model are adopted from Lindsay et al [24].The values of Henry's constant H and diffusion coefficient for the gas D g and liquid D l phases for each species are listed in table 4. The liquid diffusion coefficient for N 2 O 4 is estimated equal to 1.5 10 −9 m 2 s −1 based on its molecular weight and the values for NO and NO 2 , whereas the liquid diffusion coefficient for HNO 2 is estimated equal to that for HNO 3 [24].

Simulations set-up
The simulations are initiated from a static environment, with no gas of liquid flow (i.e.zero initial velocity through the gas and liquid phases through the whole domain), and a flat gasliquid interface.The simulations then proceed with the development of the gas jet and the induced deformation of the gasliquid interface.These purely multiphase (i.e.no species transport) stage of the simulations lead to the steady-state profiles of velocity u, pressure p, and liquid fraction α throughout the gas and liquid domains.Next, these steady-state fields are used as the initial conditions for the simulations of reactive species interphase transport.This de-coupled approach allows a drastic reduction of the computational time required to obtain time-dependent species transport results for up to 100 s of simulation time.

Multiphase flow
Steady-state velocity profiles and streamlines in the liquid and gas phases for four inlet gas flow rates, from 1.1 to 1.7 slm, are presented in figure 3. The results were attained after more than 1600 characteristic travel times, where the characteristic travel time is defined τ travel = h/u inlet .Upon reaching the unperturbed water, the transferring of the momentum from the gas flow onto the water leads to the surface deformation of the interface.Subsequently, the gas flows radially along the interface towards the water container walls.Once the gas flow reaches the wall, it turns upwards, creating vortexes that carry gas back towards the plasma jet region.Two vortexes form in the gas domain for all flow rates.With increasing flow rate, the center of the bigger recirculation area moves toward the axis of jet impingement and the secondary vortex formed near the water container wall strengthens.The results in figure 3 show that the recirculation patterns within the gas domain are similar to those obtained in the modeling of the induced flow by a corona discharge reported in [26,54].
The movement of the gas along the interface creates shear stress on the interface, inducing water motion.The water starts flowing in the same direction as the driving gas flow, radially towards the container walls.The maximum velocities within the water are located near the gas-liquid interface, close to the domain's main axis.The maximum induced liquid velocity, for the case of Q = 1.7 slm (u inlet = 11.2 m s −1 ), is 0.045 m s −1 .Upon reaching the wall, the liquid flow turns downwards, creating a recirculation region whose center is indicated by the arrows in figure 3. The center of the recirculation region moves radially towards the water container wall with the increase of gas flow rate.For the highest inlet flow rate of 1.7 slm, an additional recirculation zone in the water is formed near the bottom corner of the container.The occurrence of recirculation can have important implications for the uptake and distribution of reactive gas species, as discussed in section 4.2.
For validation of the numerical simulation results, interface cavity depths formed by different flow rates are compared with the experimental results by Winter et al [17] for a similar setup.Seven values of Q ranging from 1.1 to 1.7 slm with a step of 0.1 slm are considered.A comparison of the results is presented in figure 4(a).The error bars in the experimental results depict inherent variations in the experiments, whereas the error bars in the numerical results represent the extent along the y axis of the discretization volume enclosing the interface.The numerical results show that the dimple depth increases with the inlet gas flow rate, from 1.3 mm for Q = 1.1 slm to 3.4 mm for Q = 1.7 slm.The results in figure 4(a) show good agreement between the numerical and experimental results for the intermediate and low flow rates, and that the numerical results over-predict the dimple depth for the higher flow rates.The discrepancy with the experimental results for the higher flow rates can be attributed to the low-Reynolds number k-ε   turbulence model used.Moreover, for the flow rates Q = 1.6 and 1.7 slm, the numerical simulations present slight oscillations of the interface.These oscillations persisted after a time of 6300 τ travel , which can also be attributed to limitations of the k-ε turbulence model to describe the dissipation of the gas flow.Figure 4(b) presents the zoomed-in view of the jet impingement area for four representative flow rates.The color scheme used corresponds to that employed in figure 3. The horizontal lines and arrows depict the extent of the dimple in each case.These results depict the level of surface deformation attained with increasing flow rate.Moreover, the higher velocity magnitudes within the liquid-side of the interface indicate the region experiencing the highest shear from the gas jet.

Multiphase species transport
One advantage of numerical models is that they allow unveil the role of particular mechanisms in manners that are largely impractical or impossible experimentally, for example, by artificially activating or de-activating terms in the model equations.To discern the role of the reactivity of the species considered in the model within their assimilation within the water, multiphase species transport simulations are performed with the complete model described in section 3 (reactive) and with the same model but without the chemical reactions in table 1 (nonreactive) for the representative inlet gas flow rate of Q = 3 slm.This flow rate corresponds to the value used by Winter et al [17] of hydrogen peroxide uptake in water, even though the results in figure 4 indicate that greater discrepancies between experimental and computational results of dimple depth can be expected at flow rates above 1.6 slm.Given the results of benchmark tests in appendix A.2, the harmonic mean mixture diffusion multiphase species transport model by Haroun et al [28] is adopted in the simulations.
Following the approach described in section 3.4, first steady-state velocity, pressure, and liquid fraction distributions for Q = 3 slm are obtained.Subsequently, using the inlet concentrations listed in table 3, multiphase species transport simulations with chemical reactions de-activated (nonreactive) and activated (reactive) are run for a physical time of 1.0 s, corresponding to 3950 τ travel , leading to a relatively significant distribution of species throughout the gas and liquid domains.
The resulting distributions of species concentrations from the reactive and nonreactive flow simulations at time t = 1.0 s are presented in figure 5. Linear and logarithmic scales are  [17] for an inlet gas flow rate of 3 slm.used for the gas and liquid phases, respectively, to discern the variation of concentrations more clearly through the domains.The plots show velocity streamlines through both the gas and liquid domains.The formed streamline patterns are similar to those presented in figure 3. Nevertheless, given the higher inlet gas flow rate of 3 slm, the center of the gas recirculation region is closer to the jet region.Moreover, the depth of the dimple is significantly greater than those reported in figure 4, as expected due to the greater inlet flow rate.The greater extent of the dimple leads to a greater surface area and greater containment of reactive gas species, which affect their intake by the water.
The results in figure 5 illustrate that the reactive species generally follow the velocity streamlines in the gas, indicating the dominance of advective over diffusive transport (i.e. the term ∇ • (uC i ) dominates over ∇ • (D i ∇C i ) in equation ( 16)), as expected given the high inlet gas flow rate.Then, the species are transported into the water following the flow of the tangential advective transport at the gas-water interface together with their solubility in water quantified by the species' Henry's constant (table 4).The dissolved species are subsequently transported within the water due to the combined effects of the induced flow patterns and the species' diffusivity in the liquid phase (D l in table 4).
Since N 2 O 4 is not injected at the inlet, its concentration is zero throughout the domain for the nonreactive simulations.N 2 O 4 can only be generated by chemical reactions (R1 in table 1), which is consistent with the reactive simulation results.The reverse situation occurs for OH, which is rapidly consumed due to reactions R3-R6 in table 1 to create HNO 2 and H 2 O 2 in both phases.The comparison of the results for NO shows that it is also being consumed in reactions R3 and R6 in the jet impingement area to create HNO 2 .All NO 2 is converted to N 2 O 4 in the gas phase because of R1.
The results in figure 5 also show the drastic range of effects that chemical reactions can have on the distribution of species within the water.From one end, chemical reactions have a minor role in the uptake of H 2 O 2 , HNO 2 , NO, and NO 2 , as quantified by their comparable distribution within the water for both the nonreactive and reactive simulations.In contrast, all N 2 O 4 in water is due to its net production within the gas phase (reactions R1-R2), whereas all OH from the afterglow plasma jet is rapidly depleted.
Finally, to validate the model in regard to the uptake of chemical species into the liquid, figure 6 compares average hydrogen peroxide (H 2 O 2 ) concentrations in the water with the experimental results by Winter et al [17] at three watertreatment times, namely 40, 60, and 80 s, and for different inlet H 2 O 2 concentrations from 1 to 5 ppm with 1 ppm steps.Similar to the simulations in section 4.1, the water volume was 17.6 ml.Liquid fraction and velocity distributions are obtained from simulations of the impingement of an argon gas jet with flow rate Q = 3 slm onto water used as initial conditions for the species transport simulations.Given the results in figure 5, only hydrogen peroxide transport is modeled and without plasma-induced chemical kinetics given that its concentration within the water is relatively independent of the reactions considered.
The comparison of the modeling results with the experiments for the three times 40, 60, and 80 s in figure 6 shows that the simulation results have the same order of magnitude and are within the same range reported in [17].Overall, the obtained results show that the multiphase species transport model in section 2 is effective at describing the uptake of plasma-produced reactive species in water.

Conclusions
The interaction of an argon low-temperature APPJ with water was computationally investigated using time-dependent 3D simulations comprising turbulent gas flow and induced liquid movement, gas-water interface dynamics, multiphase species transport, and gas-and liquid-phase chemical reactions.A single-field approach based on the VoF method coupled with the CVA model is used to consistently describe the dynamics of the interface together with interfacial reactive mass transfer.Three CVA models for multiphase species transport have been implemented and evaluated.These are based on arithmetic, harmonic, and unified mixture diffusion.Validation studies for each component of the multiphase species transport model were performed and demonstrated good agreement with the reference solutions.
The model was applied to simulations of the interaction between a kINPen ® -generated argon APPJ with water 10 mm away from the device and contained in a cylindrical container, a set-up studied in [17,26].The developed model was first used to obtain turbulent multiphase fluid flow profiles and interface deformations for seven inlet flow rates ranging from 1.1 to 1.7 slm.The extent of the resulting cavity depth increased from 1.3 mm to 3.4 mm, with the inlet flow rate growing from 1.1 to 1.7 slm, and was consistent with the experimental values reported in [17].A set of six species were considered in the model, namely OH, H 2 O 2 , NO, NO 2 , N 2 O 4 , and HNO 2 , representative of those in the argon APPJ afterglow from the kINPen ® .A set of seven chemical reactions, spanning the gas-and liquid-phases, were used to describe the reactivity of the plasma-water multiphase system.
First, the same simulation set-up and conditions, but for an inlet gas flow rate of 3 slm, were used to investigate the multiphase transport of reactive species.Subsequently, simulations with and without chemical reactions for 1.0 s were carried out to discern the role of species reactivity on their within-water concentration.The modeling results showed that species injected at the inlet followed the fluid flow in gas, dissolved in the liquid, and were further advected in it by recirculation patterns.The species' uptake to the water varied depending on the Henry's constant, and gas and liquid phase specie's diffusion coefficients.The comparison of reactive and nonreactive cases showed that the effect of the reactions was significant.The interphase species transport model was further validated against experimental results of the volumeaverage concentration of hydrogen peroxide in water for different treatment times and inlet concentrations for an inlet gas flow of 3 slm.Modeling results showed reasonable agreement with the experimental findings.The obtained results show that the multiphase species transport model is effective at describing the uptake of plasma-produced reactive species in water.
The results of this study indicate that the VoF-CVA approach is suitable for the modeling of general multiphase systems, including those with complex plasma-liquid interfaces such as plasma-within-liquid or plasma interacting with bubbles or mist.

Appendix. Model validation
Three benchmark cases are used to validate the main aspects of the interphase species transport model, namely the description of interface dynamics, multiphase species transport, and multiphase chemical kinetics.These are the standing gravity wave problem, the transfer of species across the interface of two immiscible fluids, and a chemically-reacting three-species system, respectively.

A.1. Standing gravity wave
The model's capabilities to capture the dynamics of twophase flows are assessed with the standing gravity wave problem from Deshpande et al [33].The model is illustrated in figure A1(a) and consists of a fixed amount of water in a closed rectangular box.The water is subject to an initial perturbation that initiates periodic fluctuations of the air-water interface.
The problem is specified by defining as boundary conditions (along each side of the box) zero gradient for pressure p and liquid fraction α, and a slip boundary condition for velocity to suppress vorticity generation near the walls.The effects of viscosity and surface tensions are neglected.Liquid and gas phase densities are set to 1000 and 1 kg m −3 , respectively; both kinematic viscosities are set to 0 m 2 s −1 (i.e.inviscid fluids), and acceleration due to gravity is set to 10 m s −2 .The closed rectangular box is set 1.0 m wide and 1.5 m height.The spatial domain is discretized with a uniform mesh with 100 elements in each direction to capture the interface movement correctly.The initial perturbation of the interface is given by: where d is the average depth of the container, L is the container's width, and A is the wave amplitude.Their values are equal to 1.0, 1.0, and 0.05 m, respectively.An analytical expression for the frequency of the wave ω λ , where λ is the wavelength of the standing wave, can be obtained from the linearized potential flow equations by assuming a small amplitude of the oscillations [55], leading to: Given the problem conditions, the above expression leads to ω λ ≈ 7.92 s −1 .
To validate the computational multiphase species transport model, computational results of the relative height of waves y/d at x = L/2 attained with the multiphase species transport model described above were compared against analytical results and the numerical results reported in Deshpande et al [33].Comparison results are presented in figure A1(b).The analytical solution is given by: where is ω = 2π /L.The results in figure A1(b) show generally good agreement with the analytical frequency and the reference results by Deshpande et al.The small undershoots of the wave amplitude were also reported in the results in [33,56].They were attributed to initial nonlinearities of the waveform and limited numerical accuracy and deemed of secondary relevance.

A.2. Interface species transport
Two benchmark cases are used to validate the interface species transport model, namely a steady-state case and a transient case, both with known analytical solutions.

Steady-state
This benchmark test deals with species mass transfer across the planar interface of two immiscible static fluids inside an enclosed domain.The set-up is identical to the one presented in Haroun's work [28].The problem is inherently onedimensional (1D), depicting species distribution along the direction perpendicular to the interface.Nevertheless, the problem is solved in a 2D domain for validation purposes.The problem set-up is schematically presented in figure A2(a).
The spatial domain is a square with side L = 0.1 m and height equal to 2d, with the interface fixed horizontally at half the height, d.The problem is solved for Henry's constant values of 0.1, 1.0, and 10.The density and viscosity of the liquid phase are equal to 10 3 kg m −3 and 10 −6 m 2 s −1 , respectively, while the values for the gas phase are set as 1 kg m −3 and 1.48 10 −5 m 2 s −1 , respectively.The ratio of liquid and gas phase diffusivities is equal to 0.1.No acceleration due to gravity is included, and the initial concentrations in the liquid and the gas are set equal to 0 and 1 mol m −3 , respectively.
The analytical solution is given by: In the above expressions C 0 l , C 0 g , H, D l , and D g are the initial concentrations in the liquid and gas phases, Henry's constant, and diffusivities in the liquid and the gas, respectively.
Simulation results are presented in figures A2(b)-(d) for the three different values of Henry's constant.Comparison of the numerical results with the analytical solution shows good agreement for all three considered species transport models, namely the models based on arithmetic mean, harmonic mean, and the unified diffusion model.

Transient
The transient multiphase species transport case corresponds to the benchmark test described in Marschall et al [29].The model is presented in figure A3(a) and consists of a square domain with sides 2d = 10 cm initially at rest, with its lower half filled with water.Even though the problem is essentially 1D, i.e. the distribution of species only varies along the direction perpendicular to the interface, the computational model is solved in a 3D domain.
The initial concentration of species in the liquid (water) is C l (t = 0,x) = C l 0 = 1, while the gas phase (air) is devoid of species.The densities of water and air are equal to 998.2 and 1.122 kg m −3 , respectively; dynamic viscosities are set to 1.0 and 18.24 10 −3 mPa s, respectively; and diffusion coefficients are 0.2 and 1.0 cm 2 s −1, respectively.The Henry's coefficient is set to H = 3.As boundary conditions, zero-gradient is applied over all boundaries.The spatial domain is discretized with a mesh of 20 × 256 × 20 cells along the x, y, and z directions, respectively (the higher number of cells corresponds to discretization along the direction normal to the interface).
Simulation results are presented in figures A3(b)-(d) for three times.As can be observed, the numerical concentration profiles obtained with the three models agree with the analytical solution in all cases.Therefore, it is concluded that all three transport models can model transient multiphase species transport and can be used to investigate plasma-liquid interaction phenomena further.

A.3. Chemical kinetics in a three-species system
This benchmark case is used to assess the developed solver's capability to simulate multiphase species transport with chemical reactions.The problem is schematically depicted in figure A4(a).It consists of the temporal evolution of the concentrations of three species, denoted as A, B, and C, that react irreversibly according to the reactions: The computational domain consists of a cube with a side of 1 cm that is at rest and filled with water.The domain initially has only species A with a concentration of 0.1 mole m −3 and zero concentration of species B and C. The domain is discretized with 20 × 20 × 20 cells.To validate the multiphase species transport solver, the numerical simulation results are compared against those obtained with a global solver (0D) based on Matlab ® 's ode15s stiff ordinary differential equation

Figure 1 .
Figure 1.Model of an atmospheric pressure plasma jet (APPJ) impinging on water.Schematic depiction of an argon APPJ impinging onto liquid water within a container.The model encompasses gas jet dynamics, interface deformation, induced liquid transport, interface reactive species transport, and chemical kinetics in the gas and liquid phases.

Figure 2 .
Figure 2. Plasma jet impinging on water model set-up.(a) Computational domain depicting the three-dimensional discretization mesh, the gas and liquid phases, and the main geometric parameters.(b) Cross-sectional view of the spatial domain discretization using hexahedral finite volume cells and definition of domain boundaries, together with a representative depiction of the deformed gas-water interphase.

Table 1 .
Chemical reactions considered in the model.Set of gas-and liquid-phase reactions and their rate constants.The simulations consider a uniform temperature T = 300 K throughout both phases.Reaction k f,r (s −1 , m 3 mol −1 s −1 , or m 6 mol −2 s −1 ) k f,r at T = 300 K References Gas R1 2NO 2 → N 2 O 4 6.02 10 −4 (300/T)

Figure 3 .
Figure 3. Velocity magnitude and streamlines in the gas and liquid phases.Results for inlet gas flow rate Q equal to (a) 1.1 slm, (b) 1.3 slm, (c) 1.5 slm, and (d) 1.7 slm.The arrows indicate the center of induced recirculation regions within the gas and liquid domains.

Figure 4 .
Figure 4. Water dimple depth.(a) Numerical results of interface dimple depth as a function of gas inlet flow rate compared against the experimental results from Winter et al [17].(b) Depiction of the depth of the interface dimple for representative flow rates.

Figure 5 .
Figure 5.Chemical species concentration distributions.Concentrations of the six reactive species considered in the model throughout the domain for t = 1.0 s with (reactive) and without (nonreactive) chemical reactions for an inlet flow rate Q = 3 slm.

Figure 6 .
Figure 6.Average H 2 O 2 concentration in water as a function of inlet concentration and treatment time.Comparison of modeling results against the experimental results by Winter et al[17] for an inlet gas flow rate of 3 slm.

Figure A1 .
Figure A1.Standing gravity wave problem.(a) Schematic depiction of the problem and its main parameters, including boundary conditions and the initial deformation of the interface.(b) Comparison of results of the relative height of waves y/d obtained with the computational multiphase species transport model in this work against the analytical solution and the numerical results reported in[33].

Figure A2 .
Figure A2.Multiphase species transport benchmark-steady-state.(a) Problem set-up, including domain geometry, the extent of phases, and boundary conditions.Comparison of the analytical solution against the three multiphase transport models for Henry's constant equal to (b) 1.0, (c) 0.1, and (d) 10.
parameters k i are the chemical rate constants, set with values: k 1 = 0.1, and k 2 = k 3 = 0.01.

Figure A3 .
Figure A3.Multiphase species transport benchmark-transient.(a) Problem set-up, including domain geometry, the extent of phases, and boundary conditions.Comparison of the analytical solution against the solutions attained with three multiphase transport models at t equal to (b) 0.05 s, (c) 0.2 s, and (d) 0.8 s.

Figure A4 .
Figure A4.Chemical kinetics in a three-species system test.(a) Schematic depiction of the problem set-up, including initial conditions.(b) Comparison of numerical and reference results obtained with a global (zero-dimensional) solver of the temporal evolution of species concentrations.

Table 2 .
Boundary conditions used in the simulations.The subscript n denotes the outer normal to the corresponding boundary.

Table 3 .
Inlet concentrations of RONS included in the model.

Table 4 .
Henry constants and diffusion coefficients.Properties evaluated at 300 K and 1 atm.