The MEMOS-U code description of macroscopic melt dynamics in fusion devices

The MEMOS-U physics model, addressing macroscopic melt motion in large deformation and long displacement regimes, and its numerical schemes are presented. Discussion is centred on the shallow water application to the metallic melts induced by hot magnetized plasmas, where phase transitions and electromagnetic responses are pivotal. The physics of boundary conditions with their underlying assumptions are analysed and the sensitivity to experimental input uncertainties is emphasized. The JET transient tungsten melting experiment (Coenen et al 2015 Nucl. Fusion 55 023010) is simulated to illustrate the MEMOS-U predictive power and to highlight key aspects of tokamak melt dynamics.


Introduction
For next generation tokamak devices such as ITER (currently under construction) and DEMO (currently in the planning phase), the plasma-facing component (PFC) lifetime and power handling capabilities constitute major design challenges. ITER features a beryllium (Be) first wall and a full tungsten (W) divertor [1]. In the divertor, intense steady state heat loads will elevate surface temperatures close to or above the W recrystallization range, raising concerns for edge melting and full surface melting during downward vertical displacement events (VDEs), edge-localized modes (ELMs) or loss of detachment control [2][3][4]. On the Be main chamber PFCs, the principal concern is the thermal and current quench phases of unmitigated VDEs or major disruptions.
Metallic tokamak melts are subject to strong electromagnetic forces which inevitably displace molten material, deforming the armour surface topology which could harm the power-handling capabilities by altering the incident magnetic field angle. Severe melt motion could also bridge castellations if melt fills the gaps between PFC elements. The main concern is the accumulated material excavation which may ultimately necessitate replacement of damaged components. Moreover, strong surface melting will drive immediate and deep recrystallization, leading to the propagation of macrocracks. In order to minimize the risk of PFC damage, ITER aims to achieve complete ELM suppression and to deploy an effective disruption mitigation system [3,5]. Nevertheless, some PFC melting appears to be inevitable; the Be wall is already at risk during early plasmas at low plasma current [6], while melting of W PFCs in the divertor becomes a risk later in the operational stages, when the stored plasma energies are higher [3].
In order to assess the PFC lifetime and power handling capabilities, it is necessary to acquire reliable predictions for the onset of melting, the macroscopic melt motion and the final surface deformation. This is a complex, multifaceted problem which involves the temperature, momentum and electromagnetic response of the metal to the incident plasma fluxes. Moreover, the fluid domain is bounded by two dynamically evolving interfaces, namely the plasma-facing free interface and the liquid-solid interface. These not only delimit the fluid motion domain, but also exert dynamic influence on the liquid-in particular the free surface spatial alignment and local curvature.
Thus far, fully self-consistent simulations of this class of problems have been carried out in two spatial dimensions, within the assumption of rotational symmetry or local invariance in a spatial dimension. Applications concerning gasmetal arc welding [7,8], vacuum arcs [9][10][11] and unipolar arcs [12] have been modelled successfully. Under fusion-relevant conditions, 2D simulations of W melt layer instabilities have been performed that neglect phase transitions [13][14][15][16]. Such simulations require prescribed melt pools and are valid only for short timescales.
In the context of fusion reactor walls, the main objective is to predict the mass displacement and the gross surface damage inflicted over large spatial and long temporal scales. The potentially molten PFC area can range from ∼cm 2 to tens of cm 2 while relevant time scales vary from tens of milliseconds up to seconds [6,[17][18][19][20][21]. Naturally, numerical modelling of the complete problem that features the entangled heat transfer including phase transitions and the three-dimensional liquid dynamics is a formidable computational task given the above spatiotemporal scales. However, the problem can be simplified by taking advantage of key melt flow characteristics, identifiable from experiments in contemporary machines.
It has been recognized from the early experiments in the TEXTOR tokamak that the shallow water approximation is applicable, courtesy of the very large ratio between the molten surface extent and the melt depth [22]. The shallow water approximation makes the long length and timescale problem computationally feasible by reducing the dimensions of fluid flow by one [23,24]. This feature has been implemented in the original macroscopic melt motion code MEMOS-3D [22,25] which solves 3D heat transfer in combination with the shallow water equations.
More recent evidence obtained in tokamak experiments dedicated to ELM-induced W melting in the divertor (JET [17], ASDEX Upgrade [18,19]) and VDE-induced Be melting in the main chamber (JET [20]) not only confirmed the shallow aspect of the pool dimensionality, but also brought up another crucial aspect of melt dynamics. The post-mortem analysis of these experiments revealed that frozen displaced material piles-up at distances of the order of the liquid pool extent, thus challenging the tacit assumption of MEMOS-3D [17,26,27] that the heat diffusion is much faster than the characteristic fluid motion. In scenarios where melt spills onto a progressively colder solid surface, the resolidification rate determines the viscous shear force and virtually governs melt dynamics [28]. Resolidification simultaneously stabilizes the liquid flow, as confirmed by experiments where a predominant flow direction can be discerned in the frozen material without signatures of liquid disintegration [17][18][19][20]. A new model has been developed and implemented in the MEMOS-U code, a significantly upgraded version of MEMOS-3D, aiming to address this novel fluid dynamics regime that is ubiquitous in contemporary fusion devices [28] and is also expected to be realized in melting scenarios on ITER [6].
The MEMOS-U code has been demonstrated to accurately describe the overarching physical processes behind shallow melting and surface deformation. Comprehensive benchmarking activities with four state-of-the-art tokamak experiments in JET and ASDEX Upgrade revealed that MEMOS-U is capable of quantitative predictions [28,29]. The objective of the present work is to provide a complete detailed description of the physics model and the numerical schemes of MEMOS-U, especially since such a description has always been absent for the original MEMOS series of codes and is often sought after by the fusion plasma boundary research community. The shallow water fluid equations are established in the literature [23,24] and have been widely utilized in studies of oceanic flows and atmospheric flows [30,31]. Here, discussion will focus on peculiarities in the shallow water treatment of metal melts induced by hot magnetized plasmas, where gravity is negligible while phase transitions and electromagnetic responses are pivotal. Emphasis will be placed on the physics of boundary conditions, the underlying assumptions as well as the sensitivity to experimental uncertainties. Finally, aiming to highlight the most important aspects of melt dynamics under tokamak-relevant conditions, the JET ELM-driven W melting experiment [17] will be analysed in detail, demonstrating the shortcomings in the original MEMOS-3D analysis of experimental observations.

The MEMOS-U model equations
This section introduces the partitioning of domain dimensions, the key aspect of the physics approach behind MEMOS-U. The code features 3D heat transfer and 3D current density pathways, but only 2D fluid motion, the latter realized by employing the shallow water approximation where the Navier-Stokes equations are integrated across the vertical column depth, yielding two horizontal momentum equations and a continuity equation for the liquid column height h. The continuity equation then describes the surface deformation by providing the evolution of the sample surface position. The description of the free surface interface with a 2D vector field makes for extraordinary ease of interpretation and efficient numerical implementation [32,33]. Thus, this approach facilitates MEMOS-U simulations of large domains and over long times within feasible time frames.
The shallow water equations are valid when the horizontal liquid length scales are much larger than the vertical length scale. It is not possible to set a strict maximum allowable ratio between the vertical and horizontal scales, but one can consider the value 0.05 as a soft validity limit for the shallow water approximation [24]. This threshold is typically not surpassed for molten PFCs in tokamaks, where the melt pool commonly extends several centimetres while the melt layer and surface build-up generally remain at around or below a millimetre [17,19,20].
Since the shallow water equations have been rather extensively described in the literature [23,24], their complete derivation will not be repeated here. Instead, the discussion will focus on the applied physical boundary conditions, the neglected terms and the treatment of the transient boundary at the bottom of the liquid pool, i.e. the melting or re-solidification front.

The thermoelectric magnetohydrodynamic equations
The self-consistent description of the liquid metal is based on the incompressible resistive thermoelectric magnetohydrodynamic (TEMHD) equations, together with the temperature convection-diffusion equation [34]. The temperature equation is derived from conservation of energy; it includes the volumetric energy source terms from Joule and Thomson heating [34], but the negligible viscous heating term has been omitted. The fluid description is given by the two incompressible Navier-Stokes equations and includes the Lorentz force density.
The electromagnetic description is limited to the irrotational condition for the electric field and the solenoidal condition for the current density. This is a consequence of the magnetostatic approximation which allows important simplifications of Faraday's and Ohm's law, as well as of the strong tokamak magnetic fields that essentially allow Ampere's and Gauss' magnetic law to be bypassed [35]. Note that the thermoelectric current does not explicitly appear in the current density expressions, since it can be shown to be an irrotational vector for uniform material composition [34].
The emerging fluid, temperature and current density equations are the following: with (v, J, B, p, T) the fluid velocity, current density, magnetic flux density, fluid pressure, temperature and (ρ m , c p , k, µ, ρ e , S) the mass density, specific isobaric heat capacity, thermal conductivity, dynamic viscosity, electrical resistivity, absolute thermoelectric power. Note that the relative dimensions have been exaggerated for illustrative purposes.

The shallow water equations
In what follows, the shallow water equations will be derived from the Navier-Stokes equations supplemented by appropriate boundary conditions corresponding to the free surface as well as the liquid-solid interface. In order to formulate boundary conditions on the free surface, not only the liquid but also the plasma-plus-vapour ambient medium need to be implicitly treated as a (immiscible) TEMHD fluid. Ultimately, the plasma and vapour effects appear exclusively through the pressure and drag force in the stress balance conditions as well as the incident plasma heat flux and the vapour cooling flux in the boundary condition for the temperature equation. Those plasma properties are physical quantities that can be extracted from external (kinetic) models or experimental observations. The notations b 1 and b 2 are employed for the surfacenormal positions of the liquid-solid and the free surface boundaries, respectively. Figure 1 offers an illustration of their relative placement with respect to the undeformed (pristine) surface. The melt layer thickness h is defined as the distance between the two interfaces, h = b 1 − b 2 . Throughout the text, the subscript n is reserved for the inward surface normal (vertical, directed into the domain) components and the subscript t is reserved for the surface tangent (parallel) components of all vector quantities.
On the free interface, the kinematic boundary condition is applied with the surface normal velocity given by the time derivative of the surface deformation in the Lagrangian frame. Accounting for the rate of change of the interface position due to vaporization (ẋ vap ), we obtain In the application of the dynamic boundary condition, the Cauchy stress tensor is considered but the magnetic component of the Maxwell stress tensor (that results from the J × B force density) is discarded due to the negligible induced magnetic field and the non-ferromagnetic nature of both media. Stress balance normal to the surface leads to a pressure jump at the free interface that is equal to the surface tension force: where γ is the surface tension, κ is the free surface curvature and P is total ambient pressure (plasma plus metal vapour). In the shallow water applications of interest, the surfaces remain relatively flat and the surface curvature is generally very small so that surface normal tension can be neglected. Stress balance tangential to the surface leads to the boundary condition: Here the viscous plasma contribution has been converted to the external drag force density f d . Assuming uniform material composition, the surface tension gradient leads to the thermocapillary term (∂γ/∂T)∇ t T.
On the liquid-solid interface, the no-slip condition is applied at the bottom of the melt pool. The relation: forces a vertical velocity gradient which induces a shear force that counters the flow velocity. In the appendix, we discuss a more general kinematic boundary condition. The first step in the derivation of the shallow water equations involves integration of the incompressibility equation, equation (1), across the shallow melt layer between the free surface and liquid-solid interface that introduces the depth-averaged velocity U = (1/h)´b 1 b2 v t dn. Consideration of the kinematic boundary condition at the free surface and the no-slip condition at the liquid-solid interface yields a continuity equation for the surface-normal column height, or melt layer thickness h. The equation features a source or sink term accounting for melting or re-solidification (∂b 1 /∂t) and a sink term accounting for vaporization (ẋ vap ). The height equation reads as: The derivation of the shallow water equations proceeds with inspection of the surface-normal momentum equation. Orderof-magnitude estimates reveal that the inertial terms are commonly negligible, and hydrostatic balance of the liquid column can then be assumed, as if the column were at rest at all times [24]. For the regimes that are encountered in tokamak melting scenarios, the vertical force components are typically negligible, resulting in ∂p/∂n = 0 from normal momentum balance. Accounting for the dynamic boundary condition at the free surface, we acquire: A non-negligible normal force will be counteracted by a surface-normal pressure gradient, and the average pressure in the liquid column will depend on h through hydrostatic balance. The next step of the shallow water equation derivation involves the depth averaging of the tangential momentum equation across the shallow melt layer, the application of the kinematic boundary condition and combination with the column height continuity equation. Neglecting, for simplicity, the variation of the velocity from the average, we get: Explicit momentum transfer across the evolving liquid-solid boundary has been simplified, see the appendix for details. The last viscosity term can be evaluated by prescribing a parabolic profile for the depth-dependence of v t . The three unknown polynomial coefficients are then computed with the (a) tangential stress condition equation (8), (b) no-slip condition equation (9), (c) equality of the depth-averaged polynomial to the depth average velocity field. Denoting the depth average operation with brackets ⟨·⟩ and introducing the surface temperature T s , we have: Thus, the shallow water approximation reduces the unknowns from four (p, v) to three (h, U), i.e. to the melt layer thickness and depth averaged tangential velocity.

The model equations
The depth-averaged melt velocity U should also be used in the advective term of the temperature equation. The present treatment of the convective derivative is based on two simplifying assumptions that are aligned with the shallow water approximation: (a) the surface-tangential velocity is constant along the melt column depth; (b) the surface-normal velocity component is zero. Moreover, the irrotational electric field allows the definition of an auxiliary potential ψ and conversion of the current continuity condition into a differential equation for this scalar field. Overall, the MEMOS-U system of equations reads as: with σ e = ρ −1 e the electrical conductivity.
Since P, f d as well as B = B ext are external quantities and b 1 ,ẋ vap are fully determined by the temperature, there are essentially 5 unknowns (h, U, T, ψ).
The boundary condition for the height equation is formulated by specifying the inflow at the edge of the fluid domain. Since no material can flow into the domain from outside, a trivial condition emerges when the vectors U,n t are antiparallel, where the latter is the outward pointing normal to the 2D fluid domain. On the other hand, there is a free outflow, that is no boundary condition is enforced, when U,n t are parallel. Hence, we have: The boundary condition for the momentum equation is: for the advective term and shear stress due to horizontal gradients. The procedure of removing or adding occupied fluid cells is slightly involved; details of the numerical implementation will be discussed in section 3. The boundary condition for the temperature equation stems from prescribing a heat flux across the plasma-exposed surface. The modification of Fourier's law by thermoelectricity has also been included [34] to obtain: In the above;n is the inward directed unit normal to the 3D domain, q inc is the incident heat flux from the plasma-facing side, q cool is the surface cooling flux. Finally, the boundary condition for the auxiliary potential stems from current continuity and reads as [35]: where J s denotes the total current density on the sample surface. The description and the physics behind all the so far unspecified terms (q inc , q cool , J s ) will be thoroughly discussed in section 4.

The coupling between equations
The shallow water equations for incompressible 3D flows are fully analogous to a compressible flow in 2D, where the column height plays the role of the mass density [24]. Typically, in shallow water applications gravity acts as a restoring force with the corresponding equation of state given by the hydrostatic balance of the liquid column. In the present case, no such explicit relation exists between equations (10) and (11). The evolution of the column height h is determined by the 2D velocity field, but there are no restoring forces in the momentum equation. Only the vertical shear stress, which is the third rhs term of equation (11), is a strong decelerating term depending on h.
The spatiotemporal evolution of the solidification front b 1 is determined by heat transfer. Thus, equation (12) implicitly affects equation (10). Temperature diffusion dictates surface cooling and in this manner affects the momentum source terms which depend on surface temperature, for instance the thermocapillary and Lorentz force densities. Transversely, the velocity field of the molten material determines the convective term of the temperature equation and the continuity equation for h determines filled grid points to be added or removed in the process of surface deformation, determining the heat conduction domain.
Equations (10) and (11) describing liquid metal motion are two-dimensional, while equations (12) and (13) describing heat transfer and the current pathways are three-dimensional. Therefore, there is a respective partition of the numerical domains in the MEMOS-U code; the quantities (T, ψ), or equivalently (T, J), are defined on a 3D Cartesian grid, whereas the quantities (h, U, b 1 , b 2 , p, f d ) are defined on a 2D grid parallel to the exposed side, with correspondingly overlapping grid point positions.

The MEMOS-U numerical implementation
The MEMOS-U code solves the system of equations (10)-(13) on a rectilinear grid with a finite difference approach [36] The sub-steps (a)-(f) are discussed in detail below.

Temperature diffusion
The temperature equation is solved on the 3D domain with an implicit scheme, splitting the treatment into one tridiagonal matrix inversion for each dimension so that the solution requires O(N) operations with N the number of grid points. Figure 2 illustrates a sample geometry at two exposure phases-before and after some surface deformation by fluid motion. The height function b 2 determines the surface index that constitutes the boundary of the heat conduction domain where condition equation (16) is applied by utilizing ghost Here the domain is cross-sectioned along the n, t 1 -plane on a point of symmetry of the t 2 -direction where t = (t 1 , t 2 ). Note that the relative dimensions have been exaggerated and that a low resolution has been selected for illustrative purposes. point temperatures. As demonstrated in the figure, the heat conduction domain approximates the deformed surface by a step function.
The transformation between liquid and solid phases is treated with the heat integration method [37], where an algorithm keeps track of the latent heat of fusion for each individual grid point, and the temperature is kept at the melting point until phase transformation is completed. This casts the melt interface problem in a weak formulation, where each grid point can be considered as fully molten, fully solid, or undergoing a phase transition. The interface position within the grid point is unknown, thus the melt layer thickness accuracy is directly related to the grid size.

Temperature advection
Fluid motion displaces the sample's temperature profile in the fixed grid reference frame. The treatment of advective temperature derivatives is standard in computational fluid dynamics and is implemented in MEMOS-U via the finite difference discretization of the advective equation: by employing a first order upwind scheme. The typical melt layer shape is indicated in red and the characteristic replacement current pathways are depicted in a simplified manner, see [35] for numerical solutions. The alignment of the local Cartesian coordinate system is also displayed (different for each exposure geometry). In both scenarios, the melt layer displacement is directed out of the page, i.e. perpendicular ton,t.
Owing to the dimensionality mismatch, the melt velocity is applied to the whole molten material column at the corresponding surface point. As a consequence of the step function approximation of the melt layer thickness and the deformed surface, some fluid elements will be index-adjacent to solid or non-material elements and the advective derivative is set to zero for the corresponding direction. This can be discerned in the illustrative melt pool example of figure 2, where molten elements are adjacent to empty or solid grid points due to the step-wise approximation. In a highly resolved simulation, there will be few such grid points compared to the fully surrounded ones, thus introducing negligible errors.

Current density distribution
The finite difference approximation of equation (13) leads to the linear system of equations AΨ = B where Ψ is the auxiliary potential array, B the boundary value array and A the matrix operator for the spatial derivatives [35]. The system is solved with the iterative generalized minimum residual method (GMRES) [38]. Figure 3 features 2D cross-section sketches for two sample geometries which have been used on tokamaks in recent years in order to investigate melt dynamics: a sloped lamella (a) and a leading edge (b). The typical shape of the melt layer is included and simplified, but realistic current paths are drawn, see [35]. for the exact numerical solutions. The sketches illustrate the important differences in the near-surface current pathways between the leading edge and sloped samples as well as in the current density component that mainly induces melt motion through the J × B force. In the sloped lamella, the normal current density component is sought, implying a negligible difference between the current density within the melt layer and the surface current density J s . On the other hand, in the leading edge, the tangential current density component is sought which can be significantly smaller than J s depending on the bending of the current paths [35].

Column height continuity
The column height h is determined by equation (10) and is propagated with the aid of an explicit first order upwind scheme. The liquid-solid interface b 1 is determined solely by the temperature equation and incremented in the heat diffusion sub-step, described above. Thus, surface deformation b 2 = b 1 − h is also determined in this sub-step.

Momentum equation
The velocity field is updated by employing a first order upwind scheme for the advective term and an explicit representation of the momentum source terms. When a liquid metal element resolidifies, the corresponding melt velocity is set to zero. On the other hand, when a solid element melts, a distinction should be made. In case the surface element melted under the effect of the surface heat flux only, then the melt velocity begins from zero. In case the surface element melted due to material outflow from an adjacent boundary, then the melt velocity of the corresponding upstream element is assigned.

Re-evaluation of the simulation domain
The re-evaluation of the simulation domain constitutes the last sub-step of each time step propagation. For each grid point on the surface-parallel plane, the surface position b 2 is used to detect the corresponding surface-normal grid point index that constitutes the domain boundary.
The updated surface indices delimit the heat conduction domain during the next computational cycle (time step), when the heat flux boundary condition will be enforced. This procedure results in the removal or the addition of grid points for each surface-parallel index and the temperature treatment at these points needs to be considered. When elements are being discarded, no action is taken and their removal is assumed not to impact the adjacent filled grid points. When elements are added, a distinction should be made similar to that in the momentum equation. In the case in which a molten grid point is added on top of a molten column, the temperature of the previous surface element at the same position is assigned. On the other hand, in the case in which a molten grid point is added on top of a previously solid element, the surface temperature of the upstream element is assigned.

The surface current density
The only non-trivial boundary condition in the replacement current module of MEMOS-U is the surface current density J s , whose nature depends on the plasma scenario and the sample material composition. In case of steady-state driven, ELM driven and disruption thermal quench driven melting in tokamaks, the surface current is given by the local escaping thermionic current that is dictated by the temperature field and the incident plasma fluxes. In contrast, in case of melting due to disruption current quench heat fluxes, the surface current is given by the halo current which strongly depends on the plasma conditions and in general requires global modelling.

Thermionic surface current.
Within the assumption of plasma ambipolarity, the release of thermionic electrons from the surface of a grounded bulk metallic sample generates a replacement current which flows through its volume [4,35]. The associated current density is computed from the boundary value problem of equations (13) and (17) and is responsible for the dominant J × B force density acting on the molten metal [35]. In the absence of external electromagnetic fields, the nominal thermionic current is described by the Richardson-Dushman formula that reflects the classical over-the-barrier physical mechanism [39]: with W f the room temperature work function and A eff the effective, material-dependent Richardson constant. It is evident that thermionic emission drastically increases for elevated surface temperatures and/or low work functions. As a result, considering the main PFC materials on ITER, it is very important for W melting (T m = 3695 K, W f = 4.55 eV), but it is totally negligible when considering Be melting (T m = 1560 K, W f = 4.98 eV). The unimpeded thermionic current from W is known to greatly exceed the incident plasma current at the strike point position of typical Type-I ELMing H-mode AUG and JET discharges already at the melting point [40,41]. However, complications arise due to the self-consistent sheath electrostatic fields and strong external magnetic fields. Such high emission intensities are incompatible with the classical Bohm presheath structure and eventually space charge accumulation leads to the formation of a virtual cathode limiting the escaping current density to a maximum value [42,43]. This is referred to as a space charge limited sheath. In the case of oblique magnetic field lines, a fraction of the low energy thermo-electrons is recaptured by the surface during the first Larmor gyration [44,45]. This is referred to as prompt re-deposition.
Both suppression mechanisms need to be properly considered, since thermionic emission has a rather complex impact on the melt dynamics, directly influencing both liquid motion and heat transfer. The emitted current density determines the volumetric J × B force exerted on the melt layer and thus the displaced liquid volume as well as the displacement itself. Furthermore, thermionic emission constitutes an important cooling mechanism (see section 4.3) which affects heat balance and thus the quantity of liquid metal volume produced. The intricate coupling between the momentum and energy transfer aspects of thermionic emission is so important that any miscalculation of the escaping current density would result in large errors in the melt displacement estimates. The central role of thermionic emission has been recognized in MEMOS-U and led to its rigorous treatment via dedicated particle-in-cell (PIC) simulations [40,41,46]. 2D3V PIC simulations of magnetized plasma sheaths with intense thermionic emission revealed that the escaping thermionic current saturation is a global characteristic of the spacecharge limited regime, irrespective of the inclination angle of the magnetic field [46]. Extensive simulations led to the identification of an accurate and physically transparent expression that describes the dependence of the limited thermionic current on the plasma conditions, magnetic field strength and inclination angle. The semi-empirical expression was constructed on the basis of the Takamura model of unmagnetized space charge limited sheaths [43], the Chodura picture of inclined magnetized sheaths [47] and a single particle field-free approximation of prompt re-deposition. It reads as [46]: with n e the plasma density, v Te the electron thermal velocity, α the magnetic field inclination angle with respect to the surface normal. The above expression is valid for the inter-as well as intra-ELM plasmas of contemporary devices. It is accurate within 25% for non-grazing incidence (α ≳ 5 • ) and within 60% for α = 2.5 • [46]. Note though that sin 2 (2.5 • ) ≃ 0.002, which suggests a rather negligible value of J lim th . Overall, due to the strict current limitation in the space charge limited regime, the thermionic current is given by: The following remarks are pertinent to the MEMOS-U implementation. (a) The escaping thermionic current expression stems from PIC simulations which assume an infinite emitting wall with a homogeneous prescribed temperature [46]. Since equation (20) is applied to all surface elements, it is implicitly assumed that the escaping current density contributions from each surface element are independent. PIC simulations of hot spots should be carried out to confirm the validity of this convenient assumption. (b) The low temperature branch of equation (20), which is selected prior to current limitation, follows the unimpeded Richardson-Dushman expression. This means that prompt re-deposition is only considered after the virtual cathode has formed. This simplification does not lead to errors since, for oblique inclination angles, the transition to the space-charge limited regime takes place at temperatures well below the W melting point. (c) As a result of surface deformation by melt motion, the surface normal is spatiotemporally varying which modifies the escaping current density [40,41,46]. The overall sin 2 (α) angular dependence implies a rapid escaping current decrease for shallow inclinations. Thus, depending on the geometry, it may be important that the inclination angle is not assigned its static value defined for a pristine unexposed sample but its local dynamic value [29].
(d) The Child-Langmuir type expression for the escaping current density depends on the plasma density and electron temperature profiles. Since such profiles are generally not available during ELMs, this expression can be transformed into a dependence on the incident field line parallel heat flux by assuming that the heat flux spatiotemporal variations can be attributed to temperature variations, that the plasma is isothermal, and by adopting a standard value for the sheath heat transmission coefficient [48]. We emphasize that the present description of the escaping emitted current density applies to both inter-and intra-ELM periods in contemporary fusion devices as well as to inter-ELM periods in ITER. A number of complications emerge during the ITER intra-ELM periods. (a) Owing to the combination of high plasma densities with elevated electron temperatures (T e ≳ 200 eV), the contributions of secondary electron emission [49] and electron backscattering [50] start to become important along with thermionic emission. (b) Due to the elevated plasma temperatures and densities, the wall electrostatic field becomes high enough (E w > 10 7 V m −1 ) so that the Schottky effect becomes important leading to thermionic emission enhancement [39]. (c) Due to the elevated electron thermal current, a virtual cathode is formed at unrealistically high-even for melting scenarios-surface temperatures which implies that the emitted currents are only suppressed by prompt re-deposition. PIC simulations of the magnetic presheath and electrostatic sheath in the presence of all emitted electron populations should be carried out in order to acquire reliable expressions for the escaping emitted current during ITER intra-ELM periods.

Halo surface current.
In case of PFC melting that is induced during disruption current quench phases, the halo current density determines the volumetric Lorentz force. The charge continuity equation inside the metal implies that the current density gradient scales should have the same magnitude in each spatial direction in the proximity of the surface. Thus, depth-variations of the current density inside the molten material [35] can be neglected provided that the horizontal gradient scales are very large compared to the melt thickness. Then, J = J s and the current module can be circumvented with J s × B providing the Lorentz force density. Given the overlap between the plasma-wetted area and the halo region, the horizontal extent typically reaches several centimetres during disruptions [6,20], while the melt layer thickness remains below a millimetre. Hence, there is no need to accurately model the complicated current paths through the vessel that is a topic of considerable complexity.
In contrast to thermionic emission where J s = J th is calculated self-consistently by the surface temperature, incident heat flux and local inclination angle, J s = J halo constitutes an external input for the MEMOS-U code. In disruption current quench-driven Be melting in JET [28], 100 kA m −2 was employed based on empirical input [21]. In predictive Be melt simulations of a 5 MA VDE ITER scenario, 20 kA m −2 was used based on the total halo current and the wetted area [6].

The surface heat flux
In ELM-driven and disruption-driven melting events, the incident energies of the plasma particles are always less than few keV. The corresponding depth range, i.e. the total material thickness traversed by the impinging particles prior to thermalization, lies well in the nanometre range and is thus much smaller than the characteristic length scales of interest and the typical MEMOS-U mesh sizes. To be more specific; for E inc = 10 keV the electron, proton and helium depth ranges are ∼300, ∼150, ∼100 nm in Be and ∼80, ∼40, ∼20 nm in W [50][51][52]. As a result, plasma incidence does not result in volumetric heating, but in surface heating which can be accounted for only by the boundary condition equation (16).
For completeness, we note that runaway electron (RE) driven melting events are subject to volumetric heating, since when E inc = 1 MeV the electron depth ranges are ∼3 mm in Be and ∼0.4 mm in W [53]. For the relativistic initial electron energies, electronic stopping as well as nuclear stopping are important, but also internal photon and secondary electron transport are central in the determination of the energy deposition profile [54]. Heavy RE-driven melting events cannot be simulated with MEMOS-U not because of volumetric plasma heating, but because the shallow water approximation breaks down.
The incident plasma heat flux is an external input to the MEMOS-U code and is inferred either by experimental measurements or theoretical considerations. In the ELM-driven W melt experiments on JET and AUG using dedicated leading edge or sloped samples [17,19], spatiotemporal heat flux profiles are deduced from infrared diagnostic measurements of the temperature response of adjacent flat W samples. Employing the optical approximation and in accordance with the local magnetic field line inclination, the extracted heat flux is mapped to the exposed surface of interest. In some cases, 3D diffusion effects and isotropic background contributions need to be accounted for [55]. The inaccuracies of such heat flux profiles are expected to be within ∼25% [19,28,29]. Since the reference sample temperature is low, the total q inc is extracted.
In disruption driven (current quench) Be melting in JET, direct time resolved measurements are unavailable. The thermocouple measurements in the upper dump plates yield the total energy deposited onto the Be plates during the whole discharge. The heat flux is evaluated by estimating the elapsed deposition time and the wetted area [20,56]. Given the thermocouple measurements, in this case only the conducted heat flux, q inc − q cool , is extracted. Finally, in predictive disruption current quench-driven Be melt simulations for ITER, the adopted heat flux is based on theoretical estimations and simulations of the conducted power across the separatrix and of the deposition profile on the first wall [6].
The sensitivity of the MEMOS-U output to the heat flux input depends on the melting scenario. In repetitive-ELM melting of W samples, most of the incident energy is expended in raising the base sample temperature close enough to the melting point, so that superimposed ELM temperature excursions will suffice to initiate transient melting near the end of the exposure. This is a threshold phenomenon, where a minuscule (<1%) quantity of input energy is consumed for the phase transformation. Thus, small-well within experimental uncertainties-heat flux variations will significantly affect the amount of melt generated. A quantitative discussion will follow in section 7.4. On the other hand, in disruption current quench-induced melting of Be samples, a larger fraction of the input energy is expended for melting. Therefore, the amount of metallic melt produced is less sensitive to heat flux variations.

The surface cooling flux
Surface cooling consists of electron emission cooling, vapour cooling and thermal radiation cooling [57,58]: Their implementation in MEMOS-U is detailed below.

Thermionic cooling.
For any uncompensated emitted electron, not only its kinetic energy but also its binding energy (given by the work function) is released from the PFC surface. Local thermodynamic equilibrium implies that thermionic electrons follow a Maxwellian distribution with a temperature equal to the local surface temperature [39]. Hence, the cooling flux is given by [39,59]: where e denotes the elementary charge. Since W f ≫ k B T s holds even beyond the normal W boiling point, the binding energy contribution is dominant. This cooling flux plays an important role in bulk W heat balance where it can reach values in excess of 50 MW m −2 , significantly limiting the surface temperatures attained by an exposed W surface.

Vapour cooling.
Within the assumption of thermodynamic equilibrium between the exposed surface and the emitted vapour, the evaporated mass flux is described by the Knudsen formula [60]: where m a is the atomic mass and P vap is the vapour pressure whose dependence on the surface temperature is described by the exponential Antoine correlation. Evaporation plays a role in both surface deformation and cooling. As mentioned earlier, surface mass loss is translated into an interface velocityẋ vap which acts as an independent sink term in equation (10). The respective cooling flux has a similar form to the thermionic cooling flux, but with the binding energy equal to the vapourization enthalpy h vap , i.e.
Since h vap ≫ k B T s holds even beyond the normal boiling point of the materials of most current interest in tokamaks and for ITER, the binding energy contribution is dominant. For completeness, we note that h vap exhibits a very weak temperature dependence with ∼3% and 5% differences between its room temperature and normal boiling point values for W and Be, respectively [61,62].
At this point, we discuss the relevant issue of vapour shielding of the incident plasma heat flux. For high surface temperatures, saturated vapour pressure can exceed the plasma pressure (∼1 kPa for tokamak conditions). In this case, the vapour species is expected to strongly interact with the incident plasma particles, leading to a reduction of the absorbed plasma heat flux. Atomic processes such as ionization, recombination as well as scattering govern the local plasma-vapour energy transfer, while the vapour transport characteristics are also influenced by the plasma and surface parameters. In spite of the fact that the ionized vapour-plasma interplay is highly complicated, it has been clearly demonstrated that the choice of vapour shielding model has a minor effect on the absorbed heat flux [63,64]. This is a direct consequence of the exponential temperature dependence of P vap , which effectively masks any variable interaction strength. One must therefore distinguish between the anticipated accuracy in the evaluation of the absorbed heat flux and of the evaporated material loss. For the latter, inaccurate evaporation rates are expected in scenarios involving very high surface temperatures. On the other hand, the absorbed heat flux is expected to be relatively insensitive to the vapour shielding description, with the beneficial consequence of robust melt production predictions.

Thermal radiation cooling.
Surface cooling due to the emission of thermal radiation is described by the Stefan-Boltzmann law [65]: where σ SB is the Stefan-Boltzmann constant and ϵ T is the total hemispherical emissivity of the PFC which depends strongly on the surface temperature. This cooling term is significant at surface temperatures below melting and for longer time scales, being overpowered at higher temperatures by the exponential nature of thermionic emission and/or evaporation. This contribution is rather negligible for typical tokamak transient melting experiments. The relevant significance of these surface cooling terms can be better understood by plotting their magnitude as a function of the surface temperature for W and Be, see figure 4. In the case of W, thermionic cooling dominates in the rather wide neighbourhood of the liquid-solid phase transition. The emergence of current limitation close to the melting point, here calculated for a typical heat flux of q inc = 1 GW m −2 , switches the monotonic character of the q te (T s ) increase from exponential (courtesy of the Richardson-Dushman formula) to weak linear (courtesy of the kinetic cooling contribution). For W, vapour cooling remains negligible even beyond melting, since W possesses one of the lowest known vapour pressures, but due to its continuous exponential dependence, ultimately surpasses thermionic cooling for higher surface temperatures (not shown here). As mentioned above, radiative cooling provides an important contribution only prior to melting. In the case of Be, thermionic and radiative cooling are always negligible due to the low melting point, while vapour cooling becomes significant after melting.

The surface acceleration
Here we discuss the physics and importance of the secondary momentum terms, i.e. all terms other than the main volumetric Lorentz force and viscous deceleration.

Surface pressure gradient.
As discussed in section 2.2, the surface pressure is equal to the ambient total pressure P which has two contributions, namely the plasma P pl and vapour pressure P vap . The plasma pressure is an external input to MEMOS-U, whereas the vapour pressure depends on the surface temperature through the Antoine equation. For the typical gradients of <1 kPa cm −1 , ∇P pl is small compared to the J × B force density that arises in contemporary tokamak exposures and also remains negligible in the ITER scenarios considered in dedicated studies [6]. Hence, it is evident that P vap can only be important once it becomes comparable to P pl . This is realized at very high surface temperatures, well beyond melting, where the accuracy of calculations strongly depends on the choice of the vapour shielding model, see section 4.3.

Thermo-capillary convection or the Marangoni effect.
Thermodynamic arguments predict a linear decrease of the surface tension beyond the melting point ∂γ/∂T = −α 2 , as confirmed by experimental results for most liquid metals [66,67]. The resulting thermo-capillary flows are introduced through the tangential stress balance boundary condition equation (8) on the free surface and, after depth averaging, yield the fifth momentum source term on the rhs of equation (11). By simple inspection, it is evident that intense surface temperature gradients lead to an effective force on the melt layer which pulls the liquid metal from the hot zones toward the colder areas. This is the so-called Marangoni term that is calculated selfconsistently through the known temperature field. Given the spatial extent of the metal melt in the case of tokamak experiments, strong tangential gradients are rather atypical for leading edge ELM-induced melting or disruption current quenchinduced melting [6,28], thus thermo-capillary accelerations are generally not important compared to the main J × B acceleration mechanism. Nevertheless, the effect can be important in some exposure geometries, such as in the sloped W lamella experiments in AUG [29].

Plasma surface drag.
Plasma flows oriented along the surface lead to a drag force per unit area f d . After depth averaging, the respective force density is given by the sixth momentum source term in equation (11), formally also introduced through the tangential stress balance boundary condition equation (8). In contrast to the Marangoni term, the acceleration explicitly depends on the incident plasma parameters and is an external input to MEMOS-U.

The MEMOS-U material properties
The typical melt velocities found in contemporary tokamak experiments are of the order of a few m s −1 , that is orders of magnitude lower than the ∼ km s −1 thermodynamic sound speed in liquid metals at their melting point [68]. These extremely low Mach numbers allow us to circumvent an equation of state and safely assume the validity of the incompressibility condition. Moreover, given the kPa plasma pressures; the ∼100 GPa inverse compressibilities and the ∼0.1 K MPa −1 initial slopes of the melting curves of the liquid metals of interest [68][69][70][71] suggest that the thermophysical and electrochemical material properties are independent of the pressure and depend solely on the temperature.
MEMOS-U employs analytical expressions for the temperature dependence of (ρ m , c p , k, ρ e , S, µ, γ, P vap , ϵ T ) valid from room temperature up to the atmospheric pressure boiling point, and also state-of-the-art values for the quantities (T m , W f , A eff , h fus , h vap ). All expressions are based on reliable experimental data rather than theoretical calculations and great effort has been taken to ensure that extrapolations are minimized.
Most of the adopted W material properties follow the recommendations of a recent comprehensive literature review [72]. A similar-yet unpublished-investigation has been carried out for the Be material properties by one of the authors [73]. In contrast to W, the properties of polycrystalline Be are much less explored, particularly for the liquid phase. As a consequence, few data sources are available for most properties beyond the melting point. Hence, one should remain critical of the adopted values, but keep their accuracy in reasonable relation to the overall accuracy and the sensitivity of any specific modelling problem. The discussion below exemplifies this remark.
Under atmospheric pressure, pure Be undergoes a polymorphic phase transition from an hcp (alpha-phase) to a bcc solid (beta phase), that takes place at 1543 ± 5 K and is naturally accompanied by discontinuities in most material properties [74]. Since the melting phase transition from a bcc solid to a liquid takes place at 1560 ± 5 K, there is a narrow temperature window of 17 ± 10 K for the beta phase of Be. Such desired accuracy would lead to unnecessary computational cost given the uncertainty in the plasma-facing boundary conditions. Therefore, in MEMOS-U, the alpha-beta transition is ignored and the enthalpy of hcp-to-bcc transformation (6.855 kJ mol −1 ) is combined with the enthalpy of fusion (7.959 kJ mol −1 ) in a single latent heat contribution that is absorbed or released only at the melting point.
In addition, only pure polycrystalline W and Be are considered. In fact, experience with JET tokamak operation in a Be/W environment has revealed that the thickness of deposited Be layers and the extent of Be-W intermetallic compound formation in the divertor as well as the depth of beryllium oxide formation in the chamber are of the order of few micron [20,[75][76][77] and thus much thinner than the typical melt layers. There are no significant benefits in attempting to incorporate the mostly unknown BeO or Be x W y material properties.
Finally, it is also worth pointing out that, in view of planned transient and steady state melt experiments with proxy materials in AUG and DIII-D, the MEMOS-U material library has been recently expanded to include niobium, iridium, and aluminium.

The small melt displacement and deformation simplification
The initial computational domain should match the 3D geometry of the pristine sample. As the computation/exposure progresses, the sample is heated and eventually a melt pool is produced. The 2D outline of the melt delimits the fluid domain and the exposed surface of the 3D domain is enclosed by the dynamic scalar field b 2 . As surface deformation progresses driven by the Lorentz force, the 3D heat conduction domain should be re-evaluated for consistency. When molten material is displaced along the surface and excavates or amasses in different parts of the sample, both the sensible and the latent heat are displaced with the fluid. The numerical difficulties arising in the treatment of the transient 3D domain along with the energy displacement constitute the main complications in the implementation of the full model, equations (10)- (13).
These complexities can be circumvented on the basis of two assumptions for the molten surface characteristics: (a) The surface deformation is small compared to the melt depth. Within this approximation, the column height continuity equation (10) can be rewritten as: by substituting for h = b 1 − b 2 and employing h ≃ b 1 in the convective term, where b 1 is the distance between the melt interface and the initial pristine surface, as illustrated in figure 1.
The accumulated surface deformation is small so that the temperature response of the deformed surface is essentially the same as that of the pristine surface.
These assumptions drastically simplify the coupling between the fluid equations and the temperature equation, see equations (10)- (12). The fluid solver does not provide any input to the heat solver and the continuity equation for the surface deformation becomes linear. The surface deformation is stored, but it does not constitute input for the progressive heat transfer calculations, or for any other parts of the code. In this manner, surface deformation is not considered in the heat transfer domain, but the deformed surface shape can still be predicted after several accumulated events (for example in repetitive ELM scenarios) or continuous displacements (for example in longer timescale disruption current quench scenarios). Moreover, since the wetted PFC area is typically much smaller than their total area, surface melting is localized so that regions of freshly produced melt are surrounded by pristine solid or by re-solidified regions. The above simplifications lead to liquid pools which thin out toward the edges being a direct reflection of the incoming heat flux spatial profile. Liquid motion is therefore essentially arrested past the pool edges. Section 7 provides a detailed discussion of this particular aspect.
The aforementioned assumptions further reduce the computational complexity and have been implemented in the original MEMOS code employed in [17,26,27]. However, as already pointed out in the Introduction, recent cross-machine experiments provided strong evidence that contradict their validity. In particular, post-mortem analysis after the repetitive ELM-driven melting of misaligned W samples in the JET and AUG divertors [17,19], as well as the current quenchdriven melting of the Be upper dump plates in JET [20], clearly revealed that the melt displacements are large-of the order of the pool extent. The metallic melt is accelerated out of the pool and displaced on top of the adjacent solid, where the strong temperature gradient facilitates heat conduction and rapid solidification near the contact. This regime necessitates a solution of the full model described in section 2. This fact, together with the unjustified circumvention of the field equations and the lack of rigorous boundary conditions in MEMOS-3D, provided the motivation for the development and implementation of the new physics model in MEMOS-U.

Application to tokamak melting experiments
In this section, aiming to demonstrate the capabilities of the MEMOS-U code and the predictive power of the new physics model, we will present simulation results for the JET W leading edge exposures [17], the first attempt in a tokamak to provide a controlled, transient-induced melt for the purposes of model validation. In the experiment, a specially engineered lamella located in the outer, full W divertor, was deliberately misaligned, exposing a tapering leading edge of width 0.25-2.5 mm, intercepting the full heat flux flowing parallel to the magnetic field lines-see a schematic of the geometry in figure 5. The outer separatix strike point was swept across the sample for a few seconds at a time, during nine consecutive type I ELMing H-mode pulses. High resolution visible camera images taken after each pulse revealed that the majority of the surface deformation was inflicted during pulse #84779 [17]. Hence, the MEMOS-U simulations reported here will focus on this single discharge. For a complete description of the experiment and its results, see [17].
An experimental spatiotemporal heat flux is available that was reconstructed from IR thermography data [55]. The typical ELM duration was <3 ms, the average ELM frequency 30 Hz and the toroidal field 2.9 T. The maximum ELM heat flux varied within 0.5-1.5 GW m −2 and the inter-ELM heat flux varied within 50-200 MW m −2 . The experimentally extracted heat flux results in superficial melting that is incompatible with the observed surface deformation. To produce melting in accordance with the observations, the heat flux has been uniformly scaled up by a constant factor 1.10, an increase that lies well within the estimated 25% experimental uncertainty.
The replacement current module is coupled to the heat solver through the escaping thermionic current and the temperature-dependent electrical conductivity, but also coupled to the fluid solver owing to surface deformation. In the leading edge case, as a result of the exposure geometry and the dominant toroidal magnetic field, the main tangential ⟨(J × B) t ⟩ force component lies along the poloidal direction and requires the determination of the radial current density [35]. Dedicated MEMOS-U simulations of current flow in an un-deformed leading edge domain have revealed that the average radial current component inside the melt layer is ⟨J r ⟩ t ≃ J th /4. This approximation is invoked in the present simulations instead of solving the aforementioned boundary value problem in the deformed domain for each time step. This circumvention of the replacement current module results in a drastic decrease of the computational cost and is deemed appropriate on the basis that uncertainties in the experimental heat flux should outweigh the marginal accuracy gained by a fully self-consistent bulk current.
In what follows, in order to demonstrate the inapplicability of the small melt displacement and deformation simplification to leading edge exposures, we shall compare simulation results for heat transfer independent of melt motion as employed in the original MEMOS-3D code analysis of the JET experiment in [17]. and for fully coupled heat transfer-melt motion as employed in the MEMOS-U code. Hereafter these are labelled as 'the simplified model' (see section 6) and 'the full model' (see section 2), respectively. Both models will be realized with the MEMOS-U code in order to quantify this specific model aspect without interference from (a) all other differences between the MEMOS-3D and MEMOS-U models and their numerical implementation, (b) differences in the employed heat fluxes (the reconstruction of the heat flux as in [55]. employed here was not available at the time of [17]).

Response to a single ELM
Let us first consider the effect a single JET ELM whose energy is deposited onto a W sample with an assumed uniform elevated temperature. Here a 2D sample geometry suffices for the illustration. The chosen initial temperature of 3200 K is high enough so that absorption of the ELM energy will increase the sample temperature beyond the W melting point. The temperature response and the surface deformation predictions of the simplified and full model are illustrated in figure 6. A 2D sample-cross section is considered at three time instants. In both cases, the surface is mainly cooled by thermal conduction, with respectable contributions from thermionic and radiative cooling.
The top row in figure 6 features the temperature response and surface deformation resulting from the simplified model. As a result of the non-symmetric spatial profile of the ELM heat deposition, an irregular melt pool is produced that features irregular isotherms. The metallic melt reaches a depth of ∼150 µm and a horizontal extent of ∼4 cm. The ELM lasts ∼3 ms, and the induced melt pool resolidifies within 24 ms. Since the heat conduction domain is not updated to reflect the lateral motion of the molten material, the latent heat of fusion is always deposited back in the area where melting was originally induced. Consequently, after complete resolidification, the hottest point resides between 3.5 and 4.0 cm. As mentioned above, the spatiotemporal dependence of the solid-liquid interface b 1 , obtained by the heat solver, is forwarded as input to the fluid solver which then predicts the surface deformation. The induced melt layer thickness tapers off toward the pool edges owing to the spatial profile of the heat flux pulse. This is significant because: (a) the thickest central molten region is subject to a stronger acceleration for a longer duration compared to the pool edge regions, (b) the depth-dependent viscous shear force (∝1/h 2 ) is relatively weak in the pool centre and gradually increases in magnitude toward the pool edge. These suggest that the depth-averaged fluid velocity develops a smooth gradient from its maximum value at the pool centre to a zero value at the pool edges. Since the heat conduction domain is undisturbed by melt motion, this gradient is present at all times. As a result, the linearized column height continuity described by equation (26) will lead to a surface deformation that builds up smoothly across the whole downwind region of the melt pool. The resulting surface deformation is illustrated by the dashed line that is superimposed on the temperature profile.
The bottom row shows the temperature response and the surface deformation resulting from the full model. Up to 2 ms, the temperature responses of the simplified and the full model are nearly identical owing to the validity of the small melt displacement approximation, since the short elapsed time does not allow for appreciable melt motion. However, at longer times, the predictions of the two models for the temperature profile and surface deformation begin to deviate strongly. Molten material that is produced in the central pool region eventually reaches the edge of the initially induced pool and spills over. The latent and sensible heat are transported with the melt, elevating the temperature of the colder surface underneath which provides an additional conduction channel for the deposited energy. Consequently, after complete resolidification, the hottest point resides around ∼1 cm.

Response to multiple consecutive ELMs
We shall now investigate the combined effect of multiple consecutive ELM pulses by simulating the full exposure during JET pulse #84779. The evolution of the induced melt pools will be easier understood on the basis of the single melting event discussion. The melt layer thickness at the top of the leading edge is illustrated in figure 7 as a function of poloidal coordinate and exposure time.
As for figure 6, the top row of plots illustrates the predictions of the simplified model. The ELMs repeatedly strike near the centre of the sample, where they locally induce melting. The omission of poloidal motion in the heat solver implies that there is no convective heat transport (latent and sensible). Thus, the corresponding region is efficiently heated. At the beginning of the exposure, the pools completely resolidify in-between ELMs. As time progresses and the base temperature increases, melt layers begin to overlap. Melting becomes sustained and the surface temperature gradually rises well above the melting point. The surface deformation as described by equation (26) leads to a significant overestimation and the material pile-up (constrained to be uniquely within the melt pool extent) resides in a place along the lamella (∼2 cm) that is inconsistent with post-mortem observations of the lamella melt profile (see figure 9).
The bottom row gives the predictions of the full model. Due to the inclusion of deformed domains in the heat solver, the absorbed energy is transported together with the molten material toward the lower poloidal coordinates. As a consequence, the centre of the sample is less heated and the local induced melt pools re-solidify faster. Moreover, the transported melt gradually heats the lower coordinate regions, where no melting is originally induced (<2 cm), allowing the subsequent melt layers to move further along the sample prior to re-solidification.
The direct comparison of the outcome of the two models with experimental evidence demonstrates how the violated assumptions of the simplified model distort the temperature response and surface deformation profile. In fact, in order to avoid disagreement with observations, early MEMOS-3D simulations [17] needed to invoke vapour shielding despite the modest heat and energy loads, as well as the lack of direct experimental evidence in favour of significant W vapour generation. The MEMOS-U model, within the experimental heat flux uncertainties, is capable of reproducing all main experimental features without introducing ad hoc factors or vapour shielding effects.

Build-up and final surface deformation profile
The final surface deformation profile is naturally determined by the accumulation of the material displaced from the individual transient melt pools depicted in figure 7. In what follows, we shall limit the discussion to the predictions of the full MEMOS-U model. It is worth pointing out that melting is induced only by the last few tens of ELMs in the discharge, since the absorbed energy is exclusively expended in raising the sample temperature during a large fraction of the ∼1 s exposure.
The gradual build-up of the surface deformation profile is illustrated in figure 8 at the onset of transient melting. The surface profile has been plotted at five instants that are evenly spaced over ten incident ELMs. The region of material pile-up is clearly observed to progressively move further toward the lower coordinates. In the course of the exposure, the deformation builds up after repeated cycles of melting, poloidal acceleration and displacement, surface excavation, material pile-up and resolidification.
The final MEMOS-U surface deformation profile is displayed in figure 9(a), together with a post-exposure image of the JET-exposed lamella from a similar angle. The results correspond to the constant heat flux pre-factor 1.10 that well reproduces the experimental surface profile of [17]. The specific deformation metrics of interest concern; the crater position (1.7 cm from the edge), its depth (0.6 mm) and extension (2 cm), but also the melt pile-up location (0.7 cm from the edge), its height (up to 1.8 mm) and extension (1 cm).

Sensitivity to the incident heat flux uncertainties
The MEMOS-U input that is required for the simulation of the JET ELM-driven W melting exposures consists of the sample geometry, local magnetic field and incident heat flux. Since only the latter quantity is characterized by non-negligible uncertainties, we shall investigate the sensitivity of the final deformation profile to the heat flux input. In principle, the plasma surface drag force per unit area and plasma pressure also constitute input, but their effect here is insignificant, see section 4.4. In order to quantify the sensitivity to heat flux uncertainties, the final surface deformation profile is illustrated in figure 9(b) for a heat flux pre-factor of 1.15. Compared to the reference case of a 1.10 pre-factor, the displaced volume doubles from ∼10 to ∼20 mm 3 . This seemingly disproportional melting enhancement can be understood by considering the partition between the latent and sensible heat during the exposure. The total energy incident on the W sample is ∼11 kJ and the total energy released by the W sample in the form of thermionic electrons, W atoms and thermal photons is ∼1 kJ. Hence, the total energy absorbed by the sample during the exposure is ∼10 kJ. In contrast, the energy required to melt the observed ∼6 mm 3 of W is ∼30 J, calculated assuming a constant mass density of 19.25 g cm −3 and using 280 kJ kg −1 for the enthalpy of fusion [72]. Therefore, as already mentioned in section 4.2, approximately 99% of the absorbed energy is expended to raise the base sample temperature, while only a very small fraction corresponds to latent heat, which explains the emerging threshold-like behaviour.
The strong sensitivity of the displaced volume to the heat flux input casts a doubt on the value of such benchmarking activities. However, careful comparison between figures 9(a) and (b) reveals that the surface morphologies are nearly identical and the deformation profiles are approximately scaled versions of each other. Furthermore, the key deformation characteristic (molten material strongly displaced from the transient pool locations which appear as a crater upon resolidification) is well captured by the full MEMOS-U model and cannot be predicted by the simplified model featuring stagnated melt pools. In general, it is crucial to recognize that benchmarking activities should not be limited to comparisons of the highly sensitive displaced volume, but should be focused on deformation profile comparisons. It is worth pointing out that in the similar leading edge exposures performed on AUG, molten material was actually ejected from the sample edge, a salient feature that was successfully reproduced only by the MEMOS-U full model, see [28]. for details. Finally, it is important to reiterate that such strong sensitivity is a peculiarity of leading edge exposures to type-I ELMs and that other exposure scenarios exhibit much more moderate sensitivity, see the discussion in section 4.2.

Surface cooling and thermionic emission
Let us expand the discussion on heat flux sensitivity by including the effect of cooling fluxes. In figure 10, the integrated incident plasma energy and the integrated cooling energy are plotted as functions of the exposure time. Clearly, the total cooling energy released from the sample constitutes only ∼5% of the overall energy budget. As discussed in section 4.3 and observed in figure 4, the major cooling contribution stems from thermionic cooling (the maximum surface temperature  reaches 4500 K in this simulation). Therefore, a high sensitivity to the heat flux input directly translates to a high sensitivity to thermionic cooling.
The above statements highlight the intricate role of thermionic emission in the determination of the surface deformation profile of these W leading edge exposed samples. The escaping thermionic current density not only dictates the melt acceleration through the volumetric Lorentz force, but also affects melt production through thermionic cooling. In the case that the escaping emission is underestimated, the amount of melt will substantially increase, but melt acceleration will decrease resulting to smaller displacements for a given melt layer. On the other hand, an overestimation of the escaping emission will result in stronger melt acceleration but less melting due to the efficient cooling.

Summary and outlook
The MEMOS-U code simulates the melting, macroscopic liquid motion and surface deformation induced by incident plasma fluxes on metallic PFCs of magnetic confinement fusion devices. The macroscopic melt dynamics are selfconsistently described by the set of incompressible resistive TEMHD equations and temperature evolution is governed by the convection-diffusion equation. The subset of fluid equations is simplified within the shallow water approximation which transforms the three momentum equations and incompressibility condition to two momentum equations and a continuity equation for the melt layer column height. The electromagnetic field equations are simplified with the magnetostatic approximation for strong external magnetic fields which reduces Maxwell's equations and Ohm's law to a single, Poisson-like equation for an auxiliary potential that describes the bulk current pathways.
The MEMOS-U code is not only based on a transparent physical description but is also void of adjustable coefficients that could be utilized as free parameters owing to its rigorous boundary conditions and robust numerical schemes. The simulation outcome depends solely on the incident heat flux, local magnetic field and problem geometry for ELM-driven melting scenarios and also on the halo current density for melting induced by thermal plasma impact during the current quench phases of disruptions.
The physics model has been demonstrated to offer a unified description of ELM-induced W divertor melting and disruption-induced Be first wall melting, in spite of the significant differences in terms of heat loading, bulk current source and material composition [28]. A quantitative agreement with state-of-the-art fusion experiments performed in two tokamaks (JET and AUG), concerning two plasma scenarios (ELMs and VDEs) and involving two PFC geometries (sloped and leading edge), as well as two metal PFCs (W and Be) has been achieved, with the only variations allowed concerning the incident heat flux and strictly within the experimental uncertainties [28,29].
It has been demonstrated that the upgraded MEMOS-U code is equipped with the essential physical and numerical features that allow reliable simulations of significant displacements of metallic melts, i.e. of the order of the melt spatial extent, in fusion environments. These novel features have been discussed in detail in the present investigation and comprise the continuously deforming heat conduction domain, free-flow boundary condition at the melt pool edge and inclusion of convective heat transport. Concerning the latter, it has been revealed that the horizontal displacement of latent and sensible heat affects the temperature response in an intuitive, but important fashion, in which the colder adjacent material is heated by the displaced liquid and the initial induced melt area cools down. This process explains the gradual surface deformation and horizontal melt displacement under repeated transient melting scenarios (ELMy H-mode W exposures, discussed herein), or continuous melting scenarios (disrupting plasma Be exposures, discussed elsewhere [28]).
The shallow water assumption leads to a drastic reduction of the computational cost, making the MEMOS-U code a highly suitable simulation tool for parametric predictive studies and input sensitivity studies, such as investigations of the impact of heat flux input uncertainties and of the adopted nonambipolar current description on surface deformation. The natural trade-off has to do with loss of high-order accuracy relevant to the surface pinch-off in droplet ejection, instability growth, wetting dynamics, deep melting and melting of small scale or complex geometries. In the scenarios addressed here, re-solidification stabilizes the flow by controlling the melt depth and thus the viscous damping [28]. However, to reliably address the general flow stability and melt splashing, the full Navier-Stokes system needs to be employed. Simulation tools beyond the shallow water approximation are computationally heavy and cumbersome to use in parametric predictive studies. Nevertheless, the methodology is well developed, as mentioned in the Introduction.
Concerning transient W melting ITER scenarios, we reiterate that the present MEMOS-U description of the escaping emitted current is reliable for inter-, intra-ELM periods of contemporary fusion devices, where thermionic emission is dominant and becomes suppressed by space charge effects and prompt re-deposition. This expression should also be valid during inter-ELM ITER periods but is expected to be inaccurate during intra-ELM ITER periods due to role of electron-induced electron emission, the Schottky effect and the elevated incident plasma currents. An investigation has been initiated with the objective of identifying an accurate boundary condition for the escaping emitted current in this complicated regime [78], within the framework of the earlier studies [40,41,46].
After depth-averaging the momentum equation (2), and combining equation (A1) with equations (6) and (10), we acquire: There are subtleties involved with the first rhs term which describes the change of mean velocity in the column as liquid is added or removed by melting or solidification. For the no-slip condition v t | b1 = 0, this term represents the velocity change resulting from the addition or removal of stationary fluid at the bottom of the pool, without explicitly changing the parcel momentum. Two possible consequences of this term need to be considered which depend on the sign of ∂b 1 /∂t. When this quantity is positive, the melt layer thickness increases due to melting and there is a stationary liquid that is 'added' to an already moving column. The net effect is a reduction of the average velocity of the column. On the other hand, when it is negative, the melt layer thickness decreases due to solidification and the average velocity should increase since the material that is removed from the column is stationary.
It might appear that this term adequately captures the explicit effects of phase change on the melt column velocity. However, the term arises because of the prescribed parabolic velocity profile which neglects the finite viscous propagation speed in the liquid. Let us consider the case ∂b 1 /∂t < 0 in infinitesimal time steps. A small portion of stationary liquid is removed from the parcel, which will increase the column velocity average. Since now there is a slightly higher average velocity and a thinner melt layer, the parabolic velocity profile is assumed to 'equilibrate', adjusting to the new average, which implies a steeper gradient with higher maximum velocity value. This process repeats for each infinitesimal step, providing a nonphysical acceleration of the liquid. The opposite effect is true for the melting phases. It is evident that this process does not adequately describe the dynamics in the melt layer-the actual scenario needs a more detailed description of viscous propagation of melt velocity upwards in the melt layer as the liquidsolid boundary moves. In the absence of such a treatment, and in light of other model uncertainties, the term has been neglected. This is equivalent to stating that the velocity of the 'added/removed' fluid is identical to the average melt column velocity, v t | b1 = U.