A Cancellation Nanoflare Model for Solar Chromospheric and Coronal Heating II. 2D Theory and Simulations

TL;DR

This paper develops a 2D resistive MHD simulation to study magnetic flux cancellation as a mechanism for solar atmospheric heating, confirming analytical estimates and exploring atmospheric responses like hot and cool ejections.

Contribution

It advances the theory by providing 2D simulations that validate analytical energy estimates and examine the atmospheric effects of magnetic reconnection during flux cancellation.

Findings

Simulation results agree qualitatively with analytical energy estimates.

Reconnection can produce hot coronal-like ejections or cool surges depending on height.

Atmospheric responses include hot coronal loops and cool spicule-like structures.

Abstract

Recent observations at high spatial resolution have shown that magnetic flux cancellation occurs on the solar surface much more frequently than previously thought, and so this led Priest et al 2018 to propose magnetic reconnection driven by photospheric flux cancellation as a mechanism for chromospheric and coronal heating. In particular, they estimated analytically the amount of energy released as heat and the height of the energy release during flux cancellation. In the present work, we take the next step in the theory by setting up a two-dimensional resistive MHD simulation of two canceling polarities in the presence of a horizontal external field and a stratified atmosphere in order to check and improve upon the analytical estimates. Computational evaluation of the energy release during reconnection is found to be in good qualitative agreement with the analytical estimates. In addition, we go further and undertake an initial study of the atmospheric response to reconnection. We find that, during the cancellation, either hot ejections or cool ones or a combination of both hot and cool ejections can be formed, depending on the height of the reconnection location. The hot structures can have the density and temperature of coronal loops, while the cooler structures are suggestive of surges and large spicules.

Figures

Tables

Table 1. Table 1: Initial Conditions for the Simulations

Name	$B_0$ (G)	$y_N$ (Mm)
Case 1	45	7.6
Case 2	210	2.9
Case 3	300	2.2
Case 4	360	1.8
Case 5	600	0.9

Equations

$$d_0^{3D}={\left(\frac{F}{\pi B_0}\right)}^{1/2}.$$

$$v_i=\frac{\eta }{l}.$$

$$L{v}_i=l{v}_{Ai},$$

$$M_{Ai}=\frac{1}{R_{mi}^{1/2}},$$

$$\frac{dW}{dt}=\frac{v_i{B}_i^2}{\mu }L{L}_s=\frac{1}{R_{me}^{1/2}}\frac{v_{Ai}{B}_i^2}{\mu }L{L}_s,$$

$$L{v}_i=l{v}_{Ai},$$

$$\frac{dW}{dt}=\frac{4}{5}\frac{v_i{B}_i^2}{\mu }L{L}_s,$$

$$\mathbf{B}=\frac{F}{\pi }\frac{{\hat{\mathbf{r}}}_1}{r_1}-\frac{F}{\pi }\frac{{\hat{\mathbf{r}}}_2}{r_2}+B_0\hat{\mathbf{x}},$$

$${\mathbf{r}}_1=(x-d)\hat{\mathbf{x}}+y\hat{\mathbf{y}},{\mathbf{r}}_2=(x+d)\hat{\mathbf{x}}+y\hat{\mathbf{y}},$$

$$d_0=\frac{2F}{\pi B_0},$$

$${\bar{B}}_x=\frac{B_x}{B_0},\bar{d}=\frac{d}{d_0},\bar{y}=\frac{y}{d_0}.$$

$${\bar{B}}_x=\left(-\frac{\bar{d}}{{\bar{y}}^2+{\bar{d}}^2}+1\right).$$

$${\bar{y}}_N=\sqrt{\bar{d}-{\bar{d}}^2},$$

$$B_y+i{B}_x=kz,$$

$$B_y+i{B}_x=k(z^2+\frac{1}{4}{L}^2{)}^{1/2},$$

$$B_i=\frac{1}{2}kL,$$

$${\bar{B}}_x=2ϵ(1-\bar{d}),$$

$$B_x=2\sqrt{\frac{d_0-d}{d}}\left(\frac{y-y_N}{d_0}\right)B_0.$$

$$\frac{B_i}{B_0}=\sqrt{\frac{d_0}{d}-1}\frac{L}{d_0}.$$

$$v_i{B}_i=\dot{\psi },$$

$$\frac{v_i}{v_{A0}}=\frac{\dot{\psi }}{v_{A0}{B}_0}\frac{B_0}{B_i},$$

$$v_{A0}=\frac{B_0}{\sqrt{\mu {\rho }_i}}.$$

$$\psi$$

$$=\frac{2F}{\pi }{\int }_0^{\overline{y_N}}1-\frac{\bar{d}}{({\bar{y}}^2+{\bar{d}}^2)}d\bar{y}$$

$$=\frac{2F}{\pi }\left(\sqrt{\bar{d}-{\bar{d}}^2}-{\tan }^{-1}\frac{\sqrt{1-\bar{d}}}{\sqrt{\bar{d}}}\right).$$

$$v_i=v_0\frac{d_0}{L}.$$

$$S_i=2\frac{v_i{B}_i^2}{\mu }L$$

$$L=\frac{\eta v_{Ai}}{v_i^2},$$

$$\frac{L}{d_0}=\frac{1}{R_{m0}}\frac{B_i}{B_0}\frac{v_{A0}^2}{v_i^2},$$

$$\frac{L_{SP}^2}{d_0^2}=R_{m0}\frac{v_0^2}{v_{A0}^2}\frac{1}{\sqrt{d_0/d-1}}.$$

$$S_{i_{SP}}=2\frac{v_0{B}_0^2}{\mu }{d}_0\sqrt{d_0/d-1}{R}_{m0}{M}_{A0}^2,$$

Full text

Sun: coronal heating – Sun: magnetic reconnection – Sun: activity Sun: Magnetic fields –Magnetohydrodynamics (MHD) –methods: numerical

1 Introduction

The emergence of new magnetic flux from below the photosphere (Harvey & Martin, 1973) and its reconnection with the overlying magnetic field has long been recognised as being one way of heating part of the solar corona, namely, X-ray bright points (Golub et al., 1974), and of heating small flares (Heyvaerts et al., 1977). It has also been proposed as a possible source of coronal X-ray jets (Shibata et al., 1992), for which there has been a host of observational papers (e.g., Moore et al., 2010; Shimojo & Shibata, 2000) and numerical experiments (e.g., Yokoyama & Shibata, 1996; Archontis & Hood, 2010; Moreno-Insertis & Galsgaard, 2013; Syntelis et al., 2015). Indeed, it is now appreciated that reconnection can produce a mixture of hot and cold structures and that their origin can be highly subtle and complex (e.g. Hansteen et al., 2017; Nóbrega-Siverio et al., 2017, 2018).

