Flow-Area Relations in Immiscible Two-Phase Flow in Porous Media

We present a theoretical framework for immiscible incompressible two-phase flow in homogeneous porous media that connects the distribution of local fluid velocities to the average seepage velocities. By dividing the pore area along a cross-section transversal to the average flow direction up into differential areas associated with the local flow velocities, we construct a distribution function that allows us not only to re-establish existing relationships between the seepage velocities of the immiscible fluids, but also to find new relations between their higher moments. We support and demonstrate the formalism through numerical simulations using a dynamic pore-network model for immiscible two-phase flow with two- and three-dimensional pore networks. Our numerical results are in agreement with the theoretical considerations.


INTRODUCTION
When two immiscible fluids compete for the same pore space, we are dealing with immiscible two-phase flow in porous media [1]. A holy grail in porous media research is to find a proper description of immiscible two-phase flow at the continuum level, i.e., at scales where the porous medium may be treated as a continuum. Our understanding of immiscible two-phase flow at the pore level is increasing at a very high rate due to advances in experimental techniques combined with an explosive growth in computer power [2]. Still, the gap in scales between the physics at the pore level and a continuum description remains huge and the bridges that have been built so far across this gap are either complex to cross or rather rickety. To the latter class, we find the still dominating theory, first proposed by Wyckoff and Botset [3] and with an essential amendment by Leverett [4], namely relative permeability theory. The basic idea behind this theory is the following: Put yourself in the place of one of the two immiscible fluids. What does this fluid see? It sees a space in which it can flow limited by the solid matrix of the porous medium, but also by the other fluid. This reduces its mobility in the porous medium by a factor known as the relative permeability, which is a function of the how much space there is left for it. And here is the rickety part: this reduction of available space-expressed through the saturation-is the only parameter affecting the reduction factor or relative permeability. This is a very strong statement and clearly does not take into account that the distribution of immiscible fluid clusters will depend on how fast the fluids are flowing. Still, in the range of flow rates relevant for many industrial applications, this assumption works pretty well. It therefore, remains the essential work horse for practical applications.
Thermodynamically Constrained Averaging Theory (TCAT) [5][6][7][8][9] is built on the framework of relative permeability. However, it is based on a full analysis based on mechanical conservation laws, constitutive laws, e.g., for the motion of interfaces and contact lines, and on thermodynamics at the pore level. These are then scaled up using averaging theorems, which, loosely explained, consist of replacing derivatives of averages by averages of derivatives. In principle, this approach solves the up-scaling problem. However, as Gray and Miller point out in their book [9], each component of TCAT involves significant mathematical manipulations. The internal energy has contributions from the bulk liquids, the fluidfluid and fluid-matrix interfaces at the pore level. The averaging process redefines the variables describing these contributions, but does not reduce their number. This accounts for a high level of complexity.
A further development somewhat along the same lines, based on non-equilibrium thermodynamics uses Euler homogeneity, more about this later, to define the up-scaled pressure. From this, Kjelstrup et al. derive constitutive equations for the flow [10,11].
Another class of theories is based on detailed and specific assumptions concerning the physics involved. An example is Local Porosity Theory [12][13][14][15][16][17]. Another is DeProf theory which is a mechanical model combined with non-equilibrium statistical mechanics based on a classification scheme of fluid configurations at the pore level [18][19][20].
A recent work [21] has explored a new approach to immiscible two-phase flow in porous media based on elements borrowed from thermodynamics. That is, it is using the framework of thermodynamics, but without connecting it to processes involving heat. The spirit behind this approach is like that taken in Edwards and Oakshot's pseudo-thermodynamic theory of powders [22]. The approach consists in looking for general relations that transcends details of the physical processes involved. An example of such an approach in the field if immiscible two-phase flow in porous media is found in the Buckley-Leverett theory of invasion fronts [23]. The Buckley-Leverett equation is based solely on the principle of mass conservation and on the fractional flow rate being a function of the saturation. In the approach of Hansen et al. [21], equations are derived that originate from Euler homogeneity as in ordinary thermodynamics. These equations transcend the details of the physics involved in the same way that the equations of thermodynamics are universally applicable if a set of simple underlying conditions are met.
Thermodynamics is a theory that is valid on scales large enough so that the system it refers to may be regarded as a continuum. Statistical mechanics is then the theory that makes the connection between thermodynamics and the underlying atomistic picture.
It is the aim of this paper to formulate a description of immiscible two-phase flow in porous media that may form a link between the continuum-level approach of Hansen et al. [21] and the pore-level description of the problem-a sort of "statistical mechanics" from which the pseudo-thermodynamics may be derived, but which also describes the flow problem at the pore level. After defining the system and the variables involved in section 2, we will in section 3 review the pseudo-thermodynamic approach [21]. The next section 4 we introduce the central object in the paper, the differential transversal area distribution which corresponds to the Boltzmann distribution in ordinary statistical mechanics, and relate it to the pseudo-thermodynamics relations. Then follows section 5 which then moves beyond the results of the pseudo-thermodynamics by focusing on fluctuations. In section 6 we use the dynamic network simulator [24] first introduced by Aker et al. [25] and then later refined [26][27][28] to verify the relations derived in the earlier sections. There is also a second goal behind this numerical work: the dynamic network model is a model at pore level and by its use, we show how the formalism developed here connect to the flow patterns at the pore level. Finally, we draw our conclusions in section 7.

