Global MHD Simulations of the Time-dependent Corona

We describe, test, and apply a technique to incorporate full-Sun, surface flux evolution into an MHD model of the global solar corona. Requiring only maps of the evolving surface flux, our method is similar to that of Lionello et al., but we introduce two ways to correct the electric field at the lower boundary to mitigate spurious currents. We verify the accuracy of our procedures by comparing to a reference simulation, driven with known flows and electric fields. We then present a thermodynamic MHD calculation lasting one solar rotation driven by maps from the magnetic flux evolution model of Schrijver & DeRosa. The dynamic, time-dependent nature of the model corona is illustrated by examining the evolution of the open flux boundaries and forward-modeled EUV emission, which evolve in response to surface flows and the emergence and cancellation flux. Although our main goal is to present the method, we briefly investigate the relevance of this evolution to properties of the slow solar wind, examining the mapping of dipped field lines to the topological signatures of the “S-Web” and comparing charge state ratios computed in the time-dependently driven run to a steady-state equivalent. Interestingly, we find that driving on its own does not significantly improve the charge state ratios, at least in this modest resolution run that injects minimal helicity. Still, many aspects of the time-dependently driven model cannot be captured with traditional steady-state methods, and such a technique may be particularly relevant for the next generation of solar wind and coronal mass ejection models.


INTRODUCTION
Surface magnetic fields are the key observational input to coronal and solar wind models.Used as the inner boundary condition, it is the distribution, polarity, and strength of the radial magnetic field, B r , at the surface that defines the ground state of the corona (the potential field) and largely determines the topology and geometric features of the corona's magnetic field at global scales (coronal holes and streamers).It is also the evolution of the surface flux, which includes the shearing, emergence and cancellation of magnetic features that determines how energy is stored and released in the corona.As such, ground-and space-based observatories routinely measure photospheric and chromospheric magnetic observables, which are then used to create coronal models.
With the exception of the PHI instrument onboard the Solar Orbiter spacecraft (Solanki et al. 2020), current earth and lionel@predsci.comspace-based observatories provide only an Earth-centered view of the photospheric magnetic flux on the Sun.The average rotation period of the Sun as viewed from Earth, known as a Carrington Rotation (CR), is ≈ 27.3 days, so a traditional approach for constructing model boundary conditions has been to use full-sun Carrington Synoptic maps.These maps are built up by combining the data near central meridian from from full-disk magnetograms as the sun rotates.As such, changes in the average shape and structure of the corona in time have often been characterized in a variety of models on CR-to-CR basis, ranging from simple Potential Field Source-Surface (PFSS) extrapolations to full Magnetohydrodynamic (MHD) calculations.In this modeling paradigm, successive monthly 'snapshots' of the corona are computed separately, with each calculation for a given CR being fully independent of the previous case (e.g., Riley et al. 2006;Luhmann et al. 2022).
On the other hand, we know full-well that many of the most intriguing and important phenomena in the corona (e.g., solar flares, jet eruptions, coronal mass ejections, etc.) are in-herently driven by surface flux evolution at timescales from minutes to days depending on the spatial scale of interest (Webb & Howard 2012;Raouafi et al. 2016;Benz 2017, and references therein).Similarly, at the largest temporal and spatial scales (i.e.global scales) surface flux evolution is responsible for the formation and evolution of helmet streamers and pseudostreamers over days to months, which is thought to play an important role in the processing of solar wind and ejecta as smaller structures erupt or migrate across their boundaries (Higginson et al. 2017;Scott et al. 2021;Wyper et al. 2021).How these large-scale closed structures evolve in response to surface changes and subsequently interchange reconnect with nearby open fields has important implications for our understanding of switchbacks measured by Parker Solar Probe (Fisk & Kasper 2020;Zank et al. 2020;Telloni et al. 2022), as well as the formation of composition and charge-state variations in the solar wind (Zurbuchen et al. 2002;Kepko et al. 2016).
In the context of global coronal models, one way to incorporate time-evolution at the inner boundary is by leveraging the outputs of surface flux transport (SFT) models.SFT models describe the time-evolution of magnetic flux over the full-sphere by incorporating various evolutionary processes; typically differential rotation, meridional flow, supergranular diffusion, and random flux emergence.Such models can also assimilate magnetograms from available observatories to produce a continuous approximation of the state of the photospheric magnetic field, also known as a sequence of synchronic maps.Using this approach SFT models have been successful in predicting the evolution of photospheric magnetic fields (Wang et al. 1991;Worden & Harvey 2000;Schrijver & DeRosa 2003;Arge et al. 2010;Upton & Hathaway 2014).
The magnetic flux information from an SFT model can subsequently be processed to create the full-sun boundary condition of B r for the global coronal magnetic field model.With a time-sequence of synchronic maps, it then becomes possible to model successive states of the corona and heliosphere at a much smaller time-interval than one Carrington rotation.This can be done by running successive (but independent) 3D calculations at the cadence of the synchronic maps (e.g.Odstrcil et al. 2020), or by driving the model at the inner boundary using electric fields derived from the evolving sequence of B r (e.g.Weinzierl et al. 2016).The latter, 'driven' approach is physically more attractive as it allows one to capture how surface flux evolution surface flux evolution imprints a dynamical 'memory' in the system as it evolves from state to state, and allows for the build up of magnetic stresses and energy in time at the correct dynamical timescales.
On the other hand, deriving an appropriate driving electric field from a sequence of B r maps is both challenging and not fully constrained (Fisher et al. 2010;Cheung & DeRosa 2012;Yeates 2017;Lumme et al. 2017).For global coronal coronal field models, such driving was studied with the simplified magnetofrictional approach by Weinzierl et al. (2016) who focused on the importance of the non-inductive freedom in the boundary electric field.While convenient for studying energy injection and storage, magnetofrictional models fundamentally cannot capture the dynamical timescales of the system (set by the magnetoacoustic speeds), non-force free processes (such as eruptions), nor how the 3D plasma state of the corona and solar wind evolve in tandem with surface flux changes (flows, density, temperature).Another challenge for full-sun, global driving is the fact that data assimilation in SFT models, which ingests new measurements from the Earth-Sun line as the Sun rotates, will instantaneously overwrite existing flux.These imposed changes imprint both an unphysical dynamical timescale and forcing on the system as well as a floating magnetic monopole that must (typically) be corrected by some means.
In this light, we describe our efforts to develop and test a suitable boundary driving approach for global MHD models of the solar corona.In §2, using a technique similar to Lionello et al. (2013), we illustrate how a full electric field may be expressly determined from time-evolving B r maps and the flow profile from an SFT model.At the boundary, the electric field obtained through our formulation evolves the magnetic field smoothly at the code time-step, thus avoiding the instantaneous overwriting of the flux.We also discuss additional corrections that can help eliminate spurious currents at the inner boundary of the MHD model.In §3 we test and compare the various approaches on an idealized case where the true surface flows are known.
Next in §4 we demonstrate the approach in a full thermodynamic MHD calculation for one-month time of coronal evolution, which is driven by a sequence of magnetic flux maps provided by the Lockheed Evolving Surface-Flux Assimilation Model (ESFAM, Schrijver & DeRosa 2003).This particular ESFAM run does not assimilate data, but is instead designed to yield Sun-like evolution with an average strength sunspot cycle at solar minimum conditions.This provides a smoothly evolving full-sun SFT dataset without the typical artifacts of data assimilation.To characterize the results, we calculate EUV synthetic emission images, as well as maps of coronal holes locations, the squashing factor (Q, Titov 2007), dips in magnetic field lines, and fractional charge states.To evaluate the importance of time-dependent evolution, we also compare these results with those of steady-state models.Finally, we conclude in §5 by discussing the our results in the context of the solar wind as well as the relevance of these techniques for future applications.