The cancellation of photospheric magnetic flux is another common process (Martin et al., 1985) which has been proposed as a mechanism for heating X-ray bright points (Priest et al., 1994; Parnell & Priest, 1995), in which magnetic reconnection is driven in the overlying atmosphere during the approach of opposite-polarity magnetic fragments before they actually cancel. Indeed, we shall include this pre-cancellation phase in our use of the words “flux cancellation”. Photospheric flux cancellation has been shown to be associated with both hot and cool jets and also with many different examples of small-scale energy release, such as Ellerman bombs, UV bursts and IRIS bombs (e.g. Watanabe et al., 2011; Vissers et al., 2013; Peter et al., 2014; Kim et al., 2015; Vissers et al., 2015; Rutten et al., 2015; Rezaei & Beck, 2015; Nelson et al., 2016; Tian et al., 2016; Reid et al., 2016; Rutten, 2016; Nelson et al., 2017; Toriumi et al., 2017; Hong et al., 2017; Libbrecht et al., 2017; van der Voort et al., 2017).

A new achievement is the remarkable observations from the Sunrise balloon mission (Solanki et al., 2010, 2017), which have revealed images of the photospheric magnetic field at a spatial resolution of 0.15 arcsec, which is six times better than the Helioseismic Imager (HMI) on the Solar Dynamics Observatory (SDO). In particular, they show that magnetic flux is emerging and cancelling at a rate of 1100 Mx cm*-2* day*-1* (Smitha et al., 2017), which is an order of magnitude higher than previously realised. Furthermore, whereas, at the spatial resolution of HMI, coronal loops have their footpoints located in regions of uniform polarity, at Sunrise resolution the footpoints are surprisingly revealed to have mixed polarity that is cancelling at a rate of 1015 Mx s*-1* (Chitta et al., 2017b). Other examples of flux cancellation producing coronal loop brightening have been presented by Tiwari et al. (2014), Huang et al. (2018), and Chitta et al. (2018)

These observations led Priest et al. (2018) to propose reconnection driven by photospheric flux cancellation as a mechanism for heating the chromosphere and corona. They set up an analytical model for the approach and cancellation of two opposite-polarity magnetic fragments of flux $\pm F$ in the photosphere in the presence of an overlying uniform horizontal magnetic field $B_{0}$, and found that the evolution of the system depends on the value of a key parameter, called the interaction distance, which, for three-dimensional sources, may be written as

[TABLE]

Suppose the magnetic flux sources are separated by a distance $2d$. Then, when $d>d_{0}^{3D}$, the sources are not connected magnetically, but when $d=d_{0}^{3D}$ a null point (or in 3D a separator) forms in the photosphere. As the sources approach closer, such that $d<d_{0}^{3D}$, reconnection is driven and the reconnection location rises in the atmosphere to a maximum height proportional to $d_{0}^{3D}$. Thereafter, the reconnection location moves back towards the solar surface, which it reaches when the two sources come into contact and cancel ($d=0$). Thus, the maximum reconnection height can be located in the photosphere or chromosphere if $d_{0}^{3D}$ is small enough or in the corona if it is large enough.

As well as calculate the way the reconnection height depends on flux ($F$) and overlying field strength ($B_{0}$) through $d_{0}^{3D}$, Priest et al. (2018) made estimates for the energy release, and found that, for reasonable values of the parameters, the heating rate is sufficient to heat the chromosphere and corona.

In the present paper, we develop the model further by setting up a two-dimensional computational experiment for flux cancellation that has the same features as our analytical model, namely, two approaching flux sources in the presence of an overlying horizontal magnetic field, so that we can test the predictions of the analytical model. However, we add an extra feature, namely, a simple stratified atmosphere in order to understand some of the effects of stratification.

Section 2 presents some more details of the theory of reconnection in two dimensions, including Sweet-Parker reconnection, fast reconnection, and energy conversion. Then, Section 3 presents our computational model and compares it with the analytical theory, before a summary discussion is given in the final Section.

2 Theory for Energy Release driven by Photospheric Flux Cancellation in 2D

Here we make some theoretical estimates of the energy release by steady-state magnetic reconnection in two dimensions, developing the basic theory from Priest (2014) in new ways. We will start by briefly describing slow Sweet-Parker reconnection and fast reconnection, and then discussing reconnection driven by magnetic flux cancellation.

2.1 Slow Sweet-Parker Reconnection

Consider first a simple Sweet-Parker current sheet of given length $L$, depth $L_{s}$ and width $l$ situated between oppositely directed magnetic fields $B_{i}$ and $-B_{i}$ (Fig. 1a). If plasma and magnetic field is brought in from both sides at a speed $v_{i}$, then a balance between inwards advection and outwards diffusion of magnetic flux implies

[TABLE]

Furthermore, if the plasma has uniform density $\rho_{i}$, balancing the rates of inflow and outflow of mass gives

[TABLE]

where $v_{Ai}=B_{i}/\sqrt{\mu\rho_{i}}$ is the outflow speed from the current sheet, namely, the Alfvén speed based on the inflow magnetic field.

Eliminating $l$ between Equation (2) and (3) produces an expression for the dimensionless inflow speed or Alfvén Mach number ($M_{Ai}=v_{i}/v_{Ai}$), i.e., the reconnection rate, of

[TABLE]

where $R_{mi}=Lv_{Ai}/\eta$ is the external magnetic Reynolds number based on the global external length-scale ($L$) and Alfvén speed ($v_{Ai}$). For a given current-sheet length ($L$) and external magnetic field ($B_{i}$), Equation (4) thus provides the Sweet-Parker rate.

Half of the magnetic energy that comes into the reconnection region from both sides is converted into heat and the other half into kinetic energy (which can later itself dissipate viscously or through shock waves). The rate of inflow of magnetic energy from one side through an area of $LL_{s}$ is just the Poynting influx $(EH_{i}LL_{s}=EB_{i}LL_{s}/\mu)$, where the magnitude of the electric field is $E=v_{i}B_{i}$, and so the rate of conversion to heat of magnetic energy coming in from both sides of the current sheet is

