Comparison of refilling schemes in the free-surface lattice Boltzmann method

Simulating mobile liquid-gas interfaces with the free-surface lattice Boltzmann method (FSLBM) requires frequent re-initialization of fluid flow information in computational cells that convert from gas to liquid. The corresponding algorithm, here referred to as the refilling scheme, is crucial for the successful application of the FSLBM in terms of accuracy and numerical stability. This study compares five refilling schemes that extract information from the surrounding liquid and interface cells by averaging, extrapolating, or assuming one of the three different equilibrium states. Six numerical experiments were performed, covering a broad spectrum of possible scenarios. These include a standing gravity wave, a rectangular and cylindrical dam break, a Taylor bubble, a drop impact into liquid, and a bubbly plane Poiseuille flow. In some simulations, the averaging, extrapolation, and one equilibrium-based scheme were numerically unstable. Overall, the results have shown that the simplest equilibrium-based scheme should be preferred in terms of numerical stability, computational cost, accuracy, and ease of implementation.


I. INTRODUCTION
The free-surface lattice Boltzmann method (FSLBM) [1] is a numerical model for simulating free-surface flows combining the latticed Boltzmann method (LBM) for hydrodynamics simulations with the volume of fluid (VOF) approach [2] for interface tracking. It successfully simulates applications such as rising bubbles [3], waves [4], dam break scenarios [5], impacts of droplets [6], and electron-beam melting [7]. Free-surface flows relate to immiscible two-fluid flow problems in which the fluid dynamics of the lighter fluid can be neglected. Therefore, the problem reduces to a single-fluid flow with a free boundary. In this article, the lighter fluid will be called gas, and the heavier fluid will be called liquid. The Eulerian computational grid is represented by lattice cells in the LBM. In the FSLBM, each lattice cell is categorized as either gas, liquid, or interface type, with the latter separating the former. In the LBM, information about the flow field is stored in each cell in terms of particle distribution functions (PDFs).
Agreeing with the free-surface definition, in the FSLBM, valid PDFs are only available in liquid and interface cells but not in gas cells. Gas cells are converted to interface cells during the simulation because of the free interface's motion. These cells must be refilled with valid flow field information, that is, their PDFs must be reinitialized. To the authors' knowledge, no other refilling scheme but the one suggested in the original FSLBM from Körner et al. [1] has yet been tested for the FSLBM. However, there have been similar studies about moving solid obstacle cells in the LBM. There, analogously, cells are converted from solid to liquid and must be refilled [8][9][10][11][12][13][14]. Based on the schemes used for this application, five different schemes for refilling cells in the FSLBM are compared in this article.
The manuscript is structured as follows. First, the numerical foundations of the LBM and FSLBM are introduced. Then, the different refilling schemes are presented and discussed in terms of mass conservation and computational costs. The first scheme under investigation initializes the PDFs with their equilibrium constructed with the average fluid velocity and density of non-newly converted neighboring interface and liquid cells [1]. The second and third scheme extend the first one by adding a contribution of the non-equilibrium PDFs [8], or by including information about the local pressure tensor using Grad's moment system [10][11][12]15], respectively. The fourth refilling scheme initializes PDFs with a second-order extrapolation from neighboring cells' PDFs [9]. In contrast to these, in the final scheme tested here, the PDFs are initialized with the average corresponding PDFs from neighboring, non-newly converted interface or liquid cells [14]. Six numerical benchmarks then compare the refilling schemes in terms of accuracy and numerical stability. These benchmarks include a standing gravity wave, the collapse of a rectangular and cylindrical liquid column, the rise of a Taylor bubble, the impact of a droplet into a thin film of liquid, and a bubbly plane Poiseuille flow. Finally, it is concluded that for the FSLBM, the simplest equilibrium-based refilling scheme is preferable in terms of numerical stability, computational costs, accuracy, and ease of implementation.
The source code of the implementation used in this study is freely available as part of the open source C ++ software framework waLBerla [16] (https://www.walberla.net). The version of the source code used in this article is provided in the supplementary material.

II. NUMERICAL METHODS
This section introduces the foundations of the lattice Boltzmann method and its extension to free-surface flows, the free-surface lattice Boltzmann method. The section is based on Section 2 from prior articles [17,18] but repeated here for completeness.

A. Lattice Boltzmann method
The lattice Boltzmann method is a relatively modern approach for simulating computational fluid dynamics. This article only introduces its fundamental aspects. A rigorous introduction to the LBM is available in the literature [19].
The LBM is a discretization of the Boltzmann equation from kinetic gas theory. It describes the evolution of particle distribution functions on a uniformly discretized Cartesian lattice with spacing ∆x ∈ R + . The macroscopic fluid velocity is discretized with the DdQq velocity set in each cell of the lattice, with d ∈ N referring to the lattice's spatial dimension and q ∈ N referring to the number of PDFs per cell. A PDF f i (x, t) ∈ R with i ∈ {0, 1, . . . , q − 1} describes the probability that there exists a virtual fluid particle population at position x ∈ R d and time t ∈ R + traveling with lattice velocity c i ∈ ∆x/∆t {−1, 0, 1} d , where ∆t ∈ R + denotes the length of a discrete time step. The successive steps of collision, also called relaxation, and streaming, also called propagation, form the lattice Boltzmann equation. In the collision step, the collision operator Ω i (x, t) ∈ R relaxes the PDFs towards an equilibrium state f eq i (x, t), which is influenced by external forces F i (x, t) ∈ R. In the streaming step, the post-collision PDFs f i (x, t) propagate to neighboring cells. For the simulations in this article, the single relaxation time (SRT) collision operator was used, where τ > ∆t/2 is the relaxation time. The PDF's equilibrium [20] f eq can be derived from the Maxwell-Boltzmann distribution and includes the lattice weights w i ∈ R, the lattice speed of sound c 2 s ∈ R + , the macroscopic fluid density ρ ≡ ρ(x, t) ∈ R + , and the macroscopic fluid velocity u ≡ u(x, t) ∈ R d . In this article, the well-established D2Q9 and D3Q19 lattice models are used. The corresponding lattice weights are available in the literature [19]. The lattice speed of sound for these velocity sets is c 2 s = 1/3 ∆x/∆t. It relates the macroscopic fluid density ρ(x, t) and pressure p(x, t) = c 2 s ρ(x, t). The PDFs' zeroth-and first-order moments yield the fluid's density and velocity with external force F (x, t) ∈ R d . The fluid's kinematic viscosity is computed from the relaxation time τ , that is, relaxation rate ω = 1/τ . In the simulations used in this article, the gravitational force, as part of F i (x, t) in the LBM collision (1), was modeled according to Guo et al. [21] with where again, u ≡ u(x, t) was used.
The rectangular and cylindrical dam break simulations in Sections IV B and IV C were performed using a Smagorinksy-type large eddy simulation turbulence model [22,23]. Based on the user-chosen relaxation time τ 0 > ∆t/2, the collision operator's relaxation time is locally adjusted by the model with a contribution τ t (x, t) ∈ R from the turbulence viscosity The turbulence viscosity is obtained from the filtered strain rate tensor where the filtered mean momentum flux is computed from the momentum fluxes The index notation with α and β refers to the components of a vector or tensor. The moment fluxes are given by the second-order moments of the PDFs' non-equilibrium parts. The turbulence model's contribution to the relaxation time is then [22] τ t (x, t) = 1 2 where ∆x LES is the filter length and C S is the Smagorinsky constant. For the simulations in this article, these parameters were chosen ∆x LES = ∆x and C S = 0.1, as suggested by Yu et al. [23].
At solid obstacles, a no-slip boundary condition were realized using the bounce-back approach, where PDFs streaming into solid obstacle cells are reflected reversely. The PDF's original direction with index i is reversed, denoted asī, with lattice velocity cī = −c i [19].
Free-slip boundary conditions are modeled similarly with the PDFs being reflected specularly.
In the resulting lattice velocity c j , the normal velocity component of the incoming velocity c i is reversed with c j,n = −c i,n at a free-slip boundary [19].
In the remainder of this article, ∆x = 1 and ∆t = 1 are assumed as this is common practice in the LBM [19]. All quantities are denoted in the LBM unit system if not explicitly stated otherwise. The LBM reference density ρ 0 = 1 and pressure p 0 = c 2 s ρ 0 = 1/3 were set in all simulations. The relaxation time τ or relaxation rate ω specified for the numerical experiments in Section IV refer to the constant user-chosen values that the Smagorinsky turbulence model did not yet adjust.

B. Free-surface lattice Boltzmann method
The free-surface lattice Boltzmann method as presented by Körner et al. [1] is used in this article. The FSLBM extends the LBM by simulating the interface between two immiscible fluids. It assumes that the heavier fluid governs the entire flow dynamics of the system with the lighter fluid's influence being negligible. Consequently, the immiscible two-fluid flow problem reduces to a single-fluid flow with a free boundary. Therefore, the hydrodynamics of the lighter fluid are not simulated in the FSLBM. A simplification such as this is valid if the fluids' densities and viscosities differ substantially, for example as in liquid-gas flows. In what follows, the heavier fluid is called liquid, whereas the lighter fluid is called gas.
The free interface between the liquid and gas is treated as in the VOF approach [2]. A fill level ϕ(x, t) is assigned to each lattice cell, acting as an indicator that describes the affiliation to one of the phases. Cells can be of liquid (ϕ(x, t) = 1), gas (ϕ(x, t) = 0), or interface type (ϕ(x, t) ∈ (0, 1)). A sharp and closed layer of interface cells separates liquid and gas cells. Interface and liquid cells are treated like regular LBM cells, which contain PDFs and participate in the LBM collision (1) and streaming (2). In contrast, conforming with the free-surface definition, gas cells neither contain PDFs nor participate in the LBM update.
The liquid mass of each cell is determined by the cell's fill level ϕ(x, t), fluid density ρ(x, t), and volume ∆x 3 . Note that in two-dimensional simulations, the cell's volume is also given by ∆x 3 . The domain is then assumed to have an extension of a single lattice cell in the third direction. The mass flux between an interface cell and cells of other types is computed from the LBM streaming step via The simplicity of this mass flux computation is an advantage of the FSLBM when compared to non-LBM-based VOF approaches. In these methods, the advection of mass commonly requires solving a partial differential equation that describes the evolution of the mass.
In the implementation used here, interface cells are not immediately converted to gas or liquid when their fill level becomes ϕ(x, t) = 0 or ϕ(x, t) = 1, respectively. Instead, they are converted with respect to the heuristically chosen threshold ε ϕ = 10 −2 that prevents oscillatory conversions [24]. Consequently, an interface cell is converted to gas or liquid if its fill level becomes below zero with ϕ(x, t) < 0 − ε ϕ or above one with ϕ(x, t) > 1 + ε ϕ , respectively.
When an interface cell converts to gas or liquid, surrounding gas or liquid cells may convert to interface cells to maintain a closed interface layer. It is important to point out that neither liquid nor gas cells can directly convert into one another. Instead, both cell types can only convert to interface cells. The separation of liquid and gas is prioritized in case of conflicting conversions. When converting an interface cell with fill level ϕ conv (x, t) to gas or liquid, the fill level is forcefully set to ϕ(x, t) = 0 or ϕ(x, t) = 1. The resulting excess mass is distributed evenly among all surrounding interface cells to ensure mass conservation.
During a simulation, unnecessary interface cells may appear, which do not have neighboring gas or liquid cells. In the implementation used in this study, these cells are forced to fill or empty by adjusting the mass flux (15), as suggested by Thürey [25].
The cells' PDFs are not modified when cells are converted from interface to liquid or vice versa. If a cell converts from interface to gas type, the cell's PDFs need not be considered further and can therefore be invalidated. Note that this does not affect mass conservation as any excess mass (16) will be distributed accordingly. However, no valid PDF information is available when cells convert from gas to interface type. The PDFs of these cells must be initialized with one of the schemes presented in Section III.
The LBM collision (1) and streaming (2) are performed in all interface and liquid cells.
Körner et al. [1] proposed to weight the gravitational acceleration with an interface cell's fill level in the LBM collision. Conforming with the work of other authors [24,26,27], the implementation used here did not weight the gravitational force with the fill level.
The macroscopic boundary condition at the free surface is given by [26,28] p where p G (x, t) is the gas pressure, p L (x, t) is the Laplace pressure, t 1 (x, t) and t 2 (x, t) are interface-tangent vectors, and n(x, t) is the interface-normal vector. As shown by Bogner et al. [29], this macroscopic boundary condition is approximated by the LBM anti-bounce-back pressure boundary condition which Körner et al. [1] suggested to use. In this equation, u ≡ u (x, t) is the interface cell's velocity and ρ G ≡ ρ G (x, t) = p G (x, t) /c 2 s is the gas density. Other formulations of the boundary condition have been investigated in the literature [29,30]. The free-surface boundary condition (18) is applied to all PDFs streaming from the gas towards the interface as these PDFs are not available. However, in the original FSLBM [1], this boundary condition is not only used to reconstruct missing PDFs. It is also used to reconstruct some PDFs that are already available. It should be pointed out that this approach overwrites existing information about the flow field. In the implementation used in the study presented here, no information is overwritten, and only missing PDFs are reconstructed at the free boundary.
This scheme was found to be of superior accuracy [18]. Note that the free-surface boundary condition (18) must also be applied at free-slip boundaries for specularly reflected PDFs that originate from gas cells.
The gas pressure incorporates the volume pressure p V (t) and Laplace pressure p L (x, t). The volume pressure stays constant in case of atmospheric pressure or results from changes in the volume V (t) of an enclosed gas volume, that is, a bubble, according to The Laplace pressure is defined by the surface tension σ ∈ R + and the interface curvature κ(x, t) ∈ R. As suggested by Bogner et al. [31], in the simulations shown in this article, a finite difference approximation was used, wheren(x, t) ∈ R d is the normalized interface normal vector. The interface normal was computed with central finite differences according to Parker and Youngs [32]. Near obstacle cells, the computation of the normal is modified as proposed by Donath [27]. This modification narrows the access pattern of the finite differences such that obstacle cells are excluded from the computation. A bubble model algorithm is used to keep track of the bubbles' volume pressure [24,33]. This is required because bubbles might coalesce or segment during the simulation.

III. REFILLING SCHEMES
In the FSLBM, gas cells do not contain valid PDF information as described in Section II B.
Therefore, when a gas cell converts to an interface cell, its PDFs must be reinitialized. This reinitialization is commonly referred to as refilling. While refilling has not yet been studied in the context of the FSLBM, it has been investigated for moving solid obstacle cells [8][9][10][11][12][13].
Analogously to the gas cells in the FSLBM, solid cells do not carry valid PDFs, such that these PDFs must be refilled after conversion. In what follows, the schemes developed for refilling moving solid cells are introduced and adapted to the FSLBM. Then, their influence on the conservation of mass and their computational costs are briefly discussed.

A. Scheme definitions
In the FSLBM proposed by Körner et al. [1], PDFs are refilled with their equilibrium (4) according to The average local fluid densityρ ≡ρ (x, t) and velocityū ≡ū (x, t) are computed from all surrounding, non-newly created interface and liquid cells. Note that this is in contrast to the same scheme applied for solid obstacle cells, where the velocity of the solid object is used instead [8,9,13].
The EQ scheme can be extended by adding the non-equilibrium contribution from neighboring fluid cells with [8] EQ+NEQ: 1} is used to access the cell that corresponds best to the local interface normal n ≡ n(x, t), that is, for which the scalar product c i · n is the largest. Note again that only non-newly created neighboring interface and liquid cells are valid directions for c n i . Another variant extends the EQ scheme with information from the local pressure tensor [10][11][12] using Grad's moment system [15], leading to where δ αβ is the Kronecker delta. In the implementation used in this study, the velocity derivatives are approximated by second-order finite differences. If not enough neighboring fluid cells are available, the derivatives are approximated by first-order backward-or forward-finite differences.
Instead of using the PDF's equilibrium, Lallemand and Luo [9] suggested refilling based on a second-order extrapolation scheme. The PDFs are extrapolated from the lattice direction closest to the surface normal c n i by In the implementations used for this study, a corresponding lower-order extrapolation is used if the number of neighboring cells in direction c n i is not sufficient. If no neighboring cell is available in this direction, the EQ scheme is employed as a fallback. A first-order extrapolation scheme led to the same numerical instabilities as will be discussed in Section IV for the second-order scheme. These numerical instabilities were not observed with a zerothorder extrapolation, that is, copying PDFs from a direct neighboring cell. However, this zeroth-order extrapolation scheme was more inaccurate than other refilling schemes, and it even led to physically implausible results in some tests. For brevity reasons, these results are not included in this article.
The AVG scheme uses the averagef i (x, t) of the identically oriented PDFs of surrounding, non-newly created interface and liquid cells. This is denoted by and has already been used by Fang et al. [14] for solid obstacle cells. However, Fang et al.
averaged the PDFs resulting from a higher-order extrapolation scheme, including a larger neighborhood of cells. In this work, only the PDFs from direct neighboring cells are used.