The MAS MHD Model
To develop, test, and study methods for evolving the magnetic flux at the inner boundary of an MHD model, we employ the Magnetohydrodynamic Algorithm outside a Sphere (MAS) code.MAS is designed to model the global solar atmosphere from the top of the chromosphere to Earth and beyond and has been used extensively to study coronal structure (Mikić et al. 1999;Linker et al. 1999;Lionello et al. 2009;Downs et al. 2013;Mikić et al. 2018), coronal dynamics (Lionello et al. 2005(Lionello et al. , 2006;;Linker et al. 2011) and CMEs (Linker et al. 2003;Lionello et al. 2013;Török et al. 2018).MAS solves the resistive, thermodynamic MHD equations in spherical coordinates (r, θ, φ) on structured nonuniform meshes.Magnetosonic waves are treated semi-implicitly, allowing us to use large time steps for the efficient computation of long-time evolution.The semi-implicit method is not harmful for obtaining time-dependent solutions, as it introduces dispersive effects only for processes occurring on timescales shorter than or equal to the time step (Schnack et al. 1987).Given that the flows that drive evolution of flux at the surface are much smaller than the typical coronal flow speeds that set the CFL limit and time step in the model, the semi-implicit method is more than suitable for the timedependently driven calculations described here.
The present version of MAS allows for several modes of operation that govern which terms are solved for and/or added to the MHD equations.Here we use the term 'thermodynamic MHD' to indicate that MAS solves for additional transport terms that describe energy and momentum flow in the solar corona and solar wind (coronal heating, parallel thermal conduction, radiative loss, and Alfvén wave acceleration; as fully described in Appendix A of Török et al. 2018).
The latest version of the thermodynamic mode in MAS also includes a physics-based specification of the coronal heating term through a Wave-Turbulence-Driven (WTD) phenomenology.In this approach, additional equations are solved to capture the macroscopic propagation, reflection, and dissipation of low-frequency Alfvénic turbulence.The physical motivation underlying the WTD approach is that outward and reflecting Alfvén waves interact with one another, resulting in their dissipation and heating of the corona (e.g., Zank et al. 1996;Verdini & Velli 2007).This follows related works, where the general formalism for the propagation of Alfvén waves (e.g., Heinemann & Olbert 1980;Zank et al. 1996Zank et al. , 2012) ) is usually approximated to produce tractable equations for the propagation of the energy density or the amplitude of the Alfvén waves (e.g., Velli 1993;Matthaeus et al. 1999;Dmitruk et al. 2001;Cranmer et al. 2007;Verdini & Velli 2007;Breech et al. 2008;Chandran & Hollweg 2009;Usmanov et al. 2011;Jin et al. 2012;Sokolov et al. 2013;van der Holst et al. 2014;Oran et al. 2015).A full description of the MAS-WTD model and equations solved is provided in the supplementary materials of Mikić et al. (2018).
We have also recently incorporated into MAS a nonequilibrium ionization module to advance the fractional charge states of minor ions according to the model of Shen et al. (2015): (1) For an element with atomic number Z, Z F i (r, θ, φ) indicates the fraction of ion i+ (i = 0, Z) with respect to the total at a grid point: For each element, the ion fractions are coupled through the ionization, Z C i (T ), and recombination, Z R i (T ), rate coefficients derived from the CHIANTI atomic database (Dere et al. 1997;Landi et al. 2013).Here T , n e , and v are respectively the temperature, number density, and the velocity of the plasma.This module has been tested in our hydrodynamic 1D, WTD wind code (Lionello et al. 2019).A similar time-dependent 3D model of charges states of minor ions is shown in Szente et al. (2022).

