Caf\'e Latte: Spontaneous layer formation in laterally cooled double diffusive convection

In the preparation of Caf\'e Latte, spectacular layer formation can occur between the expresso shot in a glass of milk and the milk itself. Xue et al. (Nat. Commun., vol. 8, 2017, pp. 1-6) showed that the injection velocity of expresso determines the depth of coffee-milk mixture. After a while when a stable stratification forms in the mixture, the layering process can be modelled as a double diffusive convection system with a stably-stratified coffee-milk mixture cooled from the side. More specifically, we perform (two-dimensional) direct numerical simulations of laterally cooled double diffusive convection for a wide parameter range, where the convective flow is driven by a lateral temperature gradient while stabilized by a vertical concentration gradient. When the thermal driving force dominates over the stabilizing force, the flow behaves like vertical convection in which a large-scale circulation develops. However, with increasing strength of the stabilizing force, a meta-stable layered regime emerges. Initially, several vertically-stacked convection rolls develop, and these well-mixed layers are separated by sharp interfaces with large concentration gradients. The initial thickness of these emerging layers can be estimated by balancing the work exerted by thermal driving and the required potential energy to bring fluid out of its equilibrium position in the stably stratified fluid. In the layered regime, we further observe successive layer merging, and eventually only a single convection roll remains. We elucidate the following merging mechanism: As weakened circulation leads to accumulation of hot fluid adjacent to the hot sidewall, larger buoyancy forces associated with hotter fluid eventually break the layer interface. Then two layers merge into a larger layer, and circulation establishes again within the merged structure.


Introduction
Layered patterns are striking features in double diffusive convection (DDC), where the fluid density depends on two scalars with different diffusivities (Turner 1974;Huppert & Turner 1981;Schmitt 1994;Radko 2013;Garaud 2018). A typical example of layer formation is found in the Ocean, where seawater density is affected by temperature and salinity. As a result of double diffusion, thermohaline staircases are found in different regions of the Ocean, such as a salt-finger regime in (sub-)tropic regions (Schmitt 2005;Johnson & Kearney 2009;Yang et al. 2019) and a diffusive regime in high-latitude regions (Kelley et al. 2003;Timmermans et al. 2008;Sommer et al. 2013).
An intriguing daily example of layered pattern can be found in Café Latte, where the corresponding laboratory experiments have recently been conducted by Xue et al. (2017). When a shot of espresso (lower-density) is poured into a glass of milk (higher-density), the system cools down from the side since it loses heat to the ambient through the sidewall, pronounced layers form in the mixture, rather than a mixed-up solution as one may expect. Xue et al. (2017) showed that the injection velocity determines the depth of milk being mixed with expresso. After a while when a stably-stratified zone forms in the mixture, the layering process is governed by double diffusion: The temperature difference between the hotter bulk and the colder sidewall fluid layer implies a horizontal thermal driving, whereas a stabilizing vertical concentration gradient exists in the coffee-milk mixture.
Examples in other physical systems also illustrate the importance of horizontal thermal driving to layer formation in a stably-stratified fluid. For instance, when sedimenting suspensions of colloidal particles are subjected to a horizontal temperature gradient, the initially uniform suspension will also develop multiple layers (Mendenhall & Mason 1923). Also, one can observe layering when ice blocks melt in salty liquid, building up a salinity gradient (Huppert & Turner 1980).
To numerically study layer formation in DDC systems with both vertical and lateral gradients, here we pick, inspired by Café Latte, laterally cooled double diffusive convection with a concentration gradient in vertical direction. In this set-up, the temperature gradient is imposed horizontally, whereas the vertical concentration gradient is stabilizing. In pioneering experimental and theoretical work of laterally cooled DDC, Thorpe et al. (1969) showed the successive growth of layers in a stratified brine solution heated from one side. They further conducted linear stability analysis to find the onset criteria of layers. Their pioneering paper motivated further experimental and numerical work focusing on how the layers form (Chen et al. 1971;Wirtz et al. 1972;Lee & Hyun 1991). Also, salinity and heat fluxes were studied extensively in laterally cooled DDC, because it is relevant to the high-latitude Ocean being affected by melting icebergs (Jacobs et al. 1981;Huppert & Turner 1980;Gayen et al. 2016). Moreover, layer merging in DDC is also an important issue because it influences the fluxes across the layer interface (Tanny & Tsinober 1988;Chen & Chen 1997).
Previous simulations (mostly in 20th century) on DDC had severe CPU-time limitation on the parameter range and on collecting long enough time series of layer evolution. Thanks to the full temperature and velocity information obtained from present numerical simulations, the extension of the parameter space and the possibility to run very long simulations, the layer formation and properties can now be understood in much more detail.
In this work, we study laterally cooled DDC over a wide range in parameter space, namely three decades of temperature Rayleigh number Ra T and four decades of density ratio Λ. We begin with the description of the governing equations and the setup in Section 2. Then we examine the flow morphologies and show the layer formation in Section 3. In Section 4, we can estimate the thickness of the initially formed layers from an energy balance. We further elucidate the mechanism of layer merging in Section 5. Finally, conclusions are given in Section 6.