B. Effect on mass conservation
It is important to point out that the choice of the refilling scheme does not affect the system's total mass and its conservation. When a cell is refilled, it has converted from gas to interface and is initially empty with fill level ϕ(x, t) = 0. Therefore, the refilled cell's mass (14) is initially m(x, t) = 0 and is independent of the cell's refilled PDFs.
In the following time steps, the interface cell's mass then changes according to where ∆m i (x, t) is the mass flux (15). If the cell x + c i ∆t is also an interface cell, the same ∆m i (x, t) affects this interface cell's change in mass (29). If the cell x + c i ∆t is of liquid type, the corresponding ∆m i (x, t) is implicitly considered in the liquid cell's density ρ(x, t) and, therefore, in its mass, because the mass flux ∆m i (x, t) is computed directly from the balance of the streaming PDFs. Consequently, ∆m i (x, t) is present in the change of the PDFs' values. The density ρ(x, t) is obtained by taking the PDF's zeroth-order moment (5).
Since the zeroth-order moment is the sum of the cell's PDFs, ∆m i (x, t) leads to an according change of ρ(x, t).

C. Computational costs
The computational costs of the individual refilling schemes strongly depend on the implementation of the FSLBM and the test case simulated. For instance, the EQ, EQ+NEQ, and GEQ schemes require the neighboring cells' density ρ(x, t) and velocity u(x, t). However, these macroscopic quantities are not necessarily available in any LBM implementation but might have to be computed only when required. The LBM algorithm of collision (1) and streaming (2) works on PDFs but does not involve ρ(x, t) or u(x, t). Therefore, depending on the implementation, these values might be computed explicitly for the refilling schemes via the moments (5) and (6), increasing the computational costs.
As another example, the EQ+NEQ and EXT schemes involve the interface normal n.
The interface curvature κ(x, t) must be computed if the Laplace pressure (21) is relevant in a test case. The curvature computation algorithm employed in this work also uses n so that no additional computations may be required to obtain n when refilling cells.
It should also be pointed out that the computational costs are affected by the refilled cell's neighboring cells, as only non-newly created interface and liquid cells are considered by the refilling schemes.
These examples show that it is not generally possible to present the specific computational cost of each refilling scheme. Only the EQ, EQ+NEQ, and GEQ refilling schemes can be put into perspective, as the EQ+NEQ and GEQ build upon f eq i (ρ,ū) as computed by the EQ scheme but add additional computations. Therefore, the EQ+NEQ and GEQ schemes are computationally more expensive than the EQ scheme.
Note that although cell conversions appear frequently, the computational costs for refilling might be insignificant compared to the costs of other algorithmic parts in the FSLBM.
However, this likewise strongly depends on the test case under investigation.