Incorporating Magnetic Flux Maps
To drive the magnetic field evolution in MAS, we can evolve the radial component of the magnetic field at the boundary using a technique similar to that described by Lionello et al. (2013).We specify the tangential electric field at the boundary E t0 (θ, φ,t) as (3) The potential Ψ controls the evolution of the normal component of the magnetic field B r0 , We use a sequence of full-sun B r (θ, φ) maps in time to specify the evolution of Ψ using Eq.(4).Since the temporal cadence of the input maps will generally be much slower than the MHD model time step, we must interpolate B r in time at every step.Because linear interpolation of the B r maps in time implies a discontinuous ∂B r /∂t as we shift from the interval between one pair of maps to the next, one must use a higher order interpolation scheme to ensures a smooth stepto-step evolution of the electric field at the inner boundary.
Here we use a simple piece-wise cubic Hermite interpolation scheme (Fritsch & Carlson 1980), to ensure continuity in at least the first derivative and preserve monotonicity in the evolution of B r (t) at all points within in the time-interval between maps.
Next, the potential Φ(θ, φ,t) can be completely specified if the flows that led to the evolution are known, i.e., we can find Φ from the equation: In general however, the full-sun distributions of the complete flow and magnetic field vectors responsible for map-to-map B r changes are not typically available from observations or SFT models that assimilate observational data.Instead we can incorporate the known large-scale flows (differential rotation and meridional flows) to at least include their contribution to the transverse electric field.In other words, while the Φ potential solve ensures that we can evolve flux to exactly match the B r component of the maps, the magnetic fields that evolve and emerge will generally have less shear and twist than is observed (especially for ARs).On the other hand, this method has a significant advantage over more complex or localized techniques in that vector magnetic field observables and/or flow correlation tracking are not required.Furthermore, additional energization can easily be explored through the freedom in the Φ potential1 , which will be explored in future work.At this point the Ψ and Φ potentials only provide information about the transverse electric field.Our next aim is to find a full boundary electric field, E 0 , that is (1) minimally diffusive and causes (2) minimal boundary layers.A possible solution is to prescribe the condition of ideal MHD at the lower boundary (method A): where t means the component tangential to the boundary of the B 0 and E 0 fields.From that it follows that which can be regularized using a small parameter ε as This method has the disadvantage of introducing artificial reconnection at the polarity inversion line (PIL), where B r = 0.
Building upon method A, we can compute an additional correction to the boundary electric field, Ẽ0 , such that Ẽ0 • B 0 ≃ 0 everywhere by slightly modifying E t0 and E * r near the neutral line (method B).To define Ẽ0 , we use which is the flow associated with E 0 as calculated with A: We obtain Whether method A or B is used, we derive again the boundary flow, ṽ0 , using the boundary electric field expression provided by either method, where we add a small parameter ε to regularize the flow at the neutral line.
Only for method B, we also add a small resistivity value at the boundary PIL: where we dropped the subscript 0 from the notation.

Asymptotic Values at the PIL
Because each method requires some form of regularization near the PIL, it follow that we should seek a robust dimensionless formulation that is equally appropriate all types of broad/compact and weak/strong flux-distributions on the Sun.
We can rewrite Eq. ( 8) in terms of f = ε 2 /B 2 r and obtain The right-hand side (RHS) of Eq. 15 can be easily computed at any state of the simulation, regardless of dB r /dt or the full electric field itself.We seek a definition of f such that E r is bounded near the PIL and f is controlled by a dimensionless constant, unlike ε, which is in units of B. Let us therefore find a form of f so that the RHS of Eq. 15 will reach a fixed, asymptotic value when B 0 → 0. We simply define f as where A is a dimensionless constant, and 1/A sets the largest value for the LHS of Eq. 15.The new form assumed by Eq. ( 8) is thus the following: Similarly, we seek a new formulation for the boundary flow ṽ0 in Eq. ( 13) that will converge to a meaningful asymptotic value for ṽ0 , when B 0 → 0 at the PIL.We rewrite Eq. ( 13) as where g is a dimensionless quantity.We proceed heuristically by comparing ṽ0 to the local sound speed, calculating the magnitude, and dividing by c s to obtain Then we pose where D is a dimensionless constant, and insert this into Eq.( 19).The ratio simplifies as which converges to 1/D as B 0 → 0. In other words, 1/D specifies the maximum fraction of |ṽ 0 | relative to the local sound speed at the boundary.The new formulation of ṽ0 thus becomes

