Self-consistent description of the replacement current driving melt layer motion in fusion devices

The bulk replacement current density triggered by surface charge loss owing to thermionic emission leads to a volumetric Lorentz force which has been observed to drive macroscopic melt layer motion in transient tungsten melting tokamak experiments in which components of different geometries (deliberate leading edges and sloped surfaces) have been exposed to edge localized mode (ELM) pulsed heat loads in high power H-mode discharges. A self-consistent approach is formulated for the replacement current which is based on the magnetostatic limit of the resistive thermoelectric magnetohydrodynamic description of the liquid metal and results in a well-defined boundary value problem for the whole conductor. A new module is incorporated into the incompressible fluid dynamics code MEMOS-3D, which numerically solves the finite difference representation of the problem. The phenomenological approach, employed thus far to describe the replacement current, is demonstrated to be accurate for the sloped geometry but inadequate for the leading edge. MEMOS-3D simulations of very recent ASDEX-Upgrade leading edge experiments with the rigorous as well as the simplified approach are reported. For these simulations, the self-consistent approach predicts a fivefold reduction of the displaced material volume, a sevenfold reduction of the maximum peak height of displaced material and a different eroded surface morphology in comparison with the previously applied simplified approach.


Introduction
ITER will begin operation with a full tungsten (W) divertor and the design of the W monoblocks needs to ensure sufficient power handling capability as well as a divertor lifetime extending well into the nuclear phase [1,2]. In spite of the high W melting point (T m = 3695 K), the risk of W melting during plasma exposure remains one of the major concerns together with deep W recrystallization (T c = 1300−1600 K). The transient or sustained liquid metal is acted upon by various forces which displace the material and may cause large-scale surface erosion, thus invariably compromising the power handling capabilities of the monoblocks in subsequent exposures. Extensive work is currently being carried out in order to categorize the heat loading characteristics of divertor tiles with different surface shaping geometries in various ITER operation Nuclear Fusion Self-consistent description of the replacement current driving melt layer motion in fusion devices E. Thorén  The bulk replacement current density triggered by surface charge loss owing to thermionic emission leads to a volumetric Lorentz force which has been observed to drive macroscopic melt layer motion in transient tungsten melting tokamak experiments in which components of different geometries (deliberate leading edges and sloped surfaces) have been exposed to edge localized mode (ELM) pulsed heat loads in high power H-mode discharges. A self-consistent approach is formulated for the replacement current which is based on the magnetostatic limit of the resistive thermoelectric magnetohydrodynamic description of the liquid metal and results in a well-defined boundary value problem for the whole conductor. A new module is incorporated into the incompressible fluid dynamics code MEMOS-3D, which numerically solves the finite difference representation of the problem. The phenomenological approach, employed thus far to describe the replacement current, is demonstrated to be accurate for the sloped geometry but inadequate for the leading edge. MEMOS-3D simulations of very recent ASDEX-Upgrade leading edge experiments with the rigorous as well as the simplified approach are reported. For these simulations, the self-consistent approach predicts a fivefold reduction of the displaced material volume, a sevenfold reduction of the maximum peak height of displaced material and a different eroded surface morphology in comparison with the previously applied simplified approach.
Keywords: melting, tungsten, divertor, melt layer motion, MEMOS, thermionic emission (Some figures may appear in colour only in the online journal) Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. regimes [3,4]. Even in the pre-nuclear phase, while full surface melting will probably not occur provided that disruption mitigation techniques are successful, it seems unlikely that the transient melting of tile corners and edges can be avoided during uncontrolled edge-localized modes (ELMs) [3].
Transient and steady state W melting tokamak experiments, performed on contemporary devices in support of ITER divertor component shaping choices [5][6][7][8], had featured deliberately misaligned leading edges (intercepting nearly the full parallel power flux) and protruding sloped surfaces (intercepting a respectable amount of the full parallel power flux). During these exposures, macroscopic W melt motion has been primarily observed along the (poloidal) direction consistent with the volumetric Lorentz or J × B force, where B denotes the toroidal magnetic field and J the volume current density passing through the molten layer [5][6][7][8]. The origin of the current density has been attributed to replacement electrons flowing through the bulk of the conductor towards the plasmawetted surface in order to replenish locally emitted electrons [1,5], as dictated by charge continuity. At the elevated temperatures relevant for liquid W, electron emission is primarily due to thermionic emission which refers to over-the-surfacebarrier electron release by thermal excitation [9,10]. It is worth mentioning that for the high plasma temperatures experienced during ELMs, the contribution of electron induced electron emission could also be significant [11]. In fact, secondary electron emission has been suggested as the trigger for the replacement current responsible for the motion of molten beryllium observed on the main chamber surfaces of JET [12,13].
MEMOS-melt motion at surfaces-is a finite difference code that solves the heat conduction equation coupled with the incompressible Navier-Stokes equations [14,15] within the shallow water approximation [16,17]. MEMOS-3D has been previously employed for the modelling of W melt experiments carried out in TEXTOR as well as JET [5,6] and is now being used in support of the interpretation of the multiple recent exposures in the ASDEX-Upgrade (AUG) divertor [8,18,19]. For MEMOS-3D to acquire predictive power, given the prominent role of the J × B force density in driving melt layer motion, efforts should mostly concentrate on refining the description of the replacement current density responsible for this force. The task essentially consists of two problems: (1) Thermionic emission into the magnetized plasma sheath: the Richardson-Dushman formula for the thermionic current density is valid in the absence of ambient electromagnetic fields [9,10]. However, space-charge accumulation can generate repelling electrostatic fields that cause electron recapture [20] and magnetic fields can lead to prompt electron re-deposition during the first Larmor gyration [21]. Recent progress has been made in the quantification of the escaping thermionic cur rent based on collisionless particle-incell (PIC) simulations with prescribed sheath edge properties [22,23]. (2) Current flow through an electron-emitting conductor: a simplified approach has been adopted thus far, summarized in [1], where it is assumed that the radial component of the replacement current density (leading to a force in the poloidal direction as observed in tokamaks) is constant inside the molten layer with a magnitude equal to the thermionic current density. This phenomenological description neglects all details of the cur rent paths inside the W sample and lacks a solid theoretical basis. The aim of this paper is to provide such a foundation. In this work, a self-consistent approach is formulated for the determination of the replacement current density flowing through an electron-emitting, partially molten, grounded conductor. The treatment is based on the magnetostatic limit of the resistive thermoelectric magnetohydrodynamic description of the liquid metal and results in a well-defined boundary value problem for the whole conductor. To handle the numer ical task, MEMOS-3D has been equipped with a replacement current module that has been coupled with the existing routines. The adequacy of the simplified approach is discussed for different exposure geometries. MEMOS-3D simulations of ASDEX-Upgrade melting experiments are reported which demonstrate the importance of a selfconsistent description of the replacement current.