IV. NUMERICAL EXPERIMENTS
The refilling schemes introduced in Section III are compared in six numerical experiments in this section. The chosen test cases are vastly similar to those suggested in prior work [17,18]. Therefore, the corresponding description of the test cases, simulation setups, and figures are based one those from these articles but are repeated here for completeness. The numerical benchmarks include the simulation of a standing gravity wave, the collapse of a rectangular and cylindrical liquid column, the rise of a Taylor bubble, the impact of a drop into a thin film of liquid, and a bubbly plane Poiseuille flow. All simulations were performed with double-precision floating-point arithmetic.

A. Gravity wave
A gravity wave is a standing wave oscillating at the phase boundary between two immiscible fluids. Surface tension forces are neglected, and gravitational forces entirely govern the wave's flow dynamics. The analytical model [34,35] is used as reference data for assessing the simulation results.

Simulation setup
As illustrated in Figure 1, a gravity wave of wavelength L was simulated in a two- to as diffusive scaling in the LBM [19]. The system is characterized by the Reynolds number where is the angular frequency of the wave, and ν is the kinematic fluid viscosity. Owing to the gravitational acceleration g, the initial profile evolved into a standing wave oscillating around d. It was dampened because of viscous forces. The dimensionless surface elevation a * (x, t) := a(x, t)/a 0 and non-dimensionalized time t * := tω 0 were monitored at x = 0 every t * = 0.01. The simulations were performed until t * = 40, which was found to be sufficient for the wave's motion to be decayed in the simulations.