TEST OF THE CORRECTIONS
To evaluate the performance of the correction methods in §2.2, we can apply them to a solution where the true electric field that evolved the boundary is known.To this end, we apply them to time-dependent simulations driven with the magnetic flux maps extracted from the 3D MHD coronal solution of Lionello et al. (2005Lionello et al. ( , 2020)).These employed the polytropic model of Linker et al. (1999) to study the effects of differential rotation on a photospheric magnetic field distribution similar to that of Wang et al. (1996).This consists of a global dipole distribution with a bipolar active region superimposed (see also Fig. 2).After a relaxation period lasting till t = 150 code units (1 CU = 1445.87s), they turned the following surface flow on: which is 10 times the value of Wang et al. (1996).We repeated this run in the same manner, evolving the surface field directly using the prescribed E = −v × B from the driving flow (only v φ in this case), and B 0 at a given timestep.We label this run with O (original), and sample the values of B r on the solar surface with a 1 CU cadence.
Then we use the surface B r maps from O to run four simulations using Eqs.(4-5) to prescribe the electric field at the r = R ⊙ surface: the first run is uncorrected (U, i.e. neither of the methods of §2.2 is applied but the 'true' radial electric field is used, E 0 r = ωR ⊙ B 0 θ ), the second run (E) is erroneously driven with E 0 r = 0, the third run is corrected with method A, and the fourth is corrected with B.
In Fig. 1 we present the time histories of the four runs from the inception of photospheric flows at t = 150 until t = 170.Panels (a-f) show quantities integrated over the whole computational domain.With the exception of currents, the U run presents values coincident with those of the O run in all panels.The magnetic energy (Fig. 1a) of the four runs are also practically indistinguishable.However, following the introduction of differential rotation at t = 150, the kinetic energy (Fig. 1b) of the A and B driven runs is ∼ 1% smaller than the reference O run.This behavior is replicated in the time history of the magnetic energy in Fig. 1c.The integrated Lorentz force (Fig. 1d) shows a jump at (t = 150) for all runs, with the curve of the O and U runs bracketed by those of of the corrected runs A and B. The last integrated quantity |J • B| measures the deviation from ideal MHD.Similarly to (d), U and O are between A and B. The erroneous run, E has larger values of the integrated Lorentz force and parallel current.
We also examine the history of the maximum current density J on three orthogonal cut-planes in spherical coordinates: the r = R ⊙ surface encompasses the whole active region; the θ = θ 0 surface intersects the active region; the φ = φ 0 lies west of the positive polarity.When the flows are turned on at t = 150, we notice in panel (f) a sudden spike on the r = R ⊙ surface for the (O) and U run, while the enhancements for A and B are much smaller.As it was the case for the global measurements in (d) and (e), the (O) and U values soon become bracketed by those of A and B. Run E shows high values of the maximum currents on the r = R ⊙ and φ = φ 0 surfaces.The (g) and (h) panels shows that, with the exception of E, all the driven runs manage to reproduce the maximum J on θ = θ 0 and φ = φ 0 surfaces in substantial agreement with the O run.
To convey a more intuitive sense of how the corrections operate, we show full-sun surface maps of various quantities at the final state of runs O, E, A, and B in Fig. 2 (skipping U because it is essentially the same as O).First, to illustrate the shape and location of the PIL near the large bipolar region, the upper left panel shows the full-sun B r distribution.The remaining panels in the left column show E • B for the computed electric field for runs E, A, and B (E • B is zero for O by construction).As expected, E • B is largest when E r  is ignored altogether in run E and the signature is strongest near the largest fields along the the PIL.Method A reduces the maximum E • B compared run E by just over a factor of 5, while method B effectively eliminates it as intended.
The middle column of Fig. 2 illustrates how field-aligned currents build up at the boundary in each case, showing a map of log 10 (|J • B|/B 2 ).Compared to the reference case O, all runs introduce some level of additional currents but the features change from case to case.As expected, by not applying any E • B = 0 correction, run E builds up the largest boundary layer, especially near the main PIL of the two polarities.Method A attenuates this current signature by a factor of two or more (note the log 10 scaling of the min/max values), while method B attenuates this signal even further.However, right at the PIL for run B we see a slight enhancement of J • B compared to the other runs, which is the signature of the small resistivity that is active where B r → 0 that is needed to ensure the B r evolution matches the driving sequence of maps (Eq.14).
Lastly, the right panels of Fig. 2 show maps of the longitudinal component of the velocity field perpendicular to B at the inner boundary, v ⊥,φ .For run O this is simply the component of the differential rotation flow that actually moves the field eastwards, while for the remaining cases this component is determined from Eq. 22.For E, A, and B we see that the true driving flow is largely recovered everywhere except near the PIL of the flux concentrations.Methods A and B do a better job than run E near the PIL, especially in filling in the strong patches of v ⊥,φ to the southeast and northwest of the central PIL segment.That said, on the opposite side of the sun where dipolar flux-distribution is basically symmetric north south (and thus B r is not changing) the maximum strength of |v ⊥,φ | is not quite recovered.This stems from the fact that the PIL for the global dipole is very broad, such that the truncation function for determining E 0 r (Eq.16) is not quite zero even at mid-latitudes.
In summary, these tests runs illustrate the inherent limitations of recovering a suitable electric field to drive surface flux evolution from limited information (a sequence of fullsun B r maps in time) as well as the pros and cons of the various methods.When the full surface flow profile is known, the combination of E t derived from the Ψ and Φ potential solves and the true E r , which can be computed if the exact surface flow is known, reproduces the true solution (runs O vs U).When E r is not available, which will be generally true when only the full-sun evolution of B r and the average macroscopic flows are known (i.e.flux-transport models assimilating observations), it must be computed from available information (runs A vs B), or ignored (run E).
Although both correction methods A and B appear to give similar results, method B is able to completely eliminate the non-ideal part of the electric field, which may help it perform better when strong magnetic flux is carefully emerged (e.g., during the simulations of CMEs).On the other hand, the advantages of method B comes at the at the cost of an additional free parameter ( η) and concentrating the resistive boundary layer right at the PIL.In favor of smoothness and simplicity, we thus opt to use method A in our first application of this approach to a more realistic case, which is described in the following section.