Numerical method and set-ups
We consider a two-dimensional rectangular box of width W and height H. The left/right wall has high/low temperature, and there is no salinity flux through the lateral boundaries. The top/bottom wall has low/high salinity and is adiabatic to temperature. No-slip velocity boundary conditions are used on all the walls. We apply the Oberbeck-Boussinesq (OB) approximation, such that the fluid density depends linearly on tem-peratureT and a scalarS:ρ(T ,S) =ρ 0 1 − β T T −T 0 + β S S −S 0 . Here,ρ 0 , T 0 ,S 0 represent the reference density, temperature and concentration, respectively. β T and β S are the thermal and solutal expansion coefficients. The governing equations are nondimensionalized by normalizing lengths by H, velocities by the free-fall velocity U = gβ T |∆ T |H, temperatures by ∆ T (the temperature difference between the sidewalls) and concentrations by ∆ S (the concentration difference between the top and bottom plates): Here, u i are the velocity components, p the kinematic pressure, T the temperature and S the concentration, all now non-dimensional. δ iz denotes the Kronecker delta and g the gravitational acceleration. The five dimensionless control parameters are the aspect ratio Γ , the thermal Rayleigh Ra T and the Prandtl P r T number for the temperature, the Lewis number Le, and the density ratio Λ, defined as: is the concentration Rayleigh number and P r S = ν/κ S the concentration Prandtl number. ν, κ T and κ S are the kinematic viscosity, the thermal diffusivity, and the solutal diffusivity, respectively. Λ measures the relative strength of the buoyancy force induced by the stabilizing concentration difference to that induced by the destabilizing temperature difference. The three key response parameters of the system are the two scalar fluxes and the flow velocity, which are measured by the two Nusselt numbers and the Reynolds number: Here . x,t / . z,t represents the average over time and the horizontal/vertical plane. In this work, we calculate N u T by temperature gradients at the two sidewalls and N u S by concentration gradients at the top and bottom plates. u 2 is the root-mean-square value of the velocity magnitude, calculated over the entire domain.
Equations (2.1)-(2.4) are solved by a second-order finite difference scheme using a fractional-step procedure and advanced in time by a low-storage third-order Runge-Kutta scheme (Verzicco & Orlandi 1996;van der Poel et al. 2015). We use fixed aspect ratio Γ = 0.5. Our simulations cover the range 10 6 Ra T 10 9 , 10 −2 Λ 10 2 with P r T = 1 and P r S = 100 (corresponding to a Lewis number Le=100 which is large enough to demonstrate the layer formation). The large Le implies that the resolution for the concentration is more demanding than that for the temperature, and thus a multiple resolutions strategy is employed , and such strategy had been already used for DDC simulations . Specifically, for the case of Ra T = 10 9 , we use 432 2 for the base mesh and 1296 2 for the refined mesh.
We are aware of the limitations of 2D DDC simulations as compared to 3D ones. However, at least for large P r 1, qualitatively the results for 2D & 3D are very similar (van der Poel et al. 2013;Chong et al. 2020) and we aim more at elucidating the physical processes originating the layering rather than detailing a specific case. Only by restricting us to 2D, we can explore a large region of the parameter space.