SYSTEM DEFINITION
In two-phase flow, the steady state [29][30][31] is characterized by potentially strong fluctuations at the pore scale, but steady averages at the REV (Representative Elementary Volume) scale. As such they differ fundamentally from stationary states that are static at the pore scale as well. Steady states have much in common with ensembles in equilibrium statistical mechanics. They are also by implication assumed in the conventional descriptions of porous media flows that take the existence of an REV for granted.
Our REV is a block of homogeneous porous material of length L and area A. We prevent flow through the surfaces that are parallel to the L-direction which is the flow direction. The two remaining surfaces, each having an area A, act as inlet and outlet for the incompressible fluids that are injected and extracted from the REV. The porosity of the material is defined as where V p is the pore volume and V = AL is the volume of the REV. Due to the homogeneity of the porous medium, any cross section orthogonal to the axis along the L-direction will have a pore area that fluctuates around the value There is also a solid matrix area fluctuating around The homogeneity assumption consists in the fluctuations being so small that they can be ignored. There is a time averaged volumetric flow rate Q through the REV. The volumetric flow rate consists of two components, Q w and Q n , which are the volumetric flow rates of the more wetting (w for "wetting") and the less wetting (n for "non-wetting") fluids with respect to the porous medium. They are related through In the porous medium, there is a volume V w of the wetting fluid and a volume V n of the non-wetting fluid so that V p = V w + V n . We define the wetting and non-wetting saturations S w = V w /V p and S n = V n /V p , so that S w + S n = 1. We define the wetting and non-wetting transversal pore areas A w and A n as the parts of the transversal pore area A p which occupied by the wetting or the non-wetting fluids, respectively. We have that As the porous medium is homogeneous, we will find the same averages A w and A n in any cross section through the porous medium orthogonal to the flow direction. We have therefore Likewise, We define the seepage velocities, i.e., the average flow velocities in the pores, for the two immiscible fluids, v w and v n as and v n = Q n A n .
The seepage velocity associated with the total flow rate Q is defined as We may express Equation (4) in terms of the seepage velocities, v p = S w v w + S n v n .