TIME DEPENDENT WTD MODEL
We now present a time-dependent simulation of the solar corona obtained with the WTD model of Mikić et al. (2018), driven with the transverse electric field from Eqs. (4-5), and the radial electric field from method A. As input to the model, we used 720 1-hour cadence frames of the ESFAM model (Schrijver & DeRosa 2003), which corresponds to 30 days of evolution or just longer than one Carrington rotation.This is the same model that is implemented within the PFSS package in SolarSoft as the SFT (surface-flux transport) module, which routinely assimilates photospheric magnetogram to best describe the surface conditions at a given time.
For this case we use a special 'fake sun' ESFAM run designed to represent a the typical conditions during solar minimum, including large-scale surface flows, the emergence and decay of bipolar patches, and the evolution of random flux across a range of scales.For our purposes, this type of sim-ulation is ideal because we can sidestep some of the typical systematics of SFT maps that assimilate data.These include: (1) zero-point offsets for the entire map because at a given assimilation step only a portion of an active region is observed; (2) unphysical evolution driven by the sweep of the assimilation window, introducing new data instantaneously, as it moves across the Sun at the solar rotation rate.We now describe these maps in more detail.

SFT Maps
The main methodology for the formation of the SFT maps is described in Appendix A of Schrijver ( 2001), but we briefly summarize it here.The emergence pattern for bipoles is based in an automated process that randomly selects bipoles with values for their flux, location, and axial tilt based in power law distributions derived from long-term observational statistics.Once the flux has emerged-a process which is scaled to the total flux of the bipole-the polarities advect and cancel with nearby flux independently of each other.For very large bipoles (i.e., active regions), the flux is emerged as a group of smaller bipoles that aggregate to the characteristics listed above.Smaller bipoles are also emerged with similarly randomized qualities (and a correspondingly larger latitudinal distribution) in order to create a realistic background and distribution of the total flux.To process this total flux and model the cancellation on the correct time scales, there is an additional exponential flux removal term (described in Schrijver et al. 2002).
The flux in the SFT module was advected with the differential rotation profile taken from Table 1 of Komm et al. (1993a), in the case for the one-dimensional crosscorrelation analysis applied to Kitt Peak magnetogram data from 1975-1991: The above is a general functional form for differential rotation fits.Here, Ω means the sidereal rotation rate, and λ is latitude.The coefficients measured by Komm et al. (1993a) are As implemented in the SFT model, a Carrington frame of reference is used, and so the solid-body sidereal Carrington rotation rate of 14.18 deg/day is subtracted off from A. Also, in the SFT model, the rotation profile does not change with time (i.e., no torsional oscillations are present, etc.).
The meridional flow profile M is slightly more ad hoc, but still takes its cues from empirical measurements.Like the differential rotation profile, it is implemented in the SFT model data as a constant-in-time function and does not vary.The SFT model's profile was originally based on the measurements described in Komm et al. (1993b), which are of the form M = f [sin(2θ)], or more specifically, sin(2θ) is the leading term in an expression also involving sin(4θ).However, the current SFT model ignores those fourth order terms, and also includes a tapering function applied to the polar latitudes, so that the functional form looks is the following: Here, M is the meridional flow speed and θ is colatitude, with the coefficient v A set to 12.7 m/s.Note that we could switch from latitude to colatitude since sin(2lat.)= sin(2colat.).The tapering functions g only apply poleward of 40 degrees in each hemisphere: where a = 3.0.Without these tapering functions, the meridional flow will concentrate flux too close to the poles in the SFT model, in contrast to observations in which it is evident that the polar cap flux during solar minimum intervals is more spread out.Schrijver & Title (2001) felt justified in using this tapering function given the higher uncertainties of measurements of meridional flows near the poles.Note that this meridional flow profile is functionally similar to that given in Eq. (3) of Schrijver & Title (2001) but with different values of v A and a (these values can be found in Barnes et al. (2023)).