[TABLE]

after substituting for $v_{i}$ from Equation (4). The phrase “two-dimensional” can refer to a situation in which the variables are situated in three dimensions but they depend on only two of them, such as $x$ and $y$, but it can also refer to variables that exist only in two dimensions, in which case the above expression would need to be divided by $L_{s}$. In what follows it should be clear which of the two definitions is being inferred.

2.2 Fast Reconnection

Next, suppose the inflow speed is faster than $v_{Ai}/R_{me}^{1/2}$, while the inflow magnetic field ($B_{i}$) and area $L_{s}$ are the same as before (Fig. 1b). Then, three possibilities have been studied. Firstly, according to fast steady-state reconnection theory (either Petschek (1964) or Almost-Uniform (Priest & Forbes, 1986)), the reconnection region possesses a complex internal structure consisting of a central small Sweet-Parker current sheet together with four slow-mode shock waves propagating from their ends and standing in the flow. As the speed increases, the central sheet diminishes in size, while the length and inclination of the shock waves increases. Most of the energy conversion then takes place at the shock waves, with $\frac{2}{5}$ of the inflowing magnetic energy being converted to heat (rather than the $\textstyle{\frac{1}{2}}$ that is found in Sweet-Parker reconnection) and the remainder going to kinetic energy.

Secondly, fast collisionless reconnection is helped by the Hall effect, when the resistive diffusion region is replaced by an ion diffusion region and a smaller electron diffusion region. In this case, a similar fast maximum rate of reconnection as in Petschek’s mechanism results (Shay & Drake, 1998; Birn et al., 2001; Huba, 2003; Huba & Rudakov, 2004; Shay et al., 2007; Birn & Priest, 2007).

Thirdly, when the central sheet is long enough, it goes unstable to secondary tearing mode instability and a regime of impulsive bursty reconnection results, first described by Priest (1986); Lee & Fu (1986); Biskamp (1986); Forbes & Priest (1987) and later studied by Loureiro et al. (2007, 2012, 2013); Bhattacharjee et al. (2009). Reconnection is then fast but time-dependent and impulsive, although the mean rate is likely to be similar to the previous cases.

For each of the three above scenarios, Equation (2) no longer holds, but the same mass conservation relation holds as before for the reconnection region as a whole, namely,

[TABLE]

where $L$ now refers to the length of the whole reconnection region (including shock waves and central current sheet) rather than the length of just the central sheet, and variables with subscript $i$ refer to values at the inflow to that whole region. Equation (6) determines the overall width ($l$) of the complex reconnection region for a given $L$, $v_{i}$ and $v_{Ai}$. The conversion rate of inflowing energy from both sides of the current sheet then becomes

[TABLE]

where $v_{i}$ possesses any value up to a maximum of typically 0.01–0.1 of the Alfvén speed ($v_{Ai}$).

2.3 Energy Conversion during Photospheric Flux Cancellation in 2D

2.3.1 Magnetic configuration

Consider sources of positive and negative photospheric magnetic flux ($\pm F$) situated at points B $(d,0)$ and A $(-d,0)$ on the $x$-axis in a region of uniform magnetic field $B_{0}\bf{\hat{x}}$, and suppose they approach one another at speeds $\pm v_{0}$.

The resulting magnetic field (in two dimensions) above the photosphere ($y>0$) is given by

[TABLE]

where

[TABLE]

are the vector distances from the two sources to a point P($x,y$).

It is natural in Equation (8) to non-dimensionalise the magnetic field with respect to $B_{0}$ and distances with respect to the 2D version of the interaction distance (Longcope, 1998), namely,

[TABLE]

and so define

[TABLE]

Then the magnetic field on the $y$-axis becomes

[TABLE]

Consider what happens as the two sources approach each other. The evolution of the topology is similar to what happens in three dimensions, as described in detail in Sec. 2.1 of Priest et al. (2018). When the two sources are far away ($d>d_{0}$), they are not connected magnetically and two first-order null points lie on the $x$-axis between the sources (Fig. 2a). When $d=d_{0}$, a high-order null point appears at the origin (Fig. 2b). As the sources approach one another ($d<d_{0}$), the null point rises above the photosphere (Fig. 2c). The location of the null point at $y=y_{N}$ is given by

[TABLE]

and it rises along the $y$-axis to a maximum of $y=\textstyle{\frac{1}{2}}d_{0}$ when $d=\textstyle{\frac{1}{2}}d_{0}$. Then, it falls, reaching the origin when $d=0$, as shown in Fig. 3a. When $d=0$, the flux of the two sources has completely cancelled.

2.3.2 Inflowing Plasma Speed ($v_{i}$) and Magnetic Field ($B_{i}$) at the Reconnection Site

To analyze the energy release during flux cancellation, the natural parameters, for each value of the source separation ($2d$), are the magnetic diffusivity ($\eta$), the critical source half-separation distance ($d_{0}$), the flux source speed ($v_{0}={\dot{d}}=dd/dt$) and the overlying field strength ($B_{0}$). We now proceed to calculate the inflow speed ($v_{i}$) and magnetic field ($B_{i}$) to the current sheet and the sheet length ($L$) as functions of these parameters in the cases of slow reconnection and fast reconnection. The magnetic configuration driven by flux cancellation is shown in Fig. 2d.

First we consider $B_{i}$. The components of the potential magnetic field sufficiently close to a 2D X-point can be written as (Priest, 2014)

[TABLE]

where $k$ is a constant and $z=x+iy$ is the complex variable. Suppose that the configuration with a reconnecting current sheet of length $L$ is represented by

[TABLE]

such that the sheet is a cut in the complex plane between $z=\pm{\textstyle{\frac{1}{2}}}iL$. Then, putting $z=0+$ implies that

[TABLE]

which is the required expression for $B_{i}$ when the $x$-component of the field in the potential state near the null has the form $B_{x}=ky$. The value of $k$ is calculated as follows. The horizontal field $B_{x}$ near $y=y_{N}$ may be obtained by putting $y=y_{N}(1+\epsilon)$, where $\epsilon\ll 1$ in Equation (10). Keeping only the linear terms, this equation gives

[TABLE]

or

[TABLE]

This determines the value of $k$, and so our required expression becomes

[TABLE]

Next, consider $v_{i}$. This may be calculated from the rate of change (${\dot{\psi}}\equiv d\psi/dt$) of magnetic flux, since

[TABLE]

or, in dimensionless form