Flow structures at various density ratios
We begin with the qualitative description on how the flow morphology changes with increasing density ratio Λ, where Λ measures the strength of the thermal buoyancy compared to stabilization due to the stable stratification. For this study we fix Ra T = 10 9 and Le = 100. In figure 1 we show the concentration and temperature fields for both regimes which emerge: At Λ = 1, when the thermal buoyancy is dominant, there is a single large-scale circulation. From the temperature field, it can be seen that the detached hot (cold) plumes travel upwards (downwards). For the chosen aspect ratio Γ = 1/2 they travel over the distance of the entire cell height. At the same time, the concentration field is advected by the thermally-driven circulation. As a result, there is a well-mixed region formed in the bulk with nearly uniform concentration, whereas the concentration only changes sharply near the top and bottom boundary layers. This flow structure is similar to that in vertical convection (Ng et al. 2015;Shishkina & Horn 2016;Wang et al. 2019). We thus classify this and corresponding cases into the so-called quasi-VC regime.
Strikingly, different flow structures are obtained for Λ larger than a threshold value, which will be calculated in Section 4. For example, for Λ = 4 beyond this threshold, we identify the formation of a layered structure from the concentration field. The physical process of layer formation is as follows: Initially, the detached hot plumes travel upwards. Due to the restoring force caused by stably-stratified concentration field, thermal buoyancy is not strong enough to maintain the upward-moving plumes throughout the entire domain height. Therefore, thermal plumes travel horizontally towards the middle of the cell, causing a sequence of thermal streaks as seen from the temperature field. In this case, the thermal driving leads to the vertically-stacked convection rolls. Because the concentration diffuses much slower than the heat (Le = κ T /κ S = 100), a wellrecognizable layered concentration field is resulted. Within each roll, the concentration is nearly uniform due to the convective mixing. At the interface between two adjacent rolls, the concentration changes sharply. For an even larger stabilization (Λ = 7), even more layers initially form as compared to the case of Λ = 4, in accordance with our physical explanation of the formation process.

Initial layer thickness and phase space
As shown above, a series of layers will form initially in the layered regime, and the initial layer thickness decreases with increasing strength of stabilization. What sets the initial layer thickness, or equivalently the size of the localized circulation? We will derive this initial layer thickness from an energy balance. A similar energy argument was adopted in stratified Taylor-Couette (TC) flow, which is TC flow subjected to vertical linear stratification (Boubnov et al. 1995). In stratified TC, spontaneous layer formation can be observed in both the low-Re (Boubnov et al. 1995) and the high-Re (Oglethorpe et al. 2013) regimes. In order to estimate the layer thickness in stratified TC, Boubnov et al. (1995) successfully employed the balance between the work exerted by the centripetal force and the potential energy for moving the fluid parcel in the stable stratification.
Likewise, in laterally cooled DDC, the work for raising the fluid parcel in the stable stratification is done by the thermal buoyancy. As the fluid movement is driven by the horizontal temperature difference ∆ T , the work done by thermal buoyancy to raise the fluid parcel over distance h is (β T g∆ T )h. In stable linear stratification, the potential The boundary between the quasi-VC and the layered regime is given by Λ = 2 as derived by the energy balance in §4. The boundary between the layered regime and the regime without convection is Λc = 0.6Ra 1/5 T as obtained already by Thorpe et al. (1969). The colors of the points denote the number of layers observed in the early stage of the layer formation; later these layers partly merge. energy to bring a fluid parcel out of its equilibrium position with vertical displacement h is N 2 0 h 2 , where N 0 = gβ s ∆ s /H is the buoyancy frequency. Assume that all work is converted to potential energy. This balance gives We emphasize that this relationship is only valid for estimating the initial thickness because it assumes a linear stratification which is only the case during the initial stage. We now check whether equation (4.1) is a good approximation to the initial layer thickness. We first manually count the number of layers formed in the very initial stage after layer development, which is well recognizable from the snapshots. Then the average layer thickness h can be estimated by dividing the cell height H over the counted number. Figure 2(a) shows the evaluated layer thickness h /H versus 1/Λ for various Ra T . It can be seen that the data points generally follow the trend of h/H = 1/Λ. Yet, close inspection suggests that all data points are actually below the estimated line (black dashed), consistent with previous experiments which also found that the measured initial thickness is in general less than 1/Λ (Chen et al. 1971). There are two reasons why h/H is slightly smaller than 1/Λ: (i) Part of the work is not converted but dissipated which is neglected in obtaining equation (4.1). (ii) The buoyancy in general is smaller than β T g∆ T . Thus, equation (4.1) can only be seen as upper limit for the initial thickness. In addition, the layer thickness is obviously limited by the system height H, and indeed in figure 2a it can been seen that h levels off at h = H for small Λ 1 or 1/Λ 1.
We now explore the full parameter space (Λ, Ra T ) to identify when a single large-scale circulation forms, and when there are layers consisting of stacking localized circulations. In figure 2(b), the transition boundary between the quasi-VC regime and the layered regime corresponds to the case with only two initially-formed layers, i.e., Λ = 2 from equation (4.1). Upon increasing Λ, the number of layers increases but eventually the system reaches the motionless state when stabilization becomes dominant. This transition density ratio Λ c to the no convection regime was deduced previously from linear stability analysis (Thorpe et al. 1969).

