FULLY-DISCRETE, DECOUPLED, SECOND-ORDER TIME-ACCURATE AND ENERGY STABLE FINITE ELEMENT NUMERICAL SCHEME OF THE CAHN-HILLIARD BINARY SURFACTANT MODEL CONFINED IN THE HELE-SHAW CELL

. We consider the numerical approximation of the binary fluid surfactant phase-field model confined in a Hele-Shaw cell, where the system includes two coupled Cahn-Hilliard equations and Darcy equations. We develop a fully-discrete finite element scheme with some desired characteristics, including linearity, second-order time accuracy, decoupling structure, and unconditional energy stability. The scheme is constructed by combining the projection method for the Darcy equation, the quadratization approach for the nonlinear energy potential, and a decoupling method of using a trivial ODE built upon the “zero-energy-contribution” feature. The advantage of this scheme is that not only can all variables be calculated in a decoupled manner, but each equation has only constant coefficients at each time step. We strictly prove that the scheme satisfies the unconditional energy stability and give a detailed implementation process. Various numerical examples are further carried out to prove the effectiveness of the scheme, in which the benchmark Saffman-Taylor fingering instability problems in various flow regimes are simulated to verify the weakening effects of surfactant on surface tension.

characteristics caused by the amphiphilic nature of surfactant molecules is postulated. By using the gradient flow method (Cahn-Hilliard dynamics) to minimize the free energy, a governing system composed of two highly nonlinear coupled Cahn-Hilliard equations can be obtained.
To simulate the free interface motion driven by the flow field (e.g., droplets coalescence/non-merging phenomena under the shear flow, or Saffman-Taylor fingering instability in the porous medium flow, etc., see [12,14,16,24,28,31,34,43,56]), the hydrodynamic equations will be coupled with the two Cahn-Hilliard equations to obtain a full flow-coupled binary fluid surfactant model. Different fluid equations will be used in different flow regimes, such as the Navier-Stokes equations for incompressible fluids, the Darcy equations for the porous medium, etc. Obviously, due to the existence of more coupled nonlinear terms, designing numerical schemes for a flow-coupled model is far more difficult than that of a model without flow. Then, for the Darcy flow-coupled binary surfactant Cahn-Hilliard phase-field model (abbreviated as DCHS) considered in this article, a natural question is how difficult it is to design an "ideal" type scheme. One may think it is not difficult at all because there exist so many effective numerical methods that can handle the phase-field equations (e.g., the Invariant Energy Quadratization (IEQ) method [39,48,53], Scalar Auxiliary Variable (SAV) method [40,49,52,57], convex splitting method [20], etc.), and the "ideal" scheme can be easily and naturally obtained after stacking these methods with any effective numerical method of the Darcy equation.
However, the fact is quite the opposite. As far as the author knows, for the Darcy coupled with the phase-field equations, the only scheme owning the properties of "decoupling, second-order time accuracy, energy stability" (partially "ideal" type scheme due to the lack of linearity) was developed in [21,22]. The scheme uses the implicit-explicit combination method to deal with the advection and surface tension terms. Using the linear relationship between the pressure gradient and velocity, the fluid velocity in the advection term is formulated by the pressure and surface tension terms, so as to achieve a fully-decoupled scheme (see the details introduced in Rem. 3.9). However, if the same technique developed in [21,22] is applied to the DCHS model, we will quickly find that the fluid velocity in the advection term involves both phase-field variables. This means that although the velocity field can be decoupled from the two Cahn-Hilliard equations, the price to pay is that the two Cahn-Hilliard equations remain coupled. In addition, another disadvantage of this method is that the Cahn-Hilliard equation with variable coefficients must be solved at each time step, which results in a higher computational cost than just solving equations with constant coefficients (see Rem. 3.9).
Therefore, in order to solve the highly complex nonlinear DCHS model, we aim to construct an "ideal" type fully-discrete numerical scheme, where the main challenge that needs to be overcome is how to obtain a decoupling structure and second-order accuracy while keeping the linearity and maintaining the energy stability unconditionally. We achieve such a scheme by assembling several proven methods, including the finite element method for the spatial discretization, projection method for the Darcy equation, the quadratization approach for the nonlinear energy functional, and a new decoupling method to deal with the coupled terms (advection and surface tension). Inspired by the "zero-energy-contribution" property (see Rem. 2.3) satisfied by these coupled terms, the decoupling method is developed by introducing an additional nonlocal variable and designing a special ordinary differential equation (ODE), which consists of the inner product of the coupled terms with some specific functions. This ODE is trivial at the continuous level because all the terms contained in it provide a zero summation. But after discretization, it can help to obtain unconditional energy stability. Meanwhile, the nonlocal variable can decompose each discrete equation into multiple sub-equations that can be solved independently, thereby obtaining a fully-decoupled structure. Besides, the high efficiency of this scheme is also reflected in that not only can all variables be calculated in a decoupled manner, but all equations have constant coefficients at each time step. We also give rigorous proofs of the solvability and unconditional energy stability of the scheme. To demonstrate the stability and accuracy numerically, we further simulate various numerical examples, including the Saffman-Taylor fingering instability caused by the continuously injected radical/uniform flow or rotating Hele-Shaw cell, in which the weakening effects of surfactants on surface tension can be verified by the number of generated fingers.
The rest of this article is organized as follows. We first introduce the DCHS model limited in the Hele-Shaw cell in Section 2, and its law of energy dissipation is derived. In Section 3, we construct a fully-discrete finite element scheme of the "ideal" type, and describe its implementations in detail. Unconditional energy stability and solvability are also proved rigorously. In Section 4, we perform several accuracy/stability tests and implement various simulations on the Saffman-Taylor fingering instability problems to demonstrate the effectiveness of the scheme. In Section 5, some concluding remarks are given finally.