[TABLE]

where $v_{A0}$ a hybrid Alfvén speed based on the magnetic field $B_{0}$ and the density of the inflowing material, namely,

[TABLE]

In turn, ${\dot{\psi}}$ may be calculated from the reconnected flux ($\psi$), as estimated from the magnetic flux below the null point, namely,

[TABLE]

It can be seen from Fig. 3b that, as expected, the reconnected flux vanishes when $d=d_{0}$ and increases monotonically to a value of $F$ as the separation ($2d$) between the sources approaches zero.

Then, differentiating Equation (17) with respect to time determines ${\dot{\psi}}$ in terms of ${\dot{d}}=v_{0}$, and Equation (15) becomes

[TABLE]

2.3.3 Energy Release

The rate of inflow of magnetic energy from one side of a current sheet of length $L$, at speed $v_{i}$, with field strength $B_{i}$ and density $\rho_{i}$ is the Poynting flux through that surface. In 2D the surface of the current sheet will be a line of length $L$, so the Poynting influx is $EH_{i}L=EB_{i}L/\mu$. Since the electric field is $E=v_{i}B_{i}$, and the magnetic energy inflow occurs from both sides of the current sheet, the Poynting influx from both sides will be

[TABLE]

This has units of energy/time/length, since we assume here a purely 2D configuration with no depth in the third dimension. To derive the energy release, the length of the current sheet and the conversion rate to heat has to be estimated. Both will depend on the type of reconnection (Sweet-Parker or fast). For a configuration with depth $L_{S}$ in the third dimension this would be multiplied by $L_{S}$.

2.3.4 Slow Sweet-Parker Reconnection

Here we calculate the energy release for Sweet-Parker reconnection. After eliminating $l$ between the Sweet-Parker relations (Equations (2) and (3)) we find that the current sheet length is

[TABLE]

which can be nondimensionalised in terms of $d_{0}$ to give

[TABLE]

where $R_{m0}=d_{0}v_{A0}/\eta$ is the magnetic Reynolds number based on $d_{0}$ and $v_{A0}$. We then substitute for $B_{i}/B_{0}$ from Equation (13) and $v_{i}$ from Equation (18) to give

[TABLE]

where the subscript $SP$ denotes the Sweet-Parker current sheet length. Finally, by substituting in Equation (19) the values of $v_{i}$ (Equation (18)), $B_{i}$ (Equation (13)) and $L=L_{SP}$ (Equation (21)), the rate of Poynting influx becomes

[TABLE]

in terms of the Alfvén Mach number ($M_{A0}=v_{0}/v_{A0}$) based on the flux source speed $v_{0}$. Since half of the magnetic energy is converted to heat during Sweet-Parker reconnection, the energy release rate will be:

[TABLE]

which has units of energy/time/length for our 2D theory.

2.4 Fast Reconnection

We derive now the energy release for fast reconnection driven by flux cancellation. During fast reconnection, the length of the current sheet is much smaller than the Sweet-Parker one. $L$ is determined by assuming the inflow speed $v_{i}=\alpha v_{Ai}$. By writing $v_{i}=\alpha v_{Ai}=\alpha v_{A0}B_{i}/B_{0}$ and then using Equation (13) and (18), $L$ becomes

[TABLE]

Then, after substituting for $v_{i}$, $B_{i}$ and $L$ into Equation (19), we find the rate of energy inflow for fast reconnection as

[TABLE]

During fast reconnection, $\frac{2}{5}$ of the magnetic energy is converted to heat and $\frac{3}{5}$ to kinetic energy. Therefore, the rates of kinetic energy release and energy release as heat become

[TABLE]

and

[TABLE]

3 Numerical Computations

3.1 Numerical Setup

To perform the computations, we numerically solve the 2D MHD equations in Cartesian geometry using the Lare3D code (v3.2) of Arber et al. (2001). The equations in dimensionless form are:

[TABLE]

where $\rho$, $\mathbf{v}$, $\mathbf{B}$ and P are density, velocity vector, magnetic field vector and gas pressure. Gravity is $g_{0}=274$ m s*-1*. We assume a perfect gas with specific heat of $\gamma=5/3$. Viscous heating ($Q_{\mathrm{v}}$) and Joule dissipation ($Q_{\mathrm{j}}$) are included. Heat conduction ($Q_{c}$) is treated using super-time stepping (Meyer et al., 2012), similarly to Johnston et al. (2017). The reduced mass is $\mu_{m}=m_{f}m_{p}$, where $m_{p}$ is the mass of proton and $m_{f}=1.2$. $k_{B}$ is the Boltzmann constant.

The normalization is based on the photospheric values of density $\rho_{\mathrm{u}}=1.67\times 10^{-7}\ \mathrm{g}\ \mathrm{cm}^{-3}$, length $H_{\mathrm{u}}=180\ \mathrm{km}$ and magnetic field strength $B_{\mathrm{u}}=300\ \mathrm{G}$. From these we obtain temperature $T_{\mathrm{u}}=6234~{}\mathrm{K}$, pressure $P_{\mathrm{u}}=7.16\times 10^{3}\ \mathrm{erg}\ \mathrm{cm}^{-3}$, velocity $v_{\mathrm{u}}=2.1\ \mathrm{km}\ \mathrm{s}^{-1}$ and time $t_{\mathrm{u}}=86.9\ \mathrm{s}$.

The computational domain has a physical size of $x\in[-30,30]$ Mm in the horizontal direction and $y\in[0,30]$ Mm in the vertical direction, on a $2048\times 1024$ uniform grid. The photosphere is at $y=0$. To mimic the steep temperature increase from the photosphere to the corona, we assume an hyperbolic tangent profile for the atmospheric temperature

[TABLE]

where $T_{ph}=6109$ K, $T_{cor}=0.61$ MK, $y_{cor}=2.12$ Mm and $w_{tr}=0.18$ Mm. These parameters create an isothermal photospheric-chromospheric layer at $0\ \mathrm{Mm}\leq y<1.96\ \mathrm{Mm}$, a transition region at $1.96\ \mathrm{Mm}\leq y<3.3\ \mathrm{Mm}$ and an isothermal coronal at $3.3\ \mathrm{Mm}\leq y<30\ \mathrm{Mm}$.