PSEUDO-THERMODYNAMIC RELATIONS
Hansen et al. [21] derived a number of relations between the seepage velocities defined in (8)-(10) based on the volumetric flow rate being an Euler homogeneous function of order one with respect to the wetting and non-wetting transversal pore areas A w and A n . We present here a short review of the main results in that paper for completeness. The meaning of the statement that the volumetric flow rate is an Euler homogeneous function of order one is that it obeys the scaling relation where λ is a scale factor. By taking the derivative of this equation with respect to λ and then setting λ = 1, we find By dividing this expression by the transversal pore area A p and using Equations (5)- (7), we may write this equation as The two partial derivatives have the units of velocity, and Hansen et al. [21] name these velocity functions the thermodynamic velocities, We use Equations (6) and (7) and the chain rule to derive Likewise, we find We now combine these two equations with the definitions (15) and (16), and use that Q = A p v p , i.e., Equation (10), to find andv Combining the definitions (15) and (16) with Equation (14) gives which should be compared to Equation (11). We see that The seepage and thermodynamic velocities are related through andv We now calculate Using Equations (6) and (7) together with Equations (19) and (20), we transform this equation into where we have used that v p = Q/A p , i.e., Equation (10). We now use Equation (21) to calculate Compare this equation to Equation (26) and we get an analog to the Gibbs-Duhem equation, Using Equations (23) and (24), we find that the seepage velocities obey and where we have combined Equations (23) and (24) with Equation (28). By combining Equations (15), (16), (23), and (24), one finds and These two equations, (31) and (32), may be seen as a The inverse of this transformation, i.e., (v w , v n ) → (v p , v m ) are given by Equations (11) and (29), i.e., where v ′ w = dv w /dS w and v ′ n = dv n /dS w . But, what is the co-moving velocity v m physically? We first need to understand the thermodynamic velocitiesv w andv n . These are the velocities the two fluids would have had if they were miscible. Equation (26) then tells us that a change in the saturation S w leads to a change in the average seepage velocity v p which is the difference in seepage velocities of the two fluids. However, the two fluids are not miscible and they do get in each other's way. How much is dictated by the co-moving velocity through Equation (29).
From Equation (26) onwards to the end of this sections, none of the equations contain the size of the REV. If we now imagine a REV associated with each point in the porous medium, we have a continuum description. We may then add equations that transport the fluids between these points. Assuming that the fluids are incompressible, these equations are [1] where t is the time coordinate and x is the spatial coordinate, and We add the two equations and get The generalization to three dimensions is straight forward.
In order to connect the equations that now have been derived to a given porous medium, constitutive equations for v p and v m need to be supplied, linking the flow to the driving forces. These may in the simplest case be pressure gradient and saturation gradient.

DIFFERENTIAL TRANSVERSAL AREA DISTRIBUTIONS
In this section, we connect the pseudo-thermodynamic results of section 3 to the properties of an underlying ensemble distribution. This concept in the context of immiscible two-phase flow was first considered by Savani et al. [32]. Here we generalize this concept. In some sense, we introduce here a statistical mechanics from which the pseudo-thermodynamics ensue.
We define a differential transversal pore area a p = a p (S w , v) where v is a velocity such that a p dv is the pore area covered by fluid, wetting or non-wetting, that has a velocity in the range [v, v + dv]. Hence, a p -and the other differential transversal pore areas that we will proceed to construct-are statistical distributions of the pore level velocities. The new idea we are introducing is that the velocity distribution is measured in terms of transversal pore areas. This makes it possible to make the connection between the flow at the pore level and the pseudothermodynamic theory reviewed in the previous section.
We must have that where the integral runs over the entire range of negative and positive velocities since there may be local areas where the flow direction is opposite to the global flow. The total flow rate Q is given by and the see page velocity defined in Equation (10) is then given by Likewise, we define a wetting differential pore area a w and a nonwetting differential pore area a n . They have the same properties except that they are restricted to the wetting or the non-wetting fluids only. That is, we have and They relate to the wetting and non-wetting seepage velocities defined in Equations (8) and (9) as and We have that a p = a w + a n .
We now combine this equation with Equation (39) to find which is Equation (11). We have here used Equations (6) and (7). We may associate a differential area a m to the co-moving velocity v m defined in Equation (29). By using Equations (39), (42), and (43) in combination with Equation (29), we find so that  (29) and (30) to the differential transversal areas. It follows that where A m is the pore area associated with co-moving velocity v m . This is to expected as the areas A w , A n , A p and A m are ways to partition the transversal pore area A p ; and we have that A w + A n = A p + 0. This implies that there is no volumetric flow rate associated with the co-moving velocity since Lastly, we may associate differential transversal areas to the thermodynamic velocities defined in Equations (19) and (20). We use Equations (23) and (24) to find andâ n = a n − S w S n a m , where a m is given in Equation (48). The thermodynamic velocities are then given bŷ We find as expected that Summing the two differential transversal areas for the thermodynamic areas giveŝ a w +â n = a w + a n = a p .
This leads us to an important remark. The differential transversal areas are statistical velocity distributions at the pore level. We see that the differential transversal areas that are associated with the thermodynamic velocities are different from those associated with the seepage velocities. However, Equation (57) shows that the combined differential transversal area based upon the thermodynamic velocity distributions is the same as that based upon the distributions giving the seepage velocities. Hence, the two types of differential transversal areas represent a redistribution of the pore level velocities, but in such a way that A w and A n are preserved. We see the same from Equation (49) showing that A m is zero and combining this Equations (51) and (52). We see from Equation (47) that a m is only zero if a w and a n are linear in S w and S n , respectively, i.e., a w = S w b w where b w is independent of S w and a n = S n b n where b n is independent of S n . Hence, this is the condition for the thermodynamic velocities to be equal to the seepage velocities.