Governing System
Now, we give a brief introduction of the DCHS system limited to the Hele-Shaw cell. Suppose that Ω is a smooth, open, bounded, connected domain in R , = 2, 3. We adopt two phase-field variables ( , ) and ( , ), where is used to represent the density (or volume fraction) of the two fluids, and is used to represent the local concentration of surfactants, i.e., with a thin, smooth transition region with a width ( ). The mixing free energy associated with and is assumed as follows (see [28,44,45,49,56]), and 1 , 2 , , , , , are all positive parameters. The total free energy includes the hydrophilic (gradient)hydrophobic (double-well) trend of the phase-field variable where the parameter represents the width of the binary fluid interface, and the hydrophilicity (gradient)-hydrophobic (quartic polynomial ) trend of the concentration variable where is a penalty parameter. The last component ( , ) includes the coupling item between the surfactant and the fluid interface, where the parameter controls the degree of the effect of the surfactant accumulated on the fluid interface, and the term ( ≪ 1) is used to ensure that the total free energy is bounded from below, see [49,56].
Remark 2.1. There are many different options for the hydrophobic functional ( ). For example, it can be selected as the Flory-Huggins logarithmic type which reads as ( ) = ln + (1 − )ln(1 − ) (cf. [44,45,56]), double-well type (cf. [25]), or quartic polynomial type given in (2.3) (cf. [49]). The last option is used here because it can easily obtain the boundedness property (from below) of the energy potential, which is significant for the well-posedness of the model and the development of numerical algorithms.
Assuming that the fluid motion conforms to the mechanical principles in porous medium, the DCHS model reads as:

