A combined fluid-dynamic and thermodynamic model to predict the onset of rapid phase transitions in LNG spills

Globally, transport of liquefied natural gas (LNG) by ship occurs on a massive scale. The large temperature difference between LNG and water means that the LNG will boil violently if spilled onto water. This may cause a physical explosion known as a rapid phase transition (RPT). Since RPT results from a complex interplay between physical phenomena on several scales, the risk of its occurrence is difficult to estimate. There is need for a tool to quantify this risk in a straightforward manner, provided that the spill scenario and geometry are well known. In this work, we present a combined fluid-dynamic and thermodynamic model to predict the onset of delayed RPT. On the basis of the full coupled model, we derive analytical solutions for the location and time of delayed RPT in an axisymmetric steady-state spill of LNG onto water. These equations are shown to be accurate when compared to simulation results for a range of relevant parameters. The relative discrepancy between the analytic solutions and predictions from the full coupled model is within 2 % for the RPT position and within 8 % for the time of RPT. This provides a simple procedure to quantify the risk of occurrence for delayed RPT for LNG on water. Due to its modular formulation, the full coupled model can straightforwardly be extended to study RPT in other systems.


Introduction
Natural gas is a common fossil fuel used for heating, cooking, propulsion and electricity-generation across the globe. Its main component is methane (about 90 wt%), with the remainder consisting of progressively smaller amounts of the heavier alkanes, such as ethane, propane, and butane. For the purpose of long-range transportation, natural gas is cooled in large refrigeration systems to form liquefied natural gas (LNG) [1]. LNG is transported across the world's oceans in large carriers, and a single carrier may carry up to around 260 000 m 3 of LNG. Combined with an increasing trend towards both processing and usage of LNG at sea [2], this means that there are many scenarios where LNG may be accidentally spilled and come in contact with seawater.
The chain of events of a marine LNG spill consists of several steps as illustrated in Figure 1. Initially, there must be a containment breach where the containment of LNG in a tank or transfer line is broken. If the breach is above sea level, the LNG may fall towards the water surface, and when this jet impacts the water surface it would break up. If the momentum of the jet is large enough, the droplets initially penetrate the surface and become submerged in water, which forms a chaotic mixing region. Since the density of LNG is about half that of water, the droplets are buoyant and rise to the surface. Subsequently, the spill forms a pool that spreads on top of the water surface. The boiling point of LNG is about −162 • C, so the pool starts to boil while spreading. Since methane is the most volatile component, Figure 1: An illustration of an accidental spill event during loading or unloading of LNG. The flow chart shows the possible pathways of the different kinds of RPT events. The jet from a containment breach impacts the water and forms a mixing region beneath the surface. At this location, there is a known possibility of early RPT. Since the LNG density is lower than that of water, the LNG droplets rise and form a floating pool. This pool then spreads out while simultaneously changing its composition via boil-off. Eventually, this change may result in a delayed RPT. In this work, we focus on the delayed RPT. the resulting vapor is almost pure methane. This causes a gradual compositional change, which increases the relative amounts of the heavier alkanes such as ethane, propane and butane.
In most cases, the spill proceeds to complete boil-off. However, in some cases it is observed to suddenly, and seemingly at random, undergo a localized explosive vaporization. This phenomenon was discovered in the 1960s and is called a rapid phase transition (RPT), although sometimes referred to as a vapor explosion or physical explosion. However, an RPT is not an explosion in the common meaning of the word, since it does not involve combustion or other chemical reactions. It is still destructive in nature and poses a danger to both people and equipment. Its peak pressures and released mechanical energy can be large enough to displace and damage heavy equipment [3][4][5] and could cause secondary structural damage and cascading containment failures [6]. The yields of single RPT events vary greatly, and may have TNT equivalents ranging from a few grams to several kilograms, which corresponds to about 0.01-25 MJ [7][8][9][10][11]. As such, RPT is considered one of the main safety concerns in the LNG industry [4,12]. Still, the attention given to RPT risk in LNG safety reviews is highly variable, ranging from significant discussions [3,4,13,14] to little more than a brief mention [5,6,[15][16][17]. There may also be a risk of RPT occurring during transport of other cryogenic liquids such as liquid hydrogen [18,19], argon, helium, nitrogen, and oxygen [12]. There is an increasing focus on hydrogen as a clean energy carrier, and liquid hydrogen at about −250 • C has been described as a promising mode for large-scale transport [18]. This merits further development of models and theoretical predictions for the occurrence of RPT in cryogenic fluids.
As illustrated in Figure 1, there is an established distinction between two kinds of RPT events depending on when and where it occurs [3,7]: Early RPT, defined as any RPT that occurs in the mixing region at any time during the spill event; and delayed RPT, defined as any RPT that is not an early RPT. Experiments indicate that delayed RPT only occurs a considerable time after the start of the LNG spill event (on the scale of minutes) [8]. The figure also illustrates the pathways along the previously discussed chain of events for the two different kinds of RPT.
In the first half of the 1970s, researchers arrived at a general consensus for a theory of RPT [20][21][22][23][24][25] which is described in depth by Aursand [26,Ch. 1]. The main idea is that the LNG pool on water is initially within the filmboiling regime. In this case, the LNG pool is insulated from the water by a vapor film that keeps the heat flux low. When the LNG composition changes due to boil-off of the volatile components, mostly methane, the film collapses due to a change of boiling regime [27]. The temperature at which this change occurs is known as the Leidenfrost temperature. This leads to direct contact between water and LNG and a rapid increase in the heat flux. The LNG pool may then be superheated and transition through a meta-stable state until the liquid approaches the superheat limit [28]. At the superheat limit, the LNG will spontaneously vaporize by homogeneous nucleation. This is called a rapid phase transition and will create a pressure wave and potentially release considerable energy through expansion work.
Whether or not an RPT event will occur in any given spill has been notoriously difficult to predict. From extensive tests performed by Lawrence Livermore National Laboratory (LLNL) in the 1980s [3,7,8], it was found that RPT occurred in about one third of all spills. It was also observed that a single spill may lead to more than ten distinct RPT events.
Cleaver et al. [13] give a summary of the experimental data on LNG safety. A more thorough history of LNG RPT research and experimental measurements is given by Aursand [26,Ch. 1]. However, although there have been extensive testing, there is still a shortage of available data on RPT events for LNG on water. In particular, there is a lack of data that is suitable as a reference for advanced models, especially with regard to RPT events [29].
The main challenge when predicting the occurrence of RPT is to predict the sudden film-boiling collapse. We will refer to this as the triggering event. The approach depends on whether one considers early RPT (droplet boiling) or delayed RPT (pool boiling). The scope of the present paper is limited to delayed RPT. In particular, we deal with the modelling and prediction of delayed RPT events as a consequence of some containment breach.
A comprehensive review of earlier work on modelling LNG spills is provided by Webber et al. [30]. Another detailed overview is provided by Hissong [31]. The main focus of earlier works has been to model the spreading, evaporation, and dispersion. To our knowledge, there has been little efforts towards predicting the onset of delayed RPT events as a consequence of spills. One exception is Horvat [32], who presents a CFD methodology based on the homogeneous multiphase formulation for simulation of LNG spills and RPT. He uses this methodology to simulate the behaviour of an LNG spill from point of release, with spreading, until the LNG evaporates and disperses. He applies the Leidenfrost criterion to predict RPT events, but there is no further discussion or characterisation of the RPT events in particular. Instead, the main focus is to predict the evaporated LNG cloud dispersion and flammability.
In our previous work [33], we proposed a method to determine how much boil-off of methane that is necessary to trigger RPT. However, the method does not quantify when and where this criterion is reached. This requires a more elaborate simulation of the spill event, which is the purpose of the present work. In particular, a complete model for studying delayed RPT must include pool formation and spreading, heat transfer between the water and the LNG, and evaporation of LNG. A thermodynamical equation of state (EoS) is necessary to describe the thermophysical properties and solving the multicomponent liquid-vapor phase transitions.
In this work, we present a simulation model for predicting delayed RPT. The model relies on the shallow-water equations to capture how LNG spreads on water and includes models for evaporation, heat-transfer, and a state-ofthe-art EoS. The Leidenfrost criterion is used to predict the risk of RPT. The model and the numerical method used to solve it are presented in Section 2. In Section 3, we present an analysis of a continuous tank spill and derive simple, predictive models for the position and time of delayed RPT for this particular case. The methodology for predicting the RPT time and position is general and applicable to any steady-state scenario. A two-dimensional axisymmetric spill is analysed in detail and explicit estimates for RPT position and time are derived. Next, in Section 4, we consider the axisymmetric spill case with constant spill rate and show results from simulations. The results are compared to the predictions from the analytical solutions derived in Section 3. The results show that the analytical models accurately reproduce the predictions from the more advanced simulation model. Finally, in Section 5, we summarize the results and provide an outlook of possible future work.

LNG spill model
In the following, we present a model to predict the flow of LNG after a spill event. The model accounts for important mechanisms such as boiling and evaporation, as well as enrichment of the LNG mixture due to evaporation, as shown in Figure 2.
When the LNG is spilled on water and starts to spread, it evaporates. Vapor bubbles are continuously formed at the water-LNG interface and rise to the surface. The bubbles displace the liquid phase, which reduces the average density of the liquid-vapor mixture. The reduced density implies that the volume and thus the height of the LNG-air Pool-spreading model surface increases. This in turn affects the spreading rate that depends on the leading-edge height [34]. Because of the coupling between spreading rate and vapor volume, the flow model must include the dynamics of the vapor mixed into the LNG.
Throughout this section, we consider axisymmetric spill geometries where the local fluid states are depth-averaged at each position throughout the pool. Because of this, all properties are functions of time t and radius r. We take the origin r = 0 to be the center of the spill source, and t = 0 to be the start of the spill event.
We assume that the liquid mass is much larger than the mass of the vapor bubbles entrained in the LNG, and that the density of this vapor is nearly constant. We further assume that the pressure is atmospheric, p = p atm , and thermal equilibrium such that the LNG temperature is always at the bubble point, that is, T = T bub (p atm , z), where z is the LNG component mass-fraction vector which changes with time.

Pool spreading
The dynamics of an LNG pool spreading on water may be captured accurately by the two-layer shallow-water equations (SWE) [35,36]. These may be derived from the Euler equations for two liquid layers by assuming a negligible vertical velocity within each layer. It is also common to use various effective forms of the single-layer SWE modified to represent the two-layer system, see e.g. Refs. [37,38]. In this work, we use a form of the singlelayer SWE derived and analysed by Fyhn et al. [34]. This form of the SWE has the benefit of being simple while still capturing essential features of the two-layer equations for liquid-on-liquid spreading. In particular, it can be used without any additional boundary condition for the leading edge of the spill, as long as the water depth is large relative to the thickness of the LNG layer.
We employ the SWE to model the flow of an LNG liquid-vapor mixture, that is, Here m i is the liquid mass per area of component i, u is the horizontal flow velocity, h v is the vapor volume per area, m spill i is the spill source term,ṁ evap i represents the liquid mass loss due to evaporation,ḣ evap v represents the increase in the vapor phase due to evaporation,ḣ release v represents the decrease in the vapor phase due to vapor bubbles released to atmosphere, h is the thickness of the liquid-vapor mixture, and g eff = δg is the effective gravitational acceleration. δ is the buoyancy factor, where ρ w is the water density and ρ is the overall density of the liquid-vapor mixture, Here we have assumed that the vapor mass is negligible, which is a good assumption as long as there is no foaming. The liquid-vapor mixture thickness is given by where ρ l = ρ l (T bub , p atm , m i ) is the liquid LNG density. Both ρ l and T bub are calculated from the EoS, see Section 2.5. As shown by Fyhn et al. [34], the water-LNG and LNG-air interface positions are given by where the reference h w = 0 is the water level where there is no LNG.

Evaporation
LNG spilled on water immediately begins to evaporate due to the large temperature difference. Throughout this work, we assume that the absorbed heat from the water into the LNG goes to evaporation and to heating of the LNG as the boiling point is shifted. Evaporation produces vapor that passes through the LNG before it is released. The effect of evaporation is represented in the flow model (1) through the source termsṁ To estimateṁ evap i , we will assume that the composition of the evaporating gas is that of the incipient vapor phase at the bubble point. The change in mass of component i in the liquid due to evaporation is given bẏ where H is the specific enthalpy. To calculate dm i /dH numerically, we start from a saturated liquid state, add a small amount of enthalpy ∆H, then calculate the corresponding liquid-vapor state from the EoS, see Section 2.5. Based on this state, we find the change of liquid mass for each component ∆m i , which is always negative. Next, since the pool is a system at constant pressure, addition of heat causes a corresponding increase of the enthalpy. This implies that whereq is the heat flux, which is further discussed in Section 2.3.
The vapor volume generated per area is given directly from the evaporation rate, where ρ v is the vapor density as predicted by the EoS. The vapor bubbles rise through the pool due to buoyancy at an average velocity u B . We have used u B = 24 cm/s, since this was measured by Chang et al. [39] for pure methane.
The timescale for bubble migration through the pool is relatively short, so we assume that the vapor bubbles are immediately dispersed uniformly throughout the liquid-vapor mixture. In this case, the release rate of vapor volume per surface area is given byḣ Due to the short timescale of bubble migration compared to the timescale of the spill flow, a quasi-steady state is quickly established where the rate of vapor released from the top of the pool is equal to the rate of vapor added at the bottom due to evaporation, i.e.ḣ evap v =ḣ release v . In this state, the vapor volume-fraction in the pool must be We further obtain a minimal mixture density which is the same as that derived by Chang et al. [39], When the LNG film is thin, the timescale of bubble migration is small compared to the timescale of the pool spreading. We therefore assume h v = h max v in the regions where h < 1 cm.

Heat-transfer model
The heat flux between a heat source and a boiling liquid strongly depends on the temperature difference. For large temperature differences, which is typically the case for LNG boiling on water [40], the heat transfer is within the film-boiling regime. When the methane fraction in the LNG decreases due to evaporation, the temperature difference drops gradually, and at some point the boiling enters the transition boiling regime. However, as will be explained in Section 2.4, a change of boiling regime from film boiling to transition boiling is thought to be part of the mechanism for RPT triggering. We therefore restrict the heat-transfer model to the film-boiling regime, since the main goal is to investigate the onset of RPT triggering.
Standard film-boiling heat-transfer correlations, such as those by Berenson [41] and Kalinin et al. [42], tend to underpredict the heat transfer between light hydrocarbons and metallic heaters. Furthermore, mixture effects on boiling are known to be significant, even for nearly pure fluids [43,44]. Experimental results on the boiling regimes of LNG [45] and pure methane [46] produce very different boiling curves. When the heater material is changed from a metal plate to water, the heat-transfer coefficient increases by a factor of 1-4 [40,[44][45][46][47][48]. There are several factors that can contribute to this, such as increased LNG to heater area because of the uneven interface or big changes in the thermal transport properties of the heater. Further experimental investigations are necessary to reduce this significant uncertainty.
In this work, we have used a film-boiling model developed by Sciance et al. [49] for light hydrocarbons, with Figure 3: Illustration of the boiling curve for saturated LNG pool-boiling. This is a plot of boiling heat fluxq against a variable surface temperature T . Of particular interest is the Leidenfrost temperature T L , which marks the lower end of the film-boiling regime. Also shown is the typical temperature of seawater T w , which is significantly above the Leidenfrost temperature of LNG. However, during boil-off (removal of methane) the Leidenfrost point will shift to the right, eventually crossing the water temperature (film-boiling collapse).
Here c 1 and c 2 are model parameters, k is the thermal conductivity, l c is the capillary length, Ra is the two-phase Rayleigh number, ∆H is the difference in specific enthalpy between the LNG and the vapor film, T r = T bub /T c is the reduced temperature of the mixture, T c is the critical point temperature, c p is the isobaric specific heat, ∆T = T w −T bub is the LNG-water temperature difference, ∆H evap is the specific enthalpy of evaporation, g is the acceleration of gravity, ρ is the density, µ is the dynamic viscosity, and σ is the interface tension of the LNG liquid-vapor interface. The subscripts l and v indicate that a property should be evaluated for the LNG liquid and vapor phase at the bubble point, respectively. The subscript vf indicates that a property should be evaluated at the vapor-film temperature. The model parameters used by Sciance et al. [49] are c 1 = 0.369 and c 2 = 0.267. We found that c 2 = 0.3 allowed us to better fit LNG-water heat-flux data [40,[50][51][52].
The thermodynamic properties and the transport properties, i.e. the viscosity, thermal conductivity, and surface tension, are described in further detail in Section 2.5.

Prediction of RPT triggering
The fundamental triggering condition for delayed RPT is dealt with in our previous work [33]. It may be summarized by the following triggering window, where T SHL is the superheat limit temperature of LNG, T w is the temperature of seawater, and T L is the Leidenfrost temperature of LNG. This criterion is deceptively simple, as it appears to only require the evaluation of three numbers. However, both T SHL and T L depend on the LNG composition. This criterion is not satisfied for typical LNG compositions. This is due to the fact that both T L and T SHL are far below the water temperature. However, as boil-off proceeds, both T L and T SHL gradually increase while maintaining T SHL < T L . This eventually causes Eq. (15) to be satisfied, in which case RPT may be triggered. Figure 3 shows the boiling curve of LNG and indicates how boil-off shifts T L towards T w . As T SHL < T L is maintained during the boiling process, calculating T SHL is not required. Note that actual triggering of RPT is not part of the spill model. That is, when the RPT criterion is satisfied, the spill simulation is continued as if RPT did not occur. Instead, we track the amount of LNG that satisfies the triggering criterion, that is, where m = (m 1 , . . . , m n ) is the mass vector of the n LNG components. The Leidenfrost temperature T L is defined in Section 2.5.

Thermodynamics and transport properties
An EoS is used to predict the thermodynamical properties of the mixture at a given state. In this work, we use an extended corresponding states EoS [53], since it represents a good compromise between accuracy and computational speed [54]. We use the Peng-Robinson EoS [55] to calculate the shape factors and methane as a reference fluid described with the modified Benedict Webb Rubin (MBWR) EoS [56]. The accuracy of the EoS is good for mixtures with a high methane content. For the predicited densities and heat-capacities of pure methane, the error is below 2 % for gas and supercritical states, and approximately 4 % in the liquid region when compared to GERG-2008 [54].
For the selected EoS, multiple implicit thermodynamic properties must be determined numerically. Especially, the LNG bubble temperature T bub is predicted at atmospheric pressure for a given composition of n components (m 1 , . . . , m n ). At specified atmospheric pressure, specific enthalpy, and total composition, the thermodynamically stable phase distribution and temperature must be calculated. The latter follows from a pHz flash, where the solution is found iteratively through an implicit relation, see Ref. [53] for details.
The Leidenfrost temperature of a fluid depends not only on the fluid properties, but also thermal transport properties of the heater [57], geometry and interface roughness [58]. For pure fluids using the van der Waals EoS, Spiegler et al. [59] derived a simple relation that relates the Leidenfrost temperature at sub-critical pressures to the critical temperature: T L = (27/32)T c . This approach was also adopted by Aursand and Hammer [33], who showed that this simple correlation gave an excellent prediction for methane, and generally a decent upper estimate for other hydrocarbons. Due to lack of data, the same correlation was used for mixtures, without any validation. As (27/32)T c corresponds to the liquid spinodal of the van der Waals EoS when p → 0, we have chosen to correlate T L = T spinodal (m, p) in this work. The method produces a similar composition dependence as that reported by Yue and Weber [43] for binary mixtures. The liquid spinodal was calculated as the temperature where the smallest eigenvalue of the Gibbs energy Hessian matrix with respect to mole numbers become zero [28,60].
Since the liquid is assumed to always be saturated, the thermodynamical properties of the vapor phase are calculated from the EoS by finding the incipient vapor phase at the corresponding temperature and pressure.
The viscosity and thermal conductivity are calculated using the TRAPP method by Ely and Hanley [61,62], with propane as a reference fluid. The surface tension is calculated using the corresponding-state correlation presented by Poling et al. [63]. These models take density as an input, and so the EoS influences the accuracy of the transport properties. Typically, we get errors less than 5 % for the transport properties for both the liquid and vapor phase.

Numerical implementation
In the axisymmetric case, the SWE (1) may be reduced to These equations are discretized spatially using a finite-volume scheme. We employ the FORCE (first-order centered) flux [64] and the second-order MUSCL (Monotonic Upstream-Centered Scheme for Conservation Laws) reconstruction with a minmod limiter [65] in each finite volume. The solutions are advanced in time with a standard third-order three-stage strong stability-preserving Runge-Kutta method [66]. 8 The time steps are restricted by the Courant-Friedrichs-Lewy (CFL) condition with a CFL number of 1.0 for all cases such that where subscript j indicates grid cell j. In addition, we also restrict the time steps to account for the spill source term, The denominator is regularized to avoid division by zero in e.g. scenarios with no initial LNG. Unless otherwise specified, we apply zero-gradient boundary conditions for all variables at the boundaries, except that we that we specify u = 0 at r = 0.

Analytical estimates
Detailed models such as the one presented in Section 2 are useful for improving our understanding of fundamental phenomena. However, for practical applications such as safety assurance, a simplified model is desired. Such models are easier to understand and apply, require fewer input parameters, and often produce satisfactory results. In this section, we develop a method to predict the radius of RPT for steady-state spills of LNG. The method is used to derive an analytical estimate of the RPT radius and triggering time for an axisymmetric two-dimensional tank spill of LNG on water. This spill scenario could be the result of an operational failure in a loading arm during loading or off-loading of an LNG tanker to an on-shore facility.
We define the boil-off limit as, wherez i is the mass fraction of chemical component i at RPT triggering measured relative to the initially spilled mass. Thus, θ is the mass fraction of the initially spilled LNG that must evaporate before an RPT triggering can occur. This boil-off limit only depends on the spill composition and evaporation model, and can be determined in many different ways. Here we calculate it numerically from the condition T L (z) = T w , which ensures that the RPT triggering criterion from Section 2.4 is satisfied. Alternatively, it can be analytically estimated to within a few percent from the spill composition alone using the correlations derived in Ref. [33]. It is also possible to measure it experimentally by placing LNG into a container with enough water, and record how its weight changes until triggering. It is worth noting that if only methane evaporates from the LNG mixture, the boil-off limit can also be written in the simpler form θ = z 1 −z 1 , where z 1 andz 1 are the methane mass fractions at spill and triggering, respectively. Our model from Section 2 and the following analysis is applicable to any multicomponent mixture. However, for simplicity, we will henceforth consider LNG mixtures z = (z 1 , z 2 , z 3 ) that have three components: methane (CH 4 ), ethane (C 2 H 6 ), and propane (C 3 H 8 ). We consider three such mixtures, which we denote z A , z B , and z C . The assumed compositions and calculated boil-off criteria for these mixtures are listed in Table 1.

Methodology for triggering prediction
Consider a spill scenario, such as a containment breach in a large tank or a rupture in an LNG transfer line, which causes LNG to leak onto sea water. If the spill continues at a steady rate S , the flow should approach an approximately Table 1: List of initial LNG compositions and the calculated boil-off limit for RPT θ.
Contours for the corresponding boil-off fractions Γ(t). Figure 5: Sketch of the RPT prediction strategy for continuous spills. (a) Here, we illustrate a cylindrical spill source in a wide channel, resulting in 2D spreading followed by 1D spreading. If the steady-state velocity field u(r) is known, the path r(t) that a control mass takes (red line) can be calculated by integrating u(r) from a point r(0) on the spill boundary (red circle). By repeating this for each such point r(0), all paths r(t) from the source can be mapped out (gray lines). We can then also describe t(r), i.e. the time LNG spends moving from the source to any given point r.
(b) Since a roughly constant amount of methane evaporates per unit time the LNG is on sea water, the fraction Γ that has evaporated at each point r can be calculated. By comparing the resulting contours Γ n to the boil-off limit θ, we can predict the region at risk of an RPT event (yellow). steady-state pattern, as shown in Figure 5. In this limit, the LNG mass distribution m and velocity field u become nearly time-independent.
Assume that the flow patterns are not affected by evaporation, that is, that the flow velocity u is the same regardless of evaporation. Then we may replace m in a case without evaporation with (1 − Γ)m in a corresponding case with evaporation, where Γ is the fraction that has boiled off at each position. This is usually a reasonable assumption since the amount of methane that evaporates per unit area and time is roughly constant, so the methane fraction changes faster where the LNG is spread thin, i.e. far from the spill source. This is also observed by numerically solving the more advanced model, as the chemical composition remains nearly constant throughout most of the spill region (see Figure 8a). Thus, if we neglect evaporation, the steady-state m(r) and u(r) can be obtained e.g. via computational fluid dynamics or analytical approximations based on mass conservation.
We now consider an arbitrary spill geometry, and derive a general equation to predict when and where RPT events may occur. This is achieved by decoupling the model of how LNG spreads from the model of how it evaporates and eventually undergoes RPT. This is an approximation that can enable simple RPT predictions in realistic geometries. The predictions are validated numerically in Section 4.
Methane typically evaporates much faster than the other components in LNG. We therefore approximate the specific enthalpy of evaporation ∆H by the methane value ∆H 1 , which is listed in standard chemical tables. Furthermore, we assume that the heat fluxq from sea water into LNG is roughly independent of position and composition. With these simplifications, the equations in Section 2.2 yield the following methane evaporation ratė The evaporation rates of all other chemical species in the spilled LNGṁ evap i have been neglected. The enthalpy contribution of mixing methane with the LNG, is also neglected.
Given the above assumptions, we next track how a small control mass flows through the system. If the control mass starts at a point r(0) at time 0, the path r(t) it traverses is given by the steady-state velocity field u(r): This equation may be solved implicitly to find a specific position r(t) at any given time t. It may be challenging to solve an implicit equation analytically, but it should be straightforward to solve numerically for any geometry.
As mentioned previously, the methane evaporation rateṁ evap 1 is assumed to be constant. However, how fast this changes the composition of LNG depends on position, since the mass per unit area m(r) is not constant. We can relate the mass fractions of each component z i to the mass components as m i = mz i , which for methane implies that the change in composition isż 1 =ṁ evap 1 /m. Note that we defined z using the mass profile m(r) calculated in the absence of evaporation. Thus, as mass evaporates, this definition of z implies that the sum of components i z i < 1, and the "missing mass" is the boil-off fraction Γ that plays a central role in RPT prediction below.
The total amount of methane that has evaporated from the control mass can then be found via integration oḟ z 1 =ṁ evap 1 /m along the path r(t). Moreover, we can substitute in Eq. (22) forṁ evap 1 . This leaves us with: .
This is negative since methane disappears during evaporation. As only methane evaporates, the boil-off fraction is: .
As introduced in Eq. (21), the timet RPT the control mass has to travel before an RPT event occurs is defined according to the boil-off limit Γ(t RPT ) ≡ θ, which implicitly definest RPT for a given path r(t). The corresponding boundary to the region where an RPT event might occur is then given by r RPT ≡ r(t RPT ). In both cases, the path r(t) is defined by the initial position along the spill boundary r(0); to obtain the full RPT risk region, one therefore has to repeat this consideration for each point r(0) along the spill boundary.
To summarize, the region at risk of RPT can in general be determined from a known mass distribution m(r) and velocity field u(r) by solving these equations for each point r(0) along the spill boundary: .
These equations may be solved for the path r(t) and the triggering timet RPT . The locations at risk of delayed RPT events are then given by r RPT ≡ r(t RPT ). Note thatt RPT is not the time it takes from the initial spill until a delayed RPT is possible. This time is typically longer, since the initial spreading will be restricted by the water.
In order to estimate the time t RPT from the initial spill starts until a delayed RPT is possible, we use the method of characteristics. In the absence of evaporation, Eq. (1) can be written where c = g eff h,û is a unit vector along u, and u = |u| is the magnitude of u. Combining these equations, we get From this we can connect the velocity at the leading edge to the steady-state velocity distribution. The characteristics described by Eq. (29) will move faster than the shock at the leading edge. Hence, we can find an equation that relates the height and velocity at the shock with the height and velocity further back where the solution can be assumed to be better approximated by the steady state solution. Integrating Eq. (29) from t a to t b , where u[r(t b )] = u LE , gives Here we assumed that u LE =n · u LE , which from the Rankine-Hugoniot condition implies that u LE = 2g eff h LE , as shown in Ref. [34]. If we assume that u[r(t a )] and h[r(t a )] are well approximated by the steady state solutions and t b − t a is sufficiently small so that we may neglect the contribution from the integral, we can write the time from initial spill to a possible RPT as where the integration contour r(t) starts at r(0) and ends at r RPT . In general, each point r(0) along the spill boundary can give a different value of t RPT , so Eq. (31) should be minimized with respect to r(0). This result can be used to perform RPT predictions for any geometry, and can in practice be coupled to a computational fluid-dynamics simulation as a post-processing step. The idea is also summarized in Figure 5. However, we note that a central assumption behind the results above was that evaporation is not very important for the flow pattern. In geometries where evaporation effects are more dominant to the flow patterns, the above formulation may yield less accurate results. This is further discussed in Section 4.

Application to an axisymmetric 2D spill
To illustrate the methodology derived in the previous subsection, we will now apply it to an axisymmetric twodimensional spill of LNG on the ocean. The spill is assumed to occur within some radius r < r 0 from the origin, see Figure 6, while the earliest point an RPT event may occur is at some much larger radius r = r RPT . In the following, we first apply the above method to determine r RPT . We then expand the analysis and derive estimates for the time t RPT when delayed RPT events first become possible.
In the absence of evaporation, mass conservation implies that the following relationship between the net spill rate S , the radial velocity u, and the mass m for r > r 0 : This equation is solved for 1/m, which substituted into the boil-off equation (27) yieldṡ However, u dt = dr by definition, and 2πr dr = π(r 2 − r 2 0 ). This can be used to rewrite the equation in terms of the position r RPT . Moreover, since we are interested in regions r RPT r 0 , we can also neglect the r 0 contribution. This leaves us with an equation for the RPT location r RPT : This result requires only three input parameters to estimate r RPT : the heat fluxq, spill rate S , and boil-off limit θ. It also has a straightforward physical interpretation: It can be reformulated as A RPT = θS /µ, where µ =q/∆H 1 is the rate at which mass disappears per unit area. The result A RPT = πr 2 RPT is simply the surface area required for a fraction θ of the spilled material S to evaporate, in accordance with the boil-off limit. We can then easily solve for the boundary r RPT of the region A RPT . Beyond that boundary, the boil-off criterion is satisfied.
To estimate t RPT , we will find the steady-state velocity distribution u(r). By integrating the shallow-water equations (18), it can be found that the steady-state height and velocity distributions must satisfy where ε is a constant to be determined,ṁ spill is the spill rate density, andṀ spill (r) = S for r ≥ r 0 . We solve Eqs. (35a) and (35b) to obtain a polynomial equation of degree three in h, The solutions of a third-order polynomial equation can be found by Cardano's formula [67]. From the three roots, we pick the desired solution to fit the boundary conditions, that the height at r → ∞ should be smaller than that at r = 0 and that u = 0 at r = 0. We find that where

13
Next, we can determine ε by demanding that the solution is continuous at r = r 0 . This is the case if χ(r 0 ) = 0, which happens when We are now in a position to estimate t RPT from Eq. (31). If r RPT r 0 , the integral will be dominated by its large-r contributions. Taylor expanding the integrand to first order in r 0 /r gives If we now integrate this from r 0 to r RPT , we get where we have written the result in terms of a function with R = r RPT /r 0 . Note that for R → ∞, f (R) → 1. In this limit, t RPT ∼ r RPT /u ∞ , where u ∞ = √ 2 is the steady-state velocity for R → ∞, giving an intuitive "distance over speed" result for t RPT . For finite R, f (R) includes more low-R corrections.
The analytical solution we derived for h(r) can also be used to estimate the mass distribution m(r). For r > r RPT , a fraction θ of the spilled LNG has already evaporated, so the mass distribution is m(r) ≈ (1 − θ)ρh(r), where h(r) is the solution derived in the absence of evaporation and is given by Eq. (37). By integrating the resulting mass per radial shell 2πr m(r) from the risk boundary r RPT to the leading edge of the spill r LE , one can also estimate the total mass M RPT at risk of RPT. This may then be combined with the analytical results for the explosive pressure and energy derived in Ref. [33] to obtain a worst-case scenario analysis.
A similar analysis is also applicable to one-dimensional spills, which approximate spills into narrow channels. The mass conservation equation then changes to 2d u(x) m(x) = S , where d is the channel width and x the coordinate along the channel, but the derivations follow along the same lines.

Numerical results and discussions
We still consider an axisymmetric two-dimensional spill similar to what was discussed above. This spill scenario will be used to demonstrate the applicability of the LNG spill model (18) and to compare the analytical estimates of r RPT [Eq. (34)] and t RPT [Eq. (41)] with simulation results from the full coupled model.
We assume that an LNG volume V 0 = 10 m 3 with composition z A (Table 1) is released onto sea water at a constant rate S through a hole with radius r 0 = 0.1 m. The spill lasts until t s = 30 s. The values for the spill volume and duration are based on the case descriptions from Ref. [68]. The initial density of the LNG at atmospheric pressure and bubble temperature is ρ 0 = 437 kg/m 3 . The spill rate is therefore S = 146 kg/s. The spill source term is set tȯ where S i is the spill rate of component i. Otherwise, the initial conditions are all zero. The simulations are run in a domain of radius 60 m with 100 grid points/m. The analytical estimates for r RPT and t RPT from Eqs. (34) and (41) require values for θ,q, and ∆H. The spilled LNG composition corresponds to θ = 89.1 %, as listed in Table 1. The heat fluxq can be estimated from Eq. (13), we findq ≈ 6.9 × 10 4 W/m 2 . The specific enthalpy of evaporation can be found for methane from standard chemical tables, here we use ∆H evap ≈ 510 kJ/kg. With this, we calculate the estimates r RPT 17.5 m and t RPT 15.3 s.     Figure 7a, one can see an initial shock travelling at a velocity of about 1 m/s, which is consistent with a Froude number for the leading edge of √ 2. This follows naturally from the present formulation of the SWE Eq. (1), as discussed in [34], At t = 30 s, the pool is a monotonically decreasing function of r until it vanishes at about r = 20 m due to evaporation. Once the spill rate is reduced to zero, the LNG thickness starts to reduce at the origin and outwards, as shown in Figure 7b.
In the simulations, we track the flow of each component. This is illustrated in Figure 8a, where we show the radial profile of the methane fraction at the steady state t = 30 s. The corresponding profile of the Leidenfrost temperature T L is shown in Figure 8b. As discussed earlier, triggering of RPT may occur when T L > T w . This is highlighted in Figure 8b. In both of the plots in Figure 8, one can see that the analytically estimated r RPT from Eq. (34) (orange dashed line) is very close to the simulated r RPT , which is located where T L = T w . It is interesting to note that the regions where the Leidenfrost temperature is above the water temperature coincide with regions where the thickness is small, cf. Figure 7b.
As discussed in Section 2.4, we use the Leidenfrost temperature to estimate the region where an RPT may occur. We estimate the total mass of LNG that is in risk of triggering M RPT by numerically integrating the total mass per area over the area where T w > T L , cf. Eq. (16). In Figure 9, we compare M RPT to a fraction of the total available mass M LNG as a function of time. A factor of 1/100 is used to more easily compare M RPT and M LNG . When the spill stops at t = 30 s, the amount of LNG that is in risk of triggering initially remains constant. When the thickness at the origin becomes sufficiently thin, M RPT first increases, then quickly decreases when the inner pool fully evaporates at t ≈ 38 s. To show the sensitivity of the Leidenfrost temperature, we include estimates of M RPT for T w = ±10 • C. The figure shows that the time of triggering is not very sensitive to the value of the water temperature. Finally, Figure 9 also compares t RPT from the simulation with the estimate from simulated Eq. (41). The simulation predicts t RPT 14.4 s, which is slightly lower than the estimated t RPT 15.3 s from Eq. (41).
The above case represents a fully unconstrained spill geometry. The spill geometry will have a strong influence on the flow pattern, and thus also on the RPT phenomenon. To illustrate this, consider a fully constrained case where an amount M 0 of LNG is spilled into a closed container of water where the LNG is not allowed to spread. The RPT triggering criterion is reached as soon as the boil-off ratio reaches θ. This means that the maximum amount of LNG that can trigger is a fraction 1 − θ of the inital spill (11-33 % for compositions z A -z C in Table 1). The thickness at the triggering time is approximately (1 − θ)h 0 , which is generally thicker than what can be achieved in the unconstrained case where the triggering thickness is in the mm range for realistic spill rates. The energy yield of an RPT event will scale with the thickness, and thus constrained spills have a larger blast potential than unconstrained spills.
The earlier results indicate that the analytical estimates of r RPT and t RPT from Eqs. (34) and (41) are good approximations to the corresponding predictions from simulations. To further validate this, we compared results for a set of simulations with different mass compositions and spill rates. We considered the three different mass compositions listed in Table 1 and four different spill rates S ∈ {10, 100, 250, 500} (kg/s). The remaining simulation parameters were similar to the above case description. As showed in the top row of Figure 10, the analytical estimate of r RPT matches the simulated predictions very well. The relative deviation is lower than 2 % in every prediction, as seen in Figure 11. The analytic estimate of t RPT also matches well with a deviation within 8 %. The prediction of r RPT relies only on the assumption that mass does not accumulate, while the prediction of t RPT also assumes that the velocity and height distributions close to the shock are accurately approximated by the steady state solution. It is therefore not unexpected that the prediction of t RPT is less accurate than the prediction of r RPT . Relative deviation (%) Figure 11: Boxplots of the relative deviation of analytically predicted and simulation predicted r RPT and t RPT .

Conclusions
We have presented a coupled model for predicting delayed RPT from an LNG spill onto water. The model combines the shallow-water equations for capturing the spreading of LNG on water, a film-boiling heat-flux model, a fundamental triggering condition for RPT, and an extended corresponding-states EoS for thermodynamic calculations. It was used to study an axisymmetric continuous spill case with a constant spill rate. The results give insight into how a continuous spill scenario develops into a steady flow. It further highlights how evaporation alters the LNG composition during the spill. This in turn allows prediction of the regions at risk of a delayed RPT. We also consider the influence of the spreading geometry on an RPT event. In particular, the potential energy yield in a constrained-geometry RPT event is larger than in an unconstrained geometry.
We have also analysed a continuous tank spill in more detail based on the governing equations. We derived a general methodology to predict the position and time of delayed RPT that is applicable to different geometries. The method requires a steady-state solution of the tank spill and should be easy to implement with numerical codes as a post-processing analysis method. We applied the methodology to derive simple, predictive equations for the triggering time and position of a two-dimensional axisymmetric spill.
Subsequently, we compared results from numerical simulations with the analytical predictions. Notably, the analytical prediction of the RPT location is accurate to within 2 % of that predicted from the numerical simulations. The analytical prediction of the time of RPT is within 8 % of the numerical prediction.
Unfortunately, there is a lack of relevant experimental data with which to compare the predictions. In future work, it would be very useful to perform experiments on a similar continuous spill case in order to assess and validate both the simplified predictive equations and the simulation model. This would also make it possible to suggest a semiempiric form of the prediction equations for the RPT event that is based on the present analysis. It would also be interesting to further develop predictive equations for other types of spill scenarios and other cryogenic fluids.