To derive the atmospheric density, we assume the atmosphere is in hydrostatic equilibrium. We do so by numerically solving the hydrostatic equation $dP/dy=-gy$, assuming a photospheric density of $\rho_{ph}=1.67\times 10^{-7}$ g cm*-3*. The atmospheric temperature (solid black) and density (solid blue) are shown in Fig. 4. For comparison, we plot with dashed lines the temperature and density for the 1D model atmosphere (model C7) of Avrett & Loeser (2008).

We adopt an anomalous resistivity

[TABLE]

where $\eta_{0}=10^{-4}$, $\eta_{1}=10^{-3}$ and $j_{crit}=10^{-3}$. The resistivity can be anomalous away from the boundaries ($x\in[-28,28]$ Mm and $y\in[2,28]$ Mm). Elsewhere, it is uniform with $\eta_{0}=10^{-4}$. Anomalous resistivity (Yokoyama & Shibata, 1994) has been previously chosen to drive fast reconnection. Other methods (e.g. hyper-diffusion Nordlund & Stein, 1990; van Ballegooijen & Cranmer, 2008; Martínez-Sykora et al., 2011) can also be used to initiate a fast reconnection by permitting enhanced resistivity in current sheets.

The initial magnetic field is the sum of two magnetic sources and a horizontal field:

[TABLE]

where

[TABLE]

are the position vectors of the left and right sources, respectively, $d_{s}=1.8$ Mm is the distance of each source from $x=0$, and $y_{0}=-0.36$ Mm is the depth of the sources below the photosphere (the sources are outside the numerical domain). The flux of each source is $F=2.5\times 10^{11}$ Mx cm*-1*. The polarities produced at the photosphere have a maximum field strength of 2.2 kG and a size of about 1 Mm (defined as the length where $|B_{y}|>100$ G) (Fig. 5a). The flux of each polarity is $F_{m}=2.2\times 10^{11}$ Mx cm*-1*. The horizontal field has a strength of $B_{0}=45$ G.

The boundary conditions on the upper boundary are $\mathbf{v}=0$ and zero gradients for $\mathbf{B}$, $\rho$, $\epsilon$. The photospheric boundary conditions are zero gradients for $\rho$, $\epsilon$. The magnetic field at the photospheric boundary changes according to the driver. To drive the cancellation, we move the sources with a velocity of

[TABLE]

where $v_{max}$ = 1 km s${}^{-1}\,$, $t_{0}=10.1$ min and $w=1.4$ min. The positions of the sources changes according to $d(t)=d_{s}-x(t)$, where

[TABLE]

The simulation is driven by changing the magnetic field at the lower boundary ($y=0-$ (ghost cells)) using Equation (36) and $d(t)$. The half-separation ($d(t)$) of the sources (below the photosphere) as a function of time is plotted in Fig. 5b (black line). The blue lines show the positions of the polarities at $y=0$ (found by measuring the location of maximum $B_{y}$). The latter reflects the response of the photosphere to the driver.

For the parametric study of Sec. 3.4, we vary the magnetic field strength ($B_{0}$) of the atmosphere in order to vary the height of the null point. The values of $B_{0}$ and the corresponding null height at $t=0$ min are shown in Table 1.

3.2 Comparison of Theory with Simulation

In this section, we discuss our reconnection experiment driven by flux cancellation and compare its results with the theory presented in Sec. 2. For this, we shall focus on Case 1 of Table 1.

3.2.1 Brief Description of Simulation

The magnetic field at $t=0$ is shown in Fig. 6a, with a null point at $(x,y)=(0,7.6)$ Mm. As the driver is switched on, reconnection is driven at the null point due to the converging photospheric polarities. The energy released by reconnection spreads above and below the null (and shows up as a “horizontal” heated region and an underlying heated arcade in Fig. 6b). The heated material is denser than the background atmosphere (Fig. 6c).

At the photosphere, the positions of the polarities at $y=0$ ($d_{m}(t)$, blue line, Fig. 5b) do not keep following the driver after $t=37$ min (black line). At this time, the magnitude of the photospheric field has decreased to the point that $\beta>1$. As a result, the driver cannot move the overlying field anymore. The reconnection at the null follows the response of the atmospheric field to the driver and gradually stops.

The interaction distance for this simulation is $d_{0}=200$ based on the sources and $d_{0_{m}}=2F_{m}/(\pi B_{0})=173.2$ based on the photospheric polarities.

3.2.2 Comparison Methodology

In our simulation, we set the gradient of $\mathbf{B}$ to be zero at the boundaries (besides the photosphere). The reason for this is as follows. After flux cancellation, if the polarities are completely cancelled, the remaining atmospheric field ought to be a horizontal field with strength $B_{0}$. This cannot happen in a finite numerical domain, but only in a semi-infinite one. To achieve that in the simulation domain, we use a zero gradient boundary condition. This “straightens” the field lines, mimicing the effect we require.

In Sec. 2, the inflow speed is found by taking into account the rate of change of flux and conservation of flux. To compare the simulation with theory, we need to carefully calculate the fluxes inside the numerical domain. The Appendix calculates how much flux should be found inside and outside the finite numerical domain, from which we deduce a flux correcting factor $f$ (Equation A14, Fig. 13). This is used to multiply several quantities ($v_{i}\rightarrow fv_{i}$, $L_{sp}^{2}\rightarrow f^{2}L_{sp}$, $L^{2}\rightarrow fL^{2}$, $dW/dt\rightarrow f^{2}dW/dt$), since the simulation uses a finite domain rather than a semi-infinite one (as on the Sun, for which $f\rightarrow 1$).

In the simulation there is a difference between the driver (the sources below the photosphere) and the response of the photosphere to the driver (the polarities at the photosphere). Our theory uses observables (such as the separation of the polarities and the photospheric velocity) to predict the inflowing magnetic field and the energy release. To compare theory with the computational experiment, we use as “observables” two sets of values:

(i) the values of the driver (such as $d$, $v$), and

(ii) the values measured at the photosphere (or elsewhere), which mimic an actual observation.

We will refer to the latter quantities with a subscript $m$. Thus, the half-separation of the sources is $d$, whereas the half-separation of the photospheric polarities is $d_{m}$.

To compare the theory with the simulation, we first identify the current sheet. It is located along the $y$-axis and is the vertical region of increased temperature located above the apex of the arcade (Fig. 7a, orange line segment). We identify the coordinates of its lowest ($S_{l}$) and highest ($S_{h}$) point throughout the simulation. The separation of these two points is the measured length of the current sheet, $L_{m}$ (solid black line, Fig. 7e).