Analytical model
The analytical model for the gravity wave is derived by linearization of the continuity and Euler equations with a free-surface boundary condition [34]. The standing wave's amplitude is obtained assuming an inviscid fluid with zero damping, such that a D (t) = a 0 . Viscous damping is considered by [35] a D (t) = a 0 e −2νk 2 t .
The analytical model applies if k|a 0 | 1 and k|a 0 | kd [34], which is true in this study 3. Results and discussion Figure 2 shows the gravity wave's non-dimensionalized amplitude a * (0, t * ) over time t * as simulated with L = 800 and the refilling schemes from Section III. All refilling schemes but the EXT generally agreed well with the analytical model. The other refilling schemes showed only minor differences. Most notably is a slight overestimation of the wave's first positive amplitude when using the AVG scheme. The duration of a gravity wave's half-period T * /2 is shown in Figure 3. The EXT refilling scheme had the largest deviations when compared to the analytical model. All other refilling schemes did not follow a clear trend. Figure 4 shows of the gravity wave's oscillations. Although the differences are relatively small, the EQ+NEQ scheme could arguably be considered most accurate in this comparison.
The simulation results presented here have converged, as illustrated in Figure 26. As pointed out in prior work [17], the FSLBM can only predict the wave's motion sufficiently well, if the amplitude spans over at least one but preferably multiple cells. This is also shown in Figure 26, where less wave periods could be simulated when using lower computational domain resolutions.

