Boundary conditions

Periodic boundary conditions are imposed at the vertical boundaries. The boundary con-dition at the lower boundary is identical to the one used for the 3D convection model of Chapter 3. The magnetic field at the lower boundary is kept vertical.

For the following simulations, we have modified the top boundary condition to allow for the bodily transport of magnetic field through the top boundary. This is important because we do not want the emerged magnetic field to be artificially trapped in the photo-sphere. In principle, stress-free conditions on all three components of the velocity at the upper boundary, i.e.

∂v_x

∂z = ∂v_y

∂z = ∂v_z

∂z = 0, (4.1)

suffice for a smooth outflow. This condition, however, is independent of whether the mass flux through the upper boundary is appropriate or unrealistically high. In order that the mass flux be kept at appropriate levels, additional constraints must be applied. To this

boundary of our simulation domain.

The bottom surface of our fiducial layer is co-incident with the horizontal plane con-taining the cell centres of the first ghost layer above the top domain boundary. To simplify the following discussion, let us use the position of this plane as the reference height, i.e.

z = 0 at the bottom of the fiducial layer. Using this reference, the height of the cell centres in the uppermost domain layer isz = −h, whereh is the vertical grid-spacing.

We impose that each vertical column of the fiducial layer be isothermal and be in hydro-static equilibrium. This means that each column of mass in the fiducial layer follow an exponential decay with height,

ρ(x, y, z) =ρ(x, y,0)e^−z/Hˆp, z ∈[0,3 ˆH_p] (4.2) whereHˆ_p is the pressure scale height evaluated at the uppermost domain layer. The fidu-cial layer has a thickness of3 ˆHp. Our aim is to specify an appropriate density distribution ρ(x, y,0)at the bottom of the layer in order to fill the values of density in the ghost cells.

To do this, we keep track of the total mass enclosed in the fiducial layer. By imposing the stress-free condition (5.1) on the velocity everywhere within the fiducial layer, we can determine the mass fluxes through the bottom and top surfaces of the layer. The total mass in the layer is given by is the mean density at the bottom of the fiducial layer. Expression (4.4) relates the mass of the fiducial layer with the mean density ρ¯₀ at the bottom of the layer. We require that density perturbations in the fiducial layer follow the same horizontal distribution as the corresponding perturbations in the uppermost domain layer. Let us denote the mean density in the uppermost domain layer asρ¯_D. Our requirement means that

ρ(x, y,0) =

At each time step, the r.h.s. of the above expression is calculated and Eqs. (4.2) and (4.5) are used to specify the density in the ghost cells. By keeping track of Mfid, we give memory to the fiducial layer. When too much mass has crossed from the domain into the fiducial layer, ρ¯₀ becomes larger than ρ¯_D and so the stratification near the top boundary becomes top-heavy. Any further mass outflow through the top boundary is then suppressed by gravity. Conversely, in case too much mass has drained from the fiducial layer into the domain, ρ¯₀ becomes much smaller than ρ¯_D, inducing a upwards directed pressure gradient that overcomes gravity to drive an outflow from the domain into the fiducial layer.

Besides the density, we also need to specify the values of the specific internal energy

²in the ghost layers. This quantity is taken to be uniform over the fiducial layer. At a particular time stepn, we specify this quantity in the ghost cells using

²ⁿ= [1−δ]²ⁿ⁻¹ +δ¯²ⁿ_D, (4.6)

4.1 Simulation setup

where²ⁿ⁻¹ is the value from the previous time step and ¯²ⁿ_D is the horizontally-averaged value evaluated in the uppermost domain layer at time stepn. The small quantityδcan be varied to control how quickly we allow² in the ghost layers to adjust to the values in the domain. In all the simulation runs described in the following, we usedδ = 10⁻³ so that the adjustment is slow. This is desirable because we do not wish our fiducial layer to excite waves with periods comparable to the magnetohydrodynamic time scales of interest in the simulations.

The magnetic field above the upper boundary of the domain is matched to a potential field. This requires an extrapolation of the magnetic field in the uppermost domain layer into the ghost cells at each time step. The method we used for the field extrapolation is appropriate for both potential fields and more generally, linear force-free fields. For further details, we refer the reader to Appendix C.

4.1.2 Initial conditions

The background stratification in the domain is initially plane-parallel and static. The horizontally-averaged profiles of² andρ as functions of height were taken from the 3D model of Chapter 3. In the presence of convection, as is the case with the 3D model, the energy flux carried by radiation at the surface is maintained by a continuous replenishment of upwelling, high-entropy material. This process keeps the average height of the visible surface approximately constant. Since there is no convection in the initial configuration here, the cooling of the visible surface causes the visible surface to sink in the course of the simulation. In the following simulations, the visible surface sinks by between100and 300km over a duration of25minutes (the typical time scale of an emergence event). This side effect has relatively unimportant consequences compared with the alternative option, which would be a removal of the radiative heating term from the energy equation. Such a modification would eliminate what is arguably the most interesting result of this set of simulations, namely, the effect of radiative transfer on an emerging flux tube.

Att = 0, an axisymmetric magnetic flux tube was introduced into the sub-photospheric layers. The longitudinal and transverse components of the magnetic field have the form:

B_l(r) = B₀exp (−r²/R²₀), (4.7) B_θ(r) = λr

R0

B_l, (4.8)

wherer ∈[0,2R0]is the radial distance from the tube axis andR0 the radius of the tube.

Its longitudinal flux is given byΦ₀ = R

B_ldS = 0.98πR₀²B₀. λ is the twist parameter.

In the absence of the magnetic field, the divergence of the total stress tensor (viscous+ Maxwell+pressure) is zero. Since we are primarily interested in studying the buoyant rise of the flux tube, we require that this condition holds when the flux tube is introduced.

In such a case, the flux tube experiences only the gravitational force at timet = 0. Since the magnetic field contributes to the total pressure locally, the gas pressure in the flux tube must be decreased accordingly. Additionally, the magnetic tension of the transverse field (if present) exerts a force directed radially inwards, which must be balanced by an appropriately chosen pressure distribution. Having specified the internal gas pressure dis-tribution, one is still free to choose the distribution of one and only one of the following

tities constrains the distributions of the remaining two.

4.2 Simulation results

In this section, we describe in detail results from two simulation runs. In both of these runs, the flux tube hasB₀ = 8500G, R₀ = 200km and Φ₀ = 10¹⁹ Mx. The flux tubes have the same internal specific entropy as the average value for the upflows in the 3D model of Chapter 3 (s = 6.0R^?). Both are placed1.35Mm below the visible surface at t= 0. In the first run, the flux tube is untwisted (λ = 0). In the second run, the flux tube has twist (λ = 0.5). The relative density deficit is largest at the axes of the tubes, about 40%. In both cases, the density deficit averaged over the whole tube is much smaller, about5%.