Bernaise: A flexible framework for simulating two-phase electrohydrodynamic flows in complex domains

Bernaise (Binary ElectRohydrodyNAmIc SolvEr) is a flexible high-level finite element solver of two-phase electrohydrodynamic flow in complex geometries. Two-phase flow with electrolytes is relevant across a broad range of systems and scales, from 'lab-on-a-chip' devices for medical diagnostics to enhanced oil recovery at the reservoir scale. For the strongly coupled multi-physics problem, we employ a recently developed thermodynamically consistent model which combines a generalized Nernst-Planck equation for ion transport, the Poisson equation for electrostatics, the Cahn-Hilliard equation for the phase field (describing the interface separating the phases), and the Navier-Stokes equations for fluid flow. As an efficient alternative to solving the coupled system of partial differential equations in a monolithic manner, we present a linear, decoupled numerical scheme which sequentially solves the three sets of equations. The scheme is validated by comparison to limiting cases where analytical solutions are available, benchmark cases, and by the method of manufactured solution. The solver operates on unstructured meshes and is therefore well suited to handle arbitrarily shaped domains and problem set-ups where, e.g., very different resolutions are required in different parts of the domain. Bernaise is implemented in Python via the FEniCS framework, which effectively utilizes MPI and domain decomposition, and should therefore be suitable for large-scale/high-performance computing. Further, new solvers and problem set-ups can be specified and added with ease to the Bernaise framework by experienced Python users.


