Consistent and flexible thermodynamics in atmospheric models using internal energy as a thermodynamic potential. Part II: Non‐equilibrium regime

In numerical models of the atmosphere, the non‐equilibrium thermodynamic processes involving moisture are not always treated consistently—possibly leading to inconsistencies and errors in the energy budget. Therefore, a more consistent formulation of (moist) thermodynamics is important, for short‐timescale weather models and long‐timescale climate models. In Part I of this work, we derived a thermodynamically consistent framework, describing condensation, evaporation, freezing, and melting of cloud droplets, in which all thermodynamic quantities of interest were derived from an internal energy potential and with the moist thermodynamics coupled to a 2D semi‐implicit semi‐Lagrangian dynamical core. While this framework was primed to express non‐equilibrium processes, it was solved for the equilibrium regime only. Here, we follow the methods in Part I, but with the expression for the non‐equilibrium processes “turned on”, for example, allowing freezing of supercooled water or evaporation into subsaturated air. To implement the proposed approach, it is necessary to translate conventional atmospheric microphysics expressions for transfer rates of matter and entropy in and around a cloud droplet into the formalism of non‐equilibrium thermodynamics. This procedure is first derived for some simple idealised cases, beginning with liquid droplet growth by vapour diffusion, and proceeding to more complex three‐phase cases. To demonstrate the approach, we then simulate some idealised cloudy thermals, comparing the equilibrium and non‐equilibrium regimes—finding a robust decrease in the vertical velocity in the non‐equilibrium regime, as expected. Thus, this work demonstrates the feasibility of building a numerical model that includes a framework for modelling the moist non‐equilibrium thermodynamics of an atmospheric system consistently and provides a step towards this type of more consistent atmospheric modelling.