Properties of the Simulation
The ESFAM/SFT maps were given as input to the MHD WTD model (Török et al. 2018;Mikić et al. 2018).For the latter we used a nonuniform grid in r × θ × φ of 269 × 181 × 361 points extending from R ⊙ to 30R ⊙ .The smallest radial grid spacing at r = R ⊙ was ∼ 400 km; the angular resolution in θ ranged between 0.8 • at the equator and 1.7 • at the poles; the φ mesh was uniform.To dissipate structures that cannot be resolved since they are smaller than the cell size, we prescribed a uniform resistivity η corresponding to a resistive diffusion time τ R ∼ 4 × 10 2 hours, which is much lower than the value in the solar corona.The Alfvén travel time at the base of the corona (τ A = R ⊙ /V A ) for |B| = 2.205 G and n 0 = 10 8 cm −3 , which are typical reference values, is 24 minutes (Alfvén speed V A = 480 km/s), so the Lundquist number τ R /τ A was 1 × 10 5 .Also, in order to dissipate unresolved scales without substantially affecting the global solution, we introduced a uniform viscosity ν, corresponding to a viscous diffusion time τ ν such that τ A /τ ν = 0.015.We prescribed fixed chromospheric values of density and temperature at the base of the domain of n 0 = 4 × 10 12 cm −3 and T 0 = 17,500 K, respectively.These values were set to form a chromospheric "temperature plateau" that remains sufficiently large (Lionello et al. 2009) during the calculation no matter how large the heating.
For the coronal heating term, we use the same WTD model parameters as the simulation described in Boe et al. (2021Boe et al. ( , 2022)), which is a slight update to the numbers used in Mikić et al. (2018).The Poynting flux of wave energy is prescribed at the base of the corona through an amplitude of the Elsässer variable z 0 = 9.63 km/s, and we set the transverse correlation scale λ 0 = 0.02R ⊙ along with a scaling factor B 0 = 8.53 G such that λ ⊥ = λ 0 B 0 /B in the corona.Similar to Mikić et al. (2018) in adding two small exponential heating terms to heat the low corona: H 0 = 2.7 × 10 −5 erg/cm 3 /s, λ 0 = 0.03R ⊙ ; H 0 = 1.6 × 10 −8 erg/cm 3 /s, λ 0 = R ⊙ .Likewise, the wave pressure was specified from the WKB model (Lionello et al. 2009).
We started the calculation using the first of the 720 magnetic flux maps to calculate a potential field extrapolation.The plasma temperature, density, and velocity were imposed from a 1D solar wind solution that had been calculated previously.Then we advanced the MHD equations for about 80 hours to relax the system to a steady-state.After the relaxation was accomplished, we turned the surface evolution on.The differential rotation and meridional flow parameters were the same as those of Eq. ( 25), except that, since our model is in the corotating frame of reference of the Sun, A = 0.24 deg/day for us.

Coronal Evolution
The left-hand side of Fig. 3 and its associated online animation show the B r field evolution during the course of the simulation, while the right-hand side shows the open field for the same times.The evolution of B r at higher latitudes is dominated by differential rotation, and it is possible for the eye to pick up features as they are advected from west (right) to east (left), particularly in the associated animation available online.A few of the most prominent of these dynamic regions are indicated with circles on the figure.The black (left) and gold (right) rectangles appearing in all frames in the figure highlight a particularly persistent and complex region of mostly open field with a sizable embedded parasitic polarity; the whole feature drifts slowly to solar east over time due to differential rotation, and the open field in the area changes significantly during the month of simulated time.
However, at lower latitudes, differential rotation becomes less and less discernible, while emergence and dispersion of magnetic flux become dominant.The magenta circle in the same frames shows a relatively strong bipolar region which is present at the beginning of the simulation and which un- The results of these time-dependent effects on simulated emission are illustrated in Figs. 4 and 5.We use the temperature and density from the simulation to forward-model the emission in the SDO AIA 171 Å channel for this figure.The left-hand column of Fig. 4 shows a 3D projection observed from an Earth-like point of view at the same times as the preceding figures (i.e. the longitude of the observer changes with the Carrington rotation rate).To better understand the evolution of the emission, we can compute Carrington latitude-longitude maps of forward modeled AIA emission by integrating along the radial direction for all points on the model.These are shown in the right column of Fig. 4, with the same annotations as the preceding figure to aid in by-eye comparisons.Here, the bipoles appear as bright active regions, while particularly bright patches within the active regions which evolve rapidly (best seen in the associated animation) are signatures of thermal nonequilibrium within these areas of highly-stratified heating.The open fields contain somewhat lower emission than the closed field, though not as dark as would be expected from a higher-temperature emission line (such as 193 Å, shown in the following figure).
In Fig. 5 we show of the surface B r prescribed at the lower boundary, overlaid with forward-modeled SDO AIA 193 Å emission.This composite presentation shows the close relationship between the changing magnetic boundary conditions and the coronal evolution.Features visible in emission in the polar coronal holes are matched by the unipolar signatures in the underlying B r map, and are likewise advected by differential rotation within the coronal holes.The black box repeated on each panel of Fig. 5 highlights one such persistent open flux region.If only differential rotation were present, the open flux evolution would be (mostly) rigid, as observations (Timothy et al. 1975) show, and both potential (Wang et al. 1996) and MHD (Lionello et al. 2005) models indicate.However, magnetic flux evolution contributes to continuous changes in the boundary between open and closedfield region.This is particularly evident (but not limited to) the equatorial region between 270 • and 360 • longitude.Here on frame 543, flux emergence and open-field evolution are highlighted with the green box and cyan circle, which highlight an area in which emerging flux results in an expanded neighboring open field region.In contrast, the red circle that appears in each panel illustrates a case of magnetic flux decay, as a bipolar region becomes more diffuse over the month timeframe of the simulation.Because of the richness of the EUV evolution and its relationship to interchange reconnection with open and closed fields, we explore the effects of these and related coronal signatures in the companion paper analyzing this simulation (Mason et al. 2023).