I. INTRODUCTION
Two-phase flow with electrolytes is encountered in many natural and industrial settings.
Although Lippmann already in the 19th century [1,2] made the observation that an applied electric field changes the wetting behaviour of electrolyte solutions, the phenomenon of electrowetting has remained elusive. Recent decades have seen an increased theoretical and experimental interest in understanding the basic mechanisms of electrokinetic or electrohydrodynamic flow [3,4]. Progress in micro-and nanofluidics [5,6] has enabled the use electrowetting to control small amounts of fluid with very high precision (see e.g. the comprehensive reviews by Mugele and coworkers [2,7] and Nelson and Kim [8] and references therein). This yields potential applications in, e.g., "lab-on-chip" biomedical devices or microelectromechanical systems [9][10][11], membranes for harnessing blue energy [12], energy storage in fluid capacitors, and electronic displays [13][14][15][16].
It is known that electrohydrodynamic phenomena affects transport properties and energy dissipation in geological systems, as a fluid moving in a fluid-saturated porous medium sets up an electric field that counteracts the fluid motion [17][18][19]. Electrowetting may also be an important factor in enhanced oil recovery [20,21]. Here, the injection of water of a particular salinity, or "smart water" [22], is known to increase the recovery of oil from reservoirs as compared to brine [23]. Further, transport in sub-micrometer scale pores in low-permeability rocks in the Earth's crust may be driven by gradients in the electrochemical potential [24], which may have consequences for, e.g., transport of methane-water mixtures in dense rocks.
Hence, a deepened understanding of electrowetting and two-phase electrohydrodynamics would be of both geological and technological importance. While wetting phenomena (or more generally, two-phase flow) on one hand, and electrohydrodynamics on the other, remain in themselves two mature and active areas of research which both encompass a remarkably rich set of phenomena, this article is concerned with the interface between these fields. For interested readers, there are several reviews available regarding wetting phenomena [25][26][27] and electrohydrodynamics [28][29][30]. Notably, the "leaky dielectric" model originally proposed by Taylor [31] (and revisited by Melcher and Taylor [28]) to describe drop deformation, is arguably the most popular description of electrohydrodynamics, but it does not describe ionic transport and considers all dielectrics to be weak conductors. In this work, we shall employ a model that does not make such simplifications. Recently, Schnitzer and Yariv [32] showed rigorously that models of the latter type reduce to the Taylor-Melcher model in the double limit of small Debye length and strong electric fields. The simplified model may therefore have advantages in settings where those assumptions are justified, e.g., in simulations on larger scales; while the class of models considered here are more general and expected to be valid down to the smallest scale where the continuum hypothesis still holds.
Experimental and theoretical approaches [33][34][35] in two-phase electrohydrodynamic flows need to be supplemented with good numerical simulation tools. This is a challenging task, however: the two phases have different densities, viscosities and permittivities, the ions have different diffusivities and solubilites in the two phases, and moreover, the interface between the phases must be described in a consistent manner. Hence, much due to the complex physics involved, simulation of two-phase electrohydrodynamic phenomena with ionic transport is still in its infancy. It has been carried out with success e.g. in order to understand deformation of droplets due to electric fields [36][37][38], or for the purpose of controlling microfluidic devices (see e.g. [39]). Lu et al. [40] simulated and performed experiments on droplet dynamics in a Hele-Shaw cell. Notably, Walker et al. [41] simulated electrowetting with contact line pinning, and compared to experiments. In practical applications, such as in environmental remediation or oil recovery, the complex pore geometry is essential and it is therefore of interest to simulate and study electrowetting in such configurations. However, to our knowledge, there have been few numerical studies of these phenomena in the context of more complex geometries.
In this article, we introduce and describe Bernaise (Binary Electrohydrodynamic Solver), which is an open-source software/framework for simulating two-phase electrohydrodynamics. It is suitable for use in complex domains, operating on arbitrary unstructured meshes. The finite-element solver is written entirely in Python and built on top of the FEniCS framework [42], which (among other things) effectively uses the PETSc backend for scalability. FEniCS has in recent years found success in related applications, such as in high-performance simulation of turbulent flow [43], and for single-phase, steady-state electrohydrodynamic flow simulation in nanopores [44] and model fractures [45]. Since Bernaise was inspired by the Oasis solver for fluid flow [43], it is similar to the latter in both implementation and use.
In this work, we employ a phase-field model to propagate the interface between the two phases. Such diffuse interface models, as opposed to e.g. sharp interface models (see for instance [46]), assume that the fluid-fluid interface has a finite size, and have the advantage that no explicit tracking of the interface is necessary. Hence, using a phase-field model has several advantages in our setting: it takes on a natural formulation using the finite element method; in sub-micrometer scale applications, the diffuse interface and finite interface thickness present in these models might correspond to the physical interface thickness (typically nanometer scale [47]); and the diffuse interface may resolve the moving contact line conundrum [27,48]. Note that although ab initio and molecular dynamics simulation methods are in rapid growth due to the increase in computational power, and do not require explicit tracking of the interface or phenomenological boundary conditions, such methods are restricted to significantly smaller scales than continuum models are. Nevertheless, they serve as valuable tools for calibration of the continuum methods [48][49][50][51]. We note also that sharp-interface methods such as level-set [52,53] and volume-of-fluid methods [38,54,55] are viable options for simulating electrohydrodynamics, but such methods shall not be considered here.
The use of phase field models to describe multiphase flow has a long history in fluid mechanics [56]. Notably, the "Model H" of Hohenberg and Halperin [57], for two incompressible fluids with matched densities and viscosities, is based on the coupled Navier-Stokes-Cahn-Hilliard system, and was introduced to describe phase transitions of binary fluids or single-phase fluid near the critical point. Lowengrub and Truskinovsky [58] later derived a thermodynamically consistent generalization of Model H where densities and viscosities were different in the two phases, however with the numerical difficulty that the velocity field was not divergence free. To circumvent this issue, Abels et al. [59] developed a thermodynamically consistent and frame invariant phase-field model for two-phase flow, where the velocity field was divergence free, allowing for the use of more efficient numerical methods.
Lu et al. [40] proposed a phase-field model to describe electrohydrodynamics, but was restricted to flow in Hele-Shaw cells, using a Darcy equation to describe the flow between the parallel plates [60]. A phase-field approach to the leaky-dielectric model was presented by Lin et al. [61]. Using the Onsager variational principle, Campillo-Funollet et al. [62] augmented the model of Abels et al. [59] with electrodynamics, i.e. inclusion of ions, electric fields and forces. This can be seen as a more physically sound version of the model proposed by Eck et al. [63], which only contained a single "net charge" electrolyte species. A model for two-phase electrohydrodynamics was derived, with emphasis on contact line pinning, by Nochetto et al. [64], but this does not appear to be frame-invariant, as the chemical potential depends quadratically on velocity [62]. In this work, we will therefore focus on the model by Campillo-Funollet et al. [62].
There is a vast literature on the discretization and simulations of immiscible two-phase flows including phase-field models (see e.g. [46,56]), but here we focus on research which is immediately relevant concerning the discretization and implementation of the model by Campillo-Funollet et al. [62]. Grün and Klingbeil [65] discretized the model in Ref. [59] (without electrohydrodynamics) with a dual mesh formulation, using a finite volume method on the dual mesh for advection terms, and a finite element method for the rest. Based on the sharp-interface model benchmarks of Hysing et al. [66], Aland and Voigt [67] provided benchmarks of bubble dynamics comparing several formulations of phase-field models (without electrodynamics). Energy-stable numerical schemes for the same case were presented and analyzed in [68,69]. Campillo-Funollet et al. [62] provided preliminary simulations of the two-phase electrohydrodynamics model in their paper, however with a simplified formulation of the chemical potential of the solutes. A scheme for the model in [62] which decouples the Navier-Stokes equations from the Cahn-Hilliard-Poisson-Nernst-Planck problem, was presented and demonstrated by Metzger [70,71]. In the particular case of equal phasic permittivities, the Cahn-Hilliard problem could be decoupled from the Poisson-Nernst-Planck problem. Recently, a stable finite element approximation of two-phase EHD, with the simplifying assumptions of Stokes flow and no electrolytes, was proposed by Nürnberg and Tucker [72].
The main contributions of this article is to give a straightforward description of Bernaise, including the necessary background theory, an overview of the implementation, and a demonstration of its ease of use. Solving the coupled set of equations in a monolithic manner (as is done in Ref. [62] using their in-house EconDrop software) is a computationally expensive task, and we therefore propose a new linear splitting scheme which sequentially solves the phase-field, chemical transport and the fluid flow subproblems at each time step. We demonstrate the validity of the approach and numerical convergence of the proposed scheme by comparing to limiting cases where analytical solutions are available, benchmark solutions, and using the method of manufactured solution. We demonstrate how the framework can be extended by supplying user-specified problems and solvers. We believe that due to its flexibility, scalability and open-source licensing, this framework has advantages over soft-ware which to our knowledge may have some of the same functionality, such as EconDrop (in-house code of Grün and co-workers) and Comsol (proprietary). Compared to sharpinterface methods, the method employed in the current framework is automatically capable of handling topological changes and contact line motion, and the full three-dimensional (3D) capabilities allows to study more general phenomena than what can be achieved by axisymmetric formulations [38]. We expect Bernaise to be a valuable tool that may facilitate the development of microfluidic devices, as well as a deepened understanding of electrohydrodynamic phenomena in many natural or industrial settings.
The outline of this paper is as follows. In Sec. II, we introduce the sharp-interface equations describing two-phase electrohydrodynamics; then we present the thermodynamically consistent model of electrohydrodynamics by Campillo-Funollet et al. [62]. In Sec. III, we write down the variational form of the model, present the monolithic scheme, and present a linear splitting scheme for solving the full-fledged two-phase electrohydrodynamics. Sec. IV gives a brief presentation of Bernaise, and demonstrates its ease use through a minimal example. Further, we describe how Bernaise can be extended with user-specified problems and solvers. In Sec. V, we validate the approach as described in the preceding paragraph.
Finally, in Sec. VI, we apply the framework to a geologically relevant setting where dynamic electrowetting effects enter, and present full 3D simulations of droplet coalescence and breakup. Finally, in Sec. VII we draw conclusions and point to future work.
We expect the reader to have a basic familiarity with the finite element method, the Python language, and the FEniCS package. Otherwise, we refer to the tutorial by Langtangen and Logg [73].

II. MODEL
The governing equations of two-phase electrohydrodynamics can be summarized as the coupled system of two-phase flow, chemical transport (diffusion and migration), and electrostatics [62]. We will now describe the sharp-interface equations that the phase-field model should reproduce, and subsequently the phase-field model for electrohydrodynamics. For the purpose of keeping the notation short, we consider a general electrokinetic scaling of the equations. The relations between the dimensionless quantities and their physical quantities are elaborated in Appendix A.