B. Rectangular dam break
In the rectangular dam break benchmark case, a rectangular liquid column collapses and spreads at the bottom surface. This test case is regularly used as a numerical benchmark to validate free-surface flow simulations [5,36,37]. The experiments from Martin and Moyce [38] were used as reference data for the simulations in this section.

Simulation setup
The simulation setup resembled that of the reference experiments [38], and is illustrated in Figure 5.   negative y-direction, the liquid was initialized with hydrostatic pressure. Therefore, the LBM pressure at y = H was initially equal to the constant atmospheric gas pressure p V (t) = p 0 .
Wetting effects were not considered, and free-slip boundary conditions were set at all domain walls. Conforming with diffusive scaling, the relaxation rate ω = 1.9995 was kept constant for all tested computational domain resolutions. The simulations were performed using the turbulence model from Section II A with Smagorinsky constant C S = 0.1 [23]. Two dimensionless numbers describe the fluid mechanics of the system. The Galilei number with kinematic viscosity ν, relates the gravitational to viscous forces. The Bond number defines the relation between gravitational and surface tension forces. It is defined by the surface tension σ, and the density difference between the liquid and gas phase ∆ρ = ρ − ρ G .
Note that in a free-surface system, the gas phase density is assumed to be zero so that ∆ρ = ρ.
The reference experiments [38] were performed with liquid water, but the authors did not 2. Results and discussion Figure 6 compares the simulated dam break with the experimental measurements [38] in terms of the non-dimensionalized width w * (t * ) and height h * (t * ). The simulations were performed with a computational domain resolution, that is, initial dam width of W = 200. The simulations with the EXT refilling scheme became numerically unstable because the macroscopic velocity gradually increased after refilling, and eventually locally exceeded the lattice speed of sound c s , as illustrated in Figure 8. Exceeding the lattice speed of sound is often an effect of a scheme being numerically unstable in the LBM [19]. These numerical instabilities eventually lead to the collapse of the simulation. Note that these instabilities were not an immediate consequence of a certain cell being refilled with the EXT scheme.
Instead, high macroscopic velocities appeared in the later course of the simulation. All other refilling schemes produced results of similar accuracy and agreed with the trend of the experimental observations for h * (t * ). The simulated dam width w * (t * ) generally agreed better with the experimental measurements than the simulated height. However, the choice of the refilling scheme had more effect on w * (t * ). The AVG scheme was the most accurate, while the EQ and GEQ schemes were less accurate. Although the EQ+NEQ scheme also produced accurate results, it temporarily deviated more significantly from the experimental data at t * ≈ 3.7. Except for the unstable simulation with the EXT refilling scheme, the simulated dam contours at t * = 3 are visualized in Figure 7. While there are noticeable differences between the schemes, no accuracy assessment could be made due to the lack of suitable experimental reference data.
As shown in Figure 27, the simulation results presented in this section were converged in terms of computational domain resolution. A resolution equivalent to W = 50 was sufficient to reasonably agree with the experimental data.