Dipped Field Lines
Parker Solar Probe, during its first perihelion, discovered that the radial magnetic field was continuously interrupted by switchbacks on a time scale of less than a second to more than one hour (Dudok de Wit et al. 2020).Although our simulation was aimed at reproducing lower frequency phenomena and at larger length scales, we also investigated the formation of dips in magnetic field lines.First, for each cell in the computational domain we determined whether there was a dip in its neighborhood by tracing a magnetic field segment and checking whether B r changed sign along it.Then, from each point on an r = 19R ⊙ surface, we traced a field line and verified whether it passed through a cell having a dip.The resulting map is shown in Fig. 6a, where the seed points of dipped field lines are painted in white.Figure 6b has a map of s-log Q (signed logarithm of the squashing factor Q, Titov 2007) on the same surface.At low latitude, the so-called S-web generally coincides with the region where the slow solar wind is observed (Antiochos et al. 2011).A comparison between panel (a) and (b) indicates that dipped field lines have connections with the S-web.In particular, we enlarged two regions in the S-web map, the corresponding areas of the dipped field-lines map, and superimposed them in panels (c) and (d).While the field lines associated with reversal in the sign of B r lie along high Q lines, only some have footpoints aligned with the current sheet (Fig. 6c).Other dipped fieldlines cross single-polarity separatrices as in Figs.6d.In panel (e), the dipped field-line map of panel (a) is shown as a semitransparent surface with two groups of representative field lines intersecting the areas of panels (c) and (d).Likewise, panel (f) presents the equivalent point of view for the s-log Q map.Panels (g) and (h) present 3D enlargements around the (c) and (d) areas respectively.Some of the field lines in (g), which are all associated with the current sheet, are arranged in a flux-rope.On the other hand, some of the field lines in panel (h), which are all associated with a same-polarity high Q line, form a V shape, indicative of interchange reconnection (e.g., see Fig. 5c and 6c of Lionello et al. (2005)).

Charge-States in the Static and Time Dependent Corona
To understand the effects of the continuous evolution of the surface magnetic fields on the ion charge-state distributions of the corona and solar wind, we compared the results of the time-dependent model at t = 524 hrs with its corresponding steady-state model.This steady-state-corona (SSC) model was obtained by stopping all surface flows and magnetic flux evolution and letting the system relax for approximately 80 hours.We calculated the maxima of the C 06 /C 04 ratio along  magnetic field lines.This is the ratio recommended by Landi & Lepri (2015) for analysis and comparison between models and in situ data as less sensitive to photoionization.In Fig. 7 we show statistical distributions using either a linear (a) or a logarithmic scale (b) in the x axis.The solid lines are associated with the SSC model, the dotted lines with the TDC one.We distinguish four regions: loops (closedfield regions), coronal holes (open-field regions), which are mapped in panel (c), edge loops (i.e., long loops bordering open-field regions), edge corona holes (i.e., open-field lines but close to closed-field areas), which appear in panel (d).Although the loops (violet) distributions are similar, we notice that the TDC model has a C 06 /C 04 distribution in coronal holes (green) that is significantly higher for smaller ratios than that of the SSC model.This is even more evident in the distribution of field lines at the edge of coronal holes (gold), with visibly higher distribution values for C 06 /C 04 ≲ 1.

DISCUSSION AND CONCLUSION
We applied to our MHD model a technique that drives the evolution of the photospheric magnetic flux and of the surface flows to calculate the response of the solar corona.The electric field, which the technique provides, may be corrected with two methods to reduce the formation of current bound-ary layers.We tested the two corrective methods by comparing results with a simplified, reference simulation and found that either gives results in satisfactory agreement.We then used our technique with method A to calculate the evolution of the corona for a whole month, using a sequence of balanced magnetic-flux maps obtained with the SFT module.From the simulation data, we produced emission images that show a continuous reconfiguration of the corona as active regions emerge and disperse and surface flows rearrange the flux.Dipped field lines that are formed during the computation appear to be associated with the S-web.In particular, we identified a fluxrope in the current sheet, while field lines of single-polarity areas show the typical pattern of interchange reconnection.Ions fractional chargestates were evolved alongside the dynamics.The distribution of the C 06 /C 04 ion charge-state ratio appear to differ between the time-dependent model and the steady-state corona.This difference is visible in the open-field regions, for which the time-dependent model has a distribution more skewed to lower charge states than the steady-state corona, and particularly evident in the field lines close to the border of coronal holes.
Counterintuitively, it is the TDC model that has the lower ratios.The conventional wisdom for enhanced ratios in the slow-wind is that they represent hotter plasma from closed loops that has been liberated through interchange reconnection.In the TDC model this can occur either by surface flows and/or evolution or by a thermal and/or tearing instability introducing reconnection at the streamer cusp itself (e.g.Réville et al. 2020), while in the SSC model this can only occur through the latter processes.On the other hand, the TDC model is slightly more open and the streamer loops naturally have a shorter lifetime as the surface flux evolves and their geometry changes in time, leading to a lower temperature on average.It seems that in this case, with relatively modest resolution and minimal energization in the transverse field ( §2.2) that the latter process wins out.In future work we aim to study how additional shear (helicity) emerged at global scales will influence the charge-state ratios as well as the role of small-scale flux patches that inducing interchange reconnection near the coronal base as they evolve (e.g.Sterling & Moore 2020; Bale et al. 2022).
Ultimately we have demonstrated a practical technique for time-dependently driving a global coronal model from only a sequence of magnetic maps.This leads to an inherent timeevolution of features in the model corona that is not possible to capture with traditional steady-state methods.Such a capability is crucial for answering long-standing questions about how closed and open fields evolve in the corona, and how the field and plasma properties of the observed slowsolar wind are formed.Lastly, with the advent of modern SFT models that readily assimilate observations to produce time-sequences of full-sun magnetic maps, we expect simple but robust driving techniques like those introduced here to play a key role in the next generation of solar wind and CME models.

