Estimating the rate of field line braiding in the solar corona by photospheric flows

In this paper we seek to understand the timescale on which the photospheric motions on the Sun braid coronal magnetic field lines. This is a crucial ingredient for determining the viability of the braiding mechanism for explaining the high temperatures observed in the corona. We study the topological complexity induced in the coronal magnetic field, primarily using plasma motions extracted from magneto-convection simulations. This topological complexity is quantified using the field line winding, finite time topological entropy and passive scalar mixing. With these measures we contrast mixing efficiencies of the magneto-convection simulation, a benchmark flow known as a ``blinking vortex', and finally photospheric flows inferred from sequences of observed magnetograms using local correlation tracking. While the highly resolved magneto-convection simulations induce a strong degree of field line winding and finite time topological entropy, the values obtained from the observations from the plage region are around an order of magnitude smaller. This behavior is carried over to the finite time topological entropy. Nevertheless, the results suggest that the photospheric motions induce complex tangling of the coronal field on a timescale of hours.


Introduction
To understand the dynamics of solar plasmas, it has become evident that we need to study the topology of the magnetic field lines, particularly in relation to reconnection and heating (e.g., Parker 1983;Greene 1988;Lau & Finn 1990;Pevtsov et al. 1996;Craig & Sneyd 2005;Pontin et al. 2011Pontin et al. , 2016Guo et al. 2013;Low 2013;Park et al. 2013;Čemeljić & Huang 2014;Brookhart et al. 2015;Russell et al. 2015;Sun et al. 2016). For topologically nontrivial magnetic field constructions in the form of knots (Rañada & Trueba 1995;Candelaresi & Brandenburg 2011;Smiet et al. 2015), braids (Yeates et al. 2010;Wilmot-Smith et al. 2011), andlinks (Del Sordo et al. 2010), we know that the magnetic helicity, which is a manifestation of linkage, knottedness, and braiding (Moffatt 1969), imposes restrictions on the evolution of the field (Chandrasekhar & Woltjer 1958;Woltjer 1958;Arnold 1974;Taylor 1974;Ricca 2008. This is particularly pronounced in the solar atmosphere, where, due to the high magnetic Reynolds number, the helicity is conserved on dynamical timescales. Braiding or twisting of the field lines can then play a critical role in the dynamics: this braiding/ twisting can be induced either below the solar surface before the flux emerges into the atmosphere, or after the flux emerges in response to photospheric flows. Explaining the huge temperature increase from the solar surface to the corona remains one of the most enduring problems in solar physics (the "coronal heating problem"). The energy that must be supplied to the corona to maintain its million-degree temperature can be estimated by quantifying energy losses (Withbroe & Noyes 1977). A number of previous works have sought to quantify the energy injected into the coronal field by examining photospheric flows. Initial estimates were based on assessing "typical" motions of a twisting or braiding type (Berger 1993;Zirker & Cleveland 1993). More recently, Poynting flux estimates have been made based on velocity fields extracted from observations. Yeates et al. (2012) and Welsch (2015) used Fourier Local Correlation Tracking (FLCT) to estimate photospheric velocities in a plage region. Using these velocities, they obtained a Poynting flux into the corona of around 5×10 4 W m −2 . Shelyag et al. (2012) used magneto-convection simulations to study the vertical Poynting flux, identifying a dominant contribution from plasma motions in strong intergranular magnetic flux concentrations.
These above studies can give lower bounds on the energy injected into the corona, but do not provide clues to the mechanism by which this energy might be dissipated. Here we aim to assess specifically the braiding mechanism proposed by Parker (1972). A key ingredient of Parker's hypothesis is that once the magnetic field topology becomes "sufficiently complex" (i.e., the field lines become "sufficiently tangled"), current sheets spontaneously develop leading to a rapid conversion of magnetic energy to thermal energy. This hypothesis was recently put on firmer footing by Pontin & Hornig (2015), who proved that any force-free equilibrium with tangled magnetic field lines must contain thin current layers (the argument was extended to include fields close to force-free by Pontin et al. 2016). Nevertheless, the efficiency of the braiding mechanism for heating the coronal plasma relies crucially on the-so far unknown-timescales for the energy build up (the field line tangling) and energy release processes.
In this paper, we study the efficiency by which the coronal field lines are tangled/braided in a topological sense by photospheric motions. We principally address this by examining flows derived from an MHD simulation of solar convection. These are then contrasted with the results from flows inferred 1 from observations. We quantify the field line tangling using various measures, including the passive scalar spectrum, the finite time topological entropy (FTTE), and the winding number (e.g., Prior & Berger 2012;Mangalam & Prasad 2018). With these, we aim to address the questions: "How efficient are the motions on the photosphere in inducing twisting?" and "What time is necessary for the Sun to induce a complex entangled braid?" In the following, we will first explain the details of our benchmark flow (Section 2) and the numerical data from the convection simulation we use (Section 3), before we explain how we preprocess these data (Section 4) and analyze them in Section 5. We then go into the details of measuring entangling by using the winding number and the FTTE in Section 6 and conclude by comparing the simulation results with calculations using observational data in Section 7.

Benchmark Flow: The Blinking Vortex
Below we will discuss the tangling induced by the flows extracted from simulations that mimic the convective layer of the Sun. In order to quantify the efficiency of that tangling, we make use of a well-studied benchmark flow, termed a "blinking vortex" (Aref 1984). This pattern of opposing, overlapping twists (see Figure 1) has been shown to constitute a maximally efficient protocol in terms of generating topological entropy (a measure for mixing) per twist operation (Thiffeault & Finn 2006). The number of twist operations can be related to a gross winding angle in a general random flow. Thus, the blinking vortex flow is efficient in generating small scales. The flow is a standard example of a highly mixing flow in fluid dynamics (see Aref 1984;Boyland et al. 2000;Thiffeault & Finn 2006 and references therein). This motivates its choice as a benchmark. When imposed as boundary footpoint motions, this flow generates a braided magnetic field with a complex field line mapping, which has been studied in a series of previous works (e.g., Wilmot-Smith et al. 2009a, 2009bYeates et al. 2010;Pontin et al. 2011;Russell et al. 2015;Ritchie et al. 2016). A subset of the field lines form a "pig-tail" braid and, due to the equal but opposite alternating vortices, the net twist is zero.
This flow comprises alternating positive and negative vortices whose locations are offset such that they only partially overlap (see Figure 1). This motion is mostly confined within the domain [−4, 4]×[−3, 3] and alternates every eight time units by smoothly switching them on and off (all units being non-dimensional). Performing this motion for times between 0 and 48 results in highly tangled fluid particle trajectories. For additional detail about this flow see, e.g., Candelaresi et al. (2017), Equation (11).

Magneto-convection Simulations
We consider three-dimensional MHD simulations of convectively driven dynamos in a Cartesian domain (see Bushby et al. 2012;Bushby & Favier 2014, where more details can be found). Using the depth of the convective layer as the characteristic length scale in the system, the dimensions of the domain are given by   x 0 10,   y 0 10 and   z 0 1, with z=0 corresponding to the upper boundary. All quantities are assumed to be periodic in the horizontal (x and y) directions. The layer is heated from below and cooled from above, with the temperature fixed at each boundary. These upper and lower bounding surfaces are also assumed to be impermeable and stress-free, while the magnetic field is constrained to be vertical at z=0 and z=1.
The particular simulation that is considered here is a fully nonlinear dynamo calculation, with kinetic Reynolds number Re≈150 and magnetic Reynolds number » Re 400 M . We focus upon the time-evolution of the velocity u and vertical magnetic field B z at the top boundary (z = 0), extracting these quantities from 60 snapshots, each of which has a horizontal resolution of 512 2 . The time cadence of these snapshots is ca. 0.61 in dimensionless units. One time unit corresponds to the time taken for an isothermal sound wave to travel a distance of one unit across the surface of the domain. However, it is perhaps more informative to note that this time cadence is approximately one-fifth of the convective turnover time (which, based upon the rms velocity and the depth of the domain, is ca. 3 in these dimensionless units, see Bushby et al. 2012). This turnover time is also comparable to the time taken for the (subsonic) surface flows to travel the width of a typical convective cell. This "granular" timescale, which is arguably the most appropriate to consider when analyzing surface flows, is the temporal normalization that is typically adopted below.
While the underlying simulations are fully three-dimensional, we only make use of the horizontal velocities at the top boundary of the domain (which plays the role of the photosphere in the model). This invariably leads to apparent enhanced compression in areas of down-flows and expansions in areas of up-flows. Such compression effects lead to computational difficulties when we come to address field line tangling. To circumvent these difficulties, we decompose the velocity into divergent and rotational (incompressible) parts, and make use of the latter for our calculations (see Section 4).

Helmholtz-Hodge Decomposition
For the horizontal velocity from the simulations (and later for flows inferred from observations), we perform a Helmholtz-Hodge decomposition. The two-dimensional velocity is then written as the sum of three orthogonal terms: which are the incompressible (solenoidal), compressible, and harmonic terms, respectively. Here, orthogonal means ò · with the differentiable scalar fields ψ, f, and χ. The velocities u i and u c satisfy the boundary conditions = · u n 0 i and = · u n 0 c , where n is the normal vector to the boundary, while = · · u n u n h is the corresponding boundary condition on the harmonic term. These boundary conditions ensure that the decomposition is unique.
We compute the three flow components by solving the Poisson equation for the scalar fields while respecting the given boundary conditions: For the calculations described below of the induced topological complexity, we use the incompressible part u i only. With these boundary conditions, the use of this flow component ensures that no particle trajectories are lost across the boundaries of the domain. It is shown later that the contribution to the field line tangling of the flow component u c is insignificant.

Qualitative Mixing Patterns: Passive Scalar
A surface motion like the one simulated in turbulent convection leads to a mixing of the fluid. Any features that are initially stretching over large length scales will then be mapped into small scales, and vice versa, at a rate that depends on how efficient the mixing is. Allowing the field line footpoints to be transported by these flows, and assuming an ideal evolution in the corona, we can make a direct association between the tangling of these fluid particle trajectories in time and the induced tangling of the anchored coronal magnetic flux tubes.
We first examine the complexity induced by the flows in a qualitative way. By knowing the mapping of particles from positions ( ) x y , to ( ) F x y t , , at a time t, we can compute the mapping (pull-back) of a passive scalar (0-form or function) ( ) x c , similar to Candelaresi et al. (2017). This provides a visual representation of the mixing properties of the flow. For our initial large-scale signal, we use a simple gradient in x and y of the form . The mapped distribution, or equivalently the pull-back, is simply ( ( )) F c x y , . From the distribution of the mapped passive scalar (Figure 2), we can clearly see the turbulent nature of the thermo-convective simulations. We observe fine-scale structures and a high level of fluid mixing.
A natural question is whether there exists any characteristic scale in the pattern obtained, or a rather a whole spectrum. We cannot, however, judge from Figure 2 if there are any particular length scales forming. In order to seek the existence of characteristic length scales, we take the two-dimensional Fourier transform min max min max . From this, we compute the shell-integrated power spectrum of the passive scalar as with the shell width dk. Given that the simulation is turbulent, we expect no such scales to arise. This is indeed what we observe from the timedependent power spectrum of the passive scalar distribution (see Figure 3). Already, Candelaresi et al. (2017) showed that small-scale structures also form quickly for the blinking vortex benchmark flow, but in that case as well, a characteristic scale is absent.

Quantifying the Entanglement
While the passive scalar transport and associated spectra give a visual impression of the mixing, they do not encode any topological information about the trajectories. In order to measure the induced braiding, we consider two measures, one of which is the (finite time) topological entropy (Candelaresi et al. 2017) that is conceptually similar to the Lyapunov exponent, which was computed by  for closely related convective flows. However, first we calculate the winding number for field lines in the domain (Prior & Yeates 2014;Mangalam & Prasad 2018).

Winding Number
For a magnetic field configuration, it was shown (Prior & Yeates 2014) that the vertical magnetic field weighted winding 2 2 The angle between the two particles at time t is then calculated as x t 0 , 0 , arctan , 9 1 2 2 1 2 1 and its time derivative is Q = -´-- Summing over all trajectories, r 2 gives us the net winding rate of all other trajectories around the trajectory r 1 : By integrating over time and averaging over space, we obtain the averaged net winding around each trajectory, or averaged winding: This number will depend on the integration time T, and we need to circumvent this dependence when we make comparisons between different flows. Furthermore, there is an arbitrariness in identifying time units in both the simulations and the blinking vortex flow. (A similar arbitrariness exists in the details of the time profile in the "blinking vortex" flow: one could equally choose some other dependence that produces the same field line mapping between the bottom and top boundaries.) For these reasons, we need to apply a normalization for the averaged winding ω that eliminates such bias.
We choose to normalize by an averaged trajectory length in the xy-plane with respect to the typical granular size l granules : This number q is a measure of the average distance traveled by the trajectories in terms of the granular size l granules within the time T, while the time it takes to cross a granule of size l granules is t = T q granules . With this normalization, we can write the normalized winding number as . This quantity should be interpreted as follows. It measures the average winding of trajectories around a particular trajectory r 1 during one granular crossing time.
For the granular size in the simulation data, we use the characteristic convective scale given by Bushby et al. (2012), which corresponds to approximately one-quarter of the box size, i.e., 2.5 in code units. Since the blinking vortex set-up does not have a granular scale, we use the diameter of one of the vortices as a proxy for such a scale. Here that size is 2 units. Furthermore, to take account of the fact that the flow does not fill the entire domain (see Figure 1) Figure 4 for the magneto-convection simulations (at normalized time t = t 3.043 granules ), and in Figure 5 for the blinking vortex (at normalized time t = t 2.063 granules ). Both exhibit a complex pattern, although there exist much smaller scales in the pattern obtained from the simulations. Note that the times are only by chance this similar.  For a quantitative comparison of the tangling between the two flows, we can consider the maximum or the spatial average of the modulus of Ω. While both of these quantities depend on the integration time T, it is seen in the animated version of the figures that they converge as T is increased.
Selecting representative times to make a comparison, we find a mean of the absolute value of Ω for the convective simulations of 0.122 and a maximum of 0.521. For the blinking vortex, we obtain a lower value of mean W | | (computed within the region -´-[ ] [ ] 3, 3 2.5, 2.5 , corresponding to the region of significant twisting) of the absolute value of 0.0243, and a maximum of 0.2122. This region is smaller than the region of significant velocities, as peripheral particles tend to induce lower amount of twist, despite their significant velocities. This is contrary to our expectations from the known highly efficient mixing property of this flow. This is discussed further below. In order to achieve a comparable winding number for the blinking vortex motion, we would need to decrease l granules relative to the size of the vortex pattern by a factor of = 0.434 0.2122 2.045, leading to a granular size of » 2 2.045 1. We finally comment on the effect of the decomposition described in Section 4. When splitting the velocity field into incompressible and compressible parts using the Helmholtz-Hodge decomposition, we claimed that, for the winding number, contributions from u c can be neglected. We show this by computing á W ñ | ( ( ) )| r t 0 , r 1 1 -which is the spatial average of the norm of W( ( ) ) r t 0 , 1 -for the two flow components u i and u c separately as well as their sum. We clearly see that for the compressible part of the velocity u c alone, this number quickly drops close to zero (see Figure 6), while for u i and + u u i c , it levels off at a finite value. This justifies the usage of u i alone when computing the winding number.

FTTE
The FTTE is a measure of the topological complexity of a flow (Adler et al. 1965;Sattari et al. 2016;Candelaresi et al. 2017). It was shown by Newhouse & Pignataro (1993) that it can be interpreted as the exponential stretching rate of a material line γ by the flow (specifically, the maximal stretching exponent over all possible material lines γ in the domain). We can make use of this interpretation to measure the topological complexity using the FTTE, given a mapping g ( ) F of a material line γ. Through the repeated application of the mapping F, we can compute the FTTE. However, for the given velocity fields, we cannot meaningfully construct such a repeated mapping, since the flows are not time-periodic. Therefore, we make use of the mapping g ( ) F t , at time t and compute the FTTE as where l 0 is the length of the initial line γ, and l(t) is the length of the mapped line, g ( ) F t , .
In order to resolve the potentially highly entangled mapped line, we adaptively insert points on the original line γ where needed (see Candelaresi et al. 2017, for details). If we want to compare the FTTE of the blinking vortex (investigated in detail by Candelaresi et al. 2017) with the convective flows, we again need to normalize the time in some appropriate way. Therefore, as before, we scale the time t with the granule crossing time t granules for both the simulations and the blinking vortex.
We choose γ to be the horizontal line centered in y crossing the entire domain in the x-direction. (It is found that the choice of initial line makes little difference to the obtained value of h, due to the complexity of the flow.) A mapping of this initial line is shown in Figure 7. Due to the computational expense and the exponential growth of computing time with increasing time, we compute the FTTE by using only a small number of time steps from the simulations. This turns out to be more than sufficient for estimating the entropy, which requires only that we access the phase in which the material line length grows exponentially. As seen in Figure 8, a good fit is obtained even for the short period of time over which we are able to compute, allowing for a relatively accurate estimate to be obtained for h, c.f. Equation (15). For the magneto-convection simulations, we use 20 time steps and calculate a value for the FTTE of h=2.078 (see Figure 8).
For the blinking vortex mapping, we performed the calculations for the topological entropy in a previous paper (Candelaresi et al. 2017). For comparison of notation, the blinking vortex mapping corresponds to one-third of the "E3" mapping considered therein, and here we normalize the time by the granule crossing time. We find that the blinking vortex flow exhibits a level of efficient mixing which results in an FTTE of 0.4928. So, the FTTE for the simulations is higher by a factor of 4.217. This means that it takes 4.217 times longer to reach an equivalent state of entanglement in the blinking vortex flow than in the flow from the simulations. Similar to the winding number we attribute the higher value for the FTTE in the simulations to the fact that the turbulent and space filling motions of the simulations induce a velocity pattern that can be considered to be more chaotic.

Analysis of Flows Derived from Observations
We finally turn to compare the results obtained so far to those from observational data. While horizontal photospheric flows cannot be measured directly, they can be inferred by various different methods, though each has its limitations (Welsch et al. 2007). We use line-of-sight magnetogram data of the active region 10930 that is based on Hinode/SOT observations (Tsuneta et al. 2008) of Fe I at 6302 Å (see references in Yeates et al. 2012, for more details). From these magnetograms, Fisher & Welsch (2008) extracted velocity fields using local correlation tracking (November & Simon 1988). With an observational noise level of 17 G, the tracking threshold for the magnetic features was chosen to be 15 G when extracting the velocities. Subsequently, the magnetograms were binned into blocks of 2×2 pixels from a resolution of 0 16 to 0 32, which corresponds to 232.09 km (with 1″ on the Sun corresponding to 725.281 km). The observations start at 14:00 UT on 2006 December 12 and run until 02:58 UT on 2006 December 13. To reduce noise, this time series is averaged over four images, which results into a cadence of 121 s. An image of the magnetogram data is shown in Figure 9. This figure is available as an animation in the online journal.
Observational data are inherently affected by noise. This is already seen in the magnetic field. For the computed velocity field from local correlation tracking, the noise is only increased, resulting in frequent spikes of over 1000 km s −1 , while for most of the domain the velocities remain below 1 km s −1 in t granules together with linear fits for the magneto-convection simulations, the benchmark blinking vortex flow and the observational data. Note that the data points for the blinking vortex (stars) go well beyond the range of the plot, and its linear fit is performed on 13 data points. magnitude. This noise is greatly reduced by over two orders of magnitude once we remove the purely divergent part of the velocity data using the Helmholtz-Hodge decomposition as described in Section 4.
From observations, we know that the average granular size is of the order of 1 Mm (e.g., Priest 2014), which we use in our normalization for the winding number, which leads to = q 2.72 obs and t = 4.696 hr granules using the definitions in Section 6.1. For the observed velocity field, we find maxima of the normalized winding number W( ) r of the order of 0.05 (see Figure 10), while for the mean of the magnitude, we obtain 0.00994. This means that for every time the plasma travels (on average) a distance of the granular size, the average winding angle of all field lines around a given field line increases/ decreases by =  0.0517 rad 2 . 959. The spatial distribution of this averaged, normalized winding Ω is surprisingly smooth. With granules of the size of 4 pixels (1Mm), one might expect variations on the same scale. However, we observe in Figure 10 relatively homogeneous patches with the same sign of winding number over lengths of 10 Mm.
Again, with an appropriate time normalization, we can compute the FTTE for the observed flows for comparison with the previous results. Doing this, we obtain a value for the FTTE of h=1.598 (see Figure 8). We can now compute the physical time for which the observed motions would lead to a braiding of an initially vertical magnetic field equivalent to-in the sense of having the same FTTE as-the benchmark blinking vortex flow (with = l 2 granules ). The blinking vortex flow has a total logarithmic line stretching of 1.41, which is obtained within a normalized time of 2.235. With an FTTE of h=1.598, such a line stretching is obtained by the observations within 0.682 normalized times, which corresponds to the physical time of´= 0.682 4.485 hr 3.059 hr.

Magneto-convection Degraded
Comparing Figures 4 and 10, we observe significantly smaller lengths scales in the winding number plots for trajectories from the simulated flows compared to those derived from observations. In particular, the granular scales are now visible. Perhaps more importantly, the peak winding number is a factor of 10 larger for the simulated flows. One natural explanation for this could be the higher resolution of the simulated flows compared to the observed flows. To test this conjecture, we degrade the velocity data from the simulations such that granules have the same resolution as in the observations. This is done by applying a Butterworth filter in k-space such that small-scale motions are filtered out. For the simulation data, we know that a typical convective cell extends ca.2.5 code units , which translates into 128 pixels. For the observational data, we have a resolution of 232.09 km per pixel and a granular size of ca. 1 Mm, which means ca. 4 pixels per granule. Based on these numbers, we degrade the simulation data by using a Butterworth filter that filters out motions below the size of 32×32 pixels. In k-space, this translates to a value of 16, since the resolution is 512 (512/32). We can then calculate the winding number using this degraded velocity data.
Recomputing the winding number for these degraded flows we find that, contrary to the above conjecture, the values of Ω are similar to those computed for the highly resolved original velocity data (see Figure 11). These values are significantly higher than those obtained from the observational data. We conclude that the observational resolution alone is not sufficient to explain low values for the winding number, and one needs to take into account the effects from the velocity extraction method, as was suggested by Welsch et al. (2007). Note in particular that in that study-where a rate of helicity injection was compared between exact values and those obtained by using FLCT to infer the velocity field-the FLCT method was shown to underestimate the helicity injection by a factor of ∼10, consistent with what we observe here (albeit for different underlying flow fields). Similar to the non-degraded calculations, we compute a mean value for W | | of 0.128. In summary, the observed results using FLCT should be treated as inconclusive by themselves.

Conclusions
We have quantitatively and qualitatively investigated the efficiency of photospheric motions in inducing nontrivial magnetic field line topology into coronal magnetic loops. The efficiency of this tangling is crucial for evaluating the contribution of the braiding mechanism first proposed by Parker (1972) for explaining the high temperature of the solar corona. Recently, Parker's hypothesis was put on firmer footing by Pontin & Hornig (2015), who proved that any force-free equilibrium with tangled magnetic field lines must contain thin current layers (the argument was extended to include fields close to force-free by Pontin et al. 2016). They estimated that for coronal plasma parameters, significant energy release should be expected after three to six iterations of the blinking vortex pattern (where the mapping shown in Figure 5 with T = 48 corresponds to three iterations). Crucial to determining the efficiency of the braiding mechanism, then, are the timescales for injecting such tangling and for the associated energy release. Herein we have investigated the complexity induced in the coronal field by photospheric flows, and we find that (even assuming a trivial identity mapping of minimal complexity at t = 0) such levels of complexity can be induced in a matter of hours.
For the magneto-convection simulations we saw a relatively high degree of trajectory winding compared to the benchmark case of the blinking vortex flow. This is reflected in the calculation of the FTTE that shows a similar difference. We attribute this surprisingly high degree of winding to the presence of volume filling turbulent motions. Those lead to highly tangled particle trajectories, which we observe in our calculations.
We also compared the results to observed velocity fields of a solar plage region where the velocities had been extracted using the Fourier correlation tracking method. The winding number and the FTTE were significantly below the simulations and the blinking vortex benchmark. This cannot be due to the lower resolution alone (factor 32 difference). Tests with degrading the simulation data show similar values as the high-resolution simulation data. So it must be attributed to other effects, like the method of extracting velocity information from a time series of magnetograms. Welsch et al. (2007) compared different methods for extracting velocity data from photospheric magnetic field time series and found that the local correlation tracking approach used here leads to a factor 10 inferior magnetic helicity injection rate, consistent with the factor we find for the winding number. On the other hand, other approaches for velocity extraction tend to require vector magnetogram data and/or further assumptions about the magnetic field structure. A further potential source of the discrepancy is the different regions on the surface that we consider. For the observed field, we limit our calculations to a plage region that is relatively quiet. By comparing with the pigtail braid, we also calculated that the plage region induces as much mixing as the pig-tail braid in just 3.059 hr.
Our calculations shed some light into the Parker braiding problem for coronal magnetic fields. It is clear that the build up of braids and tangles leads to small-scale variations in magnetic fields that can further induce strong electric currents, reconnection, and possibly heating. Although the velocity fields do not induce any net winding globally, we observe locally a build up of net winding in local patches comparable to or larger than the granular scale. This affects the tangling of the magnetic field on a timescale of hours. The introduction of this tangling is expected in general to lead to reconnection and associated magnetic energy conversion, as was shown in magnetic braiding experiments by Wilmot-Smith et al. (2011), Pontin & Hornig (2015. S.C., D.I.P., A.R.Y., and G.H. acknowledge financial support from the UK's STFC (grant Nos. ST/N000714 and ST/N000781). P.B. would like to acknowledge the support of the Leverhulme Trust (grant No. RPG-2014-427). For the graphs, we made use of the package MATPLOTLIB (Hunter 2007).