C. Cylindrical dam break
In this section, the simulation setup and results for a cylindrical dam break are presented.
The numerical simulations resemble the laboratory experiments from Martin and Moyce [38].
This test case was chosen to evaluate whether the model's isotropy is affected by the choice of the refilling scheme.

Simulation setup
As illustrated in Figure 9, The liquid column's radius r(t) was monitored during the simulation. It was obtained by finding the distance of the liquid front to the column's initial center of symmetry. The liquid column's collapse can not be assumed to be perfectly symmetric. Consequently, r(t) was computed for every interface cell detected by a seed-fill algorithm [40]. The starting point of this algorithm was set to an arbitrary domain boundary. With this configuration, the algorithm only detected the outermost interface cells, that is, only interface cells at the It collapsed due to the gravitational acceleration g that acted in negative z-direction. spreading liquid's front. A statistical sample was used to evaluate r(t) by computing the maximum, minimum, and mean values of r(t) at every t * = 0.01. The radius r * (t) := 2r(t)/D and time t * := t 4g/D were non-dimensionalized as suggested in the reference data [38].
Conforming with the experimental data [38], the simulations were stopped at r * max (t * ) ≥ 4.33 with the non-dimensionalized maximum liquid front radius r * max (t * ). Large error bars indicate a deviation from rotational symmetry. As for the breaking dam test case in Section IV B, the EXT refilling scheme was numerically unstable. All other refilling schemes generally agreed well with the experimental data [38]. Although the error bars for the EQ+NEQ scheme indicate asymmetry, the liquid spread's front remained qualitatively symmetric as shown in Figure 11. However, tiny droplets detached from the main liquid spread. The evaluation algorithm detected these droplets as part of the liquid surge front.

Results and discussion
The EQ and GEQ schemes are approximately of equal accuracy, while being less accurate than the AVG scheme. The shape of the collapsing liquid column at time t * = 4 is visualized in Figure 11. The solid black lines indicate the initial center of symmetry. In general, all refilling schemes remained rotationally symmetric and did not move from their initial center of symmetry.
As illustrated in Figure 28,   indicates the column's initial center of symmetry. The simulation with the EXT refilling scheme was numerically unstable and is not included here. All other schemes kept their initial center of symmetry and remained rotationally symmetric. There are slight differences between the individual refilling schemes. Owing to the lack of reference data, no accuracy assessment could be made in terms of shape.

D. Taylor bubble
A Taylor bubble is a gas bubble that rises in a cylindrical tube filled with a stagnant liquid due to buoyancy forces. It has an elongated shape and a round leading edge with a length of multiple times its diameter. The simulation results were compared to the experimental data from Bugg and Saad [41].