7)
u + ( )u + ∇ + ∇ + ∇ = 0, (2.8) u is the dimensionless seepage velocity, is a positive parameter, 1 , 2 are two mobility parameters, is the dimensionless hydraulic conductivity, ( ) = 1 2 1 (1 − ) + 1 2 2 (1 + ) is the fluid viscosity [14], 1 and 2 are the viscosity for fluid I and II, respectively, is the pressure. Note that the time derivative of the seepage velocity u is retained for flows in porous medium, cf. [4,21,22,30]. The boundary conditions read as where is the unit outward normal on the boundary Ω. Note that the above boundary conditions also implies | Ω = 0. It is also possible to assume periodic boundary conditions for all variables. The initial conditions read as The total mass of local density variable and concentration variable are conserved over time since which can be obtained by integrating (2.4) and (2.6) and using the boundary conditions (2.11). Some notations are introduced here. We denote the 2 inner product of any two functions ( ) and ( ) is denoted by ( , ) = ∫︀ Ω ( ) ( ) , and the 2 norm of ( ) is denoted by ‖ ‖ 2 = ( , ). It is easy to see that the entire DCHS system (2.4)-(2.9) holds the law of energy dissipation by performing the standard derivation as follows.
Lemma 2.2. The following energy law holds for the system (2.4)-(2.9): Proof. First, we take the inner product of (2.4) with , of (2.6) with in 2 , and use integration by parts to get Second, we take the inner product of (2.5) with − , of (2.7) with − in 2 , and use integration by parts to get Third, we take the inner product of (2.8) with u in 2 , and use (2.9) to obtain (2.20) Combining the above five equations, we derive the energy law (2.14).
Remark 2.3. In the process of deriving the energy law (2.14), we find that the four nonlinear terms associated with the advection and surface tensions are cancelled out, namely, These equalities are derived by using the integration by parts and the boundary conditions for u (periodic or u · | Ω = 0). These two equalities can be regarded as the contribution of two types of nonlinear terms (advection: ∇ · (u ) and ∇ · (u ), and surface tensions: ∇ and ∇ ) to the total free energy of the system is zero. These unique "zero-energy-contribution" properties will be used to design decoupling type numerical schemes.

Numerical scheme
In this section, we aim to construct an "ideal" type fully-discrete finite element scheme to solve the DCHS system (2.4)-(2.9). Some special processing is needed to develop appropriate temporal discretizations for the challenging terms, including the advection, the surface tension, the nonlinear cubic term ( ), and the linear coupling between velocity and pressure through the divergence-free condition.