Layer merging and its mechanism
We next address the merging of the layers, which sucessively occurs as time proceeds. It obviously coincides with the number of layers decreasing monotonically with time. Previous experiments had observed an eventual single roll state after successive layer merging (Kamakura & Ozoe 1993) for certain parameters. To study the merging process in detail, we simulated the case of Ra T = 10 9 and Λ = 7 for 6 × 10 4 free fall time units. Figure 3(a) shows the time evolution of the lateral heat flux (N u T,right and N u T,lef t ) and the vertical solutal flux (N u S,top and N u S,bot ) for this prolonged run; the instantaneous concentration fields at different time instants are also shown in figure 3(b). Between the consecutive merging events, the system reaches a meta-stable state with fluxes fluctuating around an average value. However, just before the transition to another meta-stable state (i.e., layer merging), we observe spikes in the heat flux time series which are characteristic of layer merging. After the neighbouring layers have merged, the heat fluxes then reach another average value. In contrast, the solutal fluxes are less sensitive to the merging event as compared to the heat fluxes. The reason is that the solute diffuses a hundred times slower than the temperature, and thus it takes a longer time for the local concentration change to affect the top and bottom solutal fluxes.
The N u behaviour leads us to ask (i) Why are there spikes just before the merging of layers? (ii) What is the merging mechanism in the laterally cooled DDC? To answer these questions, we examine a particular merging event. Figure 4(a) shows the vertical and horizontal Reynolds number, Re z and Re x , computed by the globally-averaged horizontal and vertical velocities. Between t 1 and t 2 , Re z is almost unchanged, whereas Re x progressively decreases with time during this period. The weakened horizontal velocity first explains why there is a gradual decrease in lateral heat fluxes shown in figure 4(b).
To study the merging process in even more detailed, we visualize the temperature snapshots at the corresponding time instants t 1 and t 2 in figure 4(d). when the layers are just about to merge (t = t 2 ), we observe that the thermal streak becomes shorter than before (t = t 1 ), reflecting the weakened local circulation within a layer. Eventually, the weakened circulation leads to an accumulation of hot fluid adjacent to the heated walls. Indeed, in figure 4(c), we see that the averaged temperature over the area near the heated wall (denoted by zone A 1 ) increases gradually before t 2 . Finally, when the hot fluid has large enough potential energy to overcome the stable stratification, two nearby layers merge into a larger one at t 3 . A new circulation establishes after layer merging, it further carries a chunk of accumulated hot fluid near the heated wall to the opposite cold wall, and subsequently leads to the sharp increase in the heat fluxes. Our results demonstrate the complete process of layer merging in laterally cooled DDC, and further explain why there are spikes in the heat flux time series.

Concluding remarks
In summary, inspired by the layer formation in Café Latte, we have numerically studied laterally cooled and stably density stratified double diffusive convection which is regarded as a simplified version of Café Latte. We have clearly demonstrated the layer formation, which occurs after the initial expresso injection process. We numerically explored a large range parameters, namely for the temperature Rayleigh number 10 6 Ra T 10 9 and for the density ratio 10 −2 Λ 10 2 , with Le = 100.
Upon increasing strength of the stabilizing concentration gradient at a fixed lateral thermal driving, the study reveals the flow transition from the quasi-VC regime to the layered regime. In the quasi-VC regime, the flow structure is a large-scale circulation. However, in the layered regime, multiple localized circulations initially form which therefore leads to the layered structure in the concentration field. Based on the energy balance between the work done by thermal buoyancy and the potential energy to bring a fluid parcel out of its equilibrium position in stratification, we can estimate the initial layer thickness, obtaining that it roughly follows h/H = 1/Λ. Such a relationship also allows us to find the boundary between quasi-VC and layered regimes, which is Λ = 2.
We finally focused on merging events by running a specific case with a long time series. With sufficiently long simulations, we showed that the layered structures eventually merge into a single large-scale circulation. The merging mechanism is that the weakened circulation within a layer leads to the accumulation of hot fluid over the hot sidewall. The hot fluid parcel at some point obtained enough buoyancy to overcome the energy barrier set by the stable stratification, and it forms a new circulation of larger size. The formation of the new circulation leads to the spikes in the heat flux time series which is characteristic for layer merging.
Until now, we have only considered the cases with fixed temperature at the heated and cooled walls. However, in many circumstances, the double diffusive convection may also be subjected to time-dependent boundary conditions, for example, abrupt temperature change caused by falling icebergs, or seasonal temperature variation. Those time-dependent forcing may have a pronounced effect on the layer formation and merging, which is worthy to be studied in the future.