Simulation setup
The simulation setup is illustrated in Figure 12 and resembled that of the reference experiments [41].  where ν is the kinematic viscosity, for different computational domain resolutions. The bubble's rise velocity u was computed from the bubble's center of mass in the z-direction at time t * = 10 and t * = 15. The simulations generally agreed well with Re from the reference data [41], and there were only minor differences between the refilling schemes. Similarly, as illustrated in Figure 13 for t * = 15, the choice of the refilling scheme had almost no effect on the shape of the simulated Taylor bubble. Figures 15 to 17 compare the non-dimensionalized axial u * a = u a /u and radial u * r = u r /u velocities at the locations defined in Figure 14. As for the Taylor bubble's shape, the refilling schemes led to only small differences in the velocity profiles. Only the radial velocity u * r at a radial line at 0.111D from the Taylor bubble's front, was arguably predicted more accurately by the EQ and GEQ schemes, as shown in Figure 16. Table I and Figure 29, the simulation results shown here have sufficiently converged in terms of computational grid resolution.        Figure 16: Simulated non-dimensionalized axial velocity u * a (a) and radial velocity u * r (b) along the radial monitoring-line at 0.111D from the Taylor bubble's front (see Figure 14).

As depicted in
The comparison with experimental data [41] is drawn in terms of the non-dimensionalized radial location r * at dimensionless time t * = 15 with tube diameter D = 128 lattice cells.  Figure 17: Simulated non-dimensionalized axial velocity u * a (a) and radial velocity u * r (b) along the radial monitoring-line at −0.504D from the Taylor bubble's front (see Figure 14).
The comparison with experimental data [41] is drawn in terms of the non-dimensionalized radial location r * at dimensionless time t * = 15 with tube diameter D = 128 lattice cells.

E. Drop impact
In the fifth test case, the vertical impact of a drop into a pool of liquid was simulated.
As no quantitative experimental measurements are given for the reference experiments from Wang and Chen [42], only a qualitative comparison with photographs could be made.

Simulation setup
As illustrated in Figure 18, a spherical droplet with a diameter of D = 80 lattice cells was initialized in a three-dimensional computational domain of size 10D × 10D × 5D (x-, y-, z-direction). The droplet was initially located at the surface of a thin liquid film of height 0.5D with impact velocity U in the negative z-direction. The gravitational acceleration g also acted in the negative z-direction. Accordingly, hydrostatic pressure was initialized such that the pressure at the pool's surface was equal to the constant atmospheric volumetric gas pressure p V (t) = p 0 . The walls in the x-and y-direction were periodic, and no-slip boundary conditions were set at the top and bottom domain walls in z-direction. The relaxation rate was chosen ω = 1.989. The droplet's impact is described by the Weber number which relates inertial and surface tension forces, and by the Ohnesorge number which relates viscous to inertial and surface tension forces. These dimensionless numbers include the surface tension σ, dynamic viscosity µ, and liquid density ρ. Assuming g = 9.81 m/s 2 , and using ρ = 1200 kg/m 3 and µ = 0.022 kg/(m·s) [42], the Bond number Bo = 3.18 (36) with characteristic length D, closes the definition of the system. The nondimensionalized time t * := tU/D is offset by t * = 0.16 [6] to allow a comparison with the numerical simulations performed in this study.

Results and discussion
The simulated drop impact, that is, the splash crown formation at t * = 12, is qualitatively compared with experimental results in Figure 19. The solid black line indicates the crown's contour in a central cross-section, oriented with the normal in the x-direction. There was no scale bar provided for the photograph of the laboratory experiment [42]. Therefore, the simulations performed here could only be validated qualitatively. As in the dam break simulations in Sections IV B and IV C, the EXT refilling scheme became numerically unstable which led to too high macroscopic velocities. The EQ+NEQ scheme was subject to numerical instabilities for the same reason. Qualitatively plausible results could be obtained with all other refilling schemes. However, with the GEQ scheme, the droplets detaching from the crown's top formed thin and long threads of liquid. In contrast, in the photograph of the experiment, the detaching droplets rather form thicker and shorter liquid threads that then detach as spherical droplets. This kind of crown formation is arguably resembled best by the EQ scheme.

F. Bubbly plane Poiseuille flow
This final benchmark case is inspired by Peng et al. [8], where a particle-laden turbulent channel flow was simulated. The choice of the refilling scheme for solid obstacles affected the particle dynamics, that is, the particles' position during the simulation. In the study presented here, a similar test case is used with randomly initialized spherical bubbles rather than solid particles. The flow is force-driven between two parallel plates, also called plane Poiseuille flow.