Reformulated equivalent system and energy law
We introduce a nonlocal variable ( ) and design an ODE system for it, that reads as: It is easy to see that the system (3.1) is the same as a trivial ODE system of = 0, | =0 = 1 with the exact solution of ( ) = 1.
We further define a nonlocal variable ( ) as where ( , ) = 1 , is a constant to guarantee the radicand always be positive. We can always find such a constant since the nonlocal term in the radicand is always bounded from below. This is because it is easy to see that − 2 |∇ | 2 ≥ − 4 |∇ | 4 − 2 for some constant , and − 2 can be always bounded by the quartic polynomial ( ).
We rewrite the original system (2.4)-(2.9) using the new variables and to the following: The transformed system (3.3)-(3.10) satisfies the following initial conditions, Some detailed explanations about the reformulation are given in the remarks below.
Remark 3.1. We multiply the advection terms (∇ · (u ), ∇ · (u )) and surface tension terms ( ∇ , ∇ ) with . Since the nonlocal variable ( ) is equal to 1, the PDE system will not be changed by this modification. Meanwhile, note that the three new equations ( Remark 3.2. The mass-conserved property of the two phase-field variables , remains unchanged in the new system (3.3)-(3.10) noting that the variable is nonlocal (i.e., ∫︀ Ω ∇ · (u ) = ∫︀ Ω ∇ · (u ) ). By integrating for (3.4) and (3.6) and using integration by parts, we still derive This means that the mass of two fluids and the concentration of the surfactants are still conserved over time.
Proof. First, we take the inner product of (3.3) with , and of (3.5) with in 2 , respectively, and use integration by parts to get Second, we take the inner product of (3.4) with − , and of (3.6) with − in 2 , respectively, and use integration by parts to get By multiplying (3.7) with 2 , we obtain Third, we take the inner product of (3.8) with u in 2 , and use (3.9) to obtain By multiplying (3.10) with , we obtain (3.21) Combining (3.15)-(3.21), we derive the law of energy disspation (3.14), that is now related to the new variables , .

Fully-discrete finite element numerical scheme
We develop the fully-discrete finite element scheme for the modified model (3.3)-(3.10) in this section. We first formulate the PDE system (3.3)-(3.10) to the weak form. Some Hilbert spaces are introduced as follows: The weak formulation of the system (3.3)-(3.10) reads as: find ( , , , , u, Before the numerical scheme, we introduce a few finite dimensional discrete subspaces. Assuming that the polygonal/polyhedral domain Ω is discretized by a conforming and shape regular triangulation/tetrahedron mesh ℎ that is composed by open disjoint elements such thatΩ = ⋃︀ ∈ ℎ¯. We use to denote the space of polynomials of total degree at most ( ≥ 1) and define the following finite element spaces: (3.31) . Using the second-order backward differentiation formula for the time derivative, we construct the fully discrete scheme to solve the system (3.23)-(3.30) as follows: find +1 and 1 and 2 are two positive stabilization parameters. Several remarks are in order.
Remark 3.4. The computations of the second-order scheme (3.32)-(3.40) require all values of = 1 . In practice, we can obtain these values by constructing a similar first-order scheme based on the backward Euler method as long as we set = = 2, = 0, * = 0 , ∀ . Meanwhile, similar to the continuous case, the massconserved property of +1 ℎ and +1 ℎ still holds, which can be proved by taking Θ ℎ = ℎ = 1 in (3.32) and (3.34). We can derive ∫︀ Remark 3.5. A second-order pressure-correction scheme is used to decouple the computation of the pressure from that of the velocity.ũ +1 is the intermediate velocity and u +1 is the final velocity field that satisfies the divergence free condition. u +1 ℎ satisfies the discrete divergence free condition that can be deduced by taking the 2 inner product of (3.40) with ∇ ℎ where ℎ ∈ ℎ , that is where (3.39) is used.
Remark 3.6. When the system has very high stiffness issues caused by the model parameters or other conditions, while some numerical methods are formally unconditionally energy stable, but exceedingly small time steps are needed to achieve reasonable accuracy, see the stabilized-IEQ/SAV methods in [9,39,47,55]. To fix such an inherent deficiency, a commonly used effective way is to add one or more extra linear stabilization terms with the corresponding temporal order (cf. the second-order term associated with 1 in (3.33) and 2 in (3.35)). The scale of the splitting errors caused by this term are about 1 2 2 (·) and 2 2 2 (·), which are actually consistent with the error caused by the second-order extrapolated nonlinear term ( ) and ( ). In Section 4, we present numerical evidence to show that this stabilizer is effective to improve the energy stability while using large time steps in Figure 3. Similar linear stabilization techniques had been widely used in the numerical scheme for solving the phase-field type models, e.g., the methods of linear stabilization, IEQ, SAV, convex-splitting methods, etc., see [8,9,11,38,39,47,50,51,54,55].
) is actually a temporal second-order approximation of the term tot ( , , u, , ) at = +1 . Since for any smooth variable with time, we always have the following heuristic approximations as We now show that the scheme (3.32)-(3.40) is unconditionally energy stable. We will use the following three identities repeatedly: Theorem 3.8. The following discrete energy dissipation law holds for the fully-discrete scheme (3.32)-(3.40) that reads as Taking the 2 inner product of the above equality with u +1 ℎ , we derive where (3.42) is used. By using (3.51) and (3.52), we deduce (3.53) We rewrite (3.40) as Taking the 2 inner product of the above equation with itself and multiply the result with 3 2 , we derive We rewrite (3.40) as By taking the 2 inner product of the above equation with 3u +1 ℎ , using (3.42) and (3.44), we obtain (3.57) and and We multiply (3.36) with 2 +1 to obtain (3.63) Hence, by combining (3.57)-(3.63) and using (3.45), (3.46), we arrive at (3.64) Finally, by dropping the positive terms in { } of (3.64), we obtain (3.47).

Decoupled implementation
In this subsection, we develop the decoupled implementation process for the scheme (3.32)-(3.40) in which we make full use of the nonlocal property of the auxiliary variable .
Step 5: we update u +1 ℎ and +1 ℎ from (3.39) and (3.40). From the above-detailed implementation process, it can be seen that the calculation of all unknown variables are completely decoupled. At each time step, the total cost only includes the computations of several elliptic equations. The decoupling of all equations and the characteristic of having only constant coefficients means very efficient practical calculations. Remark 3.9. For the sake of completeness, here we present the fully-decoupled scheme developed in [22] for solving the Darcy coupled Cahn-Hilliard equation. For simplicity, only the first-order time marching scheme is given here, since the second-order version follows the same line. We only discretize related terms, while other irrelevant terms remain unchanged. The scheme in [22] reads as Note that the advection ∇ · (u ) and surface tension ∇ are discretized using the traditional explicit-implicit combination method. Due to the special format of the Darcy equations, by using the explicit linear relationship between u +1 and other items, the decoupling structure of the above scheme can be obtained. More precisely, one can rewrite (3.91) as Then, by replacing u +1 in the scheme (3.90) using (3.92), the scheme for the Cahn-Hilliard equation is formulated as This scheme does achieve a full decoupling structure, i.e., the computation of +1 is independent of the velocity field u +1 . However, the paid price is that the phase-field equation is then equipped with variable coefficients (i.e., the coefficients of +1 ) at each time step, which increases the practical computational cost. Moreover, it is worth noting that if we apply similar techniques to the DCHS model studied here, (3.92) will also include ∇ +1 . This means the two Cahn-Hilliard equations must be coupled together. Therefore, the decoupling technique given in [22] can not actually obtain the expected full decoupling structure for the DCHS model.

Numerical simulation
In this section, we investigate the accuracy, energy stability, and effectiveness of the proposed scheme (3.32)-(3.40), numerically, in which all the numerical simulation are implemented using the FreeFem++ [23].

Accuracy and stability test
In this example, we perform several convergence and stability tests in 2D of the finite element scheme (3.32)-(3.40), referred to as DSAV for short. We set the computational domain as Ω = [0, 2] 2 . We use = ℓ = 2, that is, 2 finite elements for , , , , , and 1 for u.
We first verify the convergence order of the scheme DSAV in space and time, where we impose some suitable source terms so that the following solutions satisfy the system (2.   To verify the temporal convergence order, the linear relation between the spatial resolution and temporal resolution is assumed such that the temporal error is dominant. In Figure 1a, the errors measured in 2 norm  between the numerical solution and the exact solution at = 1 are plotted, where we vary different time step size and the spatial grid size is set as ℎ = . It can be observed that the scheme DSAV presents second-order temporal accuracy for all variables. To verify the spatial convergence order, in Figure 1b, we plot the errors measured in various norms at = 1 for various mesh size ℎ, in which we choose sufficiently small ( = 1 − 5) so that the errors are dominated by the spatial discretization error. We can see that the third-order convergence rates are followed by the pressure , , in the 2 norm, while the second-order convergence rates are observed for the u in the 2 norm, in the 1 norm. These results are in full agreements with the theoretical expectation of accuracy for the adopted finite element space. In Fig. 2, we test the convergence rate for low stiffness parameters = = 1, = 0.001. It can be seen all variables follow the second-order convergence order.
We further perform the stability test by plotting the time evolution of the total free energy where we set the initial conditions to to be two circles with different radii, that read as follows (shown in the small inset figure in Fig. 3b), We still use the model parameters given in (4.2) and grid size ℎ = 2 256 . In Figure 3a, the evolution curves of the total free energy (3.48) computed by the scheme DSAV are shown, where different time step sizes are used. All energy curves computed show monotonic decays, thereby verifying its unconditional stability. When the time step is large (e.g., ≥ 0.1/2 2 ), the obtained energy curves display large deviations from others, which means that the larger the time step, the larger the error. When the time step is small (e.g., ≤ 0.1/2 3 ), the obtained energy curves are overlapped, which means that the smaller the time step, the more accurate the computation result. To illustrate the effectiveness of the stabilizations ( 1 and 2 ), for comparison, in the small inset figure of Figure 3a, we append the energy curves computed by DSAV but with 1 = 2 = 0. We find that, with the absence of these two stabilizers, the energy blows up quickly even with very tiny time steps (e.g., = 0.1/2 13 ), which shows these two stabilizers are very effective in improving stability. In Figure 3b, we plot the time evolution of the total free energy in the original form (2.15) and the discrete form (3.48) using the same time step size = 0.1/2 3 . Both energy curves overlap very well. We also attach the profiles of at various times and find that the coarsening effect causes the small circle to gradually absorb into the large circle.

Saffman-Taylor fingering instability
In this subsection, we simulate the Saffman-Taylor fingering pattern instability problem, which demonstrates the formation and evolution of elaborate patterned structures, as one of the most in-depth benchmark research problems in fluid dynamic systems. It considers the development of interfacial instabilities when a fluid with low viscosity displaces another fluid of higher viscosity between the narrowly spaced plates of a Hele-Shaw cell. Several widely studied versions of this problem are the so-called continuous injection (radial or uniform) [5,7,13,35], or rotating Hele-Shaw cell [2,6,35,43]. Below we carry out numerical simulations on these three problems respectively. In each case, we compare the number of fingers formed to verify the weakening effects of surface tension by surfactants.

Radial injection of a less viscous fluid
The first numerical example considers a radial injection of low-viscosity fluid to displace high-viscosity fluid. Initially, except for a small circular area in the center of the domain occupied by the low viscosity fluid, the entire domain is filled with a fluid with high viscosity. The two fluids are immiscible, and the low-viscosity fluid intrudes from the center of the domain at a constant injection rate. Physical and numerical experiments show that as the size of the fluid interface increases outward, fingers are formed, spread out and move in the direction of separation from each other, and finally get very complex branching patterns, see [1,3,7,13,17,29,35,37,43].
We implement numerical simulations in 2D with the computational domain Ω = [0, 2 ] 2 . To represent the injection flow, we adopt the Gaussian method designed in [43] where a potential radial velocity u is imposed in the momentum equation and the Cahn-Hilliard equation for . Namely, the two corresponding equations are modified as , is the injection strength, 0 is the radius of the circular injection region, is a small quantity such that + ̸ = 0. We set the initial conditions as follows:  where 1 is the (high) viscosity of the displaced fluid, 2 is the (low) viscosity of the injected fluid. The profile of the initial condition 0 is shown in Figure 4a. We change the coupling parameter to investigate the effect of the surfactant on finger formation. We first set the coupling parameter = 0 to study the way the fingers are formed without surfactants. In Figure 4, snapshots of the phase-field variable at different times are drawn, in which a grayscale image is used to obtain a clearer view. It can be seen that when the fluid with lower viscosity is continuously injected into the domain, and the size of the droplet increases to a certain size, some short and thick fingers are formed (12 fingers marked in the last subfigure of Fig. 4).
We further apply the surfactant effect by setting = 0.004 and show snapshots of at different times in Figure 5a. Unlike the case without surfactants, the number of finger-like patterns finally obtained not only far exceeds (19 fingers, shown in the last subfigure of Fig. 4), but the formed finger shows more slender and longer structure, which indicates that the applied surfactant effectively weakens the surface tension. In Figure 5b, the profile of at = 2.5 is drawn, where one can see that the surfactant concentration is high around the fluid interface due to the coupling term between and . In Figure 5c and d, we plot the pressure and the velocity field at = 2.5 where the contour line { = 0} is appended. It can be seen that the velocity field presents a radial pattern in the vicinity of each finger.
In Figure 6, by using various coupling parameter , we show the fingering pattern at the same time instance = 3.5, where the larger (stronger surfactant action) results in a slender finger structure and an increase in the number of fingers, which means that the surfactant has a stronger weakening effect on surface tension. For example, when = 0.001, 11 fingers are formed; when = 0.002, 14 fingers are formed; when = 0.003, 16 fingers are formed; and when = 0.004, 19 fingers are formed. These numerical simulations are qualitatively consistent with the physical experiments/numerical simulations in [1,3,5,7,17,29,43].

Uniform injection of a less viscous fluid
The second example is the injection problem similar to the previous example. We still inject low-viscosity fluid, but the injection location is changed to the bottom of the computational domain. The initial conditions are  First, by setting = 0 to remove the surfactant effect, in Figure 7, we plot the profiles of at various times using the grayscale image. When the interfacial layer rises to a certain height, fingering structures (8 fingers) appear, and they finally touch the upper wall (at = 4.5). Then, we apply the surfactant effect by setting = 0.004. In Figure 8a, snapshots of the phase-field variable at different times are drawn, in which we can see the significant influence of the surfactant on the surface tension, that is, more slender fingers (15 fingers) are produced over time. The surfactant concentration , pressure and velocity field u at = 2 are plotted in In Figure 8b Figure 9, where we can see that a large number of 3D fingers grow vertically upwards over time. These fingering patterns in 2D and 3D are well consistent with the numerical/physical experiments in [15,32,33,[35][36][37]46], qualitatively.
We set Ω = [0, 2 ] 2 and implement 2D simulations first. The initial conditions of are still from (4.6) with = 1.3 + 0.01rand( ), and the profile of 0 is plotted in Figure 10a using the grayscale image. The model  Similar to the previous injection example, we still perform the simulation without the influence of surfactant first, and then study the pattern differences caused by surfactants. In Figure 10, for the case of no surfactant ( = 0), snapshots of the phase-field variable at different times are drawn using the grayscale image. We observe that over time, many slender fingers are formed, some of which break into satellite droplets or more branched fingers. We further apply the surfactant effect by setting = 0.004. In Figure 11, the phase-field variable and the concentration variable are plotted at different times. It can be seen that the number of fingers obtained far exceeds the number of fingers without surfactant, which again illustrates the effect of surfactants on weakening the surface tension. In Figure 12, by using various coupling parameter , we show the fingering pattern at the same time instance = 3, where the larger causes a larger number of fingers with more elongated and slender shapes.    where 0 = 0 = , 0 = 2 25 ,˜= 1 + 0.01rand( ). We set the grid size ℎ = 4 800 and the time step size = 0.001, and all other order parameters are still from (4.12). In Figure 13, the isosurfaces of { = 0} at different times are plotted. We observe that that plenty of 3D fingers are formed over time, which are consistent with the 2D simulations qualitatively. Similar dynamics were also observed experimentally/numerically in [2, 3, 6, 10].

Concluding remarks
We design a decoupled fully-discrete finite element scheme to solve the highly complex nonlinear Darcy flow-coupled Cahn-Hilliard phase-field model of binary surfactants. The scheme is constructed based on a combination of a variety of effective numerical methods, including the finite element method, projection method, quadratization method, as well as a new decoupling technique by designing a special type of ODE. While maintaining the second-order time accuracy and linearity, it is natural to obtain unconditional energy stability and practically realized decoupling structure. The detailed actual realization, solvability and stability are given rigorously. By simulating plenty of 2D and 3D numerical examples, including the fingering instability caused by radial/uniform injection and rotating Hele-Shaw cell, we numerically prove the effectiveness of the developed scheme. To the best of the author's knowledge, the developed scheme is the first "ideal" type fully-discrete scheme for the Darcy flow coupling binary surfactant phase-field model.