A. Sharp-interface equations
In the following, we present each equation of the physical (sharp-interface) model. With validity down to the nanometer scale, the fluid flow is described by the incompressible Navier-Stokes equations, augmented by some additional force terms due to electrochemistry: Here, ρ i is the density of phase i, v is the velocity field, µ i is the dynamic viscosity of phase i, p(x, t) is the pressure field [74], c j (x, t) is the concentration of solute species j, and g c j is the associated electrochemical potential. The form of the right hand side of Eq. (1) is somewhat unconventional (and relies on a specific interpretation of the pressure), but has numerical advantages over other formulations as it avoids, e.g., pressure build-up in the electrical double layers [75].
The transport of the concentration field of species i is governed by the conservative (advection-diffusion-migration) equation: where K ij is the diffusivity of species j in phase i. The electrochemical potential is in general given by where α (c) = ∂α/∂c(c), and α(c) is a convex function describing the chemical free energy, β ij is a parameter describing the solubility of species j in phase i, z j is the charge if solute species j, and V is the electric potential. Eq. (3) With this choice of α, the solubility parameter β ij can be interpreted as related to a reference concentration c ref,i j , through the relation This gives a chemical energy j (see also [76]).
Since the dynamics of the electric field is much faster than that of charge transport, we can safely assume electrostatic conditions (i.e., neglect magnetic fields). This amounts to solving the Poisson problem (Gauss' law): Here, ε i is the electrical permittivity of phase i, and ρ e = j z j c j is the total charge density.
In the absence of advection, for the case of two symmetric charges, and under certain boundary conditions, Eqs. (3)-(7) lead to the simpler Poisson-Boltzmann equation (see Appendix B).

Fluid-fluid interface conditions
It is necessary to define jump conditions over the interface between the two fluids. We denote the jump in a physical quantity χ across the interface by [χ] + − , and the unit vector n int normal to the interface.
Firstly, due to incompressibility, the velocity field must be continuous: The electrochemical potential must be continuous across the interface, Due to conservation of the electrolytes, the flux of ion species j into the interface must equal the flux out of the interface, and the normal flux of the electric displacement field D = − i ∇V , and the electric potential, should be continuous (since by assumption, no free charge is located betweeen the fluids): Finally, interfacial stress balance yields the condition where σ is the surface tension, κ is the curvature, and E = −∇V is the electric field.
Moreover, we have defined the shorthand symmetric (vector) gradient, Further, all gradient terms have been absorbed into the pressure. Note that Eq. (12) leads to a modified Young-Laplace law in equilibrium, which include Maxwell stresses.

Boundary conditions
There are a range of applicable boundary conditions for two-phase electrohydrodynamics.
Here, we briefly discuss a few viable options. In the following, we letn be a unit normal vector pointing out of the domain, andt be a tangent vector to the boundary.
For the velocity, it is customary to use the no-slip condition u = 0 at the solid boundary.
Alternatively, the Navier slip condition, which may be of use for modelling moving contact lines [50], could be used: where γ is a slip parameter, and the subscripts denote tangential (t) or normal (n) components. The slip length µ/γ is typically of nanometer scale and dependent on the materials in question. However, since the implementation of such conditions may become slightly involved, we omit it in the following.
With regards to the electrolytes, it is natural to specify either a prescribed concentration at the boundary, c i = c 0 , or a no-flux condition out of the domain, For the electric potential, it is natural to prescribe either the Dirichlet condition V =V , or a prescribed surface charge σ e (x),n

B. Phase-field formulation
In order to track the interface between the phases, we introduce an order parameter field φ which attains the values ±1 respectively in the two phases, and interpolates between the two across a diffuse interface of thickness . In the sharp-interface limit → 0, the equations should reproduce the correct physics, and reduce to the model above, including the interface conditions. A thermodynamically consistent phase-field model which reduces to this formulation was proposed in Ref. [62]: Here, φ is the phase field, and it takes the value φ = −1 in phase i = 1, and the value φ = 1 in phase i = 2. Eq. (19) governs the conservative evolution of the phase field, wherein the diffusion term is controlled by the phase field mobility M (φ). Here, ρ, µ, ε, K j depend on which phase they are in, and are considered slave variables of the phase field φ. Across the interface these quantities interpolate between the values in the two phases: These averages are all weighted arithmetically, although other options are available. For example, Tomar et al. [54] found that, in the case of a level-set method with smoothly interpolated phase properties, using a weighted harmonic mean gave more accurate computation of the electric field. However, López-Herrera et al. [55] found no indication that the harmonic mean was superior when free charges were present, and hence we adopt for simplicity and computational performance the arithmetic mean, although it remains unsettled which mean would yield the most accurate result.
Further, the chemical potential of species c j is given by where we, for dilute solutions, may model α(c) = c(log c − 1) to obtain consistency with the standard Nernst-Planck equation. Further, we use a weighted arithmetic mean for the solubility parameters β j : which, under the assumption of dilute solutions and with the interpretation (6), corresponds to a weighted geometric mean for the reference concentrations: In analogy with g c j being the chemical potential of species c j , we denote g φ as the chemical potential of the phase field φ. It is given by: The free energy functional f of the phase field is defined by where where σ is the surface tension, is the interface thickness, and W (φ) is a double well potential. Here, we use W (φ) = (1 − φ 2 ) 2 /4. With this free energy, we obtain We will assume this form throughout.

Phase field mobility
Given a proper definition of the phase-field mobility M (φ), the phase-field model should reduce to the sharp-interface model given in the previous section. As discussed at length in Ref. [62], the two following ways are viable options: Here M 0 is a constant, and (·) + = max(·, 0). Other formulations of M are possible; some of these will in the limit of vanishing interface width reduce to a sharp-interface model where the interface velocity does not equal the fluid velocity [59,62].

Boundary conditions
Some of the interface conditions from the sharp-interface model carry over to the phase field model, but in addition, some new conditions must be specified for the phase field. Here we give a brief summary. We assume that the boundary of the domain Ω, ∂Ω, can be divided into an inlet part ∂Ω in , an outlet part ∂Ω out , and a wall part ∂Ω wall . We shall primarily discuss the latter here.
For the velocity field, we assume the no-slip condition v(x, t) = 0 for x ∈ ∂Ω wall .
Alternatively, a no-flux condition and a slip law could have been used; in particular, a generalized Navier boundary condition (GNBC) has been shown to hold yield a consistent description of the contact line motion [48,49]. However, to limit the scope, the moving contact line paradox will in this work be overcome by interface diffusion.
With regards to the flow problem, the pressure gauge needs to be fixed. To this end, the pressure could be fixed somewhere on the boundary, or the pressure nullspace could be removed.
For the concentrations c j , we may use a prescribed concentration, or the no-flux condition For the electric potential, we use either the Dirichlet condition V =V (which is reasonable at either inlet or outlet), or in the presence of charged (or neutral) boundaries, the condition similar to the sharp-interface condition. Note that σ e (x) is prescribed and can vary over the boundary.
We assume that the no-flux conditons hold on the phase field chemical potential, For the phase field itself, a general dynamic wetting boundary condition can be expressed as [77]: where θ e is the equilibrium contact angle, τ w is a relaxation parameter, and f w (φ) = (2 + 3φ − φ 3 )/4 interpolates smoothly between 0 (at φ = −1) and 1 (at φ = 1). In this work, we limit ourselves to studying fixed contact angles, i.e. considering Eq. (38) with τ w = 0. For a GNBC, the phase-field boundary condition (38) must be modelled consistently with the slip condition on the velocity [48].

III. DISCRETIZATION
For solving the equations of two-phase EHD, i.e. the model consisting of Eqs. (17)- (21), there are four operations that must be performed: 1. Propagate the phase field φ.
2. Propagate the chemical species concentrations c i .

Update the electric potential V
4. Propagate the velocity v and pressure p.
The whole system of equations could in principle be solved simultaneously using implicit Euler discretization in time and e.g. Newton's method to solve the nonlinear system. However, in order to simulate larger systems faster, it is preferable to use a splitting scheme to solve for each field sequentially. One such splitting scheme was outlined in [70], based on the energy-stable scheme without electrochemistry as developed by [68,69]. However, that scheme did not take into account that the electric permittivities in the two fluids may differ, and when they do, the phase field and the electrochemistry computations become coupled through the electric field [71]. We will here discuss two strategies for solving the coupled problem of two-phase electrohydrodynamics. First, we present the fully monolithic, non-linear scheme, and secondly, we propose a new, fully practical linear operator splitting scheme. As we are not aware of any splitting schemes that are second-order accurate in time for the case of unmatched densities, we shall constrain our discussion to first-order in time schemes.
In the forthcoming, we will denote the inner product of any two scalar, vector, or tensor fields A, B by (A, B). Further, we consider a discrete time step τ , and denote the (first-order) discrete time derivative by The equations are discretized on the domain Ω ⊂ R d , d = 2, 3, with the no-slip boundary Γ.
Since we do not consider explicitly in-and outlet boundary conditions in this work, we will omit this possible part of the domain for the sake of brevity.
We define the following finite element subspaces:

A. Monolithic scheme
Here we give the fully implicit scheme that follows from a naïve implicit Euler discretizion of the model (17)- (21), and supplemented by Eq. (31).
is given. The scheme can then be summarized by the following.
for all test functions (u, q, ψ, and the shorthands Note that Eqs. (46) constitute a fully coupled non-linear system and the equations must thus be solved simultaneously, preferably using a Newton method. This results in a large system matrix which must be assembled and solved iteratively, and for which there are in general no suitable preconditioners available. On the other hand, the scheme is fully implicit and hence expected to be fairly robust with regards to e.g. time step size. There are in general several options for constructing the linearized variational form to be used in a Newton scheme.

B. A linear splitting scheme
Now, we introduce a linear operator splitting scheme. This scheme splits between the processes of phase-field transport, chemical transport under an electric field, and hydrodynamic flow, such that the equations governing each of these processes are solved separately.
We have also used the discretization of Eq. 38 where we have used the linearization HereJ k c j is a linear approximation of the diffusive chemical flux J c j = K j (φ)c j ∇g c j . For conciseness, we here constrain our analysis to ideal chemical solutions, i.e. we assume a common chemical energy function on the form α(c) = c(ln c − 1). To this end, we approximate the flux by: for all test functions (u, q) ∈ V h × P h . Here, we have used the following approximation of the advective momentum:m Note that the terms in (54a) involving ∂ − τ ρ k + ∇ ·m k−1 , which is a discrete approximation of ∂ t ρ + ∇ · m = 0, is included to satisfy a discrete energy dissipation law [78] (i.e., to improve stability). This step requires solving for the velocity and pressure in a coupled manner. This has the advantage that it yields accurate computation of the pressure, but the drawback that it is computationally challenging to precondition and solve, related to the Babuska-Brezzi (BB) condition (see e.g. [79]). Alternatively, it might be worthwhile to further split the fluid flow step into the following three substeps, at the cost of some lost accuracy [80].
• Pressure correction step: Find p k ∈ P h such that for all q ∈ P h , we have • Velocity correction step: which we solve by explicitly imposing the Dirichlet boundary condition u k = 0 on Γ.
We will in this paper refer to the coupled solution of the fluid flow equations, unless stated otherwise. Specifically, the fractional-step fluid flow scheme will only be demonstrated in the full 3D simulations in Sec. VI B.
The scheme presented above consists in sequentially solving three decoupled subproblems (or five decoupled subproblems for the fractional-step fluid flow alternative). The subproblems are all linear, and hence attainable for specialized linear solvers which could improve the efficiency. We note that the splitting introduces an error of order τ , i.e. the same as the scheme itself. Moreover, our scheme does not preserve the same energy dissipation law on the discrete level, that the original model does on the continuous level. We are currently not aware of any scheme for two-phase electrohydrodynamics with this property, apart from the fully implicit scheme presented in the previous section.

IV. BERNAISE
We have now introduced the governing equations and two strategies for solving them.
Now, we will introduce the Bernaise package, and describe an implementation of a generic simulation problem and a generic solver in this framework. For a complete description of the software, we refer to the online Git repository [84]. The work presented herein refers to version 1.0 of Bernaise, which is compatible with version 2017.2.0 of FEniCS [42].

A. Python package
Bernaise is designed as a Python package, and the main structure of the package is shown in Fig. 1. The package contains two main submodules, problems and solvers. As suggested by the name, the problems submodule contains scripts where problem-specific geometries (or meshes), physical parameters, boundary conditions, initial states, etc., are specified. We Further, a default solver compatible with a given problem is specified in the problem, but this setting can-along with most other settings specified in a problem-be overridden by providing an additional keyword to the main script call (see below  The role of the main module sauce.py is to allocate the required variables to run a simulation, to import routines from the specified problem and solver, to iterate the solver in time, and to output and store data at appropriate times. Hence, the main module works as a general interface to problems and solvers. This is enabled by overloading a series of functions, such that problem-and solver-specific functions are defined within the problem and solver, respectively. The structure of sauce.py is by choice similar to the NSfracStep.py script in the Oasis solver [43]; both in order to appeal to overlapping user bases, and to keep the code readable and consistent with and similar to common FEniCS examples. However, an additional layer of abstraction in e.g. setting up functions and function spaces is necessary in order to handle a flexible number of subproblems and subspaces, depending on e.g. whether phase field, electrochemistry or flow is disabled, or whether we are running with a monolithic or operator splitting scheme. To keep the Bernaise code as readable and easily maintainable as possible, we have consciously avoided uneccessary abstraction. Only the boundary conditions (found in common/bcs.py are implemented as classes.

B. The problems submodule
The basic user typically interacts with Bernaise by implementing a problem to be solved. This is accessible to Bernaise when put in the subfolder problems. The implementation consists in overloading a certain set of functions; all of which are listed in the problems/ init .py file in the problems folder. The mandatory functions that must be overloaded for each problem are: • mesh: defines the geometry. Equivalent to the mesh function in Oasis [43].
• problem: sets up all parameters to be overloaded, including defining solutes and types of finite elements. The default parameters are defined in the problems/ init .py file.
• create bcs: sets all subdomains, and defines boundary conditions (including pointwise boundary condtions, such as pressure pinning). The boundary conditions are more thoroughly explained below.
Further, there are functions that may be overloaded.
• constrained domain: set if the boundary is to be considered periodic.
• start hook: hook called before the temporal loop.
• tstep hook: hook called at each time step in the loop.
• end hook: hook called at the end of the program.
• rhs source: explicit source terms to be added to the right hand side of given fields; used e.g. in the method of manufactured solution.
Note here the use of three hooks that are called during the course of a simulation. These are useful for outputting certain quantities during a simulation, e.g. the flux through a cross section, or total charge in the domain. The start hook could also be used to call a steadystate solver to initialize the system closer to equilibrium, e.g. a solver that solves only the electrochemistry subproblem such that we do not have to resolve the very fast time scale of the initial charge equilibration.
In Listing 1, we show an implementation of the problems function, which sets the necessary parameters that are required for the charged droplet case to run. Here, the solutes array (which defines the solutes), contains only one species, but it can in principle contain arbitrarily many. # Default parameters to be loaded unless starting from checkpoint. parameters = dict( solver="basic", # Solver to be used.
folder="results_charged_droplet", # Folder to store results in. Listing 2: The initialize function for the charged droplet case.

C. The solvers submodule
Advanced users may develop solvers that can be placed in the solvers subdirectory.
In the same way as with the problems submodule, a solver implementation constists of overloading a range of functions which are defined in solvers/ init .py.
• get subproblems: Returns a dictionary (dict) of the subproblems which the solver splits the problem into. This dictionary has points to the name of the fields and the elements (specified in problem) which the subspace is made up of.
• setup: Sets up the FEniCS solvers for each subproblem.
• solve: Defines the routines for solving the finite element problems, which are called at every time step.
• update: Defines the routines for assigning updated values to fields, which are called at the end of every time step.
The module solvers/basicnewton.py implements the monolithic scheme, while the module solvers/basic.py implements the segregated solver [85]. The problem is split up into the subproblems corresponding to whether we have a monolothic or segragated solver in the function get subproblems. Within the setup function, the variational forms are defined, and the solver routines are initialized. The latter are eventually called in the solve routine at every time step. Note that the element types are defined within the problem, and that the solvers in general can be applied for higher-order spatial accuracy without further ado.
The task of get subproblems is simply to link the subproblem to the element specification.
In Listing 3, we show how the get subproblems function is implemented in the basic solver. As can be readily seen, the function formally splits the problem into the three subproblems NS, PF, and EC. The other functions (such as setup) are somewhat more involved, but can be found at the Git repository [84].
Note that the implementations of the solvers presented above are sought to be short and humanly readable, and therefore quite straightforwardly implemented. There are several ways to improve the efficiency (and hence scalability) of a solver, at the cost of lost intuitiveness [43]. -Charged: A charged boundary-a Neumann conditon intended for use with the electric potential V .
-Open: An open boundary-a Neumann condition is applied.
We note that when a no-flux condition is to be applied, no specific boundary condition class needs to be supplied, since the boundary term in the variational form then disappears (in particular when considering conservative PDEs).
As • description: routine called when a question mark is added to the end of the method name during a call from the terminal, meant to obtain a description of the method without having to inspect the code.
• method: the routine that performs the desired analysis.
The implementation hinges on the TimeSeries class (located in utilities/TimeSeries.py), which efficiently imports the XDMF/HDF5 data files and the parameter files produced by a Bernaise simulation. Several plotting routines are implemented in utilities/plot.py, and these are extensively used in various analysis methods.

V. VALIDATION
With the aim of using Bernaise for quantitative purposes, it is essential to establish that the schemes presented in the above converges to the correct solution-in two senses: • The numerical schemes should converge to the correct solution of the phase-field model.
• The solution of the phase-field model should converge to the correct sharp-interface equations [86].
Unless otherwise stated, we mean by convergence that the error in all fields χ should behave like, where · h is an L 2 norm, χ is the simulated field, χ e is the exact solution, h is the mesh size, τ is the time step, k h is the order of spatial convergence, k τ is the order of temporal convergence (k τ = 1 in this work), and C h and C τ are constants.
In the following, we present convergence test in three cases. Firstly, in the limiting case of a stable bulk intrusion without electrochemistry, an analytical solution is available to test against. Secondly, using the method of manufactured solution, convergence of the full two-phase EHD problem to an augmented Taylor-Green vortex is shown. Thirdly, we show convergence towards a highly resolved reference solution for an electrically driven charged droplet.
We note that the aim of Bernaise is to solve coupled multi-physics problems, and while the solvers may contain subtle errors, they may be negligible for many applications, and dominant only in limiting cases. In addition to testing the whole, coupled multi-physics problem of two-phase EHD, a proper testing should also consider simplified settings where fewer physical mechanisms are involved simultaneously. A brief discussion of testing and such reduced models is given in Appendix C. In this section, we show the convergence of the schemes in a few relevant cases, which we believe represent the efficacy of our approach.
Tests of simplified-physics problems are found in the GitHub repository [84].

A. Stable bulk intrusion
A case where an analytic solution is available, is the stable intrusion of one fluid into another, in the absence of electrolytes and electric fields. A schematic view of the initial set-up is shown in Fig. 2 Due to the Galilean invariance, we expect the velocity field to be uniformly equal to the inlet and outlet velocities, i.e. v(x, t) = v 0x . The exact analytical solution for the phase field is given by for which we shall consider the error norm. Note that the only parameters this analytical solution depends on are the initial position of the interface x 0 , the injection velocity v 0 , and the interface width . We consider the parameters ρ 1 = ρ 2 = 1000, µ 1 = 100, µ 2 = 1, σ = 2.45, = 0.03, M (φ) = M 0 = 2 · 10 −5 , x 0 = 1, L x = 5, L y = 1 and v 0 = 0.1. Fig. 3 shows the convergence to the analytical solution with regards to temporal resolution. The order of convergence is consistent with the order of the scheme, indicating that the scheme is appreciable at least in the lack of electrostatic interactions.  The scheme is seen to converge at the theoretical rate, ∼ h 2 .
B. Method of manufactured solution: a two-phase electrohydrodynamic Taylor-

Green vortex
Having established convergence in the practically one-dimensional case, we now consider a slightly more involved setting where we use the method of manufactured solution to obtain a quasi-analytical test case.
The Taylor-Green vortex is a standard benchmark problem in computational fluid dy- The fields are given by c ± = c 0 (1 ± cos x cos y C(t)), (59d) Here, the time-dependent coefficients are given by where U 0 , C 0 and Φ 0 are scalars, and where     In Fig. 6, we plot the L 2 errors of the same fields as in Fig. 5, but as a function of the time step τ . In the simulations plotted here, we used a fine grid resolution with h = 2π/256 to rule out the contribution of spatial discretization to the error, cf. Eq. (57). Clearly, first order convergence is achieved for sufficient refinement, for all fields including the pressure.

C. Droplet motion driven by an electric field
We now consider a charged droplet moving due to an imposed electric field; a problem for which there is no analytical solution available. However, by comparing to a highly resolved numerical solution, convergence for the fully coupled two-phase electrohydrodynamic problem can be verified. This problem has already been partly presented in the above, and is implemented in problems/charged droplet.py. A sketch showing the initial state is shown in Fig. 7. We consider an initially circular droplet, where a positive charge concentration is initiated as a Gaussian distribution, with variance δ 2 c , in the middle of the droplet. In this set-up, we consider only a single, positive species. The total amount of solute, i.e. integrated in the surrounding fluid, and the diffusivity here is very low here to prevent leakage. The droplet is accelerated by the electric field towards the right, before it is slowed down due to viscous effects upon approaching the wall.
With regards to reproducing the sharp-interface equations, we consider now the case of reducing the interface thickness → 0. To this end, we keep the ratio h/τ between mesh size and time step fixed, and further we keep the interface thickness proportional to h.
The latter spans roughly 3-4 elements. Since the interface thickness changes, an important parameter in the phase-field model changes, which couples back to the equations, and thus the L 2 norm does not necessarily constitute a proper convergence measure. We therefore resort to using the picture norm or contour of the droplet as a measure, i.e. the zero-level set of the phase field φ = 0. In particular, we will consider two observables: circumference and the center of mass (along x) of the droplet, as a function of resolution. A similar approach was taken for the case of phase-field models without electrodynamics by Aland  The 'water' droplet contains positive ions and is driven by the electric field set up between the high potential on the left wall and the grounded right wall.
[67] who compared their benchmarks to sharp interface results by Hysing et al. [66].
The resolutions used in our simulations are given in Table III. In order not to have to adjust the phase field mobility when refining, whilst still expecting to retrieve the sharpinterface model in the limit → 0, we choose the phase field mobility given by (33b). All parameters for the phasic quantities are given in Table IV, while the remaining parameters are given in Table V. From these parameters, using the unit scaling adopted in this paper, we find an approximate Debye length λ D = ε/(2z 2 c R ) 1/(2 · 10) 0.2 (see Section B 2 in the Appendix for this expression), since we can approximate the order of magnitude of c R < C/(πR 2 ) = 10/(π · 0.25 2 ) for a moderate screening.
In Fig. 8, we show the contour of the driven droplet at two time instances t = 4 and t = 8, and compare increasing resolution (simultaneously in space, time and interface thickness).
Qualitatively inspecting the contours by eye, the droplet shapes seem to converge to a well defined shape with increasing resolution at both time instances. However, qualitive comparison is clearly not enough to assess the convergence. As in Refs. [66,67], we define three observables:   The effect of the four resolutions given in Table III is shown. The legend shown in the figure refers to both spatial (h) and temporal resolution (τ ).
• Drift velocity: Similarly as above, the velocity at which the droplet is driven is mea- • Circularity: Defined as the ratio of the circumference of the area-equivalent circle to the droplet circumference, .
The circumference and the integrals are computed by the post-processing method geometry in time which is built into Bernaise.
for a given quantity q. We can compute an empirical convergence rate of this norm, for two successive resolutions (h i+1 > h i ). Here we shall consider the L 2 error norm in time, i.e. p = 2, and in practice we compute the integrals in time by cubic spline interpolation of measurement points saved at every 5 time steps. There is no exact solution, or reference high-resolution sharp-interface solution available for this set-up. However, if we now assume that the finest resolution is the exact solution, and use this as the reference field in Eq. (68), we can compute error norms and convergence rates. These values are reported in Table VI.
The computed convergence rates increase for all three observables and reach 1.6-1.7 with increasing resolution, indicating also quantitatively a convergence that is in agreement with the anticipated convergence rate. Considering Eq. (57), from the temporal discretization, we expect k 2 1, and from the spatial k 2 2. Depending on which term contributes most to the error, we will measure either of these rates. The values measured here indicate that both terms may be comparable in magnitude; however if we instead of using directly the finest solution as reference, extrapolated the trajectories further, we would presumptively  have achieved lower convergence rates. This might indicate that the convergence error is eventually dominated by the temporal discretization, cf. Eq. (57).

VI. APPLICATIONS
A. Oil extrusion from a dead-end pore Here, we present a demonstration of the method in a potential geophysical application.
We consider a shear flow of one phase ("water") over a dead-end pore which is initially filled with a second phase ("oil"). The water phase contains initially a uniform concentration of positive and negative ions, c ± | t=0 = c 0 , and the water-oil interface is modelled to be impermeable. The simulation of the dead-end pore is carried out to preliminarily assess the hypothesis that electrowetting could be responsible for the increased expelling of oil in low-salinity enhanced oil recovery. The problem set-up is schematically shown in Fig. 10. The phasic parameters used in the simulations are given in Table VII, and the remaining parameters are given in Table VIII. This problem is implemented in the file problems/snoevsen.py.
To investigate the effect of including electrostatic interactions, we show in Fig. 11   The geometry is specified by the two lengths L x , L y , and the radius R used to define the dead-end pore in the center of the channel by a circle and a circular smoothed inlet. The roman numerals indicate the phase, along with the tone of gray. The darker phase is the oil-like phase (I), and the lighter one is the water-like phase (II).
For the uncharged case, the frames that are shown are almost indistinguishable. In fact, the main difference is the numerical noise of the total charge, which is due to roundoff errors of machine precision. The initial dynamics of the oil plug interface, which is to equilibrate with the neutral contact angle and the shear flow, mainly happens before the first frame presented; compare Figs. 10 and 11(a).
A markedly different behavior is displayed in the right column, Figs. 11(b), 11(d), and 11(f), where a uniform surface charge density is enforced the walls at the simulation start, t = 0. Here, we see first that two tongues are intruding on both sides of the droplet, which push the droplet out into the center of the dead-end pore. The process is continued, as shown in the second frame, and finalized, as shown in the third frame, with the complete release of the droplet as the two tongues meet at the bottom of the dead-end pore, cutting Finally, to demonstrate the ability of Bernaise to simulate 3D configurations, we present simulations of two oppositely charged droplets that coalesce. In order to achieve this efficiently, a fully iterative solver was implemented. The solver consists of a fractional step version of the basic solver, in the sense that within the fluid flow step, it splits between the velocity and pressure computations, as shown in Eqs. (56a), (56b), and (56c). The splitting introduces a weak compressibility which suffices to stabilize the problem [80] (with respect to the BB condition) and thus we can use P 1 finite elements also for the velocity.
The combination of fewer degrees of freedom and the applicability of iterative linear solvers imparts significant speed-up compared to coupled solvers, which is of paramount importance for 3D simulations. This yields advantages over solvers which rely on a mixed-element formulation of the hydrodynamic subproblem [71]. The detailed analysis of the fractional step solver will be published in a separate paper, but the implementation can be found in solvers/fracstep.py. For solving the linear systems iteratively, we use an algebraic multigrid (AMG) preconditioner and a generalized minimal residual (GMRES) linear solver for the electrochemical and the pressure correction step; Jacobi preconditioner (Jacobi) and a stabilized bi-conjugate gradient method (BiCGStab) for the velocity prediction, and Jacobi and GMRES for the velocity correction. For the phase field we use Jacobi and a conjugate gradient method.
To prevent leakage of ions out of the two coalescing droplets, a weighted geometric mean was used for the diffusivities: The average concentration of the respective ion species within each droplet is c 0 , such that the total charge in the system is zero, and the initial spread (standard deviation) of the Gaussian distribution is R/3. A potential V 0 is set on the top plane at z = L z and the bottom plane at z = 0 is taken to be grounded. We assume no-slip and no-flux conditions on all boundaries, except for the electrostatic potential V at the top and bottom planes, and the fluid is taken to be in a quiescent state at the initial time t = 0. The phasic parameters used in the simulations are given in Table IX, and the remaining parameters are given in Table X. The problem is implemented in the file problems/charged droplets 3D.py.   we stress that the division of labour between the steps is highly problem-dependent, and in particular, the electrochemical subproblem is susceptible to how far into the non-linear regime we are (see e.g., [45]).

VII. DISCUSSION AND CONCLUSION
We have in this work presented Bernaise, a flexible open-source framework for simulating two-phase electrohydrodynamics in complex geometries using a phase-field model. The solver is written in its entirety in Python, and is built on top of the FEniCS/DOLFIN framework [42,88] for solving partial differential equations using the finite element method on unstructured meshes. FEniCS in turn interfaces to, e.g., scalable state-of-the art linear solvers through its PETSc backend [89]. We have proposed a linear operator-splitting scheme to solve the coupled non-linear equations of two-phase electrohydrodynamics. In contrast to solving the equations directly in a monolithic manner, the scheme sequentially solves non-trivial settings. Firstly, to test the applicability of the code in a complicated geometry, and to illuminate the effects of dynamic electrowetting, we simulated a shear flow of water containing an electrolyte over a dead-end pore initially filled with oil. This problem is relevant from a geophysical standpoint, and exemplifies the potential of the method to simulate the dynamics of the interaction between two-phase flow and electric double layers.
Secondly, the ability of the framework to simulate three-dimensional configurations was demonstrated using a fully iterative version of the operator-splitting scheme, by simulating the coalescence and subsequent breakup of two oppositely charged droplets in an electric field. The parallel scalability of the latter solver was tested on in-house computing facilities.
The results presented herein underpin our aim that Bernaise can become a valuable tool both within the micro-and nanofluidics community and within geophysical simulation.
There are several possible avenues for further development and use of Bernaise. With regards to computational effort, the linear operator-splitting scheme constitutes a major computational improvemnt over a corresponding monolithic scheme. For the resulting smaller and simpler subproblems, more specialized linear solvers and preconditioners can be used.
However, the implementation of the schemes are still not fully optimized, as in many cases it is not strictly necessary to reassemble entire system matrices (multiple times) at every time step. Using ideas e.g. from Ref. [43] on how to effectively preassemble system matrices in FEniCS, one could achieve an implementation that is to a larger extent dominated by the backend linear solvers. However, as the phase field is updated at every time step, there may be less to gain in performance than what was the case in the latter reference.
With regard to solving the Navier-Stokes equations, the solvers considered herein either rely on a coupled solution of the (the basic and basicnewton solvers) or a fractional step approach that splits between the computations of velocity and pressure (the fracstep solver that was considered in Sec. VI B). Using direct linear solvers, the coupled solvers yield accurate prediction of the pressure and can be expected to be more robust. However, direct solvers have numerical disadvantages when it comes to scalability, and Krylov solvers require specifically tailored preconditioners to achieve robust convergence. An avenue for further research is to refine the fracstep solver and develop decoupled energy-stable schemes for this problem, which seems possible by building on literature on similar systems [68-71, 76, 78]. The implementation of such enhanced schemes in Bernaise is straighforward, as demonstrated in this paper.
A clear enhancement of Bernaise would be adaptivity, both in time and space. Adaptivity in time should be implemented such that time step is variable and controlled by the globally largest propagation velocity (in any field), and a Courant number of choice. Adaptivity in space is presently only supported as a one-way operation. Adaptive mesh refinement is already used in the mesh initialization phase in many of the implemented problems. However, mesh coarsening has currently limited support in FEniCS and to the authors' knowledge there are no concrete plans of adding support for this. Hence, Bernaise lacks an adaptive mesh functionality, but this could be implemented in an ad hoc manner with some code restructuring.
In this article, we have not considered any direct dependence of the contact angle (i.e. the surface energies) on an applied electric field. However, the contact angle on scales below the Debye length is generally thought to be unaffected, albeit on scales larger than the insulator thickness, an apparent contact angle forms [90,91]. Using the full two-phase electrohydrodynamic model presented herein, effective contact angle dependencies upon the zeta potential could be measured and used in simulations of more macroscopic models; i.e. models admissible on scales where the electrical double layers are not fully resolved.
This would result in a modified contact angle energy that would be enforced as a boundary condition in a phase field model [92].
Physically, several extensions of the model could be included in the simulation framework.
Surfactants may influence the dynamics of droplets and interfaces, and could be included as in e.g. the model by Teigen et al. [93]. The model in its current form further assumes that we are concerned with dilute solutions (i.e., ideal gas law for the concentration), and hence more complicated electrochemistry could to some extent be incorporated into the chemical free energy α(c). for solute transport, which in physical units can be written as where k B is Boltzmann's constant, T is the temperature, and q e is the elementary charge.
The Poisson equation is in physical units given by where the net charge is given by ρ e = q e j c j . The Navier-Stokes equations are given by the usual Continuity of the normal stress across the interface between the phases can be formulated where p is a pressure which has absorbed some extra gradient terms. We introduce now dimensionless versions of all physical variables, and indicate the dimensionless versions by a tilde. Further, all reference values are marked with an asterisk. Hence, we lett = t/t * , ρ = ρ/ρ * ,ũ = u/u * ,p = p/p * ,μ = µ/µ * ,c j = c j /c * ,Ṽ = V /V * ,D ± = D ± /D * ,˜ = r / * , andσ = σ/σ * . All spatial dimensions are scaled by a reference linear size R * , such that x = x/R * . The electrostatic potential V is scaled by a thermal voltage, The other reference values are given by [91] This constitutes an invertible set of relations between the physical and dimensionless variables. In particular, adopting the dimensionless variables and subsequently dropping the tildes, results in the set of equations (A1) to (A5) with q e = k B T = 1 and 0 r → . This is essentially the scaling adopted in this paper.

Appendix B: Poisson-Boltzmann equation for two phases
Here, we derive a generalized Poisson-Boltzmann equation for the case of two phases, valid in equilibrium. We are here considering the steady state of the sharp interface equations.
Considering Eq. (3) with ∂ t = 0 and v = 0, taking the inner product of it with g c j , and integrating over the domain Ω, we obtain Ω K ij c j |∇g j | 2 dΩ = ∂Γ K ij c j g c jn · ∇g c j dΓ = 0, where the last equality holds, since at equilibrium the fluxes must vanish at the boundary (and hence also in the bulk). Since c j is positive, g c j may not vary. Hence, the electrochemical potential associated with electrolyte j must satisfy: where C j is a constant. We assume that one of the two phases is connected to a reservoir far away, such that here β ij = β R , c j = c R and V = V R . Evaluating Eq. (B2) at the reservoir, we have By defining χ(·) as the inverse function of α (·), we may combine Eqs. (B2) and (B3) and invert with respect to c j : Hence, by Eq. (7), we obtain a closed equation for V : with the above boundary conditions at the reservoir. The interface condition between the phases is [V ] + − = 0, i.e. continuity in V , and the boundary condition at the reservoir is V = V R . Next, we consider some special cases of this equation.

Standard Poisson-Boltzmann
With two symmetric electrolytes, j ∈ {±}, z ± = ±z, β ij = β i , V R = V and the ideal gas chemical potential, we have that α (c) = ln c, χ(a) = e a , and we obtain from Eq. (B5): where we have defined a phase-dependent Debye length λ i,D = ε i e −β R +β i /(2z 2 c R ). Now, Eq. (B4) yields that the concentration is retrieved by When |zV | 1, we may expand Eq. (B6) to the first order to obtain the linearized Poisson-Boltzmann equation: In principle, we can also expand Eq. (B7): so that the total charge density is given by

Net charge
Now we consider the single "net charge" model which was proposed in Ref. [63] and used in the simulations of Ref. [62]. (Note that these papers redefined the diffusivity to absorb the net charge, effectively K ij c j → K; but this does not have consequences in the forthcoming.) Here, we have only one species c 1 = ρ e with charge z, and α (c) = λc, such that χ(a) = λ −1 a.
Hence this form of the chemical potential yields an exact solution. Note however, that there is in principle no mechanism controlling the sign of ρ e , which is natural given that it should here signify a net charge.
The total charge is Q = which approaches the (negative) applied surface charge in the limit of infinite domain, L → ∞, as it should.
Appendix C: Some considerations on testing and applicability of the framework To simplify the complexity of the problem, the scheme can be reduced to describe settings where fewer physical mechanisms are present simultaneously.
• The very simplest is pure single-phase flow, containing only point 4 from the list in Sec. III.
• Slightly more demanding is single-phase flow with transport of a tracer dye (in the absence of electric charges and fields), using points 2 and 4.
• More demanding, pure two-phase flow (with unmatched densities and viscosities), where only points 1 and 4 from the list above enter.
• In the absence of electric charges and external electric potential, two-phase flow with passive transport of a tracer dye can be modelled, i.e. using the points 1, 2 and 4.
• Time-dependent single-phase EHD can be modelled using points 2, 3 and 4.
• There is also a subtle case where all concentrations c i = 0, but the electric field acts as a force on the interface of the two-phase flow only due to the jump in permittivity ε. This includes points 1, 3 and 4.
• Finally, the full-fledged two-phase EHD includes all the points in the list in Sec. III. The color scale is relative to the maximum deviation from neutral charge for an entire simulation; therefore the neutral simulations display numerical noise (which is of the order of machine precision). In the left column the surface charge is zero, and in the right column, a uniform surface charge density σ e = −10 is set. The rows show snapshots at different times t.