Simulation setup
A three-dimensional domain of size 2L × L × L (x-, y-, z-direction) with a channel width of L = 100 lattice cells was filled with liquid, as illustrated in Figure 23. As shown in Figure 24, there were 381 randomly distributed spherical bubbles with a diameter of 0.1L in the channel, leading to a gas volume fraction of approximately 0.1. The bubbles were arranged so that their center was not closer than 0.05L to the domain wall. The random distribution was chosen once and kept the same for all simulations. Therefore, the simulation with any refilling scheme started from an identical initial situation. The domain's walls were periodic in the x-and y-direction and set to no-slip in the z-direction. With the force F acting in the x-direction, the fluid velocity profile took a parabolic shape in the z-direction with zero-velocity at the no-slip domain walls. This velocity profile is commonly referred to as plane Poiseuille flow and analytically given by [43] where µ is the dynamic viscosity of the liquid. The setup is defined by the Morton number Mo = 10 −5 (37) and by the Reynolds number Re = 10 4 (38) with characteristic length L and analytical maximum velocity u max = u(0.5L). The relaxation rate was chosen ω = 1.989, and the time t was non-dimensionalized with t * = t u max /L.

Results and discussion
As for the drop impact test case in Section IV E, the simulations with the EXT and EQ+NEQ refilling schemes were numerically unstable. The simulation results for the AVG, EQ, and GEQ schemes at t * = 4 are shown in Figure 25. The bubbles gathered and coalesced in the center of the domain, where the velocity was the highest. Although all simulations started from the same initial situation, the refilling schemes led to noticeable differences in the bubble dynamics, that is, the bubbles' positions. This observation agrees with those made by Peng et al. [8], where a turbulent particle-laden channel flow was simulated with different refilling schemes for solid obstacles. However, due to the lack of missing reference data from experiments, the refilling schemes' accuracy could not be assessed in this test case.

V. CONCLUSIONS
This study has compared different refilling schemes for the free-surface lattice Boltzmann method [1] (FSLBM). In the FSLBM, it is distinguished between cells belonging to the heavier fluid, lighter fluid, and the interface located between them. These cells are here referred to as liquid, gas, and interface cells, respectively. The gas phase is neglected and the interface is treated as a free surface. Consequently, gas cells neither participate in the LBM flow simulation, nor carry valid information about the flow field. In the LBM, such information is stored in terms of particle distribution functions (PDFs) in each lattice cell.
Because of the free interface's motion, gas cells regularly convert to interface cells. As the hydrodynamic LBM simulations are performed in interface cells, those cells' PDFs must be initialized with valid information during the conversion. This initialization of PDFs is commonly referred to as refilling. The first refilling scheme under investigation was the one suggested in the original FSLBM as introduced by Körner et al. [1]. In this model, PDFs are initialized according to their equilibrium (EQ), which is constructed using the average density and velocity from the neighboring, non-newly converted interface and liquid cells. This scheme was extended by adding the contribution of the neighboring cells' non-equilibrium PDFs (EQ+NEQ) [8], or by including information about the local pressure tensor using Grad's moment system (GEQ) [10][11][12]15]. Additionally, the PDFs could also be extrapolated (EXT) from neighboring cells' PDFs [9] or were taken as the average (AVG) of neighboring, non-newly converted interface and liquid cells' PDFs [14].
These schemes' accuracy and stability properties were investigated in six numerical experiments, with reference data for five of them as either analytical models or laboratory measurements from the literature. In the experiments conducted here, the EXT and EQ+NEQ schemes often led to numerical instabilities. These instabilities were caused by the lattice velocity exceeding the lattice speed of sound. The AVG refilling scheme was also unstable in the cylindrical dam break test case in one of the computational domain resolutions used in the convergence study. In contrast, the EQ, and GEQ schemes were numerically stable in all simulations performed here. Although, the AVG scheme was more accurate than the EQ and GEQ schemes in the dam break test cases, it slightly overestimated the gravity wave's amplitude in the first period. The EQ, and GEQ schemes' simulation results hardly differed when compared in terms of the quantitative reference data available in the literature. Nevertheless, qualitative differences between these schemes could be observed in the dam break, drop impact, and bubbly Poiseuille flow test cases. Because of lacking appropriate reference data, a final accuracy comparison could only be made vaguely based on visual comparison. In the drop impact benchmark, the EQ scheme arguably seemed to be favorable over the GEQ scheme when compared qualitatively. Additionally the GEQ scheme is computationally more expensive than the EQ scheme. In summary, for the numerical simulations performed here, the EQ scheme should be preferred in the FSLBM with respect to ease of implementation, computational costs, numerical stability, and accuracy.

SUPPLEMENTARY MATERIAL
The following supplementary material is available as part of the online article: An archive of the C++ source code used in this study. It is part of the software framework waLBerla [16] (version used here: https://i10git.cs.fau.de/walberla/walberla/ -/tree/01a28162ae1aacf7b96152c9f886ce54cc7f53ff). The ready-to-run simulation setups for all numerical experiments performed in this article are included in the directory apps/showcases/FreeSurface.