INTRODUCTION
There is growing interest within the atmospheric modelling community in improving the formulation of weather and climate models to permit a consistent treatment of the thermodynamics and to enable a closed energy budget to be diagnosed (e.g., Lauritzen et al., 2018;Staniforth and White, 2019;Eldred et al., 2022;Lauritzen et al., 2022;Staniforth, 2022). In this article we build on our previous work exploring the feasibility of using a thermodynamic potential to obtain a consistent formulation of the moist thermodynamics, extending it to include non-equilibrium processes.
In Thuburn (2017), a 2D semi-implicit semi-Lagrangian model of atmospheric flow was developed in which the equation of state and other thermodynamic properties of moist air were expressed in terms of the Gibbs potential and solved numerically for a number of test cases. This method has the following benefits.
• Consistency-Any approximations made remain consistent with the laws of thermodynamics. • Self-consistency-All branches of numerical code share the same thermodynamic potential; this prevents different parts of the code making approximations to the thermodynamics in different ways. • Flexibility-Modifying or swapping out the equation of state is more straightforward.
For more details and examples see Thuburn (2017). We note that, while the idea of using a thermodynamic potential for this purpose may be novel in atmospheric physics, this idea has been used in the field of oceanography for some time. Guidance on the development of this idea may be found in the Thermodynamic Equation of Seawater 2010 (TEOS-10) 1 manual (IOC and IAPSO, 2010) and Feistel (2003;, Feistel et al. (2010). Bowen and Thuburn (2022) (Part I) extended the model developed previously in Thuburn (2017). The first aim of these extensions was to include ice. However, using the Gibbs potential proved unsuitable because of an ambiguity at the triple point. A way to fix this was instead to derive all thermodynamic quantities of interest from an internal energy potential.
The second aim was to construct a framework that was primed to express non-equilibrium processes in a thermodynamically consistent way. Indeed, this was possible, but required the phenomenological equations to be expressed in terms of a resistivity matrix (see Section 2). In this way, it was possible to recover the 1 http://www.teos-10.org equilibrium regime smoothly by setting the resistivities to zero.
An argument can be made (e.g., Guggenheim (1933)) that all "natural" thermodynamic processes arise through disequilibrium and so must evolve irreversibly, producing entropy. Reversible evolution through a sequence of equilibrium states is an idealisation, and can only ever be an approximation to reality. Nevertheless, in our atmospheric system, when thermodynamic processes occur faster than the fluid transport or external forcing, the equilibrium regime is an excellent approximation. In this regime, all of the thermodynamic processes (e.g., phase changes and heat transport) within every material parcel occur instantaneously-implying that no internal entropy is generated. However, in a realistic atmospheric system some processes happen more slowly than the fluid transport (for example, freezing of supercooled water), so we need non-equilibrium equations to capture these slow processes.
Considering this previous work, our aims here are twofold.
1. Translate a simple (two-phase) microphysical case from the typical expression in (standard) microphysics into standard non-equilibrium thermodynamics formalism. However, the task is non-trivial for a number of reasons, as outlined in Section 1.2. The final form must be compatible with the framework developed in Bowen and Thuburn (2022). A similar procedure must then be repeated for the other two-phase case and more complicated three-phase microphysical cases. 2. Confirm that the resulting thermodynamic formulation can be coupled to a typical numerical method used for atmospheric dynamics, and examine what effect including the non-equilibrium processes has on the model dynamics.
Before addressing each of these aims, we should discuss the importance of the non-equilibrium processes and the various difficulties involved in modelling atmospheric systems using the framework of Bowen and Thuburn (2022).

Non-equilibrium processes, entropy, and dynamics
Non-equilibrium processes produce entropy. Although entropy production rates have been investigated and outlined at the theoretical level for some time (Onsager, 1931a;de Groot and Mazur, 1984), this has not necessarily extended into concrete applications beyond small scales. In modelling the climate, for example, the entropy is not always given proper attention. Neglecting proper treatment of the entropy can lead to inconsistencies and possible errors in the global energy budget (e.g., Johnson, 1997;Woollings and Thuburn, 2006;Thuburn, 2017).
There have been a number of attempts to estimate the different entropy production rates on a global scale. Many works quantify the average production rates in different parts of the atmosphere by considering the entropy flux in/out of the system. This has been used to gather estimates for the planetary, transfer, and material rates in works such as Lesins, 1990;Stephens and O'Brien, 1993;Li and Chylek, 1994;Pelkowski, 1994;Wu and Liu, 2010;Lembo et al., 2019;Gibbins and Haigh, 2020. The entropy production has also been investigated in some general circulation models (Goody, 2000;Lucarini et al., 2011;Lucarini et al., 2014;Laliberte et al., 2015), with particular focus on poleward transport of energy. The difficulties involved in its description have been thoroughly explored in Lucarini et al., 2014;Bannon and Lee, 2017;Gibbins and Haigh, 2020. Detailed estimates of the material entropy production have also been attempted on the smaller scales, as in Pauluis and Held, 2002a;2002b;Romps, 2008;Volk and Pauluis, 2010;Pauluis and Dias, 2012;Singh and O'Gorman, 2016, though there are some difficulties in doing so. For the material entropy production to be accurate, every non-equilibrium process needs to be described accurately. This includes processes such as heat transport, turbulence, frictional dissipation, and irreversible phase changes-to name but a few.
In this article, we are particularly interested in the non-equilibrium processes involving moisture and how these affect the dynamics. Moisture is involved in some of the most significant contributions to the entropy production (Pauluis and Held, 2002a;2002b); the two dominant sources of entropy production are the irreversible phase changes of water and diffusion of water vapour. Water is also found to have a strong interaction with the frictional dissipation entropy terms, a significant contribution coming from hydrometeor fallout (Pauluis et al., 2000;Rennó, 2001;Makarieva et al., 2013). Interestingly, the entropy production associated with moisture turns out to be more important than dry frictional dissipation or turbulent dissipation.
The non-equilibrium microphysical processes (and the implicit entropy production) define the mass and energy transfers between the various species in an atmospheric parcel (dry air, vapour, liquid, ice). The finite timescale of some of these processes can have a significant impact on the dynamics. For example, the finite timescale for freezing of cloud droplets may lead to delayed release of latent heat and an attenuated cloud ascent (e.g., Pauluis and Held, 2002b;Kurowski et al., 2018); the finite timescale for evaporation of rain below deep precipitating cumulus clouds may lead to the production of cold pools that affect the subsequent development of convection strongly (e.g., Williams and Renno (1993)).
Neglect of non-equilibrium effects can also impact the consistency of model diagnostics. For example, the finite timescale for evaporation of hydrometeors into subsaturated air and for heat transfer between hydrometeors and the surrounding air means that the temperature of precipitation when it reaches the surface is different from the temperature at which condensation occurred and also different from the surface air temperature (e.g., Anderson et al. (1998)). Neglecting the heat flux associated with these temperature differences can make it difficult to diagnose a closed energy budget in weather and climate models (e.g., Lauritzen et al. (2022)).
There are of course other non-equilibrium effects to consider that influence the dynamics, such as those associated with frictional dissipation and the seeding of clouds. The model developed here does not include all the processes mentioned above, but the following approach does permit their inclusion.

Discussion of the difficulties involved
In Bowen and Thuburn (2022), the suggested framework to capture the non-equilibrium processes (based upon the formalism of Onsager (1931a)) involved defining a number of thermodynamic "forces" F, and conjugate thermodynamic "transfer rates" 2 J. The forces and transfer rates are conjugate in the sense that F T J gives the material rate of (internal) entropy production. The transfer rates are related to the forces by phenomenological equations, that is, Here, L is the conductivity matrix, while R = L −1 is the resistivity matrix. We have some freedom over which transfer rates we include in J: for example, energy or entropy transfer rates. Any choice for J determines the conjugate forces F, such that F T J remains the internal entropy production. The entries in L must ultimately be determined by physical measurement. The second law of thermodynamics implies that L and R are positive semidefinite. The Onsager reciprocal relations (Onsager, 1931a) imply that L and R are symmetric.
2 In non-equilibrium thermodynamics literature, the "transfer rates" are more commonly known as "fluxes". To avoid confusion with movement of material, we use the term "transfer rates" instead.
Conveniently, setting R = 0 reduces the system to the equilibrium regime, implying that the internal thermodynamic processes equilibrate instantaneously. Essentially, this enforces a balance of thermodynamic forces, that is, For the non-equilibrium regime, R is non-zero, giving a finite rate of approach toward equilibrium. The information needed to specify the coefficients R is the same information used in conventional modelling of cloud microphysics. Our main task in Section 3 and Appendix B will be to derive the coefficients R for diffusional growth and evaporation, and for freezing and melting of cloud drops, from the conventional formulas. This effort is made non-trivial for a number of reasons.
• Conventional atmospheric formulas are most closely related to conductivities, but we require resistivities to be consistent with the framework of Bowen and Thuburn (2022). • R must be symmetric.
• Some parts of the system are in equilibrium (e.g., dry air and water vapour have the same temperature), hence some portions of R must contain zeros. • For consistency with our 2D model dynamical core, we want to predict specific humidities and specific entropies for the different water phases. Thus, eventually we want J to comprise mass and entropy transfer rates (with conjugate forces F and resistivity matrix R). However, to obtain R it is convenient first to determine the mass and energy fluxesJ (with conjugate forcesL and resistivity matrixR); a transformation is then required to obtain the mass and entropy transfer rates in order to be consistent with Bowen and Thuburn (2022). • For the three-phase microphysical case of a liquid core surrounded by an ice shell, we need to model the thermodynamic profile across the ice shell in order to determine R. • We want the system to behave smoothly at transitions between different microphysical cases.
To give some perspective on one of the difficulties involved, we consider a motivating example of liquid particle growth by vapour diffusion in a moist atmospheric parcel. A typical practice to determine the particle growth is to use differences between the vapour density "far enough" away from the particle surface (far-field) and the equilibrium vapour density at the particle surface (Rogers and Yau, 1989;Mason, 2010;Pruppacher and Klett, 2011). To represent this, we often see something of the form (assuming multiple identical particles) where the growth rate of the liquid mass m l has dependence on the following: r s , particle radius; N, number density of particles; V, volume of air parcel; D v , vapour diffusivity; ∞ , vapour density in the far field; v e , equilibrium vapour density at the particle surface. The thermodynamic force driving the diffusive mass growth is simply this difference in vapour densities. The driving thermodynamic forces used in Bowen and Thuburn (2022) (within F) are Δg and ΔT, where g is the specific Gibbs potential and T the temperature. We may find an analogous expression to Equation (2) (ignoring cross-effects) to be where q l is the mass fraction of liquid, L v vv is some conductivity coefficient, and T d is the temperature of the dry air. To relate L v vv to D v , it is reasonable to suggest that we could make a comparison of Equation (3) with Equation (2). However, this is not as straightforward as it first seems. While the force Δg is related to Δ , there is some non-obvious transformation or approximation that must be performed to match the thermodynamic forces properly. This gives some indication that the determination of some of the conductivity coefficients requires some careful thought.

Summary of sections
In this article we further the development of a consistent non-equilibrium formulation for moist thermodynamics, based on the internal energy potential, for use in atmospheric models. We illustrate how the resistivity matrix can be calculated for some typical microphysical processes, and we demonstrate the feasibility of the approach by coupling the thermodynamics to a typical (albeit two-dimensional) atmospheric dynamical core.
To give a necessary overview of the prior theory, Section 2 reviews the important features of the non-equilibrium framework we will be using. Notes on the coupling of the non-equilibrium system to a semi-implicit semi-Lagrangian dynamical core are given. We stress that this section is intended to provide a brief overview of the theory only; readers are advised to refer to Bowen and Thuburn (2022) for the prior non-equilibrium theory.
Section 3 develops a procedure to link formulas describing an idealised microphysical case to the expressions provided by the non-equilibrium framework. This allows the coefficients of the resistivity matrix (required for numerical solution) to be related to known quantities. This is by no means a trivial process, and as such will be presented fairly explicitly. Part of the motivation here is that more microphysical cases may then be built with increasing complexity in further work. Section 4 presents some example results comparing idealised cloudy thermals for equilibrium and non-equilibrium regimes. The initial conditions are given, and simulations of the regimes are compared with one another. The results demonstrate the feasibility of our proposed approach for coupling non-equilibrium thermodynamics to a typical dynamical core, and confirm that the non-equilibrium processes affect the dynamics in the expected way.
Section 5 reviews the aims of this article and how successfully each has been achieved. Some comments on the numerics are made, comparing the computational expense of this model with Bowen and Thuburn (2022). The assumptions/limitations of the microphysical representation are then collated and minor improvements suggested. Lastly, some general extensions are suggested for the next steps in developing the approach. This section concludes our work.
Appendix A lists the various internal energy potentials used. The derivations of some more complicated microphysics cases are detailed in Appendix B. Appendix C contains calculation of cloud condensate particle radii for both pure and mixed-phase condensate. For reference, the thermodynamic constants used are listed in Table 1.

Constant
Description Value Specific heat capacity of dry air, at constant volume Specific heat capacity of vapour, at constant volume Specific latent heat of vaporisation, at T = 0 and p = 0

NON-EQUILIBRIUM FRAMEWORK AND NUMERICAL SOLVER
In our companion article, Part I (Bowen and Thuburn [2022]), a framework was developed to represent a thermodynamically consistent non-equilibrium system-although it was solved for the equilibrium regime only. It is helpful to highlight some of the important points in order to understand this article on a relatively separate basis. We will first review the derivation of the non-equilibrium framework, and then review the process of solving it numerically. For the full details, see Bowen and Thuburn (2022).

Non-equilibrium framework
As in Part I, we want to build a non-equilibrium framework that provides prognostic equations for the individual water mass fractions q i and specific entropies i . To build this framework, we start by assuming conservation of water and specific internal energy following each Lagrangian parcel, that is, for = {v, l, f}, = {d, v, l, f}; the species of dry air, vapour, liquid, and ice are denoted by d, v, l, f, respectively. Here, using superscripts to indicate the species of concern, q i is the mass fraction, e i is the specific internal energy, p i is the pressure, i is the specific volume, andQ (x i ) is the external heating/cooling. Henceforth, we will omit the "specific" prefix; assume quantities to be in specific form unless stated otherwise. As well as water conservation, we also observe dry mass conservation, where  d represents the production of dry air. After some manipulation (see Bowen and Thuburn (2022) for details), it is possible to use Equation (5), along with Equations (4) and (6), to find D Dt where  i = ∑ j  ij is a mass transfer rate and  i = ∑ j  ij is an energy transfer rate, for i, j ∈ . The mass transfer rate  i is the incoming mass transfer rate minus the outgoing masstransfer rate to species i; the energy transfer rate  i is the incoming energy minus the outgoing energy to species i. This implies Next, we require equations to determine  i and  i . These should be given by a set of phenomenological equations of the form J = LF as in the writings of (Onsager, 1931a;1931b;de Groot and Mazur, 1984). As mentioned previously, the form J = LF is not unique; the selection of J, L, F is chosen to be convenient for the problem.
To identify the forces F conjugate to the transfer rates where Q (x) is the total external entropy source. The change in entropy is simply a product between the internal thermodynamic forces and transfer rates, as well as some external entropy flux. The internal contribution F T J must be positive as per the second law, but the external contribution Q (x) may take either sign. We may then choose J and (using Equation (10)) find F. Thus, we consider the total entropy evolution obtained from Equation (8), This can be rewritten in terms of differences of thermodynamic properties, using Equation (9), The differences (in parentheses on the right-hand side (RHS) of Equation (12)) clearly correspond to thermodynamic forcesL for the mass and energy fluxesJ = Since we prefer to work with mass and entropy-transfer rates instead of mass and energy-transfer rates, we rearrange Equation (12) as where Δ ij (X) = X i − X j , for a general thermodynamic property X, and we have substituted the entropy transfer rate Note that, while Equation (9) is true, Hence, we have found the F conjugate to J. Then, it is convenient to write Equations (7) and (8) where (with the evolution of q v and q d d given implicitly by mass and energy conservation). Finally, we may take Equation (19) and substitute J = LF to find a neatly packaged form, Alongside Equation (10), this closes the thermodynamics. Anticipating that for the full equilibrium regime we will require L −1 = R = 0, Equation (22) is rewritten as Together, Equations (10) and (23) provide two convenient expressions to work with. Further argument on the preference of using R over L can be found in Bowen and Thuburn (2022).

Coupling to fluid dynamics, and non-negativity constraints
One of our aims is to demonstrate a numerical solution to the above thermodynamic equations coupled to the fluid dynamical equations. The governing equations for a 2D system (identical to the one in Bowen and Thuburn (2022)) are Here, u = (u, w) is the 2D velocity vector, Φ is the geopotential, and q is the total specific water. Subscripts in Equation (30) indicate a partial derivative. Coriolis terms are neglected, as are forcing and dissipation terms.
This system comprises 15 equations to be solved for 15 unknowns. These unknowns are the following: specific volumes of gaseous species d , v ; specific entropies d , v , l , f ; Lagrange multipliers associated with liquid and ice mass fractions l , f ; water mass fractions q v , q l , q f ; pressure p; density ; horizontal velocity u; vertical velocity w. An inequality-constrained method is used to prevent any negative mass fractions. Equation (29) and the term (̇XC) T in Equation (28) serve to constrain the thermodynamics, ensuring non-negative mass fractions. In these parts of the system, C ( l , f ) T ; C represents the constraints to apply, represents the switches that turn on/off the constraints, and represents the Lagrange multipliers. The underlying principles of precisely how the non-negativity constraints are applied are detailed more fully in Bowen and Thuburn (2022), with similarities to a method found in Gupta and Hauser (2007).

Numerical solution
In summary, the procedure we follow to prepare the system for numerical solution is the following: 1. discretize the continuous system of equations in a semi-implicit semi-Lagrangian scheme, which provides a set of non-linear spatially coupled equations to be solved at each time step; 2. solve this system using a quasi-Newton method, which involves linearizing the system of equations; 3. reduce the linearized equations to a typical Helmholtz-type problem, in terms of the pressure increments p ′ .
Once this is done, we seek to solve numerically. At each model time-step, we take a number (typically four or five) of Newton iterations. At each Newton iteration, we do the following: 1. use the residuals in the governing equations to construct the right hand side (RHS) of the Helmholtz problem; 2. solve the Helmholtz problem for the pressure increments p ′ ; 3. back-substitute for the increments of the other unknowns, taking fractional increments where necessary to satisfy constraints-this serves to reduce the residuals of the discretized system; 4. repeat iterations of the constrained iterative solver until desired convergence.
For full guidance on the procedure, see Bowen and Thuburn (2022).

THE MICROPHYSICS
In Bowen and Thuburn (2022), the 2D system was solved numerically for the equilibrium regime only, even though the system was set up with a full non-equilibrium framework. Here, we will develop the theory further to capture some non-equilibrium processes properly. Essentially, this section will demonstrate the following: 1. expression of conventional microphysics formulas in the form expected in the theory of non-equilibrium thermodynamics; 2. extraction of some typical L or R matrix entries to demonstrate the coupling of consistent non-equilibrium thermodynamics to a semi-implicit semi-Lagrangian dynamical core.
For simplicity, we will consider six basic microphysical cases.
• Liquid condensate: The vapour in the gaseous mixture may condense on to suspended liquid particles, or liquid may evaporate into the gaseous mixture. • Ice condensate: Vapour may deposit directly on to ice particles, or ice may sublimate into the gaseous mixture. • Ice-core, liquid-shelled condensate: Ice particles melt from the surface inwards, with a liquid shell fully encapsulating the ice core.

• Ice-core, circulating liquid-shelled condensate:
This is a special case of the ice-core liquid-shelled condensate. At a certain thickness, internal circulations develop within the liquid shells, breaking down thermodynamic gradients. • Liquid-core, ice-shelled condensate: Liquid particles undergo freezing, an ice shell forms and grows inwards towards the centre of each particle. • No condensate: All water is present in the vapour state.
If this is the case, we must still find typical conductivity coefficients in order to determine whether condensate should form.
All of these cases are discussed in cloud physics textbooks (e.g., Pruppacher and Klett, 2011;Khain and Pinsky, 2018); our aim here is to provide a translation of these cases into the non-equilibrium thermodynamics formalism.
We will first work through the simplest non-trivial case of liquid condensate, detail the microphysical processes, and pin down all resistivity coefficients within R, ready for usage in the semi-implicit semi-Lagrangian numerical scheme outlined in Section 2.3.
This section does not cover every thermodynamic process relevant to atmospheric physics, only the mass and energy changes associated with the simple idealised microphysical case of liquid droplet growth. To keep focus, some other microphysics cases modelled are provided in Appendix B. Note that some microphysics effects have not been included, such as solute and surface tension, precipitation, or droplet collisions. Nevertheless, the approach taken here should indicate how other processes such as these might be included in future work.

Liquid condensate
First then, let us build a typical case: a liquid particle undergoing condensation/evaporation in a parcel of well-mixed moist air. We make the assumption that the particle is thermodynamically homogeneous and spherically symmetric. Figure 1 illustrates this particle, along with the associated heat and mass fluxes. Assume, however-for any given fluid parcel-that instead of just one particle we have a number of identical liquid particles suspended in some well-mixed moist air. Let us define the changes in mass in such a way as to be consistent with Equation (2): where M is the total mass of the air parcel under consideration and  l is a specific mass-transfer rate to the liquid within the air parcel. The mass budgets may then be written in terms of diffusive fluxes acting within the thin boundary layer surrounding each liquid particle (layer I in Figure 1), that is, where r is the radius from the centre of a particle, N is the number density of cloud condensation nuclei (CCN) in the parcel, V is the total parcel volume, and J v is a diffusive F I G U R E 1 An idealised liquid particle suspended in a well-mixed, moist gaseous parcel. Layer I: evaporation/condensation layer (not to scale) mass flux of vapour. Essentially, the change in total mass of the liquid is given by the diffusive mass flux acting across the thin boundary layer surrounding the liquid particles. As this diffusive flux acts across the thin boundary layer surrounding the surface, it is scaled with the surface area. For conservation of mass, we must also have Note that we have still defined a mass budget for the ice production, even though there is none. Next, it is sensible to define the energy budgets for the various species. We consider the thin boundary layer surrounding the liquid particle (see Figure 1) and the thermodynamic gradients with associated thermodynamic fluxes acting across this thin layer. The energy budgets can then be defined in terms of mass and heat fluxes, as where  i is an energytransfer rate to species i, h i is a specific enthalpy, J v q ′ is a conductive heat flux, J v q is a heat flux, and we have taken advantage of the relation to re-express the internal energy change of liquid in terms of a heat flux instead of a conductive heat flux.
To be explicit, the internal energy change of vapour ( v in Equation (37)) is given by the enthalpy transfer associated with the mass flux of water vapour across the thin boundary layer surrounding the liquid particles (see Figure 1); the internal energy change of the liquid ( l in Equation (38)) is made up of a heat transfer and an enthalpy transfer associated with the mass flux of water vapour;  f in Equation (40) is zero, as there is no ice present. Note that  v +  l +  f does not add up to zero: the remainder is the energy transfer given to the fraction of dry air.
Our next goal is to link these mass and heat fluxes to what we commonly see in literature. Once this is done, we should be able to compare with the phenomenological equations, and relate the conductivity coefficients to known quantities.
First, we consider the mass flux. Diffusional mass flux (near a liquid particle) is usually given (Pruppacher and Klett, 2011) as where D v is the vapour diffusivity in air (assumed to be constant), v is the density of vapour, p v is the vapour partial pressure, R v is the gas constant of vapour, and T is a reference temperature, taken to be equal to the dry air temperature. Second, we look for a form of the heat flux. The conductive heat flux is given as where K a is the heat conductivity of air (also assumed to be constant). This leads to the heat flux Equations (42) and (46) are now of the right form to relate to the phenomenological equations.
To describe the mass and heat fluxes using phenomenological equations, we assume a quasi-steady state, implying non-divergent fluxes in the thin boundary layer (following de Groot and Mazur (1984)), and write where we have assumed T d = T v ; v and v q are constants, and L v vv , L v vq , L v qv , L v qq are some conductivity coefficients we must relate to known quantities. The fluxes are acting across the thin layer surrounding the liquid particle, hence the flux scales with the inverse surface area, that is, 1∕4 r 2 .
To pin down the conductivity coefficients, first neglect p d ∕ r to obtain the approximation and rewrite the mass and heat fluxes as Then it is possible to perform comparisons with Equations (42) and (46). This implies and hence where we achieve the expected symmetry: L v qv = L v vq . These coefficients may then be used in Equations (47) and (48). Now that we have pinned down the conductivity coefficients, we need to rewrite the mass budgets (Equations 34 and 36) and energy budgets (Equations 37, 39, and 40). Currently, the mass and energy budgets are in a continuous form, but to be in the form required for Equation (13) we need to re-express the forces in Equations (47) and (48) to be in terms of discrete differences.
To do this, we take Equations (47) and (48) and integrate up radially, between the particle surface r s and sufficiently far from the particle surface r ∞ . This results in the discontinuous expressions where Δ ∞s (X) indicates a difference between a thermodynamic property in the far-field and at the particle surface, that is, Δ ∞s (X) = X ∞ − X s . Calculation of r s is given in Appendix C.
Furthermore, to eventually be consistent with the forces in Equation (13), we must approximate terms involving g d . We again neglect gradients of p d to make the approximation where we have replaced T s with T l and T ∞ with T v . The approximation of Equation (60) allows re-expression of the thermodynamic force where we have assumed g v = g l at r = r s . Therefore we may rewrite the mass and heat fluxes, Equations (56) and (57), as These expressions are now consistent with Equation (13). Hence, the mass budgets (Equations 34 and 36) and energy budgets (Equations 37, 39, and 40) become where, to consider multiple particles in a moist atmospheric parcel, we have also made the substitution = M∕V. Alternatively, we may write Here,J l comprises mass and energy-transfer rates andL l the conjugate forces. The superscript "l" onJ l ,L l , and L l 3×3 indicates that we are dealing with the liquid condensate case. As introduced in Equation (13), the tilde has been used to indicate forces, transfer rates, and the conductivity matrix in terms of mass and energy. Eventually we will construct corresponding the quantities for mass and entropy (no tilde), as required for our numerical framework.
From Onsager (1931a) we expectL l 3×3 to be symmetric, yet in its current form this is not the case. Fortunately, however, it is possible to establish the expected symmetry via some basic transforms. As we will eventually have Δ vd (1∕T) = 0, the coefficients in the second column can take arbitrary values of our choosing. Considering this, L l 3×3 can indeed be made symmetric, that is, Equation (70) becomes where T is used to enforce thermodynamic equilibrium between T d and T v ; we must eventually let T → ∞, but we will defer taking this limit until after we have inverted L to obtain R. This infinite conductance effectively enforces an equilibrium between T d and T v . To eventually allow for a smooth transition to the equilibrium regime of R = 0, we must invertL l 3×3 . In doing so, and taking the limit as T → ∞, where, to help write this in a more compact form, we have made the substitutions This inversion thus provides us with However, at this pointR l 3×3 ,J l ,F l are subsets ofR l ,J,F, respectively. To findR l ,J,F, we consider the following which is simply a form of Equation (28)-now including the non-negativity constraints. We could then write the resistivity matrix as where the question marks indicate coefficients not provided by the above derivation. Next, we consider what the remaining unknowns of R l should be. Ideally, we want to ensure a smooth transition between this microphysical case and the other cases (derived in Appendix B). Suppose this microphysical case is currently active, then the non-negativity constraints acting on the ice ( q f ) must be switched on-by definition.
Hence, as long asR l remains positive semi-definite, we are free to define elements of the second and fifth rows and columns and not affect the solution. Therefore, we propose that we can write the resistivity matrix as This ensures a smooth transition betweenR l and the resistivity matrixR m derived later in Appendix B.2.
At this stage we are not quite in the final form required. To be consistent with Equation (23), we must make transforms toJ,F, andR l , to obtain J, F, and R l . To do this, we consider the relationship between transfer rates  i ,  i ,  i , given in Equation (15). Thus, we define a transformation matrix which allows us to relate J toJ by Next, by using the identities we find that the forces are related by Therefore, the substitution of Equations (79) and (82) into Equation (75) provides R l , where The idea-now that we have found R l -is that R l can be substituted for R in Equation (28) of our 2D model. This R is then used in the discretization and subsequent numerical solution. However, this is just one microphysical case; see Appendix B for the derivation of the other microphysics cases.

Transitioning between microphysical cases
Now that we have defined all of the different microphysical cases of our model, we want to consider what will happen in the transition between cases. In each parcel (at each model grid-point), we only have one microphysical case active at any one time (this is a simplification we are making at this stage in developing the approach, not an inherent limitation and not a property of the real world); when another becomes "active", the R matrix must be replaced appropriately. As such, we are interested in whether the transition between cases is continuous or not.
To aid in visualising how the transition between microphysical cases may take place, we refer to Figure 2 . It is immediately apparent that, for the cases within the shaded grey area, an identical resistivity matrix R m is used-transitioning between these is continuous or not.
It is straightforward to identify the character of the transitions (a), (b), and (c): (a) when freezing starts, the radius of the ice core r c is equal to r s , and G 2 = 0 (G 2 defined in Equation (B52)), so R is continuous across the transition; (b) similarly, when melting starts R is continuous; (c) when the liquid-shell starts to circulate, the resistance due to liquid resistivity drops to 0, hence there is a discontinuity in R; The temperature T f must remain at T 0 , so T s must drop to T 0 .
Transition (d) is more complicated. As the liquid core freezes, G 2 becomes large (G 2 → (r s − r CCN ) ∕r CCN , where r CCN is the radius of the condensation nuclei), but at the same time w 0 becomes small: To understand what is going on here, it is helpful to rewrite Equation (B72): F I G U R E 2 Possible pathways for transitioning between the microphysical cases. Arrows indicate how the microphysical cases are allowed to transition into one another.R where w 2 0 terms have been approximated to zero and G 1 substituted for −w 0 G 2 . As the liquid core freezes, G 2 becomes large, which works to enforce (considering row 1 and row 4 ) that is, there must be zero heat flux out of the liquid (cooling and freezing). This is as expected, as there must be zero heat flux out of liquid when the liquid disappears. It must also be the case (implied by Equation (86)) that the G 1 terms in row 2 and row 5 become negligible as r c → r CCN . Thus,  f and  f should vary smoothly as the liquid core disappears. Another perspective on this is that w s becomes close to 1 as r c → r CCN ; thus T s → T f , so the coupling to the gas phase changes smoothly as the liquid disappears. Therefore, the transition (d) is almost continuous.
In conclusion, transitions (c) and (d) are not continuous. However, we find that this does not cause any issues for the numerical solver. Now that we have defined a number of microphysical cases, we need a way of choosing which is currently active. We are assuming that all condensate particles within any parcel are identical, which implies that only one microphysical case can be active at any one time.
In the numerical scheme, for each time-step, the solver iterates towards a solution. At each iteration (and at each grid-point), it is permissible that the solver may switch to use a different microphysical case. We make use of the algorithm outlined in Algorithm 1. This provides a method to switch between the microphysical cases on the fly.

EXAMPLE RESULTS
In this section, some example results from our 2D model are presented. Although the test case is simple and somewhat artificial, it demonstrates that the proposed consistent non-equilibrium thermodynamic formulation can indeed be coupled to an atmospheric dynamical core based on typical numerical methods. Comparison with the equilibrium case (Bowen and Thuburn, 2022)  the sensitivity to "switching on" non-equilibrium effects is qualitatively as expected.
To compare with the results of Bowen and Thuburn (2022), we set up an idealised thermal in the same way. A 2D vertical slice is initialised with a total domain size of 20 km × 10 km, a grid resolution of 384 × 192, and a time-step of Δt = 10 s. An approximate atmospheric profile is prescribed: the initial pressure and densities are set up in such a way as to ensure hydrostatic balance; the boundary layer has a depth of 1,000 m, with a constant potential temperature of 300 K; above the inversion we have d ∕dz = 3.3 × 10 −3 K ⋅ m −1 ; the inversion has a change in of 0.2 K over a depth of 100 m with a cubic fit smoothly connecting the potential temperature profiles above and below. Moisture is then added. The boundary layer is adjusted to have a constant water mass fraction of 0.0125; above the inversion, a relative humidity of 0.7 is set. The initial vertical profiles for potential temperature and water mass fraction are shown in more detail in Bowen and Thuburn (2022).
Similarly to Bryan and Fritsch (2002), we add in a buoyancy perturbation centred about x = 10 km, z = 0 km, with no perturbation to pressure. This perturbation is identical to the one prescribed in Bowen and Thuburn (2022). The buoyancy perturbation is prescribed as follows: where is the acceleration due to gravity, and with x c = 10 km, x r = 2 km, z c = 0 km, z r = 0.5 km. From the same initial conditions, two model runs are performed for comparison: an equilibrium run and a non-equilibrium run. In the atmosphere, typical CCN concentration (N) ranges from approx 10 8 -10 9 m −3 (Twomey and Wojciechowski, 1969 ); we choose N to be at this lower value, to maximise the differences between the two simulations. Typical values for r CCN range from 10 −2 -10 μm (Khain and Pinsky, 2018 ); we choose r CCN = 0.2 μm . Note that a monotone limiter is used on the advected quantities (w, , q, X) to avoid spurious oscillations. We first follow the simulations to t = 1000 s, to show the differences in the equivalent potential temperature and vertical velocity w. The results of this are shown in Figure 3. Both regimes have a similar form at this time, allowing for reasonable comparison. Figure 3 a shows the equivalent potential temperature of the rising thermal for the equilibrium regime. At this time, the thermals in both regimes are still more buoyant than their surroundings and are rising vertically. The water vapour is changing phase to liquid as the thermal is rising-no ice has been produced. While the potential temperatures of both the equilibrium and non-equilibrium regimes are similar geometrically, we see the differences in Figure 3 b-the thermal of the equilibrium regime is more buoyant. We note that an absolute maximum difference of ≈ 7.38 K is found around the peaks of the thermals; however, much of this difference may be attributed to the small vertical offset between the thermals. Figure 3 c compares the vertical velocities zonally averaged across a 2D column-centred at x = 10 km with a width of approximately 2 km . Notably, there is already a robust difference between the equilibrium and non-equilibrium runs. As we expect, the resistance to producing liquid in the non-equilibrium regime means that the latent heat is released more slowly and hence the rising thermal is noticeably less buoyant, leading to a decreased vertical velocity.
Perhaps this robust difference is simply due to choosing N to be on the lower end of what is typical. If we increase N to 10 9 m −3 (therefore reducing thermodynamic resistance) and again look at the vertical velocity (Figure 4), we see that-while closer to the equilibrium regime-there is still a noticeable difference. The non-equilibrium regimes both have an attenuated vertical velocity when compared with the equilibrium regime.
As there is a lag in the condensation, this leads to a reduction in latent heat release, reduced buoyancy, and
reduced vertical velocity. This gives some indication (as in Pauluis and Held (2002a;2002b) that the irreversible phase changes involving water contribute to a decreased kinetic energy in convection. As well as N, there is also the sensitivity to changes in r CCN to consider. Simulations were run with various values of r CCN (within the typical range of 10 −2 -10 μm) and outputs compared (not shown); while some minor structural differences were found, on the whole the thermals were very similar.
To give another reassurance that the non-equilibrium processes are in effect, we plot the current condensate state (for both regimes) at t = 1400 s -demonstrated in Figure 5a,b. N is reverted to 10 8 m −3 . Both thermals are still rising and undergoing phase changes, this is long before either of the thermals reaches its neutral buoyancy. At this time, as the thermals rise they are also cooling-producing ice. In the equilibrium regime (Figure 5a), the phase changes equilibrate instantly, meaning that the current phase state is thermodynamically stable. Notice that the mixed-phase regions are small, only present in thin regions where the conditions are at the triple-point. In the non-equilibrium regime (Figure 5b), the phase changes are happening on a finite timescale. We find somewhat larger mixed-phase regions, with elongated structures reaching vertically through the cloud. This is due to the non-zero resistivities, allowing the mixed-phase states to exist for longer even though they are thermodynamically unstable. This existence of mixed-phase regions is not uncommon; in the right conditions, clouds are known to have portions of supercooled liquid water present (see Rauber and Tokay (1991) and references therein).
To verify further that the non-equilibrium model cloud is evolving in a non-equilibrium manner, we present Figure 5c. We consider the absolute value of the two thermodynamic forces Δ vl (g), Δ vf (g) across the condensate region, and take the natural logarithm of the maximum of these values. As expected, the entire cloud is out of thermodynamic equilibrium. On comparison with the vertical velocity (not shown), the regions of maximum

F I G U R E 6
Equilibrium regime, t = 2250 s. Contours of ice mass fractions, contour levels at 1 × 10 −3 . Note that, while liquid is allowed, none is present as condensate is colder than T 0 . disequilibrium are where there is maximum ascent/descent of the cloud and rapid phase changes.
Finally, we follow the development of the thermals to t = 2250 s; the equilibrium regime in Figure 6 and the non-equilibrium one in Figure 7. This is at a time when both thermals have reached (roughly) their altitude of neutral buoyancy. There are differences in structure of

F I G U R E 7
Non-equilibrium regime, t = 2250 s. Contours of ice mass fractions, contour levels at 1 × 10 −3 . Note that, while liquid is allowed, none is present as condensate is colder than T 0 . the cloud formed, but this is more likely due to turbulent effects.
To give an idea of how realistic our cloud composition is (for the non-equilibrium regime only), we look at the distribution of sizes of cloud condensate across the entire model domain. This is demonstrated in Figure 8. Interestingly, the peak radius is lower than we might F I G U R E 8 Distribution of cloud particle radii across the entire model domain, at t = 2250 s. expect for a typical cloud, and there is a sharp cut-off at around 20 μm. Given that our representation of condensate growth is highly idealised, a deviation from common observations is expected (e.g., see Curry and Webster (2008); Igel and van den Heever (2017)). While our distribution of r s is not wholly unrealistic, cloud particles do grow much larger than this. The difference here is that our model of droplet mass growth is entirely due to vapour diffusion. Droplet collision, coalescence, and condensate surface effects would need to be included to achieve a more realistic distribution. However, this extension is non-trivial, and beyond the scope of this article.

CONCLUSION AND DISCUSSION
First of all, let us summarise the goals we stated at the beginning of this article.
1. To extend the work of Bowen and Thuburn (2022) to the non-equilibrium case, specifically to outline a procedure to provide resistivity coefficients for use in the numerical solver. 2. To confirm that the resulting non-equilibrium thermodynamic formulation can be coupled to a typical numerical method used for atmospheric dynamics, and examine what effect including the non-equilibrium processes has on the model dynamics.
For the first of these, we have indeed managed to create a procedure to link some basic physical processes to the particular form of equations the model required. This was done by first examining a simple case (liquid particles in moist air suspension), then writing out mass and heat fluxes involved in the problem and relating these to the phenomenological equations-in doing so, it is possible to relate the resistivity coefficients to known quantities. Then the mass and heat fluxes can be linked to the mass and energy budgets and transformed into the final form required. A similar procedure was then carried out for the other microphysical cases. We found that, while the more complex microphysical cases do require more thought, various symmetries can be used to simplify the problem.
For the second of these aims, we have demonstrated that it is feasible to build a numerical model that couples the dynamical equations to a consistent framework for modelling the moist non-equilibrium thermodynamics. Despite the additional complexity of the thermodynamics, a semi-implicit time-stepping scheme with a quasi-Newton solver can be used and leads to a Helmholtz problem of the usual form for the pressure increments at every solver iteration.
The inclusion of these non-equilibrium processes affects the dynamics in the expected way. The finite rates of condensation and freezing lead to (amongst other things) a slower release of latent heat, which has a direct effect on the buoyancy and hence the dynamics. As the behaviour of the results is as expected, this gives some confirmation of the correctness of our implementation.
We should also comment on the numerics. Like Bowen and Thuburn (2022), a constant five iterations were used throughout. The extra complication of switching between microphysical cases did not appear to require any additional iterations. Note that the full five iterations are only necessary at locations where both liquid and ice are present simultaneously. To make the model cheaper to run, we could allow the numerical scheme to vary the number of iterations taken.

Next steps and final remarks
In demonstrating the feasibility of our approach, we made a number of simplifying assumptions in order to derive expressions for R. Nevertheless, our framework captures the non-equilibrium thermodynamics in a very general way; thus, we do not anticipate any serious barriers to extending the approach to include some of the more complicated microphysical processes needed in operational weather and climate models. One set of assumptions is that, within each air parcel or model grid box, all condensate particles are identical, and hence only one microphysical case is active at any one time. Provided we can take the thermodynamic forces F to be independent of drop size (or can model the dependence of F on drop size), then it appears possible to allow for a more realistic distribution of drop sizes and microphysical cases by averaging Equation (22) over the assumed or modelled drop-size distribution as well as over the different active microphysical cases. Due to the relatively simple dependence of R on the drop radius r s , this extension appears potentially tractable, though care will be needed in converting between L and R because of the equilibrium between T v and T d .
For very small droplets, diffusional growth is strongly modified by surface tension and solute effects (Pruppacher and Klett, 2011), which we have neglected here. Such effects are well understood, both from a cloud physics perspective and from a thermodynamics perspective (e.g., de Groot and Mazur (1984)), so it should be a straightforward exercise to derive a modified R that includes them.
Provided any colliding and merging drops have the same temperature, processes such as autoconversion and accretion affect only the drop-size distribution and so can be included without modifying the thermodynamics. However, if we wish to account for temperature differences between colliding drops of different sizes (or between liquid and ice particles), then a complete treatment would require the resulting entropy source to be included.
Another simplifying assumption made here is the neglect of precipitation. A common way to model the fallout of hydrometeors is to assume that they fall at their terminal velocity relative to the surrounding air (e.g., Lipps and Hemler (1982)). To accommodate such fallout in the present framework, Equations (7) and (8) would need to be reinterpreted for the hydrometeors as rates of change following the falling hydrometeors; care would then be needed in writing the total water Equation (4) and the total entropy Equation (10). Nevertheless, this extension appears tractable.
A further assumption made here is that the predicted quantities are distributed uniformly within each air parcel or model grid box. For some purposes it is useful to allow for subgrid variability within a grid box, particularly for the case in which only a fraction of the air in the grid box is saturated. For the case of the thermodynamic equilibrium regime, the relevant quantities can be computed by assuming a joint Gaussian subgrid distribution of total water and potential temperature (e.g., Sommeria and Deardorff (1977)). It would be interesting to explore whether a similar approach could be applied to the non-equilibrium regime within our proposed framework.
All of these possible extensions would build on ideas already available in the literature.
The growing impetus within the atmospheric modelling community to develop models based on a more consistent formulation of moist thermodynamics has motivated the present work, showing that an approach based on the internal energy potential can accommodate non-equilibrium as well as equilibrium thermodynamics. Further extensions of the approach along the lines outlined above could develop it from a research tool into a key component of operational weather and climate models.

ACKNOWLEDGEMENT
We express appreciation for the constructive reviews by Maarten Ambaum and an anonymous reviewer.

FUNDING INFORMATION
This work was funded by the Natural Environment Research Council through the GW4+ Doctoral Training Partnership (Grant NE/L002434/1) and as part of the ParaCon programme (Grants NE/N013123/1 and NE/T003863/1).

APPENDIX A. EXAMPLES OF INTERNAL ENERGY POTENTIALS
For reference purposes, here we list the internal energy potentials used. See Bowen and Thuburn (2022) for an explanation of their derivation.
An expression for the internal energy potential of dry air is given as For water vapour, the internal energy potential is given as (A2) For liquid, the internal energy potential is given as For ice, the internal energy potential is given as

APPENDIX B. MICROPHYSICS CASES
Here we provide the other microphysics cases implemented in our model, building on the simple case of liquid particle growth presented in Section 3.

B.1 Ice condensate
In this case, we picture a number of identical idealised ice particles suspended in a parcel of well-mixed moist air. The same assumption is imposed as before, that the ice particles are thermodynamically homogeneous and spherically symmetric- Figure B1 illustrates this. It is near-identical to Figure 1. albeit this time with an ice particle replacing the liquid particle. Our goal (as before) is to find a system of equations to effectively represent the thermodynamic processes, to match with Equation (13).
A procedure analogous to the one for the liquid condensate case could be carried out and would give us the set of equations required. However, we may take a shortcut by taking notice of the symmetries of the problem. In the liquid condensate case, all of the mass and heat fluxes take place within the thin layer surrounding the liquid particle. Here we expect something of a similar form: the mass and heat fluxes will be like before, but produce  f and  f , instead of  l and  l . The driving forces will also be modified: Δ vl (g∕T) is replaced with Δ vf (g∕T) , and Δ ld (1∕T) is replaced with Δ fd (1∕T) . With this in mind, it is possible to write an analogous form to Equation (74), F I G U R E B1 An idealised ice particle suspended in a well-mixed, moist gaseous parcel. Layer I: evaporation/condensation layer (not to scale).
which is simply a form of Equation (28)-again including the non-negativity constraints. Then, we propose that we can write the resistivity matrix as Similarly to the steps followed in Equations (76) and (77), we have chosen elements in the first and fourth rows and columns to ensure consistency with the other microphysics cases. Specifically, the motivation for the choices of these is to be consistent with the resistivity matrix defined for the ice-core, circulating liquid-shelled condensate case defined in Appendix B.2; this also implies consistency with the previously found resistivity matrixR l for the liquid condensate case (Section 3.1). In so doing, this ensures a smooth transition between these microphysical cases. If the ice condensate case is active, then the non-negativity constraints acting on the liquid ( q l ) must be switched on-by definition. Thus, the additional resistivities will not affect the solution. Now that we have determined theR f matrix, it is ready to undergo the final transformations to find R f and subsequently use it in the numerical solver.
B.2 Ice-core, circulating liquid-shelled condensate Suppose that an ice particle is melting. As it does so it will be enveloped by a liquid shell. In the simplest case, we picture a liquid shell where the fluid is circulating in such a way that any thermodynamic gradients are immediately broken down. The liquid shell must then be thermodynamically homogenous, and in equilibrium with the ice-core- Figure B2 illustrates this. The internal circulations usually take place when the liquid shell has sufficient depth relative to the total mixed particle radius. Rasmussen and Pruppacher (1982) state that circulation in the shell begins when the shell depth ≳10% of the total particle radius.
To build the equations we require, we can take advantage of the equilibrium condition between the ice core and liquid shell: where X represents a general thermodynamic property (g∕T or 1∕T). Alongside row 1 in Equation (B1), this implies that we must have both To find the remaining expressions we require, a similar procedure may be carried out on row 3 in Equation (B1). This results in the resistivity matrix F I G U R E B2 An idealised ice core encapsulated by a circulating liquid shell which is set withiñ where this time the non-negativity constraints are inactive. Once this is done,R m must undergo the final transforms to find R m . Note that at this moment we have three cases with identical resistivity matrices, This is a useful property, as it ensures that any direct transition between these microphysical cases is smooth.

B.3 Ice-core, liquid-shelled condensate
In this case, an ice particle begins to melt and is encapsulated by a liquid shell. If the liquid shell is not of sufficient depth, then there will be negligible internal circulations; as such, a temperature gradient within the liquid shell must develop. Due to the temperature gradient, this case is necessarily more complex than the previous one. We cannot build the resistivity matrix in the same way as before, and must complete a more complex derivation.
In this subsection, we will do the following.
1. Use Equation (B11) to build some similar expressions for total mass and total energy transfers to the mixed particle. Due to thermodynamic gradients in the liquid shell, the thermodynamic forces must be written in terms of properties at the ice-core liquid-shell interface, the liquid-shell surface, and the far-field. 2. Construct a temperature profile within the liquid shell to find an average T l and weighting coefficients. 3. Rewrite the thermodynamic forces to be in terms of average thermodynamic properties. 4. Find a general equation relating the heating profile within the liquid shell to thermodynamic fluxes. 5. Use boundary conditions to pin down the remaining equations, and constructR flJ =F .
6. Ensure the resulting resistivity matrix is symmetric.
As mentioned, we cannot simply use symmetries-similarly to previous cases-to construct R fl directly; however, we can use the previous case to help us build some relationships between the thermodynamic transfer rates and forces. From row 1 and row 2 in Equation (B11) we infer the analogous equation where this time the driving thermodynamic force is between the thermodynamic property g v ∕T in the far-field and at the particle surface. We must find some way to relate this surface property to the average within the shell for two reasons: the thermodynamic property at the surface will be different from the average shell property due to the thermodynamic gradient, and the driving force in Equation (23) only considers the average thermodynamic properties. Let T 0 be the temperature of the ice core and the freezing interface and T s be the temperature at the particle surface. Within the liquid shell, we expect the temperature profile to equilibrate quickly, hence we assume a steady state:r (r 2T r ) = 0, wherer = r∕r s is a non-dimensional radial coordinate, and whereT is a non-dimensional temperature. Solving Equation (B14) for the boundary conditions, we find the nondimensional temperature profile, wherer c is the non-dimensional radius of the ice core. From this temperature profile, it is then possible to work out the non-dimensional mean temperature in the liquid shell by calculating the volume average: where we have integrated radially using the boundary conditionsT = 0 atr =r c andT = 1 atr = 1. We then substitute into the temperature profile Equation (B16) to find an expression for the mean temperature T l , in terms of T s and T 0 , as T l = w 0 T 0 + w s T s .
Here, the weighting coefficients w 0 and w s are w 0 = 2r 2 c +r c 2 (r 2 c +r c + 1 ) and w s =r c + 2 2 (r 2 c +r c + 1 ) , (B19) where w 0 + w s = 1 as expected. Now that we have an expression relating the mean shell temperature to the surface and core temperatures, we can re-express the thermodynamic forces to avoid involving surface properties. To help with this, we may use Equation (B18) and write three useful relations: From these, letting X indicate thermodynamic properties 1∕T or g∕T, Similarly, where Δ s0 (X) indicates the difference of thermodynamic property X at the particle surface and the freezing interface, that is, Δ s0 (X) = X s − X 0 . As the extrapolations in Equations (B24) and (B28) are only taking place within the liquid shell, where the pressure is constant, then to a good approximation the variations in 1∕T and g∕T are both proportional to variations in T. Subsequently, thermodynamic forces involving a property at the surface can always be rewritten in terms of the average shell property. Using Equation (B26), we may rewrite Equation (B13) to avoid an expression in terms of surface and far-field properties: . (B31) Next, consider row 4 and row 5 in Equation (B10). We must have the analogous expression where again we have rewritten the thermodynamic force to avoid using surface and far-field properties. Alongside Equations (B31) and (B33), we also need an equation that captures the heating profile across the liquid shell; specifically, how the temperature difference Δ s0 (1∕T) drives a combination of heating the liquid shell and melting the ice core. The general idea here is to first build a general expression for the heat conduction through the liquid shell and then relate all of the quantities to the heat-flux profile.
For this we consider Figure B3. We may start with an expression for the heat conduction through the liquid shell: Although we assumed nondivergent heat fluxes earlier, here we must allow the liquid to warm. Hence, let us assume there is some continuous heating profile, where  0 (r) is some as yet unknown function representing the shape of the heating profile and A is an unknown constant. While it is possible instead to write this where A is absorbed into  0 (r), separating A out like this simplifies the following derivation. Note that at this point we are not specifying the shape of this heating profile. We integrate Equation (B35) and find r 2 J q ′ = A 1 (r) + B,