We measure the inflowing magnetic field strength, velocity and density in the following manner. At both sides of the current sheet, we identify the regions that are parallel to the current sheet and at a distance of $\Delta x=0.2$ Mm away from it (line segments AB and CD, Fig. 7a). This distance is selected so that the current density there is at least an order of magnitude lower than the one inside the current sheet. We measure the values of these inflowing quantities as the average values of their magnitudes along both AB and CD, and so find the average $B_{i_{m}}$, $v_{i_{m}}$ and $\rho_{i_{m}}$.

The total inflow of Poynting flux ($S_{i_{m}}$) into the current sheet is measured by taking into account the Poynting flux along both AC and CD.

3.2.3 Current Sheet Length, Inflowing Magnetic Field and Inflow Velocity

Fig. 7b (solid line) shows the inflow magnetic field strength, $B_{i_{m}}$. The dashed line is the predicted $B_{i}$ using Equation (13) with $(d,d_{0},L_{m})$. The dashed-dot line is $B_{i}$ using Equation (13) with $(d_{m},d_{0_{m}},L_{m})$. Comparing both approaches, we see that the theory is in good agreement with the simulation. Indeed, the second approach, where we take into account only the response of the simulation to the driver, is in better agreement.

The simulation’s inflow velocity ($v_{i_{m}}$) is plotted in Fig. 7c (solid line). The dashed line shows $v_{i}$ using Equation (18) (times $f$) with $(v_{0},d_{0},L_{m})$. The dashed-dot line is $v_{i}$ using Equation (18) (times $f$) with $(v_{0_{m}},d_{0_{m}},L_{m})$. Again, both agree well with the simulation.

Before estimating the length of the current sheet for fast reconnection (Equation (24)), we need to measure two quantities: $\alpha$, which is the Alfvén Mach number of the inflow and $v_{A0}$ (Equation (16)). In Fig. 7d, we plot the Alfvén Mach number of the inflow. Between $t=10$ min and $40$ min, when the cancellation occurs, it has an average value of $\alpha_{m}=0.05$. This value of $\alpha$ is typical for fast reconnection (Priest, 2014). For the hybrid Alfvén speed we use $v_{{A0}_{m}}=B_{0}/\sqrt{\mu\rho_{i_{m}}}$.

The length of the simulation’s current sheet ($L_{m}$) is plotted in Fig. 7e (solid line). The dashed line is $L$ using Equation (24) (times $f^{2}$) with $(\alpha_{m},d,d_{0},v_{0},v_{{A0}_{m}})$. The dashed-dot line is $L$ using Equation (24) (times $f^{2}$) with $(\alpha_{m},d_{m},d_{0_{m}},v_{0_{m}},v_{{A0}_{m}})$. Both approaches show the theory to be in good agreement with the simulation.

In Fig. 7f, we plot the quantities of panel e, and overplot the length of the current sheet assuming Sweet-Parker reconnection (triple-dot dashed line, with $L_{SP}$ calculated from Equation (21) (times $f$) using $(d_{0_{m}},v_{0_{m}},v_{{A0}_{m}})$). The predicted current sheet for an assumption of slow Sweet-Parker reconnection is longer than the simulated one by an order of magnitude, and so we deduce that fast reconnection with an inflow Alfvén Mach speed of 0.05 well describes the simulation.

3.2.4 Energy Release

We study energy release only for fast reconnection and first focus on the Poynting flux inflow ($S_{i_{m}}$), which is plotted in Fig. 8a (solid lines). The dashed curve is $S_{i}$ from Equation (25) (times $f^{2}$) based on ($v_{0}$, $v_{{A0}_{m}}$, $B_{0}$, $d_{0}$, $d$). The dashed-dot curve is $S_{i}$ based on ($v_{0_{m}}$, $v_{{A0}_{m}}$, $B_{0}$, $d_{0_{m}}$, $d_{m}$). Both agree well with the simulation. The approach of using only “measured” values is in excellent agreement with the simulation results.

Next, we consider the conversion of Poynting flux to kinetic and thermal energy, for which we calculate the energy integral terms:

[TABLE]

The curve $C$ is ABCDA in Fig. 7a and the surface $A$ is its area. The term $\int_{C}{\bf E}\times{\bf B}\cdot d{\bf C}=\int_{A}{\bf\nabla}\cdot\mathbf{S}\ dA$, where $\mathbf{S}$ is the Poynting vector, is the rate of energy inflow and can be compared with Equation (25). The term $-\int_{A}\eta\mathbf{j}^{2}dA$ is the rate of energy converted to joule heating during reconnection, which can be compared with Equation (27). The term $\int_{A}\mathbf{j}\cdot\left(\mathbf{v}\times\mathbf{B}\right)dA$ is the rate of energy converted to kinetic energy and can be compared with Equation (26). The term $-\int_{A}[1/(2\mu)]\partial B^{2}/\partial t\ dA$ is negligible.

We first check whether the energy conversion rates of the simulation agree with those of fast reconnection, i.e. whether $\frac{3}{5}$ of the total Poynting influx is converted to kinetic energy and $\frac{2}{5}$ is converted to Joule heating. If the conversion rates are such, then we should find in the simulation that

[TABLE]

We plot these terms in Fig. 8b, from which it can be seen that the left (solid line) and right (dashed line) terms are in agreement. Furthermore, we examine if

[TABLE]

These terms are plotted in Fig. 8c, from which again the left (solid line) and right (dashed line) terms are in agreement. So, indeed the energy release in the simulation agrees with the rates predicted by fast reconnection.

We now compare the energy release from the simulation with the theoretical predictions. The kinetic energy release rate is calculated from Equation (26) based on ($v_{0_{m}}$, $v_{{A0}_{m}}$, $B_{0}$, $d_{0_{m}}$, $d_{m}$) and is plotted in Fig. 8b (dashed-dot line). This is in fact just the dashed-dot line of Fig. 8a multiplied by 0.6. Next, we calculate the total rate of conversion of energy to heat from Equation (26) based on ($v_{0_{m}}$, $v_{{A0}_{m}}$, $B_{0}$, $d_{0_{m}}$, $d_{m}$) and plot it Fig. 8c (dashed-dot line). In both cases, the theoretical predictions are in excellent agreement with the simulation.

3.3 Atmospheric Response