MOMENTS AND FLUCTUATIONS
We define the qth moment of the seepage velocity distribution as By using Equation (44) where we have defined and We may work out the moments of the co-moving velocity are given by where we have used (47). The thermodynamic velocity moments may be defined as in a similar manner as the moments of the seepage velocities, (60) and (61),v and we findv q p =v q w S w +v q n S n , where we have used Equations (52) and (55). We may Fourier transform a p , a w , and a n , and 2πã n (ω) = A n e ivω n = ∞ −∞ dv e ivω a n .
From Equation (44) we find and e ivω p = S w e ivω w + S n e ivω n .
We write exp(ivω) p as a cumulant expansion, where C k p is the kth cumulant. We define the wetting and nonwetting velocity cumulants C k w and C k n in the same way. We also write exp(ivω) p as a moment expansion By expanding the cumulant expression in Equation (71) and equating each power in iω with the corresponding one in Equation (72), then repeating this for the wetting and nonwetting cumulants, and lastly combining them through Equation (70), we find for the term proportional to (iω) 2 , Noting that C 1 p = v p , C 1 w = v w , and C 1 n = v n and using that v 2 p = C 2 p , v 2 w = C 2 w , and v 2 n = C 2 n , we find from this equation We may follow this procedure for any of the cumulants. The corresponding equation between the second cumulants of the thermodynamic velocities is