Self-consistent approach for the replacement current density
Let us consider a grounded bulk metal sample that is locally emitting valence electrons to its surroundings. Unless the ambient somehow compensates for this charge loss, current conservation requires that the released electrons are replenished by electrons flowing from the external circuit. Microscopically, electron emission leads to the generation of electronic holes as well as local electrostatic fields. Bound electrons in the vicinity of the hole are swiftly attracted to fill the valence states that are vacant not only due to emission, but also due to thermal excitation. This disturbance propagates along the conductor generating a current which passes through the volume of the body. The microphysics of the replenishing electrons has been known to play a significant role in determining local energetic aspects during thermionic and field emission (Nottingham effect) [9,[24][25][26]. In the disciplines of particle emission and surface physics, these electrons are typically coined as replacement electrons [9,[24][25][26]. In what follows, we will therefore employ the term replacement current to describe bulk currents triggered by charge loss due to electron emission.

Description of the liquid metal
The macroscopic motion of liquid metals is governed by the set of incompressible resistive thermoelectric magneto hydrodynamic (TEMHD) equations together with the convectiondiffusion equation for the temperature [27]: where (v, J, B, E, T, p) denote the fluid velocity field, current density, magnetic flux density, electric field strength, fluid temper ature, scalar pressure and where (ρ m , µ, c p , k, σ e , S, µ 0 ) denote the mass density, dynamic viscosity, specific isobaric heat capacity, thermal conductivity, electrical conductivity, absolute thermoelectric power, free space permeability. We point out that all thermophysical properties of tungsten are functions of the local temperature. Note also that viscous dissipation terms have been neglected in the temperature equation. Thermoelectric phenomena manifest themselves through an additional term in Ohm's law (−σ e S∇T representing the Seebeck effect) and an additional term in the temperature equation (−∇ · [STJ] representing the Thomson effect). Thermoelectricity has been retained, but it can expected to be mostly important when different materials are in contact that are characterized by noticeably different absolute thermoelectric power S [27,28]. The above set of equations is Galilean invariant and is thus accurate for (v/c) 2 1 [29]. For the fusion applications of interest, the liquid metal characteristic speeds also satisfy v/c 1. Within this magnetostatic approximation [30], Faraday's law leads to the irrotational condition for the electric field and the induced current densities can be neglected in Ohm's law. Moreover, Ampere's law can be straightforwardly converted to the solenoidal condition for the current density. Furthermore, by eliminating the electric field in the temperature equation using Ohm's law and employing ∇S = (∂S/∂T) ∇T , the standard Joule and Thomson heating terms emerge. Overall, after introducing the reciprocal electrical resistivity ρ e = 1/σ e for convenience, we end up with Finally, the electrostatic condition and Ohm's law can be combined to yield an equation for the rotation of the current density. We first need to introduce the position vector inside the liquid metal interior r with respect to an arbitrary frame of reference. In the case of uniform material composition where S(T, r) ≡ S[T(r)], we have ∇ × (S∇T) = 0, which implies that thermoelectric effects do not influence the replacement current [27]. We obtain We point out that the total magnetic field is given by the external magnetic field and the internal magnetic field selfconsistently generated by the replacement current. Within the magnetostatic approximation, the latter is described by For the large tokamak magnetic fields, it can be assumed that |B int | |B ext | leading to B = B ext everywhere inside the liquid metal.
Even when neglecting material boundary motion, the field subset, see equations (4)- (6) or alternatively equations (5) and (7), remains explicitly coupled to the temperature equation of the fluid subset, see equations (1)-(3), through Joule heating, Thomson heating and the temperature dependence of the electrical resistivity. In the presence of thermionic emission, there is also implicit coupling through the temperature dependence of the boundary condition for the replacement current. As aforementioned, in previous simulation efforts an oversimplified picture was adopted for the determination of the current density [1]. In particular, the field subset was neglected and the current density was considered constant throughout the molten layer and equal to the impeded thermionic current density. The validity of this approach depends on details of the exposure geometry and will be discussed in section 2.6.

Description of the whole sample
In the magnetostatic approximation, both the liquid and the solid part of the sample are essentially described by the same field subset. Nevertheless, this does not imply that equations (5) and (7) can be applied inside the whole volume of interest owing to the electrical conductivity discontinuity at the melting point. The discontinuity necessitates separate application of equations (5) and (7) in the solid and liquid volumes as well as solution matching through the boundary conditions ρ e,l J t,l = ρ e,s J t,s , where the subscripts (t, n) denote the transverse and normal components with respect to the moving melting front, the subscripts (l, s) denote the liquid and solid side of the discontinuity. The numerical solution of the exact problem would manifoldly increase the computational cost, requiring the exact identification of the moving liquid-solid boundary. This can be circumvented by defining an artificial phase change temperature range, within which the electrical resistivity is assumed to linearly vary between its solid and liquid phase value [31]. Physically, this implies that the discontinuous bending of the current density field lines imposed by the boundary conditions will now be smoothed out in space, which leads to large local but small global numerical errors. The main advantage stemming from the utilization of an apparent electrical resistivity, ρ app e (T), lies in the application of equations (5) and (7) in the entire volume. The MEMOS-3D ρ app e implementation for W is illustrated in figure 1(a).

Equation for the auxiliary potential
It is now straightforward to convert the field subset into a single equation to be solved in the entire volume V of the sample. By defining the auxiliary potential ρ app e J = −∇ψ, equation (7) is automatically satisfied and equation (5) In the case of uniform material composition we have ρ app e (T, r) ≡ ρ app e [T(r)]. By considering that spatial gradients of the resistivity are solely caused by its temperature dependence and expanding the divergence, we obtain Naturally, when neglecting the temperature dependence of the electrical resistivity, the whole problem collapses to the Laplace equation ∇ 2 ψ = 0. The additional term expresses the fact that the current density will choose the least resistive path possible that is consistent with the boundary conditions. The material-dependent pre-factor ∂/∂T ln (ρ app e ) is plotted in figure 1(b). At this point, we should emphasize that, owing to the Seebeck effect, the auxiliary potential ψ is not equal to the conventional electrostatic potential E = −∇φ, as originally noted by Shercliff [27]. The equality ψ ≡ φ is only true in the absence of temperature gradients.

Boundary conditions
We decompose the surface boundary ∂V of the sample into the plasma-wetted (S w ), plasma-shadowed (S s ) and vesselconnected (S g ) parts with S w S s S g = ∂V . Considering first the plasma-wetted part; the current continuity equation, equation (5), leads to the boundary condition J n = J em , where J em denotes the current density associated with the emitted electron flux under the ambipolarity assumption for the incident electron and ion fluxes to the grounded sample. This expression constitutes a Neumann boundary condition for the auxiliary potential. In addition, after neglecting electroninduced electron emission processes (valid for the ELM parameters considered herein [11]), we obtain J em = J th .
Considering next the plasma-shadowed part; we simply have J n = 0, which constitutes another Neumann condition for the potential. Finally, considering the vessel-connected part; the auxiliary potential can be considered to be equal to an arbitrary constant, chosen to be zero (grounded vessel and nearly uniform temperature), which constitutes a Dirichlet boundary condition. Overall, we are dealing with the mixed boundary conditions where the r-dependence of the thermionic current density J th and of the apparent resistivity ρ app e stems from their T(r) dependence.
The boundary value problem consists of the partial differential equation (10) complemented by the mixed boundary conditions described by equations (12)- (14). As a consequence of the positivity of the electrical resistivity, the problem is wellposed. This can be rigorously shown by considering the volume integral V (1/ρ app e ) (∇ψ) 2 dV > 0 and following the standard methodology developed for the Poisson equation [33]. Note that J th is the only source term present in the boundary value problem. In its absence, the problem has a trivial solution and there is no current flowing through the conductor. This justifies the use of the term replacement current also from a macroscopic point of view. Finally, we emphasize that J th , or J em in general, is introduced as an input to the MEMOS-3D code and should be determined from external models.

Numerical implementation
The boundary value problem is solved in the new replacement current module of MEMOS-3D. The finite difference approx imation is employed for the spatial derivatives of each Instead of employing a direct inversion algorithm for Ψ = A −1 B, the generalized minimal residual method (GMRES) [34] is utilized. In the GMRES method, the solution Ψ n is approximated by minimizing the Euclidean norm of the residual R k = AΨ k − B, where Ψ k belongs to the Krylov subspace K k = span{V, AV, ..., A k−1 V} at iteration k and V = R 0 /||R 0 || is the normalized residual corresponding to the initial guess Ψ 0 . To speed up the convergence, the preconditioned system MAΨ = MB is solved, where M ∼ A −1 is found by the incomplete lower upper factorization of A. The matrix MA should have a smaller spread of eigenvalues compared to A and a satisfactory solution to the system can be found within fewer iterations [35,36].
Finally, we point out that, since ρ app e varies only weakly between successive heat conduction time steps, the matrix A should remain approximately the same over significant time intervals. For a static grid, this implies that the same preconditioner matrix M can be employed for several time steps, drastically reducing the overall computation time.

Comparison to the simplified approach
Here we demonstrate the necessity for the self-consistent approach in the case of leading edge exposures and the adequacy of the simplified approach in the case of sloped exposures. We shall discuss specific ASDEX-Upgrade experiments, but the conclusions reached will be of generic validity. In these experiments, W samples were exposed in the outer divertor target during Type-I ELMing H-mode discharges [8,19]. Two special samples were misaligned by ∼1 mm upwards with respect to the adjacent plate: the leading edge sample intersecting the magnetic field at near-perpendicular incidence α ∼ 90 • and the sloped sample at an angle α ∼ 17.5 • . Transient W melting due to the ELM heat pulses was achieved for both samples in repeated prolonged exposures.
The sloped and leading edge exposure geometries are depicted in figures 2(a) and (b), where the respective plasma wetted areas receiving the most intense heat fluxes are highlighted. The local cartesian coordinate systems implemented in MEMOS-3D are also displayed. In order to facilitate the discussion on the force density components acting on the liquid W, representative sketches of the molten layers and the current pathways of each sample are depicted in figures 2(c) and (d). In both types of exposure, the exper imentally observed melt layer motion was directed mainly towards the private flux region (−ẑ) and was attributed to the J × B force density  ((a), (b)) Sketch of the geometry of the W samples exposed in ASDEX-Upgrade. The plasma wetted areas which receive the most intense heat fluxes (where the melt layer emerges) are indicated in red, whereas the vessel connected areas (which serve as ground reference) are indicated in blue. ((c), (d)) Zoom in sketch of the samples in the x −ŷ plane, taken at an arbitrary z coordinate. The typical shape of the molten layer is indicated (shaded in pink) and the characteristic paths of the replacement current vector are depicted, as revealed from simulations. The MEMOS-3D local coordinate system (cartesian) is also displayed. Notice the rotation of 15 • around ẑ between the two sample geometries. Note that for the leading edge exposure: x is oriented roughly along the toroidal, ŷ roughly along the radial and ẑ roughly along the poloidal direction. generated by the ŷ -component of the bulk cur rent density [1,8].
Sloped exposure. The force density along ẑ becomes f z −J y B x , since the J x B y contribution is negligible courtesy of B x B y . The liquid pool is very shallow in the normal direction (∼100 µm along ŷ) but extends much further tangentially (1-2 cm along ẑ,x). Moreover, the tangential gradients of J th are small, since, for the 17.8 • B-field inclination, the space-charge limited regime is reached for surface temperatures below the W melting point [23]. Given the above, it is reasonable to assume that the tangential variations of the replacement current are rather weak inside the molten layer. Then, equation (5) leads to ∂J y /∂y ∼ 0 and thus J y J th inside the liquid layer, which is the assumption of the simplified approach. This is not valid near the molten pool edges, where J x , J z can have significant variations as depicted in figure 2(c). Naturally such transition regions only constitute a small part of the entire liquid volume and J y J th or f z −J th B x is deemed to be a sufficiently accurate assumption for the sloped surface.
Leading edge exposure. The force density along ẑ is again approximated by f z −J y B x . However, as illustrated in figure 2(d), both the tangential and normal current components exhibit strong local variations. In particular, J y is dominant in the lowest liquid region but becomes nearly zero at the upper corner, hence its spatial dependency cannot be ignored. More importantly, the exact boundary condition of equation (12) requires that J x = J th and is violated by the simplified approach. Previous MEMOS modelling of leading edge exposures [6,18] made use of the simplified assumption that J y = J th , which should lead to a significant overestimation of the main driving force. For the leading edge, the force density averaged over the melt layer height f z (y, z) −(1/h) h 0 J y (x, y, z)B x dx, as introduced in the shallow water version of the Navier-Stokes equations, requires the input of J y (x,y,z) from the solution of the magnetostatic boundary value problem. The difference between the simplified and self-consistent approaches is quantified in the following section.

Reference exposure and simulated scenarios
Six identical sequential H-mode AUG discharges were dedicated to the ELM-induced melting of the W leading edge [8]. The discharge parameters were B t = −2.5 T, I p = 0.8 MA, P = 7.5 MW(NBI), 2 MW(ECRH). For a detailed description of the experiments, the reader is directed to [8]. Here, we shall focus on pulse #33504. The full spatiotemporal heat flux data from infra-red camera measurements (corrected for artifacts due to reflections and brehmsstrahlung) have been employed instead of the simpler coherently averaged data [18]. The heat flux received by the W sample has been deduced by employing the optical approximation (namely without taking into account Larmor orbits of the incoming plasma ions) and assuming toroidal symmetry [8]. The heat flux incident on the leading edge of the sample is nearly equal to the parallel power flux with intra-ELM values of ∼1 GW m −2 and inter-ELM values of ∼50 MW m −2 . The typical ELM duration is ∼3 ms and the average inter-ELM period is ∼11 ms, corresponding to an average ELM frequency of ∼70 Hz. The outer strike point was shifted down to its stationary position on the sample at t = 3.3 s and maintained there until t = 3.8 s, when it was shifted back upwards.
As already mentioned, the escaping thermionic current is described by the Richardson-Dushman formula only in absence of ambient electromagnetic fields [9,10]. The emergence of the space-charge limited regime [20] for the AUG leading edge exposures has been recently investigated by dedicated PIC simulations and an analytical expression describing the limited escaping current as a function of the thermal plasma electron current has been proposed [23]. Under reasonable additional assumptions, this expression has been transformed into a dependence on the heat flux incident at each point of the plasma-wetted surface [18]. Combining with the Richardson-Dushman formula, we obtain 4 where A eff 60 A cm −2 K −2 and W f = 4.55 eV, T s denotes the surface temperature and k B the Boltzmann constant. The limiting temperature T lim s (q ) is the minimum surface temperature for which a virtual cathode forms. Its value depends on the incident heat flux and can be found by the continuity of the escaping thermionic current as a function of the surface temper ature. Equation (15) is employed for the current emitted from surface elements of the leading edge, whose normal is oriented nearly parallel to the magnetic field (α ∼ 87.2 • ). The current emitted from the flat top surface of the sample is assumed to be zero since, at such shallow B-field angles (α ∼ 2.8 • ), the escaping thermionic current should be drastically suppressed by prompt re-deposition and thus become negligible compared to the leading edge current [22,23].
We point out that currently in MEMOS-3D, heat conduction is decoupled from fluid motion (the advective term in the temperature equation is neglected) and that the Thomson heating term is not considered due to the lack of reliable data for the absolute W thermoelectric power at temper atures exceeding 2000 K [37,38]. We also note two simplifications employed in the present simulations: (i) Melt layer motion is only considered along the ẑ direction. MEMOS-3D predicts small but non-negligible Marangoni flows along the ŷ direction caused by sharp surface temperature gradients near the exposed edge. This motion has been excluded. Our conclusions will not be affected, since we aim to elucidate the behavior of the ẑ force component acting on the melt layer under different replacement current descriptions. (ii) 4 Note a misprint in [18], where the saturation limit was erroneously given Simplified grounding of the sample. It has been assumed that the entire lower sample surface is vessel-connected, while the actual grounded region in the experiment consisted of a small circular contact ( 3 mm) located at the bottom of the sample (see figure 2). The current density distribution inside the molten layer is not sensitive to the position and extent of the vessel-connected part of the sample, as long as these parts are placed sufficiently far away from the surface that provides the non-trivial boundary condition, i.e. the thermionically emitting surface. This is a consequence of the smoothing properties of the elliptic partial differential equation. By considering the small extent of the melt layer (∼0.1 × 1 mm 2 ) compared to the distance to the electrical ground (∼10 × 10 mm 2 ), this simplification is justifiable.
Three simulations have been performed which follow different approaches for the replacement current: • Scenario 1, simplified approach. This is the description previously employed in MEMOS-3D and utilized in the modelling of the TEXTOR steady state and JET transient melting experiments [5,6], which assumes that the component of the replacement current density responsible for the J × B force density is constant inside the molten layer with a magnitude equal to the thermionic current density [1]. • Scenario 2, full self-consistent approach. The replacement current density is given by the solution of the boundary value problem involving the partial differential equation ∇ · 1/ρ app e ∇ψ = 0. The problem is solved for each time-step of the fluid solver. • Scenario 3, simplified self-consistent approach. The replacement current density is given by the solution of the boundary value problem involving the Laplace equation ∇ 2 ψ = 0, which corresponds to the heuristic limit of temperature-independent resistivity. The problem is again solved for each time-step of the fluid solver. Figure 3 illustrates the current field-lines in a sample x −ŷ cross section at 3 ms after the onset of an ELM, when the melt layer thickness is ∼150 µm. The emission is limited everywhere at the constant value of J th ∼ 13 MA m −2 . The contour plots illustrate the temperature (a), and the y−component of the replacement current (b) close to the molten corner. The cross section is taken roughly at the centre of the sample. Prior to the analysis of the replacement current results, some characteristics of the emitted current from the molten layer should be mentioned. The wetted surface typically melts near the ELM pulse maximum and remains liquid until nearly half of the inter-ELM period [18]. The limited emitted current density scales with the incoming parallel power flux, see equation (15). For the given plasma parameters, saturation is reached for T s ∼ 3950 K during intra-ELM periods and the limiting value is ∼15 MA m −2 , whereas during inter-ELM periods saturation is reached for T s ∼ T m = 3695 K and the limiting value is ∼5 MA m −2 . This implies that for T s > 3950 K, emission is space-charge limited even during the ELM peak, while in-between ELMs the molten surface will nearly always lie in the saturated regime. The value of the limited current is determined by the power flux density, which depends on the time and z coordinate, but does not vary with the y coordinate. Therefore, a typical situation encountered during a melt event is that the emitted current is limited and constant with respect to y, i.e. ∂J th /∂y = 0. The solution of the full self-consistent problem features an abrupt change of J y at the liquid-solid interface, which is clearly visible in figure 3(c). This drop stems from the artificial smoothing of the local discontinuity of the electrical resistivity, discussed in section 2.2 and illustrated in figure 1. Clearly, in the absence of the apparent resistivity, the drop will also be discontinuous. In the solution of the simplified selfconsistent problem, see figure 3(d), the drop does not exist and the average J y is visibly higher. The current pathways for scenarios 2, 3 have been overlaid to highlight the difference in the bending of the streamlines. In particular, the temperature dependent resistivity facilitates the prompt alignment of the current density with the leading edge surface normal and thus leads to a reduction of the J y component within the molten layer and thus also of the resulting J × B force density.

Numerical results for the replacement current
It is apparent that J y varies weakly across the melt layer for the specific snapshot in figure 3. The relevant J y spatial gradients are actually similar for all times, originating from the distribution of the emitted current density and not its absolute value. Since the melt layer thickness never exceeds 200 µm, it can be concluded that, in general, J y does not strongly vary within the melt layer height. This justifies the application of the shallow water averaging procedure (employed in MEMOS for the pure mechanical Navier-Stokes terms), also for the electromagnetic J × B term.

Numerical results for the force density and melt displacement
The pertinent question concerns how the aforementioned differences in the replacement current propagate up to the average force density and the melt displacement profile. In table 1, three relevant quantities have been listed that enable a quantitative comparison between the three scenarios; the average force density at the specific time and for the z coordinate of figure 3, the maximum peak height of the final surface profile and the volume of displaced material at the end of the exposure 5 . As expected, Scenario 1 leads to the most pronounced liquid motion, whereas Scenario 2 leads to the most restricted liquid motion. Compared to the simplified approach, the full self-consistent approach results in five times less displaced material volume and a seven times smaller maximum peak height, emphasizing the importance of a self-consistent description of the replacement current. On the other hand, dropping the temperature-dependence of the resistivity from the self-consistent approach results in appreciable but not drastic differences.
The resulting leading edge surface profiles are illustrated in figure 4, for Scenarios 1 and 2. They refer to the final surface topography corresponding to the end of the exposure. Apart from the difference in the magnitude of surface deformation, the two melt profiles are also characterized by different morph ology. For Scenario 1, there is a monotonic increase of the melt displacement towards the edge (y = 0). For Scenarios 2 and 3, the melt displacement acquires a maximum at y 5 mm. This is a consequence of the self-consistent behaviour of the replacement current. To be more exact, at y = 0 the current field lines have completely bent, becoming parallel to the emitting surface normal across the whole molten layer (which is directed along x) in order to satisfy the exact boundary condition, which implies that J y (y = 0, x, z) = 0, ∀x, z, i.e. there is no local J × B force and no local melt acceleration. Weak emission from the flat top surface (totally neglected here) would create a small J y (y = 0, x, z) component, which would still be insufficient to drive local melt motion owing to the strong prompt re-deposition of thermionic electrons.

Summary and conclusions
The bulk replacement current triggered by surface charge loss owing to thermionic emission is responsible for the J × B force density that drives macroscopic melt layer motion in the transient W melting experiments which have been conducted in recent years on tokamaks in support of ITER divertor design (a  decisions. In the absence of a rigorous approach, it has been thus far assumed that the radial component of the replacement current density (leading to a force in the poloidal direction consistent with the melt displacement observed) is constant inside the molten layer with a magnitude equal to the escaping thermionic current density. In this work, starting from the resistive thermoelectric magnetohydrodynamic description of the liquid metal and invoking the magnetostatic approximation, we have formulated a self-consistent approach for the determination of the replacement current density flowing through the whole W sample (liquid and solid). The approach reduces to a partial differential equation for an auxiliary scalar potential that is complemented by well-defined mixed boundary conditions. The finite difference representation of this problem results in a large linear system of equations, which can be efficiently solved by an iterative procedure. A replacement current module has been developed which numerically solves the boundary value problem and has been incorporated in the structure of the MEMOS-3D melt motion code. The self-consistent replacement current and its impact on gross surface erosion have been evaluated for relevant geometries and boundary conditions in simulations of ASDEX-Upgrade experiments performed with the upgraded MEMOS-3D. In addition, it has been verified numerically that the current density varies weakly across the melt layer thickness, justifying the extension of the shallow water averaging procedure to the J × B force density. For leading edge exposures, the simplified approach fails to satisfy the exact boundary condition at the emitting surface and the selfconsistent approach needs to be employed. The latter leads to a fivefold reduction of the displaced material volume, a sevenfold reduction of the maximum peak height as well as a different eroded surface morphology. Consequently, earlier simulation efforts of the JET leading edge experiments would need to be revisited [6]. For sloped exposures, an investigation of the geometry of the configuration revealed that the simplified approach suffices for the description of the replacement current.
We emphasize that the objectives of this work have been to formulate a self-consistent approach for the replacement current density and to demonstrate its importance in tokamak experiments by simulating very recent ASDEX-Upgrade exposures. A detailed comparison with the post-mortem observed melt patterns would require the modelling of multiple exposures, which lies out of the scope of the present investigation. EUROfusion WP MST1. The views and opinions expressed herein do not necessarily reflect those of the European Commission and of the ITER Organization. PT and SR would also like to acknowledge the financial support of the Swedish Research Council.