Figure 1 .
Figure 1.Time histories of the simulations respectively labeled 0, U, E, A, and B in §3.Volume integrals of the magnetic energy (a), of the kinetic energy (b), of the thermal energy (c), of |J × B| (d), and of |J • B| (e).Maximum current density J evaluated on the surfaces r = R ⊙ (f), θ = θ 0 (g), and φ = φ 0 (h).The locations of θ 0 and φ 0 are shown in the upper left panel Fig. 2.

Figure 2 .
Figure 2. 2D maps of surface quantities at the final step (t = 170 CU) for test cases O, E, A, and B. The upper left panel shows B r at this time for reference and the positions of the θ 0 and φ 0 cut planes used in Fig. 1 are indicated with the dotted gray lines.The remaining left panels show E • B for each case, where E is the driving electric field.The middle column shows log 10 (|J • B|/B 2 ) for all four runs.The right column shows the φ component of the driving velocity that is perpendicular to B (v ⊥,φ ).For each map, the minimum and maximum of the plotted quantity at the surface are shown.

Figure 3 .
Figure 3. Evolution of the photospheric magnetic field and the open flux boundaries over one month in the time-evolving MHD model.The left hand frames show B r , derived from the ESFAM flux transport model, as a latitude-longitude map.These are the boundary maps for the calculation.The right-hand frames show open/closed boundaries (coronal holes) in the same format, with dark red indicating outward (positive) polarity, dark blue showing inward (negative) polarity, and white indicating closed-field regions.Five time instances are shown, from top to bottom: t = 0, t ≃ 180 hrs, t ≃ 360 hrs, t ≃ 540 hrs, and t ≃ 720 hrs.The box on every frame (black on the left, gold on the right for visibility) highlights a persistent open field region, while the magenta circle near the center of each frame indicates a bipole which undergoes decay.Other rectangles and circles denote more transient open field and flux emergence regions, respectively.A 7-second animation of this figure is available online (www.predsci.com/corona/tdc/animations/RL_2023_Fig3.mp4)showing both maps for the duration of the simulation time.

Figure 4 .
Figure 4.The emission in the AIA 171 Å channel over one month in the time-evolving MHD model.The left hand frames show the Sun from the point of view of an observer on Earth, as the star completes a full rotation around its axis.The right-hand frames show a projection of the emission as a latitude-longitude map; the annotations are analogous to those from the previous figure.Five time instances are shown, corresponding to those of Fig. 3: (a) t = 0 (b) t ≃ 180 hrs (c) t ≃ 360 hrs (d) t ≃ 540 hrs (e) t ≃ 720 hrs.A 7-second animation of this figure is available online (www.predsci.com/corona/tdc/animations/RL_2023_Fig4.mp4),showing both visualization styles for the duration of the simulation time.

Figure 5 .
Figure 5.Time series of time-dependent model evolution over one month, presented as latitude-longitude maps of B r overlaid with forwardmodeled log-scaled synthetic SDO AIA 193 Å emission from the temperature and density of the model (Lionello et al. 2009).The black boxes indicate the same persistent open field region throughout the simulation; green boxes illustrate several other open field areas.The red circles indicate a persistent bipole, which decays during the simulation.The cyan circle shows a bipole which emerges during the course of the simulation.A 7-second animation of this figure is available online (www.predsci.com/corona/tdc/animations/RL_2023_Fig5.mp4),showing this map evolving over the full simulation time.

Figure 6 .
Figure 6.Dipped magnetic field lines (field lines that have a reversal of the sign of B r ) at t = 524 hrs.(a) Footpoints of dipped field lines at r = 19R ⊙ as latitude-longitude map.(b) Latitude-longitude map of s-log Q at the same height.(c) Enlargement and superposition of a region of (a) and (b) where the dipped field-lines are associated with the current sheet.(d) Enlargement and superposition of another region of (a) and (b) where the dipped field lines are associated with a single-polarity separatrix line.(e) The r = 19R ⊙ semi-transparent surface colored as in (a) with dipped field-lines footpoints and marked with the (c) and (d) regions.Representative field lines are visible and the solar surface is at the center.(f) The same as (e) but colored as in (b) with s-log Q. (g) 3D enlargement of the (c) region with magnetic field lines.The solar surface is at the bottom right.(h) 3D enlargement of the (d) region with magnetic field lines.

Figure 7 .
Figure 7. (a) Statistical distributions of the maximum value of the C 06 /C 04 ratio evaluated along magnetic field lines in the time-dependent corona and steady-state corona models.We show distributions for loops (i.e., closed-field regions), coronal holes (open-field regions), loops at the edge of open-field areas, and edges of coronal holes.(b) The same as (a) with a logarithmic scale for the axes.(c) The areas over which the loops (violet) and coronal-holes (green) distributions are defined.(d) The same as (c) for the edge loops (cyan) and edge coronal-holes (gold) distributions.
This work was supported by the NASA Heliophysics Living With a Star Science and Strategic Capabilities programs (grant numbers 80NSSC20K0192, 80NSSC22K1021, and 80NSSC22K0893), the NASA Heliophysics System Observatory Connect program (grant number 80NSSC20K1285), and the NSF PREEVENTS program (grant ICER1854790).Computing resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center and by the Expanse supercomputer at the San Diego Supercomputing Center through the NSF AC-CESS and XSEDE programs.MLD acknowledges support from NASA under contract NNG04EA00C as part of the AIA science team.