NUMERICAL OBSERVATIONS
The relations presented in sections 4, 5 provide the bridge between the velocity distributions at the pore level and the pseudo-thermodynamic theory outlined in section 3. In order to test these relations, and to show how they may be used, we use a dynamic pore network simulator [24]. In pore network modeling, the porous medium is represented by a network of pores which transport two immiscible fluids. The pore-network model we consider here can be applied to regular networks such as a regular lattice with an artificial disorder as well as to irregular networks such as a reconstructed network from real samples. The flow of the two immiscible fluids is described in this model by keeping the track of all interface positions with time. This approach of pore network modeling was first introduced by Aker et al. [25] for drainage displacements in a regular network. Over the last two decades, new mechanisms have been developed to extend the model for the steady-state flow as well as for irregular networks. A detailed description of this model in its most recent form can be found in Gjennestad et al. [26,27] and Sinha et al. [28] and we therefore describe it here only briefly.
The porous medium is represented by a network of links that are connected at nodes. All the pore space in this model is assigned to the links and, hence, the nodes do not contain any volume, they only represent the positions where the links meet. The flow rate q j inside any link j of the network at any instant of time for fully developed viscous flow is obtained by [33,34], where p j is the pressure drop across link, l j is the link length and g j is the link mobility which depends on the cross section of the link. The viscosity term µ j is the saturationweighted viscosity of the fluids inside the link given by µ j = s j,w µ w + s j,n µ n where µ w and µ n are the wetting and nonwetting viscosities and s j,w and s j,n are the wetting and nonwetting fluid saturations inside the link, respectively. The term p c,j corresponds to the sum of all the interfacial pressures inside the jth link. A pore typically consists of two wider pore bodies connected by a narrow pore throat. We model this by using hour-glass shaped links. The variation of the interfacial pressure with the interface position for such a link is modeled by [34], where r j is the average radius of the link and x ∈ [0, l j ] is the position of the interface inside the link. Here γ is the surface tension between the fluids and θ is the contact angle between the interface and the pore wall. These two Equations (77) and (76), together with the Kirchhoff relations, that is, the sum of the net volume flux at every node at each time step will be zero, provide a set of linear equations. We solve these equations with conjugate gradient solver [35] to calculate the local flow rates. All the interfaces are then advanced accordingly with small time steps. In order to achieve steadystate flow, we apply periodic boundary conditions in the direction of flow. We construct a diamond lattice with 64 × 64 links in two dimensions (2D) with link lengths l j = 1 mm for each link. Disorder is introduced by choosing the link radii r j randomly from a uniform distribution in the range 0.1 mm and 0.4 mm. We use 10 different realizations of such network for our simulations in 2D. In three dimensions (3D), we use a network reconstructed from a 1.8 × 1.8 × 1.8 mm 3 sample of Berea sandstone that contains 2, 274 links and 1, 163 nodes [36]. Simulations are performed under constant pressure drop P across the network. The simulations are continued to the steady state which is defined by the global measurable quantities, such as the fractional flow or the total flow rate Q fluctuate around a steady average. In the steady state, we calculate the seepage velocities averaged over time. First we use direct measurements, where we measure the global flow rates (Q, Q w , and Q n ) and the pore areas (A p , A w , and A n ) through any cross section orthogonal to the applied pressure drop and then use Equations (8)-(10) to calculate the seepage velocities. Next, we perform the measurements using the differential pore areas (a p , a w , and a n ) and calculate the seepage velocities using description given in section 4. We then compare the results from the two measurements and calculate the co-moving velocities. We then verify the relation between the seepage velocities and their higher moments.
For the direct measurements, imagine a cross section at any place of the network orthogonal to the overall direction of flow. For the regular diamond lattice in 2D, all the links have the same length. Different moments of the seepage velocities can therefore FIGURE 1 | Verification of the relations (11), (45), and (59) between the steady-state seepage velocities v p , v w , and v n , and their higher moments for the 2D regular network. The top row represents the direct approach of measurements using Equations (79), (80), and (81). The bottom row corresponds to the velocities measured from the differential area distributions defined in Equations (58), (60), and (61). v p has a unit mm/s. Subsequently, for q = 2 and 3, the units for v q p will be mm 2 /s and mm 3 /s, respectively.
and v q n = j q j a j q a j S n,j j a j S n,j , where a j is the projection of the pore area of the jth link on the cross sectional plane. Here, all links have the same angle α = 45 • with the direction of the overall flow. However, in case of the irregular network in 3D, Equations (79)-(81) need to be modified as the links have different lengths and orientations. In such case, the one can calculate the seepage velocities by [28], v q p = j q j a j q a j l x,j j a j l x,j , and v q n = j q j a j q a j S n,j l x,j j a j S n,j l x,j , where l x,j = l j cos α j is the projection of the link length (l j ) to the direction of the overall flow. If we consider every link having the same length l j = l and same orientations α j = α in these equations, we retrieve the Equations (79)-(81). For the first moment (q = 1), the velocities v q p , v q w and v q n are equivalent to Q/A p , Q/A w , and Q/A n , respectively, in both 2D and 3D.
For the second approach, we construct the distribution of differential transversal pore areas a p , a w , and a n such that a p dv, a w dv, and a n dv express the transversal pore areas for the total, wetting and non-wetting fluids within the velocity range from v to v + dv, so that they satisfy Equations (38), (40), and (41). We therefore have, where j runs over all the sites satisfying the condition: v < v j < v + dv, v j being the local velocity of link j. In case of the 2D lattice, l x,j s are same for any j and given by l x,j = l/ √ 2. With these, different moments of the seepage velocities are then calculated using Equations (58), (60), and (61), respectively.
For any saturation, the seepage velocities and their higher moments should follow the relations (11), (45), and (59). We plot our numerical measurements in Figures 1, 2 for 2D and 3D, respectively. The upper row in each figure corresponds to the direct measurements and the lower row correspond to the FIGURE 4 | Measurement of the co-moving velocity (v m ) and its higher moments for the 2D network. The top row corresponds to the calculations using Equations (29) and (30) with the direct measurements. The bottom row shows the measurements of v q m using the differential area distributions with Equation (46) and compared with the direct measurements where higher fluctuations are observed. v m has a unit mm/s. Subsequently, for q = 2 and 3, the units for v q m will be mm 2 /s and mm 3 /s, respectively.  (29) and (30), and the bottom row corresponds to the measurement from the differential area distributions using Equation (46). Here, larger fluctuations in the results calculated with the differential pore area are observed compared to the 2D network. v m has a unit mm/s. Subsequently, for q = 2 and 3, the units for v q m will be mm 2 /s and mm 3 /s, respectively. measurements from the differential area distribution. A good agreement with the relations can be observed for first as well as for the higher moments for both the networks.
Next we measure the fluctuations in the seepage velocities which obey Equation (74). Numerically, v 2 p , v 2 w , and v 2 n are calculated from the knowledge of the 1st and 2nd moments by, In Figure 3, we plot these fluctuations for the two networks to compare with Equation (74) and good agreements is observed. There are some deviations in the results for the Berea network, since the results in 3D is based on only one network configuration whereas the results for 2D are averaged over 10 different configurations. Finally, we verify the relations between seepage velocities and their higher moments while varying the fluid saturation as given by the Equations (29), (30), and (62). For this, we first calculated the co-moving velocity (v m ) and its higher moments from Equations (29) and (30) where we used the values of the seepage velocities measured with the direct approach. This is shown in the top rows of Figures 4, 5 for 2D and 3D, respectively, which show good agreements with Equations (29) and (30). We then compare these values of v q m with the measurements from the differential transversal areas using Equation (46). For this, we first constructed the histogram for the differential pore area a m corresponding to the co-moving velocity from Equation (47) where we have used the variations of a p , a w , and a n with the saturation S w . For this purpose, we have considered 21 different values of saturations within 0 and 1 with an interval of 0.05. We then integrate a m from −∞ to ∞, weighted by the velocity and normalized by the total pore area to obtain the desired comoving velocity with Equation (46). These results are plotted in the bottom row of Figures 4, 5 where they are compared with the results from direct measurements. The data points roughly follow the diagonal straight line showing satisfactory agreement with the theoretical formulations. However, we observe deviations in the results that is higher compared to the direct measurements. We believe this is due to the numerical errors that added up from several steps in the calculation such as the binning techniques while measuring the distributions, taking the derivatives and calculating the integrals. Moreover, the fluctuations for 3D are much higher compared to 2D, which is due to the lack of averaging over different samples as we have already mentioned earlier.

