Thermodynamically Consistent Models for Coupled Bulk and Surface Dynamics

We present a constructive paradigm to derive thermodynamically consistent models coupling the bulk and surface dynamics hierarchically following the generalized Onsager principle. In the model, the bulk and surface thermodynamical variables are allowed to be different and the free energy of the model comprises the bulk, surface, and coupling energy, which can be weakly or strongly non-local. We illustrate the paradigm using a phase field model for binary materials and show that the model includes the existing thermodynamically consistent ones for the binary material system in the literature as special cases. In addition, we present a set of such phase field models for a few selected mobility operators and free energies to show how boundary dynamics impart changes to bulk dynamics and vice verse. As an example, we show numerically how reactive transport on the boundary impacts the dynamics in the bulk using a reactive transport model for binary reactive fluids by adopting a structure-preserving algorithm to solve the model equations in a rectangular domain.


Introduction
Thermodynamically consistent models refer to the ones derived from thermodynamical laws and principles. In particular, they obey the second law of thermodynamics. The generalized Onsager principle (GOP) is a protocol in which the Onsager linear response theory combined with the equilibrium maximum entropy principle is extended to describe thermodynamically consistent models comprising both reversible and irreversible processes [1][2][3][4]. The generalized Onsager principle warrants the second law of thermodynamics in the form of the Clausius-Duhem inequality and has proven to be an effective theoretical framework for developing thermodynamically consistent models at various time and length scales [3,[5][6][7][8]. In the past, the generalized Onsager principle and the equivalent thermodynamical second law has been primarily used to derive dynamical equations in the bulk while boundary conditions are largely trivialized by assuming adiabatic, static boundary conditions, or periodic boundary conditions. In many problems where the scale in the bulk matches the scale in the boundary in terms of dynamics, the dynamics at the boundary can no longer be overlooked or trivialized. For material systems confined in a compact domain with active boundary dynamics, materials properties in the bulk and those on the boundary surface may differ significantly. As a result, different order parameters may need to be introduced to best describe their respective dynamics. This is an issue that has not been paid much attention in the literature. In this study, we revisit these issues and present a general, hierarchical framework to derive thermodynamically consistent models for material systems with coupled dynamics in the bulk and at the boundary using the constructive, generalized Onsager principle [1]. We illustrate the gist of the approach by deriving thermodynamically consistent phase field models and consistent boundary conditions for binary materials in a smooth domain owing to the popularity of phase field models in the literature. However, the theoretical framework is by no means limited to the phase field models.
Phase field modeling is one of the powerful and versatile modeling paradigms for studying multi-phasic materials in domains with complex interfacial geometries and complex interfacial phenomena between distinct phases (immiscible mixtures) or even miscible mixtures [9][10][11]. It is especially useful and effective when handling dynamical phase boundaries in multi-phasic materials involving topological changes compared to other methods, such as front tracking methods, level-set methods, volume-of-fluid methods, etc. [11][12][13]. By design, it is for diffuse interfaces with certain interface thicknesses in which complex interfacial dynamics prevail. A quality phase field model should be able to capture well-known sharp interface conditions (e.g., Gibbs-Thomson condition) in the case of vanishing thickness of the interface in immiscible mixtures [14,15]. This requires one to be mindful when deriving the free energy of the system in the phase field model so that thermodynamical laws and principles are followed faithfully. Notice that the advantage of the phase field model in dealing with diffuse interfaces may also limit its applicability to resolve sharp interfaces in light fluids [16]. The Allen-Cahn and Cahn-Hilliard equations are two classical phase field models used to describe dynamics in the multi-phasic systems and both of them share the same chemical potential if the free energy functionals are the same. However, the Allen-Cahn equation does not preserve the conservation of the bulk volume while the Cahn-Hilliard equation does [17][18][19]. In life science, materials science, and many engineering fields, there is an abundance of multi-phasic material systems with diffuse interfaces, which have kept multi-phase field models popular and practical [12,[20][21][22].
There are quite a number of phase field models in the literature today. However, not all are thermodynamically consistent. Even for the thermodynamically consistent ones, the proper boundary conditions are not well studied until recently. In particular, when one considers surface dynamics, an issue emerges right away that is what is the relation between the phase field variable in the bulk and the one on the boundary? By assuming these two are the same at the boundary, a series of studies have examined the issue of thermodynamically consistency [23][24][25][26][27]. By not assuming that they coincide on the boundary, Knopf et al. derived a set of boundary conditions for the Cahn-Hilliard model (known as the Knopf-Lam model) in 2020 [28] and for non-local models (known as the Knopf-Signori model) in 2021 [29]. The latter studies opened up a new venue for one to examine the thermodynamically consistency in the bulk and surface dynamics holistically. To make the paper more readable, we list some of the existing models in the literature first. In Section 3, we will show that all these existing models are the limiting cases of the general model we proposed.
• Jing-Wang model [24]: We firstly define the free energy of the system as follows where φ is the order parameter, ∇ s is the surface gradient. The corresponding governing equation system is given by where M b and M s are mobilities, α and β are parameters, µ b and µ s are chemical potentials, and µ g is the conjugate variable, which are defined in Section 2.
• Knopf-Lam model [28]: We define the free energy as follows where φ and ψ are two order parameters to describe the materials in the bulk and on the boundary surface, respectively. , σ, and K are model parameters, H(ψ) is a function of ψ. The governing equations in the bulk and on the boundary are given by • Liu-Wu model [27]: By setting H(ψ) = ψ and K → 0, H(ψ) → φ in the Knopf-Lam model, the Liu-Wu model is obtained which are given by • Knopf-Signori model [29]: Once the bulk and surface free energies are non-local as in Section 3.3 below, the nonlocal dynamics are given by In this study, we focus on the demonstration of the paradigm through the derivation of thermodynamically consistent phase field models with coupled surface and bulk dynamics of potentially distinct phase field variables, which include the derivation of the transport equation in the fixed bulk domain, as well as the consistent one on its boundary. We will present the models with free surface domains and boundaries in a sequel. We stress that it is important to study multi-phase material systems using a thermodynamically "correct" model since it not only gives one a comprehensive description of the "correct" physics for the material system but also gives one a well-posed mathematical system to analyze and compute. Speaking of a thermodynamically "correct" model, we insist that the model must be at least thermodynamically consistent and obeys known physical laws. This humble criterion would perhaps disqualify a host of existing phase field models. In addition, we notice that most of the studies on phase field models are concentrated on equations in the bulk with static or periodic boundary conditions at fixed boundaries, where the boundary contributions to thermodynamical consistency are trivialized.
Given the recent technological advances in materials science and engineering and the abundance of natural phenomena where boundaries regulate the internal dynamics, boundaries of a material-confining device can no longer be treated as passive. They can be made with distinctive properties to interact or even control the material within the device, especially, in thin films [30][31][32]. For instance, the newly discovered boundary effect to the existence of blue phases in cholesteric liquid crystals in microscales across a quite large temperature range is one of the prominent examples [33][34][35]. In life science, dynamics in cell membranes play a significant role in determining cell behavior. These require one to derive a model for the material system confined in a compact domain to take into account the potential dynamical contribution from the boundary. There has been a surge in activities in this direction on phase field models and reaction-diffusion models recently [36][37][38][39][40][41][42].
In this paper, the idea is exemplified through a derivation of a phase field model. There in fact exists an underlying paradigm for deriving much wider classes of thermodynamically consistent models in fixed domains with smooth boundaries, whose thermodynamical variables in the bulk may differ from the ones at the boundary and the free energy depends on gradients of the variables up to the second order or strongly non-local. In this approach, we begin with the prescribed total free energy comprising the bulk free energy, surface free energy and free energy that couples the bulk and surface order parameters or phase field variables and calculate the energy dissipation rate of the total free energy. The fact that order parameters in the bulk and ones at the boundary can be different in the continuum level description of the thermodynamical system is common. For example, the order parameter in the bulk in the phase field model is the volume or mass fraction, while the one at the boundary is the area fraction. Their definitions are different so that they may not be the same while evaluated at the boundary. Thus, it is quite common that the surface dynamics and the dynamics in the bulk are described by different thermodynamical variables, dictated by the respective materials properties of the bulk and the surface. In this case, in addition to the bulk and surface energy, an additional energy, defined as the coupling energy, must be introduce to formulate the continuum theory.
In the derivation, we first apply the constructive, generalized Onsager principle to the bulk energy dissipation rate to arrive at constitutive equation system that warrants the bulk energy dissipation. We note that the mobility operator stipulated in the bulk constitutive equation may affect the boundary energy transport. Secondly, after accounting for the additional surface energy dissipation due to the non-local mobility operator in the bulk, we apply the generalized Onsager principle subsequently to the surface energy contribution to the total energy dissipation. This yields the constitutive equation for surface transport at the boundary. Combining the results from both steps, we arrive at coupled bulk and surface transport equation systems that yield a coupled, thermodynamically consistent model. We remark that this paradigm applies to any thermodynamical system. To derive hydrodynamical systems, however, additional conservation laws must be enforced in addition to the constitutive equations, which we will defer to another paper.
The paper is organized as follows. In Section 2, we present the paradigm using a phase field model for a binary system as an example and discuss its various limits. In Section 3, we discuss some special cases, the strongly non-local free energy with dynamic boundary conditions, and a reactive transport system. In Section 4, we demonstrate the effect of dynamic boundary conditions on the solution in the bulk using energy-dissipation rate preserving numerical simulations. We summarize the results in Section 5.

Thermodynamically Consistent Phase Field Models with Consistent Dynamic Boundary Conditions
We illustrate the general framework for deriving transport equations in the bulk and consistent dynamic equations on the boundary for a thermodynamical model that yields a negative energy dissipation or a positive entropy production rate, using a scalar phase field model for a binary material system with free energy consisting of up to second order spatial derivatives of the phase field variable. Then, we elucidate the path for extending it to the more general free energy functional including the non-local free energy for general thermodynamical models. We will show how existing phase field models of this kind are deducible from the general framework. We focus on the derivation in the isothermal case here so that the free energy is the proper thermodynamical potential to work with.

Generalized Onsager Principle
The classical Onsager linear response theory on which the Onsager principle for dissipative systems is based provides a viable way to calculate dissipative forces in relaxation dynamics in an irreversible non-equilibrium process [2,4,43,44]. In the general setting, the linear response theory states that given a chemical potential in an isothermal system, the generalized flux φ t is proportional to the generalized force or chemical potential µ [17][18][19] as follows where M is the mobility operator. For dissipative systems where dynamics are irreversible, the additional Onsager reciprocal relation dictates that M is symmetric; for conservative systems where dynamics are reversible, M is antisymmetric [1]. In general, M comprises the symmetric part M s and antisymmetric part M a : M = M s + M a . We note that when M is a differential operator, such as in the Cahn-Hilliard equation system, the symmetric property of M is also determined by the boundary conditions of the system as well. For a system where inertia is non-negligible and there coexist irreversible and reversible dynamics in the non-equilibrium process, we extend the force balance equation to a generalized Onsager principle [1,24] whereμ is the general chemical potential, ρφ tt represents the inertia force, and ρ is a measure of mass. We next use the generalized Onsager principle to derive the general phase field model along with its consistent boundary conditions for a binary material system. For convenience, we omit( •) over variable (•) in the following.

Models with the Free Energy up to Second Spatial Derivatives
Let the bulk free energy in a fixed material domain Ω be given by where e b is the energy density per unit volume. Especially, the kinetic energy in the bulk to account for the inertia effect in the system is included in e b , which is related to φ t and is usually be chosen as ρ 2 φ 2 t . We consider a binary material system with a boundary that may have its distinctive properties than the bulk and possesses its own surface energy of derivatives up to the second order in space where ψ(x, t) is the surface phase field variable, not necessarily the same as φ(x, t) confined to the surface, e s is the surface energy density per unit area including the kinetic energy on the surface ρ s 2 ψ 2 t and ∇ s is the surface gradient operator over smooth boundary ∂Ω as those in [45,46].
In the following, we will use the volume and area fraction to illustrate why we use two phase variables or order parameters in bulk and surface, respectively. Actually, the phase field variables can include the mass fraction or density fraction in place of the volume fraction. We define the effective volume of the two particles as v 1 and v 2 and assume there are N 1 and N 2 particles, respectively, in a material volume δV = N 1 v 1 + N 2 v 2 . Analogously, we define the effective areas of the two particles' cross-section in the surface as s 1 and s 2 and assume there are N s 1 and N s 2 particles in a material surface δS = N s 1 s 1 + N s 2 s 2 . Then, the bulk volume fraction φ and the surface area fraction ψ for the first material can be defined as follows: Clearly, there are no reasons to believe . For example, we can choose the volume fraction of the first material component as the order parameter in the bulk while choosing the area fraction of the second material component as the order parameter at the boundary. So, φ| ∂Ω and ψ can be two different phase field variables. So we introduce coupling energy in the model to account for the free energy due to the discrepancy between the two physical quantities where φ = φ(x, t)| ∂Ω and e c is the coupling energy density per unit area. For example, the coupling free energy can be used to describe different boundary surface properties, such as attractive or repulsive, between the material components in different regions of the boundary [29]. The total free energy of the system is given by We calculate the time rate of change of the free energy as follows, assuming domain Ω is fixed, where H is the mean curvature of the boundary, n is the unit external normal of ∂Ω, the bulk chemical potential µ b , surface chemical potentials µ s , µ c , and the conjugate variable µ g are given, respectively, by Remark 1. (i) There exist two surface chemical potentials, one corresponding to φ is denoted as µ c , while the other corresponding to ψ is denoted as µ s at the boundary. The surface chemical potential, µ c , includes contributions from the coupling free energy as well as that from the bulk energy confined to the boundary, while µ s only includes contributions from the surface and coupling energies. (ii) The mean curvature shows up in the surface chemical potentials, indicating that curvature of the boundary affects the surface dynamics. (iii) More surface terms can appear if the free energy density function depends on higher order spatial derivatives. We adopt the Einstein notation for tensors here, denote tensor product of vector n and v as nv = n i v j , use one dot · to represent inner product n · v = n i v i and two dots: to represent contraction of two second order tensor A : where n, v are vectors, A, C are second order tensors.

Dynamics in the Bulk
We apply the generalized Onsager principle firstly to the bulk integral in (15) to obtain the transport equation for φ in the bulk, Ω, where M b is the mobility operator and M −1 b is the friction operator which is positive semidefinite to ensure energy dissipation. We limit the mobility operator in the following form in this study where M (1) b ≥ 0 is a scalar function of φ and M (2) b ∈ R 3×3 is a semi-definite positive matrix which can be a function of φ as well. If M b is also a scalar function of φ, then a special case ∇ · M (2) b · ∇ = ∇ · (M (2) b ∇) can be obtained. We note that the derivation applies to a more general mobility operator with high order derivatives as well, which we will not pursue in this study. The presence of spatial derivatives in the mobility indicates the non-local interaction is accounted for in friction operator M −1 b . This is shown in the form of pseudo-differential operators. With this, the energy dissipation rate reduces to We remark that n · M (2) b · ∇µ b is the flux of the volume fraction across the boundary. This physical quantity is determined by the balance between the surface and bulk chemical potential. We next derive consistent boundary conditions for this model.

Dynamics on the Boundary
We recognize that the boundary energy flux density is a quadratic form and then apply the Onsager principle the second time to the energy flux density to establish a dynamical constitutive equation at the boundary: where f m = n · M (2) b · ∇µ b is the inward volume flux, φ t , ψ t , ∇ n φ t are time rate of changes of three quantities, identified as generalized fluxes, and M 4×4 ≥ 0 is the surface mobility operator, a 4 × 4 matrix or second-order tensor. M 4×4 ≥ 0 means that its symmetric part is semi-positive definite. Then, which indicates the system is dissipative. The linear response relation in (20) links the generalized fluxes with the generalized forces which gives a very general boundary surface dynamical system. We next give two examples of the mobility operator to make contact with the existing thermodynamically consistent phase field models in the literature.
Example 1 (First type dynamic boundary conditions (DBCs): a purely dissipative boundary condition). We specify a symmetric mobility operator as follows where M c , M s , M g ≥ 0 are three semi-definite positive operators, α ≥ 0 is a friction coefficient, β, γ, θ are weight coefficients which can take arbitrary values. This constitutive equation establishes a balance between the inward volume flux at the boundary and the generalized chemical potential difference between the bulk and the surface: it assumes the inward volume flux is proportional to the difference between the chemical potential in the bulk and the weighted one at the boundary. When the weighted surface energy is higher than the bulk energy confined to the boundary, the volume flux is inward; otherwise, the volume flux flows outward. In either case, the total energy dissipates since M 4×4 ≥ 0.
The governing equation together with the DBCs in this model is summarized as follows The corresponding energy dissipation rate is given by The surface transport equation of φ can be rewritten into an alternative and more suggestive form on boundary ∂Ω as follows where the inward volume fraction flux is stipulated as follows at the boundary Notice that β, γ, θ can be of any numerical values and if µ g = 0, we set θ = 0 and drop n · ∇φ t from the constitutive equation in (20). (25) indicates that so long as the volume fraction flux is balanced by (26), the relaxation dynamics of the three variables φ, ψ, n · ∇φ at the boundary are dictated by the differences between the corresponding chemical potentials and the weighted inward volume fraction flux, where the weights are delicately tied to the volume fraction flux at the boundary. This gives one a great deal of flexibility to fine-tune the volume fraction flux at the boundary to control the bulk dynamics. So, there is no surprise that this model includes many existing thermodynamically consistent phase field models with DBCs in the literature [23,24,[27][28][29].
Example 2 (Second DBCs: a dissipative and transportive boundary condition). In the second example, we specify the mobility operator in another form with a non-zero antisymmetric component: where the antisymmetric component contributes to transport dynamics at the boundary, in addition to the dissipative part given by the positive semi-definite operator in M 4×4 . The antisymmetric mobility component represents an energy exchange between the bulk and the boundary without inducing any dissipation.
The governing system of equations together with the dynamic boundary conditions in this model is summarized as follows The corresponding energy dissipation rate is given by Notice that the mobility matrix has an antisymmetric component that does not contribute to the energy dissipation. In fact, if we change the antisymmetric part into a symmetric one by negating signs of the non-zero off-diagonal entries in the antisymmetric operator, we recover the model given in the previous example. This set of boundary conditions has the following interpretation: the time rate of changes of the variables at the boundary is proportional to the differences between the surface chemical potentials and weighted bulk chemical potentials. The dynamics are thus different from the previous example and the energy dissipation rate is altered as well.
The two examples of dynamic boundary conditions are derived from two different considerations of mobility operators, as well as effective chemical potentials, which contribute to distinct energy dissipation mechanisms, following the generalized Onsager principle under a unified assumption that the boundary volume fraction flux is proportional to the difference of the bulk chemical potential confined at the boundary and weighted surface energy. In the first case, the time rate of change of the volume fraction is proportional to the difference between the surface chemical potential and the outward volume flux. As a result, the distinctive surface energy dissipation rate is directly linked to the magnitude of the volume fraction flux across the boundary surface. In the second case, the time rate of change of the volume fraction is proportional to the difference between the surface chemical potential and the weighted bulk chemical potential confined to the boundary. Consequentially, the distinctive surface energy dissipation rate is measured by the bulk chemical potential confined to the surface, and the other surface chemical potentials. Two different dissipative mechanisms define two different dynamical models at the boundary. There are more cases that one can elaborate by specifying specific form of operator M 4×4 , which we will not enumerate here. is semi-definite positive. There can be many thermodynamically consistent boundary conditions that are compatible to the given bulk transport equation. However, the dissipative property of the boundary transport equation system depends not only on the mobility operator, but also on the geometry of the boundary as well. Now let us examine the special case where M (2) b = 0, the energy dissipation rate in (19) is given by

Effect of Mobilities in the Bulk and on the Surface
One choice of the surface dynamics is given by which is the limiting case for (23) and (28) To make the energy dissipation rate non-positive, suitable boundary conditions for the unknowns must be taken into consideration. However, this is not our focus here. For the mobility operators at the boundary, we consider the following forms, analogously to the bulk, c , M (2) s and M (2) g are 3 × 3 positive semi-definite matrices. Then, Whether or not the energy dissipation rate at the boundary is non-positive depends on the last term in all three lines in (34), which are linearly proportional to the mean curvature, surface chemical potential, and the corresponding surface flux of each variable.
Having discussed the energy dissipative property of the model, let us investigate how the volume fraction variable evolves dynamically by examining special cases. If s ≥ 0 and M (2) g ≥ 0 are scalar functions of φ, the last term in each line of (34) vanishes and the energy dissipation rate at the boundary is non-positive due to n · ∇ s µ c = n · ∇ s µ s = n · ∇ s µ g = 0. Of course, H = 0 is also a sufficient condition for non-positive energy dissipation rates.
It follows from (23) This indicates a β−weighted volume fraction in the bulk and the area fraction over the surface is conserved under this dynamic boundary condition within the Cahn-Hilliard surface dynamics. Since β is not limited to be positive, this result simply reveals a delicate balance between the volume fraction in the bulk and that on the surface must be maintained during the dynamical process at all times. Similarly, it also follows from (23) that c · ∇ s µ c ]ds, (37) If M This replicates an analogous delicate balance between the bulk volume fraction and the surface one that is preserved in the dynamical process in the Cahn-Hilliard surface dynamics, where γ is that delicate weight parameter with potentially an arbitrary value. Model (23) and (28) give two fairly general phase field models with two different dynamic boundary conditions, where the surface transport equations of phase variables at the boundary set the two models apart. In the first one, the cross-boundary volume fraction flux contributes directly to the energy dissipation on the surface; while in the second, it is the bulk chemical potential limited to the boundary that contributes to the energy dissipation on the surface directly. We next examine various limiting cases (23) and (28), respectively, and show that they reduce to the thermodynamically consistent phase field models in the literature. In fact, when the cross-boundary volume fraction flux (inward or outward flux) vanishes, i.e., α = ∞ and γ = β = θ = 0, the two types of dynamic boundary conditions are identical.

Reduction to Limiting Cases
We examine three limiting cases of the model to make contact with the existing thermodynamically consistent phase field models and then present a new reactive transport limiting model with reactive dynamic boundary conditions.

The Jing-Wang Model [24]
We choose e c = 1 2K (φ − ψ) 2 , where K > 0 is a penalizing parameter, and define a symmetric mobility operator as follows The governing equation system together with the boundary conditions in this model is given as follows When K → 0, ψ = φ| ∂Ω and the governing equation system reduces to The corresponding energy dissipation rate is given by Setting θ = 0, this system reduces to the general model with the first type dynamic boundary condition derived in [24] by the authors. Similarly, we deduce the general model with the second type dynamic boundary condition in [24] by selecting the mobility operator corresponding to the second type boundary condition alluded to in the previous section.

The Knopf-Lam Model [28] and the Liu-Wu Model [27]
We define the free energy densities as follows The corresponding energy dissipation rate is given by The transport equation in the bulk is given by We choose such a relationship between the general fluxes and general forces at the surface as follows where µ c = n · ∇φ − 1 Then the energy dissipation rate is given by In the limit κ 1 → ∞, κ 2 → 0, the governing equations in the bulk and on the boundary are given by The general model reduces to the Knopf-Lam model in [28]. Moreover, when H(ψ) = ψ and K → 0, H(ψ) → φ. The energy dissipation rate in (44) reduces to The thermodynamically consistent model reduces to the Liu-Wu model in [27]

Non-Local Models including the Knopf-Signori Model [29]
We consider phase field models with a non-local bulk free energy [47,48] given by where J( x ) is the interaction kernel and f is the free energy density for the bulk. This form of free energy is perhaps more generic than the one that depends on spatial derivatives of the phase variable. We call this a strongly non-local free energy. The bulk chemical potential is given by where a(x) = Ω J( x − y )dy. Likewise, we consider the surface energy and coupling energy given, respectively, by where g(φ) is the surface energy density per area and h(φ, ψ) is the coupling energy density per area. The surface chemical potentials are obtained as follows where a S (x) = ∂Ω K( x − y )ds y , a c (x) = ∂Ω L( x − y )ds y . The total free energy is then given by We calculate the time rate of change of the free energy as follows We apply the generalized Onsager principle to the bulk term to arrive at where M b is the mobility operator. For We specify boundary dynamic equations by using the generalized Onsager principle as follows where M 3×3 is the non-negative definite boundary mobility operator. The energy dissipation rate is given by When setting the general model reduces to the non-local model with dynamic boundary conditions in [29] in the limit, κ 3 → ∞.

Reactive Transport Equation in a Binary Polymeric System
Finally, we present a new phase field model for a binary reactive polymeric system and assume the two polymer chains are long enough so that they can be viewed as ideal chains. We also assume there are reactions in both the bulk and the surface described as follows, where k + i , k − i , i = b, s are the forward and backward reaction rates in the bulk and on the surface, respectively. The reaction rates in the bulk and surface may be different due to different pressure, exposure to catalysts, enzymes, light, etc.
The corresponding bulk free energy of the reactive system is given by and the surface and coupling free energies are, respectively, given by where φ, ψ are two order parameters to describe the volume fraction and the area fraction of the polymer segments in the bulk and the surface, respectively, b is the scaled Kuhn length of the polymer segment, Q 1 , Q 2 are the dimensionless partition functions for the polymer segments, n i , i = 1, 2, are the polymerization index for the two polymer chains, and χ, χ s measure the mixing energy due to the volume exclusion effects between segments in the bulk and surface, respectively.
Using the generalized Onsager principle, we obtain the governing equation with the dynamic boundary condition of the first kind: where k b and k s are two dimensionless parameters measuring the reaction rates in the bulk and the surface, respectively. M s , readers are referred to our forthcoming paper [49]. The total energy of the system is dissipative. In the following, we briefly derive the model.
The dynamical description for the reactive transport in the binary system in bulk can be written as where Φ = (φ 1 , φ 2 ) T is the volume fraction vector, φ i , i = 1, 2 represent the volume fraction for the ith material component, Λ = Λ (1) − Λ (2) ∇ 2 , Λ (1) is the mobility for reaction and −Λ (2) ∇ 2 is the mobility for transport. The bulk free energy is given by In this model, we assume n 1 = n 2 = n. For reaction dynamics, the mobility operator is given by For transport dynamics, we set The bulk transport equation is then given by ∂ ∂t Substituting 1 into the governing equation in the bulk, and using incompressible condition φ 1 + φ 2 = 1, we obtain L = (λ 11 +λ 12 )µ 1 +(λ 12 +λ 22 )µ 2 λ 11 +2λ 12 +λ 22 where the second equation is in the form: For the reactive binary system, we could introduce different order parameters besides φ = φ 1 . For example, φ = φ 2 or φ = φ 1 − φ 2 in the bulk are all acceptable. Analogously, we obtain the mobilities/dynamics on the boundary. The choice of the order parameters on the boundary may be different from the one in the bulk. For example, φ s = (φ s,1 − φ s,2 )| ∂Ω , ψ s = ψ s,1 , where ψ s,1,2 are surface area fraction for the two polymeric components, respectively.
Recall the lattice model for polymers melts/solutions/blends [50], where one divides the domain into small lattices and uses the outermost lattices in the domain to represent the boundary. φ(x) at these outermost lattice can be identified as cφ(s), s ∈ ∂Ω, where c is a parameter to make sure φ(s) represents the area fraction while φ(x) is the volume fraction. We also add one layer of lattices at the outside of the domain and introduce ψ(s) to this layer. Although φ(s) and ψ(s) are in different layers, both of them can be viewed as representing quantities "on the surface", when the thickness of the outermost lattices in the domain and outside is small enough. For any lattice we discussed above, we assume they are infinitesimal in macroscopic scale so that the "area fraction" and the volume fraction in the infinite thin layer limit can indeed be different.
There have been quite a number of papers investigating the reaction-diffusion phenomenon with dynamic boundary conditions in the last two decades in both experiments and modeling [51][52][53][54][55]. However, some of them are not thermodynamically consistent and the concept of the free energy used in the models is not clearly defined. In this paper, we use a simple binary reactive transport to illustrate how thermodynamically consistent models can be derived following the generalized Onsager principle. For more general reactive transport systems, readers are referred to our forthcoming paper [49].

Numerical Results for a Binary Reactive System
We take the binary reactive transport system as our numerical example to showcase how boundary dynamics impact bulk ones in a compact domain in 2-D. In the rectangular domain, the four sides are labeled as Γ 1 , Γ 2 , Γ 3 , and Γ 4 , respectively. We use the same free energies as those in the reactive transport system for the binary polymer system presented in section Section 3.4. The coupled dynamic equations in the bulk and at the boundary are given by the following where M ln(Q 2 ψ)−ln(Q 1 (1−ψ)) , k b , k s , Q 1 , Q 2 are assumed dimensionless parameters. In this model, we allow dynamic, reactive boundary conditions at Γ 1 , and homogeneous Neumann boundary conditions at the other sides of the boundary. We adopt H(ψ) = ψ in the simulation and note that k b and k s parameterize the reactive rates in the bulk and boundary, respectively.
In the numerical treatment of the coupled PDE system, we use the energy quadratization technique coupled with the Crank-Nicolson method in time with time step δt = 1 × 10 −5 , and the second order finite difference method on staggered grids in space Ω = [0, 1] 2 with spatial mesh size 1/256 in both x and y direction to obtain a thermodynamically consistent numerical algorithm. The numerical algorithm guarantees that the total energy dissipates in time. For more details about the algorithm please refer to our previous papers [7,56,57]. We use the following initial condition in the simulations Note that the Robin boundary condition on Γ 1 reduces to the homogeneous Neumann boundary condition when α → ∞.
The numerical examples presented next are intended to showcase the effect of boundary reactive dynamics on the bulk dynamics at a set of parameter values. To explore the solution behavior of the model, a comprehensive study is needed which we will defer to another study. In the first simulation, we show the time evolution of spinodal decomposition of the binary polymer system with homogeneous Neumann boundary conditions at all four sides in Figure 1a-c. In the second simulation, we set the homogeneous Neumann boundary condition at Γ i , i = 2, 3, 4 and the dynamic boundary condition at Γ 1 with α = β = γ = 1, k s = 0, η = 100 and present the solution in Figure 1d-f. In the third one, we use α = β = γ = 1, k s = 5 × 10 −2 , η = 100 with the same set of boundary conditions as in the second simulation and show the solution in Figure 1g-i. The first case decouples the bulk dynamics from the surface ones. The second case includes the coupling between the surface dynamics and the bulk dynamics without surface reaction. The third case includes reactive surface dynamics in dynamic coupling. Comparing the three simulations, we find that the dynamic boundary condition on the surface affects the bulk dynamics noticeably due to the non-negligible volume fraction flux (inward or outward) across the surface. The enhanced surface reaction on the boundary accelerate the merging of the droplets in the bulk in Figure 1g-i.
In this reactive-transport model, the volume fraction flux dictates boundary dynamics, as well as bulk ones near the boundary. We then depict the time evolution of the bulk volume fraction and bulk free energy in the three cases, respectively, in Figure 2. We notice that the bulk volume fraction with the homogeneous Neumann boundary condition is constant with time as expected, whereas the bulk volume fractions in the other two cases increase with time. The increasing bulk volume fractions in case 2 and case 3 lead to two different patterns in Figure 1d,i. In addition, the bulk free energies with dynamic boundary conditions dissipate faster than that with the homogeneous Neumann boundary condition. The enhanced reaction process accelerates the energy dissipation rate in the bulk.  Figure 1. It is obvious that the bulk volume fraction with homogeneous Neumann boundary conditions keeps as a constant. However, no matter whether or not there is the reaction on the boundary, the bulk volume fraction increases with time when the boundary conditions are dynamic. Compared with the DBCs without reaction, the bulk volume fraction in the DBCs with reaction increases faster than that without reaction. The bulk free energy in each case decays in time as expected. Due to the existence of chemical reaction, the bulk free energy in the third simulation (with surface reaction) decreases faster than that in the second one (without surface reaction).
Notice that the coupling energy at the boundary is a new physical quantity in this formulation and an intriguing one that determines the difference between φ(s) and ψ(s) and possibly influences transport dynamics in the bulk. Next, we investigate the effect of the coupling free energy together with the boundary reactive dynamics on the bulk dynamics at the presence of the inward/outward flux. The coupling free energy density adopted in the simulations is given by which measures the deviation between the two surface-bound variables. We present the case where the inward/outward flux is present (α = β = γ = 1) in Figure 3. If η = 1 × 10 2 , representing a strong coupling, ψ(s) can significantly affect φ(s) as shown in Figure 3a.
When η = 0 in Figure 3b, the coupling energy vanishes so does the boundary dynamic effect of ψ(s). We also show the effect of surface reaction rate k s on bulk patterns in Figure 3b,d with η = 1 × 10 2 and η = 0, respectively. Figure 3a,b plot the numerical results without the surface reaction (i.e., k s = 0), while Figure 3c,d do the ones with the surface reaction (k s = 5 × 10 −2 ). Comparing Figure 2a,c (or Figure 2b,d), we find that the magnitude of the chemical reaction at the boundary, measured by the reactive rate parameter, affects the bulk dynamics near the reactive boundary significantly. We then examine time evolution of the bulk volume fraction for these four test cases in Figure 4a. We observe that the bulk volume fraction increases with time in test 1, 3, and 4, which explain the mechanism of merging droplets or lamellar in Figure 3a,c,d. In Figure 4a, the bulk volume fraction in test 2 decreases with time, which indicates the bulk volume fraction flows from the bulk to the boundary so that there are less droplets near the boundary in Figure 3b. We also examine time evolution of the total free energies in the tests in Figure 4b, which show that our numerical algorithm used is energy-dissipative and warrant the second law of thermodynamics numerically.   Figure 2. We find the bulk volume fractions increase with time in tests 1, 3, and 4, which are related to the merging droplets or lamellar in Figure 3a,b,d. Since there is a leakage of bulk volume fraction in test 2, we do not find any large droplets or lamellar in Figure 3b. The total free energy in each test decays in time as expected.
To explore the effects of η and k s further, we employ another set of α, β, and γ values next. Firstly, we set α = 0.1, β = γ = 10 and plot the snapshots of the numerical solution of φ at T = 8 with respect to several combinations of (η, k s ) in Figure 5. Secondly, we set α = 10, β = γ = 0.1 and plot the snapshots of the numerical solution of φ at T = 8 with respect to the same combination of (η, k s ) in Figure 6. There are almost no differences among the snapshots in Figure 5 with respect to different η and k s . Coincidentally, there are also no significant differences among the snapshots in Figure 6 with respect to different η and k s either. Finally, we show the time evolution of the bulk volume fractions in Figure 7. The results show the bulk volume fractions in tests 5-8 increase with time while the bulk volume fractions in tests 9-12 decrease with time. The inward/outward flux is influenced by the magnitudes of α, β, γ, and k s , η, which shows a slight effect on the pattern before T = 8 in the eight tests. . We find the bulk volume fractions increase with time in tests 5-8, which are related to the merging droplets or lamellar in Figure 5. Since there are leakages of bulk volume fractions in tests 9-12, we find there are lamellar formed by the other phase in Figure 6.

Conclusions
Guided by the generalized Onsager principle, we have demonstrated a paradigm for deriving coupled, thermodynamically consistent dynamic models for the bulk and the boundary systematically using a phase field model for a binary fluid system as an example. All the existing, thermodynamically consistent phase field models for binary fluid systems accounting for dynamic boundary conditions are shown to be limits of this general phase field model. Then, we use a structure-preserving algorithm that we developed to numerically solve the phase field model for a binary fluid model with a reactive boundary to illustrate how boundary dynamics, including reactive dynamics and coupling energy, impact the bulk dynamics. This general framework anatomize the thermodynamical models with respect to both the bulk and boundary equations, laying a groundwork for analysing and designing efficient structure-preserving numerical approximations to these models in the future. Funding: Xiaobo Jing's research is partially supported by NSFC awards #12147165 and NSAF-U1930402 to CSRC. Qi Wang's research is partially supported by NSF awards OIA-1655740.