In this section we briefly discuss the atmospheric response to reconnection driven by flux cancellation. First, we study the time evolution of one individual case. Then, we vary the height of the null at $t=0$ by changing the strength of the external horizontal field and studying five cases with different $B_{0}$. The values of $B_{0}$ and the corresponding $y_{N}$ are shown in Table 1. In Fig. 9, we plot the $y_{N}$ values (vertical lines) and the initial temperature stratification (solid line), in order to better visualize the initial location of the null point relative to the corona, transition region, chromosphere and photosphere.

We focus on the time evolution of the temperature and density for case 2 (Fig. 10). The null point is initially located at the base of the corona. As reconnection starts, hot material is ejected along the post-reconnection field lines (panels a1, b1). A hot “loop” of 1.8 MK and density of $2\times 10^{-14}$ g cm*-3* is formed above the reconnection site. Below the null, the top of the arcade is heated to 2.6 MK (panels a2, b2). As the polarities converge the null height decreases, as predicted from the theory. When the null point reaches the base of the transition region and below, dense, cool plasma is ejected along the reconnected field lines (Fig. 10, panels a3, b3). Due to the higher density of the region, the resulting heat released from the reconnection cannot raise the plasma temperature to millions of Kelvin. As the process continues, a cooler and denser ejection is formed, with temperature of $0.05-0.12$ MK, and a density of $2\times 10^{-12}-2\times 10^{-13}$ g cm*-3*. It propagates with velocity up to $105$ km s${}^{-1}\,$, extending from transition region to coronal heights (panels a4, b4).

In Fig. 11 we plot the temperature (first column) and density (second column) for cases $2-5$ at $t=40$ min, while Case 1 is shown in Fig. 6b,c. An important qualitative difference appears between the cases. When the null point is initially located in the corona, both a hot and a cool plasma region develop above the null during the cancellation. When the initial null point is placed at progressively lower heights (top to bottom row), the amount of hot material decreases, while the cool material increases. Eventually, for a null point placed at the chromosphere (bottom row), the resulting post reconnection plasma does not have a high-temperature component. In this case, the region above the null contains

(i) a very cool component of photospheric or chromospheric material (around 6300 K) which is “slingshotted” upwards from the tension of the reconnected lines with speed of 10-20 km s${}^{-1}\,$and

(ii) a cool plasma component which is heated by reconnection to around $0.01-0.03$ MK.

In Fig. 12 we plot the time evolution of the maximum velocity of the hot ($T>1$ MK) and cool ($T<0.2$ MK) plasma components for the cases shown in Fig. 11. The maximum velocities of the hot (cool) plasma ejections from (a) to (d) are 100 km s${}^{-1}\,$(105 km s${}^{-1}\,$), 79 km s${}^{-1}\,$(87 km s${}^{-1}\,$), 35 km s${}^{-1}\,$(70 km s${}^{-1}\,$) and 0 km s${}^{-1}\,$(57 km s${}^{-1}\,$) respectively. Notice also that the hot and cool components are produced with a time delay, as shown previously in Fig. 10. The time difference between the acceleration of the hot and cold material decreases as the null point is situated lower. In cases 2 and 3, the hot material appears first and the cold material after. For case 4, the hot and cold ejections are almost co-temporal.

3.4 Atmospheric Response

4 Discussion

In Priest et al. (2018), we proposed magnetic reconnection driven by photospheric flux cancellation as a mechanism for energizing coronal loops and heating the chromosphere. We also derived analytical expressions that predict the energy release during reconnection. In the present work, we begin to numerically validate our theory by developing the theory in 2D and comparing it with computations of two converging polarities inside a stratified atmosphere containing a background horizontal field. As the polarities converge, reconnection is driven at the null point.

To compare the theory with simulations, we evaluated several quantities from the simulations. For example, we calculated the velocity of approach of the opposite polarities in two ways. One was to use the values that correspond to the simulation’s driver and the other was to measure the response of the photosphere to that driver.

We found excellent agreement between theory and simulation, especially with the second approach. The response to the driver is to initiate motions in the numerical domain which lead to reconnection. It is found that the theory agrees well with the system’s response to the driver, which is encouraging, since it shows that that our theory could be used to derive estimates of the energy released during flux cancellation from solar observations, as observations measure the photospheric and atmospheric response, without knowledge of the sub-photospheric conditions driving the cancellation. We conclude, based on our 2D computational experiments, that the energy released during photospheric cancellation can be accurately estimated from a knowledge of the converging velocity, the separation and strengths of the converging fluxes, the strength of the background magnetic field, and the density and Alfvén Mach number of the material flowing into the current sheet.

The promising results from these 2D simulations suggest that our analytical estimates can indeed be used to predict energy release. Ideal observational candidates for such a comparison in future include many cases where photospheric flux cancellation is associated with small-scale energy release, such as Ellerman bombs, UV bursts and IRIS bombs or the energy injected into coronal loops due to flux cancellation at their feet.

We have also presented an initial study of the atmospheric response to reconnection (for more sophisticated and realistic simulations, see, e.g., Danilovic et al. (2017); Hansteen et al. (2017); Nóbrega-Siverio et al. (2018)). The maximum height of the null point in 2D is $\textstyle{\frac{1}{2}}d_{0}$. As the polarities converge, the null point moves up to its maximum height and then down towards the photosphere. The atmospheric response during photospheric cancellation is as follows. When the null point is located initially at a coronal height, a hot “loop” (around $1-2$ MK) can be formed above the reconnection region. Cooler material is ejected along the reconnected field lines when the null point is located at the base of the transition region or below, with a velocity of 95 km s*-1* and a temperature of $0.05-0.1$ MK. These ejections occur with a time difference, which is smaller when the initial null point height is lower. However, if the null point is initially located below the base of the transition region, no hot material is ejected, and we only find the formation of a cooler ejection. Thus, the location of energy release is crucial for the type of plasma structure that is created. The hot structures that we find have temperatures and densities similar to those of a coronal loop, whereas the cool structures have values that are reminiscent of surges or larger spicules.

Note that, if only part of the photospheric flux cancels, the null point stops moving towards the photosphere at some intermediate height. Then, reconnection occurs only between the initial and final height of the null, which produces a shorter, less energetic burst of energy release.