SUMMARY
The aim of this paper is to provide the link between the pseudothermodynamic theory at the continuum level developed in Hansen et al. [21] (see section 3) and the velocities occurring at the pore level during immiscible two-phase flow in porous media. This link is provided by defining the differential transversal pore areas defined in section 4, which essentially correspond to the statistical distributions of velocities at the pore level. The central quantities are the velocity differential transversal pore area a p , the wetting fluid differential velocity transversal pore area a w , the non-wetting fluid velocity differential transversal pore area a n , and the co-moving velocity differential transversal pore area a m . We also consider the thermodynamic velocity differential transversal pore areasâ w andâ n . The relations found by Hansen et al. [21] for the average seepage velocities, the co-moving velocity and the thermodynamic velocities are generalized to the differential transversal areas here. In the following section 5, the relations are generalized to higher moments of the velocity distributions.
The theoretical derivations are then in section 6 validated by numerical simulations. We used dynamic pore-network modeling where an interface-tracking model is used to simulate steady-state two-phase flow. We used both regular pore networks and an irregular pore network reconstructed from a Berea sandstone for our simulations. By measuring the seepage velocities from the differential area distributions and comparing them with the direct measurements, we validated the essential predictions from the earlier theoretical sections.
Both Hansen et al. [21] and the present paper are to be seen as installments toward a theory for immiscible flow in porous media at the continuum scale. The structure of this theory reflects that found in thermodynamics: A set of general relations between the macroscopic variables based on energy conservation (i.e., the Gibbs relation) and Euler homogeneity. These general equations then have to be complemented by an equation of state which introduces the specifics of the system at hand. In the immiscible two-phase flow theory we are presenting here, Euler homogeneity and mass conservation provide the general equations that transcend the specifics of the porous medium. These general equations then have to be complemented by the constitutive equations for v p and v m , which provide the specifics of the porous medium.
The resulting set of equations may then be solved for structured porous media where the structure are associated with length scales larger than that set by the REV. This is e.g., seen in the explicit appearance of the porosity φ in Equations (34)- (36).
An open question, though, is what happens when there is non-trivial structure in the porous medium all the way from the pore scale to the continuum scale, see [37] and [38]-or when the saturation of the system is at a critical value, see [36]. The fundamental Euler scaling assumption (12) would then need to be modified, and with it, all the ensuing equations.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
SR did the numerical simulations and analysis. SS developed the codes and performed the 3D simulations. AH developed the theory.