There is an important difference between 2D and 3D simulations. In 2D, the magnetic field of a source falls off with distance like $1/r$, whereas in 3D, the fall-off behaves like $1/r^{2}$. As a result, the interaction distance in 2D ($d_{0}=2F/(\pi B_{0})$) is larger than its 3D value ($d_{0}=\sqrt{F/(\pi B_{0})})$ (see Priest et al., 2018), which produces a higher location for the null point in 2D than in 3D for a given polarity separation distance, photospheric flux and background field. Thus, in our 2D simulations, in order to place the null at a particular height in the stratified atmosphere, we adopt a stronger background field than would be needed in 3D. During the reconnection, this higher background field produces a larger Poynting influx into the current sheet in 2D than in 3D. The result is that more energy is converted into heat and kinetic energy. Therefore, we leave a detailed discussion of temperature and density distributions along reconnected field lines for a future 3D experiment. In 3D, the energy release may well accelerate the cooler plasma to form shorter structures than in 2D.

In our model we have assumed a horizontal external field in order to be able to make a direct comparison with our analytical theory. An oblique external field would have several extra effects. Firstly, it would enhance plasma draining along field lines, changing the maximum length and density of the heated plasma structures, an effect that would be stronger for the cooler ejections. Secondly, the null point would move sideways, as well as vertically. This could affect the width of the structures and possibly produce “thread-like” ejections. Thirdly, after reconnection, the flows above the null point will have up and down components, instead of being mainly horizontal. The resulting magnetic “loop” would have its footpoints rooted in the photosphere and the ejection would be dominated by a single inclined upflow (together with a much shorter downflow), rather than consisting of two opposite directly horizontal flows (e.g., Fig. 10). Jet-like structures have been observed at the feet of coronal loops (Chitta et al., 2017a, b) which could be related to the upflows we expect in the oblique field. However, we do not expect the energy release to change drastically. In Sec. 2, we derived the rate of heating by assuming it is half the total inflow of Poynting flux. For an oblique field, the flux function $\psi$ would be different, but, during flux cancellation, the same amount of flux will be cancelled, irrespective of the orientation of the external field. Consequently, the same Poynting influx into the current sheet would be produced over the same time scale. Thus, the energy release rate should not be significantly different. We aim to check this numerically in future.

In this work, we have positively validated our analytical theory using 2D numerical computations. This suggests that nanoflares driven by magnetic flux cancellation can indeed be an important mechanism for heating the chromosphere and corona, as proposed in Priest et al. (2018), which is built upon recent observational findings. In future, we aim to extend our model in several ways to make it more realistic and to consider more cases. In particular, we shall study oblique external fields in order to determine in more detail the ways in which chromospheric and coronal loops may be heated by reconnection at their footpoints. We shall also set up a fully three-dimensional computation in order to study the extent and implications of our theory and to deduce in more detail the atmospheric response to energy release.

L.P.C. received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 707837. This research has made use of NASA’s Astrophysics Data System. The authors are most grateful for invaluable discussions with Hardi Peter, Clare Parnell and Alan Hood.

Appendix A Flux Correction Factor

At $x=0$, the flux contained between the heights $\bar{y}_{1}$ and $\bar{y}_{2}$ is $\psi=\int_{\bar{y}_{1}}^{\bar{y}_{2}}\bar{B}_{x}d\bar{y}$, which for our magnetic field (Equation (36)) becomes:

[TABLE]

Thus, the total flux from the sources to the upper boundary ($y_{max}$) is

[TABLE]

which depends on $\bar{d}$ and $\bar{y}_{max}$. For a semi-infinite domain, such as considered in Sec. 2, $\bar{y}_{max}\to\infty$ and so the flux becomes $F$, which is independent of $\bar{d}$. However, the simulation has a finite region, and so the dependence of flux on $\bar{d}$ and $\bar{y}_{max}$ has to be taken into account in order to compare with theory.

Consider the fluxes above and below the null point. The flux from the sources ($\bar{y}_{0}$) to the null ($\bar{y}^{\prime}_{N}=\sqrt{\bar{d}-\bar{d}^{2}}+\bar{y}_{0})$ is:

[TABLE]

whereas the flux from the null point to the upper boundary of the numerical domain is

[TABLE]

As the sources cancel (and $\bar{d}$ decreases from 1 to 0), the flux below the null point changes from 0 to $F$, resulting in a total cancelled flux of

[TABLE]

The flux above the null up to $\bar{y}_{max}$ changes by

[TABLE]

which becomes $\Delta\psi_{\bar{y}_{N}^{\prime}}^{\infty}=-F$ as $\bar{y}_{max}\to\infty$. Therefore, in a semi-infinite domain, when moving the sources from $\bar{d}=1$ to $\bar{d}=0$, there is flux balance between the fluxes below and above the null. However, for a finite $\bar{y}_{max}$, there is extra flux above $\bar{y}>\bar{y}_{max}$ which we do not take into account. Thus, in a finite domain, it is not possible to fully cancel the two magnetic sources, to give a configuration with a uniform horizontal field, since a flux of $|F|$ would be cancelled below the null while adding less than $-|F|$ above the null.

In the simulation, the rates of change of flux from the sources up to the null and from the null up to $\bar{y}_{max}$ are

[TABLE]

and

[TABLE]

where

[TABLE]

The flux outside the finite domain changes at a rate $\dot{\psi}_{\bar{y}_{max}}^{\infty}=-\dot{\psi}_{f}$, and the rate of change of flux added to the region above the null is $-|v_{i}b_{i}|$. Therefore, below the null, from Equation (A8), the rate of change has to be:

[TABLE]

For a semi-infinite domain, $\dot{\psi}_{f}\to 0$, and therefore $\dot{\psi}_{\bar{y}_{N}^{\prime}}^{\bar{y}_{max}}=-\dot{\psi}_{\bar{y}_{0}}^{\bar{y}_{N}^{\prime}}$.

This affects the theory in the following way. In Sec. 2, the inflow speed was found using the conservation of flux and $B_{i}$:

[TABLE]

To compare theory with simulation, we must use $\dot{\psi}$ from Equation (A10) to give

[TABLE]

or

[TABLE]

where

[TABLE]

$f$ is a flux correction factor, which is plotted in Fig. 13 and which modifies several of the previous expressions, namely, changing $L_{sp}^{2}\rightarrow f^{2}L_{sp}$, $L^{2}\rightarrow fL^{2}$, ${dW}/{dt}\rightarrow f^{2}{dW}/{dt}$. For a semi-infinite domain ($y_{max}\to\infty$), $f\rightarrow 1$ and we recover the theory of Sec. 2 with $v_{i}\rightarrow v_{0}{d_{0}}/{L}$.