Your Federal Quarterly Tax Payments are due April 15th

# Wu.ST.2007.USTC by xiagong0815

VIEWS: 7 PAGES: 65

• pg 1
```									     Coupling the Photosphere and Corona
using a Data Driven 3D MHD Model

S. T. Wu and A. H. Wang
1Center   for Space Plasma & Aeronomic Research
The University of Alabama in Huntsville
Huntsville, Alabama 35899 USA

Presentation at SEC/NWS/NOAA, Boulder, CO, February 2, 2007
Table of Contents
I.     Motivation
II.    MHD Model and Boundary Conditions
III.   Numerical Simulation Procedures
IV.    Results
V.     Concluding Remarks
I. Motivation
Most modeling work of the initiation and propagation of solar disturbances
impose a theoretical perturbation to mimic the observations. Now, there
are many valuable high resolution and reasonably high cadence
observations on the source surface. To be relevant to space weather
operations, modelers need to take these observations to drive the models.
To incorporate the observations at lower boundary is not trivial. We need
to develop a set of self-consistent, time-dependent, magnetohydrodynamic
(MHD) boundary conditions. In this presentation, we present a scheme
using the “method of characteristics” to achieve this goal. To illustrate this
procedure, we have studied the evolution of active regions and their
consequences to produce the solar eruptions.
II. MHD Model and Boundary Conditions
This model consists of two parts:
1.    Global configuration - Photosphere
2.    Local Configuration – Corona Coupling

1

2
II.1. Mathematical Model

∂ρ
Continuity:          + ∇ ⋅ (ρ u ) = 0
∂t

Momentum: ρ ⎡ ∂ u + u ⋅ ∇u ⎤ = −∇p + 1 (∇ × B )× B + F − 2 ρω × u− ρωo × (ωo × r ) + Ψ
⎢ ∂t           ⎥        4π
g            o
⎣              ⎦                            Coriolis force
Inertial centrifugal force   Viscous effect

Viscous dissipative                                  [                ]
Ψ = − ∇(μ t ∇ ⋅ u ) + μ t ∇ 2u + ∇(∇ ⋅ u ) + 2[(∇μ t ) ⋅ ∇ ]u + [(∇μ t ) × (∇ × u )]
2
function:                      3

∂p                                              ⎡       μ          2⎤
+ u ⋅ ∇p + γ p∇ ⋅ u = (γ − 1)∇ ⋅ Q + (γ − 1) ⎢η J 2 + t (∇ ⋅ u ) ⎥
Energy:            ∂t                                              ⎣        2          ⎦

∂B
Induction:                ∂t
(   ) (          )
= ∇ × u × B + λ ∇ × B + η∇ 2 B + S
where,

ρ: the plasma density                 ω : the angular velocity of solar
differential rotation
u:   the plasma flow velocity         Fg : gravitational force

p: the plasma thermal pressure
∇⋅ Q: heat conduction
B : the magnetic induction vector
μt: turbulent viscosity
J : electric current
λ: coefficient of the cyclonic
S:   source function                     turbulence

γ: specific heat ratio                 η: magnetic diffusivity
Data Driven Three-Dimensional Magnetohydrodynamic Model

Using observations, a data driven 3-D (MHD) model is developed to investigate the
photosphere-corona coupling which leads to quantifying the criteria for “initiation” of solar
eruptive phenomena. The important innovation of this model is its ability to incorporate
the solar interior (convective zone) dynamics by inputting magnetic measurements (i.e.
emerging and submerging magnetic flux) together with the differential rotation, meridional
flow, coriolis force, effective diffusion, and cyclonic turbulence effects to drive the model.
The time-dependent projected characteristic boundary conditions are fully
implemented. Thus, the effects caused by sub-photospheric dynamics (i.e. convective
zone dynamics) can be reflected through the surface boundary conditions in a new way.
Hence, this model is the first MHD model to include the challenge of the coupling
between the solar interior and corona, because the photospheric measurements
already exhibits the effects of convective zone dynamics.
II.2. Boundary Conditions

A. Global Configuration (Spherical Coordinates)

In the case of global configuration the spherical coordinate system is
used. The computational domain is
0 ≤ φ ≤ 360°, -90° ≤ θ ≤ +90° and 1Rs ≤ r ≤ 30Rs.
Grids: 76 × 32 × 32
II.2. Boundary Conditions

B. Local Configuration (Rectangular Coordinate)

In the case of local configuration, the set of governing equations
has been casted in the rectangular coordinates. The computational
domain includes six planes (i.e. four sides, top and bottom planes)
for rectangular coordinates. The boundary conditions used for the
four sides are linear extrapolation, and the top boundary is non-
reflective. The surface boundary conditions are time-dependent
boundary conditions which are derived from the method of
projected characteristics (Wu and Wang, J. Comp. Methods in
Applied Mech. & Engineering, 64, 257, 1987)

A detailed description of the boundary conditions can be found at
Ap. J., 652, 800, November 2006.
Surface Boundary Conditions
In order to accommodate the observations at the bottom boundary, the evolutionary
boundary conditions must be used, thus, the method of “projected” characteristics are used
for the derivation of such boundary conditions. For a three-dimensional
magnetohydrodynamic problem, the governing equations have been cast in a vectoral form,

∂U       ∂U     ∂U     ∂U
= −A     −B     −C     +S
∂t      ∂x1    ∂x2    ∂x3
where U and S and are column vectors consisting of primary physical quantities such as,
density, temperature and three components of velocity and magnetic field. To obtain these
characteristics on the boundary, we have chosen the unit normal on the boundary is along
certain coordinates, e.g. along the radial/z direction. Thus, the characteristics along the
projected normal will be found in the r-t or z-t plane. For example, at z-t plane for a small
time difference Δt the projected characteristics given by:
dz
= λi           i=1,2,…8
dt

can be approximated by straight lines.
The Projected Characteristics Straight Lines
t

P
(p+1)Δt

pΔ
t                      L1     L2   L3   L4,5          L6        L7     L8

x
(l-1,m,n)                  (l,m,n)                       (l+1,m,n)

where λ i = (u z , u z , u z + VA , u z − VA , u z + V f , u z − V f , u z + Vs , u z − Vs ) are eigenvalues. The projected characteristics
passing the point P at the time (p+1)Δt and the spatial location (l,m,n) then intersect the pΔt axis
at L1, L2,L3,L4,5,L6,L7 and L8. For f > VA > Vs > u z > 0
V                                         , the projected characteristic PL1can be
identified u z + V f
λ 5 = with                             λ3 2 with
, PL= u z + VA                              , PL3 = u z + Vs
λ 7 with                   λ1, = u z
, PL4,52 with          ,
λ 6 with
PL8 = u z − Vs                        λ 4 = u z, − VA                     = uz
PLλ 6with− V f                  , PL8 with                     .
7
With their left-eigenvectors, the projected normal characteristic equations (i.e. compatibility
equations) are :

∂U           ∂U          ∂U       ∂U
ηj       + λ jη j    = −η j A    − ηjB    + η jS              j=1,2,…8.
∂t          ∂z          ∂x       ∂y

They are as follows:
∂ρ ∂p
dz                       a2     −   = − a 2 u ⋅ ∇ρ + u ⋅ ∇p
= uz                       ∂t ∂t
dt
∂Bz                  ∂u   ∂u y       ∂u      ∂u
= −u ⋅ ∇Bz − Bz ( x +      ) + Bx z + B y z
∂t                   ∂x   ∂y         ∂x      ∂y

dz                                   ∂u x         ∂u y         ∂B         ∂B y
= u z + VA             − B y Bz        + Bx Bz      + B yV A x − BxV A      = A+
dt                                    ∂t           ∂t           ∂t         ∂t

∂u x         ∂u y         ∂B         ∂B y
dz
= u z − VA             + B y Bz        − Bx Bz      + B yV A x − BxV A      = A−
dt                                    ∂t           ∂t           ∂t         ∂t

∂u x              ∂u y                  2 ∂u              2 ∂p         ∂B           ∂B y
dz
= uz + V f             − B x B zV f        − B y B zV f      + ρV f (V f2 − V A ) z + (V f2 − V A ) + B xV f2 x + B yV f2      = B+
dt                                        ∂t                ∂t                      ∂t                ∂t         ∂t           ∂t

∂u x              ∂u y                  2 ∂u              2 ∂p         ∂B           ∂B y
dz
= uz − V f             + B x B zV f        + B y B zV f      − ρV f (V f2 − V A ) z + (V f2 − V A ) + B xV f2 x + B yV f2      = B−
dt                                        ∂t                ∂t                      ∂t                ∂t         ∂t           ∂t

dz                                       ∂u x             ∂u y                2 ∂u             2 ∂p        ∂B          ∂B y
= u z + Vs             + B x B zV s        + B y B zVs      − ρVs (Vs2 − V A ) z − (Vs2 − V A ) − B xVs2 x − B yVs2      = C+
dt                                        ∂t               ∂t                    ∂t               ∂t        ∂t          ∂t

∂u x             ∂u y                2 ∂u             2 ∂p        ∂B          ∂B y
dz
= u z − Vs             − B x B zV s        − B y B zVs      + ρVs (Vs2 − V A ) z − (Vs2 − V A ) − B xVs2 x − B yVs2      = C−
dt                                        ∂t               ∂t                    ∂t               ∂t        ∂t          ∂t
For the bottom boundary, if the eigenvalue (i.e. wave speed) is negative, that means the
boundary condition will be affected by computational domain, and we have to use the
compatibility equation to determine the physical parameter. The number of compatibility
equations that have to use is determined by the number of negative eigenvalues.
As an example, if            uz > 0   and            u z ≤ V A , Vs , V f ,   then only   u z − VA ,   u z − Vs , u z − V f ,   will be negative. In
this case, three out of eight physical parameters have to be determined by three compatibility
equations, and five could be given or using compatibility equations to solve. If relaxing the
condition above to only                    u z ≤ V A , Vs , V f      , then      uz   (note that this is two degenerated eigenvalues),
u z − VA ,   u z − Vs ,    uz − V f ,   could be negative. In this case, five out of eight physical parameters
have to be determined by five compatibility equations, and three could be given or using
compatibility equations to solve.
For first case, expressions which describe the physical parameters of pressure, density, the
components of velocity, and magnetic field vary with time on the boundary are:
∂ p Vs B− + V f C −
2         2
=
∂ t 2V A V f2 − Vs2
2
(                )
p
γ
= const .
ρ

=
( (
∂u x B y Vs V A − Vs B− − V f V f − V A C −
2     2            2     2
)   +
Bx A−     (                ) )
∂t                          (
2 Bz Bx + B y VsV f V f − Vs
2    2         2     2
)           (
2 Bz ( Bx + B y
2     2
)                   (                   )
∂u y
=
( (                       )
Bx Vs V A − Vs2 B− − V f V f2 − V A C −
2                         2
(                ) )− B A               y       −
∂t                          (
2
2 Bz Bx      +    2
By      )V Vs   f   (V   2
f   − Vs2   )    2 B (( B + B )
z
2
x
2
y

=
(
∂u z − Vs B− + V f C−                       )
∂t   2ρVsV f V f2 − Vs2 (                  )
=
((
∂Bx Bx V A − Vs B− − V f − V A C −
2     2       2     2
+
)    B y A−  (                   ) )
∂t    2V A V f − Vs Bx + B y
2  2     2   2
(2
)(
2V A Bx + B y
2    2
)                  (                    )
∂B y
=
((
B y V A − Vs2 B− − V f2 − V A C −
2
)   2
(            ) )− B A                     x −
∂t                    2
2V A     (
V f2   − Vs2       ) (B    2
x   +B )   2V (B + B )
2
y              A
2
x
2
y

∂B z       ∂u     ∂u      ∂B       ∂u y     ∂u       ∂B y
= − Bz x + Bx z − u x z − B y      + By z − u y
∂t         ∂x     ∂x      ∂x       ∂y       ∂y       ∂y
where the coefficients A_, B_, and C_ are given below.

⎛       ∂u         ∂u y         ∂B         ∂B y ⎞
A− = −(u z − V A )⎜ B y Bz x − Bx Bz
⎜                       + B yV A x − BxV A      ⎟
⎝        ∂z         ∂z           ∂z         ∂z ⎟⎠

(
+ Bx B yV A − u x B y Bz   ) ∂∂uxx − (Bx2VA − Bx By u x ) ∂∂uxy − ByρBz ∂p
∂x

2 ∂B y            ∂B y
(          )
2
Bz 2                                         ∂Bx B y Bz ∂Bz
−   Bx + B y      + BxV Au x      − u x B yV A     −
ρ           ∂x              ∂x                ∂x   ρ ∂x

(
+ B yV A − u y B y Bz
2
) ∂∂uy − (B B V
x
x   y A   − Bx B z u y   ) ∂∂uyy + BxρBz ∂p
∂y

∂B y Bx Bz ∂Bz
(     2 ∂Bx
)           ∂B                     2
Bz 2
+     Bx + B y     − u y B yVA x + u y BxVA     +
ρ           ∂y             ∂y            ∂y   ρ ∂y
∂u y
(         ⎛
)
B− = − u z − V f ⎜ Bx BzV f
⎜
∂u x
∂z
+ B y BzV f
∂z
2 ∂u z
+ ρV f V f2 − V A
∂z
(               )
⎝

∂Bx           ∂B y ⎞
(
+ V f2 − V A2          ) ∂p + B V
∂z
x       f
2

∂z
+ B yV f2
∂z ⎟
⎟
⎠

(                         2
(
− Bx BzV f u x + a 2ρ V f2 − V A + B yV f2
2
)               ) ∂∂ux x
(
− B yV f u x Bz + BxV f              ) ∂∂uxy

⎛ Bx BzV f
(
+ ρ u xV f V f2 − V A2         ) ∂∂ux      z
−⎜
⎜ ρ        + u x V f2 − V A2(               )⎞ ∂p − u B V
⎟
⎟ ∂x        x   x       f
2   ∂Bx
∂x
⎝                                            ⎠

∂B y                      ∂Bz                             ∂u x
− u x B yV f2
∂x
− BxV f3
∂x
(
+ Bx B yV f2 − u y Bx BzV f
∂y
)

(                               (
− u y B y BzV f + a 2ρ V f2 − V A + Bx V f2
2    2
)                 ) ∂∂uyy
(
+ ρV f V f2 − V A u y
2
)        ∂u z
∂y

⎛ B y B zV f                                                                                     ∂B y
−⎜
⎜ ρ                        (
− u y V f2 − V A2                )⎞ ∂p − u
⎟
⎟ ∂y          y   B xV f2
∂B x
∂y
− u y B yV f2
∂y
⎝                                             ⎠

∂Bz
− B yV f3
∂y
+ ρ gV f V f2 − V A2 (               )
∂u y
⎛
C− = (u z − Vs ) ⎜ Bx BzVs x + B y BzVs
⎜
∂u
∂z            ∂z
2 ∂u z
+ ρVs Vs2 − V A
∂z
(         )
⎝

⎛ ∂Bx        ∂B y ⎞
(
+ Vs2 − V A2   ) ∂p + V
∂z
s
2
⎜ Bx
⎜    ∂z
+ By
∂z ⎟
⎟
⎝                 ⎠

∂u y
( (                     )
+ a 2ρ Vs2 − V A + B yVs2 + Bx B yVs u x
2     2
) ∂∂ux    x
+ B yVs (u x Bz + BxVs )
∂x

(
− ρ u xVs Vs2 − V A
2
) ∂∂uxz     ⎛
(             )
B B V ⎞ ∂p
⎟ ∂x + u x BxVs ∂x
+ ⎜ u x Vs2 − V A + x z s ⎟
⎜
2
ρ ⎠
2 ∂Bx

⎝

∂B y
+ u x B yVs2
∂x
+ BxVs3
∂Bz
∂x
(
− Bx B yVs2 − u y Bx BzVs
∂u x
∂y
)

( (                  )
+ a 2ρ Vs2 − V A + Bx Vs2 + u y B y BzVs
2    2
) ∂∂uy y
(
− ρVs Vs2 − V A u y
2
)     ∂u z
∂y

∂B y
⎛ B y BzVs
+⎜
⎜ ρ                   (
+ u y Vs2 − V A
2
)⎞ ∂p + u B V
⎟
⎟ ∂y       y   x s
2     ∂Bx
∂y
+ u y B yVs2
∂y
⎝                                        ⎠

∂Bz
+ B yVs3
∂y
− ρ gVs Vs2 − V A2(            )
Bz
Alfvén speed                VA =
4πρ

(
1⎧ 2
) ((    )                    )      ⎫
2                1/ 2

Fast MHD wave speed   V f2 =    ⎨ a +b + a +b
2   2   2
− 4a 2V A2          ⎬
2⎩                                              ⎭

(
1⎧ 2
) ((    )                    )      ⎫
2                    1/ 2
Slow MHD wave speed   Vs2 =     ⎨ a +b − a +b
2   2   2
− 4a 2V A2          ⎬
2⎩                                              ⎭

B x2 + B y + B z2
2
with                  b=
4πρ

Sound speed            a = γ RT
II.3. Numerical Methods

The numerical method we used is simple TVD Lax-Friedrichs formulation. This scheme
achieves the second order accuracy both temporally and spatially. To achieve second
order temporal accuracy, the Hancock predictor step and corrector step are used.

Predictor Step:
⎧ F (U in j ,k + 1 δU in ) − F (U in j ,k − 1 δU in ) ⎫
⎪                                                     ⎪      ⎧ G (U in j ,k + 1 δU n ) − G (U in j ,k − 1 δU n ) ⎫
⎪                                               j ⎪
U in+,2k   = U in j ,k   − Δt ⎨                                                     ⎬ − Δt ⎨
1
,        2                 ,        2                        ,        2    j           ,        2
⎬
Δx1                                                         Δx2
,j           ,
⎪
⎩                                                     ⎪
⎭      ⎪
⎩                                                   ⎪
⎭

⎧ H (U in j ,k + 1 δU k ) − H (U in j ,k − 1 δU k ) ⎫ Δt
⎪
n                         n
⎪
− Δt ⎨                                                   ⎬+
,                         ,
2                         2
S (U in j ,k )
Δx3
,
⎪
⎩                                                   ⎪ 2
⎭

Corrector Step:

U in+,1k = U iT j ,k + 1 (Φ iF − Φ iF ) + 1 (Φ G+ − Φ G− ) + 1 (Φ k + − Φ k − ) +                   Δt
S (U in+,k )
2
H       H                                       2
,j         ,        2     +      −
1
2     2
1
2
j      j    1
2 2
1
2
1
2
1
2       2          ,j

U iT j ,k = U in j ,k + Δt[ 1 ( Fi + − Fi− ) + 1 (G LR − G LR ) + 1 ( H kLR − H kLR )]
,           ,            2
LR
1
2
LR
2
1
2
j+     j−  1
2  2
1
2
+       −    1
2
1
2
Powell’s Corrective Terms
Powell discovered that including ∇⋅B corrective terms and the corresponding
characteristic divergence wave, can stabilize the solution for the TVD type algorithms.
In our equations, the source terms include the following corrective terms:

Global                               Local
⎡        0        ⎤                   ⎡        0       ⎤
⎢                 ⎥                   ⎢                ⎥
⎢ − Br ∇ ⋅ B ⎥                        ⎢ − Bx ∇ ⋅ B ⎥
⎢ − B ∇⋅B ⎥                           ⎢ − B ∇⋅B ⎥
⎢       θ         ⎥                   ⎢       y        ⎥
⎢ − Bϕ∇ ⋅ B ⎥                         ⎢ − B ∇⋅B ⎥
S′ = ⎢                 ⎥              S′ = ⎢                ⎥
z

⎢ − ur ∇ ⋅ B ⎥        or              ⎢ − u x∇ ⋅ B ⎥
⎢                 ⎥                   ⎢                ⎥
⎢ − uθ∇ ⋅ B ⎥                         ⎢ − u y∇ ⋅ B ⎥
⎢                 ⎥                   ⎢                ⎥
⎢ − u ϕ∇ ⋅ B ⎥                        ⎢ − uz∇ ⋅ B ⎥
⎢                 ⎥                   ⎢                ⎥
⎣− (u ⋅ B )(∇ ⋅ B)⎦                   ⎣− (u ⋅ B)(∇ ⋅ B)⎦
Trial Values

Iteration                                     MHD
Equilibrium
F                                   T                     Solution
Error Check

Initial State

Drivers:
Bottom Bdy Conditions                Magnetic field
measurements
Differential rotation ω0
Meridinal flow
Predict Step

Top & Side Bdy Conditions

Evolution
Bottom Bdy Conditions               Solution

Correction Step
Numerical Code Flow Chart
Top & Side Bdy Conditions
Computational Flow Chart for the 3D MHD
Artificial Dissipation                                            code where “F” and “T” represent the “false”
and “true”, respectively. Note, the upper box
F                               T                                        represents the code to compute the
Time =TSAVE                    Save Data
equilibrium solution and the lower box is for
computing the evolutionary solution.
Time                      Terminating Code
≥TSTOP
III. Procedures for Carrying Out the Numerical Simulation

III.1.      Establishing the Observational MHD equilibrium state
To obtain the magnetohydrodynamic equilibrium state to represent the observational state based on
the observed quantities as the initial state for the evolutionary study comprises three steps;

(a) to construct a trial 3D magnetic field configuration using pre-event magnetograms,

(b) since there are no density distribution measurements on the photosphere yet. We could input
any trail values for density and temperature, according to the relaxation technique (Steinolfson,
et al. 1982; Roache, 1998) the equilibrium state will be obtained. However, if the trail values are
close to the equilibrium state, then the final equilibrium solution would be rapidly obtained. Thus,
we choose the trail values for density distribution at the photospheric level is directly proportional
to the absolute value of the magnitude of the transverse field and decreases exponentially with
z
−
Bx + B y
2     2

the scale height; thus ρ ( x, y, z,0) = ρ o
Hg
e
2
Bo                    where ρo and Bo are the constant reference with
⎛ RT ⎞
values Hg as the    scale height, ⎜ o ⎟ .
⎜ g ⎟
⎝ o ⎠
In the case where Bx and By is equal to zero (i.e. purely radial field), the value of density is
chosen to be the average of the four neighboring points. This assumption does bear some
physical meaning at the level of the chromosphere/corona interface. That is, usually the
observed brighter feature on the surface corresponds with the higher density that also
corresponds with the closed magnetic field (Pneuman and Orrall, 1986). It is also worth noting
that the conditions (a) and (b) are the trial values.
III.1. Establishing the Observational MHD equilibrium state (cont.)

(c) by inputting these trial values of (a) and (b) into the MHD model described in section 2.1
without differential rotation and meridional effects to allow the atmosphere to relax to a
magnetohydrodynamic quasi-equilibrium state. This procedure is indicated on the top
portion of the numerical code flow chart.
III.2. Evolution of the active region using photospheric measurements

To evolve the active region corona, we will input the photospheric measurements according to time-
dependent bottom boundary conditions. For example, the magnetic field variation at the bottom boundary
is described by a set of magnetograms. If the time cadence of the set is the same or smaller than the time
step of the numerical code (CFL condition) then,

Bbottom ( x, y, t ) = Bob (x, y, t )
If the time cadence is larger than the time step of the numerical code, then between two consecutive
magnetogram values, for each numerical time step, we assume the linear variation from one magnetogram
to the next. For example, the magnetogram cadence is p min, and the time step is q sec, then
Bob ( x, y, t n +1 ) − Bob ( x, y, t n )
Bbottom (x, y, t ) = Bob (x, y, t n ) + k
60 p
, k = 1,2,...
60 p q                                    q

where   Bob    is the measured magnetic fields, and subscript “n” represents the cadence of the

magnetogram. The magnetograms are deprojected and ambiguity is resolved using Metcalf methods
(Metcalf et al. 2006). For other accessible measurements (i.e. density, velocity, and temperature) a similar
expression the equation above applies.
Survey of Transport Coefficients

Before we can carry out the simulation study, we need to know two
coefficients; effective diffusivity (η) and cyclonic turbulence (λ). There
are no precise theories, observations nor laboratory experiments to
determine these coefficients. However, there are some previous
works which have discussed the choice of these two coefficients. For
example, η = 160 – 300 km2 s-1 given by Parker, (1979); Leighton’s
value of η is 800 – 1600 km2 s-1 (1964); DeVore, et al (1985) selected
η = 300 km2 s-1 for their study. Wang (1988) derived a value of η
being 100 – 150 km2 s-1 on the basis of observation of sunspot decay.
We notice that there is a wide range of values for the effective
diffusivity. The value of cyclonic turbulence is chosen according to the
scale law (λ ~ η/L), given by Parker, (1979) where L is the
characteristic length of the sunspot, it is chosen to be 6,000 km for this
study and η is 200 km2 s-1.
IV. Results
IV.1. Global Coronal Model
MHD Equilibrium State
• Solar Wind
• Density
• Temperature
Pre-Initial Field
Observed Radial Field (CR2000)
60

40

20

Latitude
0

-20

-40

-60

100           200         300
Longitude

Potential Radial Field
60

40

20
Latitude

0

-20

-40

-60

100           200         300
Longitude
Steady-State
Radial Magnetic Field
60

40

20

Latitude
0

-20

-40

-60

100            200          300
Longitude

Transverse Magnetic Field

50
Latitude

0

-50

0    100           200          300
Longitude
Steady-State
Density Distribution
60

40

20

Latitude
0

-20

-40

-60

100           200           300
Longitude

Transverse Velocity Field

50
Latitude

0

-50

0   100           200           300
Longitude
Radial Velocity Distribution With Radii

Radial Velocity
500

400                   Coronal Hole
Velocity(km/s)

300

200

Coronal Hole Boundary
100

0
1      2       3      4           5            6
Solar Radii
IV. Results cont.. (Rectangular Coordinate Code)

Initiation Processes: Active Region Evolution
Evolutionary simulation of an Active Region (AR)

To carry out this simulation study, we have chosen the
SOHO/MDI magnetic field measurements of NOAA/AR 8100
from October 29, 11:15 UT to Nov 3, 15:59 UT, 1997 for this
study. The SOHO/MDI field measurements of the active region
have a resolution of ~ 2 arc sec with 198 × 198 pixels with a
cadence of ~ 96 min. In order to assure the computational grids
are compatible with the measurements, the computation domain
are 99 in longitudinal direction (x), 99 in latitudinal direction (y)
and 99 in height (z), respectively. To match the data with the
grids, we have taken a four point average of the pixels inside the
domain. On the boundary we have taken a two point average
from the measurements. At the four corners, the measurements
are used.
The simulated initial (MHD Equilibrium) state of AR8100 at various times
Bz Contours & Transverse Bt                                            Bz Contours & Transverse Vt
10

20                                                                     20
5

Y (3.95arcsecs)

Y (3.95arcsecs)
15                                                                     15
0

-1.0

-1.0
-4.0                                                                             -4.0
10                                                                     10
2.0.5                                                                   2.0.5

-1.0

-1.0
-5

6.0

6.0
5                                                                      5

0                                                                      0                                                                     -10
0            10            20       30                                 0             10            20       30
X (3.95arcsecs)                                                         X (3.95arcsecs)
(a)                                                                     (b)

Bt & Density Contours(109cm-3)                                              Plasma Beta Contours
0.     1.00
0.80.9              0.8
9
0.8
20                                                                     20                                    0.6
0.75
Y (3.95arcsecs)

Y (3.95arcsecs)

0.
4
15                                                                     15
0.50

5
0.0
10                                                                     10                                          0.1
0.05
0.25
5                                                                      5

.8
0.4
0.6
0.9 0
0.4

0.2
0.1
0.6
0                                                                      0                                                                    0.00
0            10            20       30                                 0             10            20       30
X (3.95arcsecs)                                                         X (3.95arcsecs)
(c)                                                                     (d)
The simulated initial state of AR8100 at 14:27 UT 1997, Oct, 31; (a) the transverse magnetic field vector (5G≤⏐Bt⏐ ≤ 6G) and contours of the line-of-sight magnetic field (Bz) with the solid
lines and broken lines representing the positive and negative polarity, respectively. The color bar on the upper right side indicates the strength of the line-of-sight magnetic field
(-10G≤⏐Bz⏐ ≤ 10G) contours, (b) the transverse velocity (maximum is 1.9 km s-1) and Bz contours, (c) the density contours at surface with transverse magnetic field, and (d) the plasma
beta (β =16πnkT/B2) distribution on the surface. The color bar at the lower right side is for both density and β contours.
Simulated three-dimensional magnetic field configuration (at magnetohydrodynamic equilibrium) of
AR8100 at 14:27 UT, Oct 31, 1997.
Magnetic Field Evolution

Time = 14:27UT                                                     Time = 16:03UT
10

20                                                                  20

Y (3.95arcsecs)

Y (3.95arcsecs)
15                                                                  15
-4.0

0
-1.0

-1.
-4.0
10                                                                  10
2.0.5                                                           0 0.5

-1.0
2.                                         5

-1.0
6.0
6.0
5                                                                   5

0                                                                   0
0            10            20       30                              0          10            20       30
X (3.95arcsecs)                                                   X (3.95arcsecs)

0

Time = 17:39UT                                                     Time = 19:12UT

20                                                                  20
Y (3.95arcsecs)

Y (3.95arcsecs)
15                                                                  15                                               -5
-1.0
-1.0
0.5

10                                                                  10    0.5
2.0                                                                 2.0

-1.0
-1.0

5                                                                   5
6.0

6.0
0                                                                   0                                              -10
0            10            20       30                              0          10            20       30
X (3.95arcsecs)                                                   X (3.95arcsecs)

The simulated evolution of magnetic field at 14:27 UT, 16:03 UT, 17:39 UT and 19:12 UT Oct 31, 1997. The representation is similar to Figure 1. The color bar on the right-hand side
indicates the strength of LOS magnetic field. The white arrows represent the non-potential transverse magnetic field, and black arrows represent the potential transverse magnetic field.
Transverse Velocity Evolution (Recovering Photospheric Velocity using MHD model)
Time = 14:27UT                                                    Time = 16:03UT
10

20                                                                  20

Y (3.95arcsecs)

Y (3.95arcsecs)
15                                                                  15
-4.0

0
-1.0

-1.
-4.0
10                                                                  10
2.0.5                                                          0 0.5

-1.0
2.                                                5

-1.0
6.0
6.0
5                                                                  5

0                                                                  0
0            10            20       30                             0          10            20       30
X (3.95arcsecs)                                                  X (3.95arcsecs)

0

Time = 17:39UT                                                    Time = 19:12UT

20                                                                  20
Y (3.95arcsecs)

Y (3.95arcsecs)
15                                                                  15                                                    -5
-1.0
-1.0
0.5

10                                                                  10   0.5
2.0                                                                2.0

-1.0
-1.0

5                                                                  5
6.0

6.0
0                                                                  0                                                   -10
0            10            20       30                             0          10            20       30
X (3.95arcsecs)                                                  X (3.95arcsecs)

The simulated evolution of surface transverse velocity vector (Vt), and the contours of the line-of-sight magnetic field for AR 8100 on Oct 31, 1997; (a) at 14:27 UT with 0.0002 kms-1
≤ |Vt| ≤ 1.9 km s-1 (b) at 16:03 UT with 0.0018 kms-1 ≤ |Vt| ≤ 3.7 kms-1, (c) at 17:39 UT with 0.0088 kms-1 ≤ |Vt| ≤ 5.0 kms-1 and (d) at 19:12 UT with 0.0164 kms-1 ≤ |Vt| ≤ 7.1 kms-1.
Corresponding Vertical Velocity Evolution at the surface

Time = 14:27UT                                                         Time = 16:03UT
0.3

20                                                                      20

Y (3.98arcsecs)

Y (3.98arcsecs)
15                                                                      15
-4.0
0.2

0
-1.0

-1.
-4.0
10                                                                      10
2.0.5                                                               0 0.5

-1.0
2.

-1.0
6.0
6.0
5                                                                       5
0.1
0                                                                       0
0            10           20            30                              0          10            20         30
X (3.98arcsecs)                                                       X (3.98arcsecs)

0.0
Time = 17:39UT                                                         Time = 19:12UT

20                                                                      20                                                -0.1
Y (3.98arcsecs)

Y (3.98arcsecs)
15                                                                      15                  -1.0
-1.0
0.5

10                                                                      10    0.5                                         -0.2
2.0                                                                     2.0

-1.0
-1.0

5                                                                       5
6.0

6.0
0                                                                       0                                                -0.3
0            10           20            30                              0          10            20         30
X (3.98arcsecs)                                                       X (3.98arcsecs)

The evolution of the vertical velocity (km s-1) (color coded on the right-hand side) of AR 8100 at the same four times as Figure 4.
The magnetic energy (1022
erg/km2) across the low
boundary of the AR8100 at
four times (14:27UT,
16:03UT, 17:39UT and
19:12UT).
Simulated energy flux through the photosphere for the Active Region AR8100.
4

Total energy flux
Energy flux through the photosphere                                                          Surface flow effect
⎛ dE ⎞
⎜     ⎟      =
1
dt ⎠ photo 4π
(       )       1
∫photo Bt ⋅ Vt Bn dS + 4π ∫photo Bt ⋅ Vn dS
2
3    Emerging flux effect
⎝
where the Bt and Vt are the transverse magnetic
field and velocity, respectively. Vn is the normal

dE/dt (1029erg/s)
velocity. On the right side of the equation above,
the first term represents the surface flow effect and
the second term is due to the emerging flux. Helicity
flux across the photosphere is determined by                                            2
⎛ dH ⎞
⎜    ⎟       = −2 ∫    (        )
A p ⋅ Vt Bn dS + 2∫    (        )
A p ⋅ Bt Vn dS
⎝ dt ⎠ photo       photo                    photo

where Ap is the vector potential to be determined by
using LOS field where Bn and Bt are the measured
quantity at the photosphere and Vt and Vn as a                                          1
function of time are the products of the MHD model.

0
0    1          2          3   4
Time(96min)
B. Variation of Non-Potential Magnetic Parameters for AR8100, AR8210
and AR 10646
III.2. Simulation of non-potential magnetic field parameters and CME
initiation.

Using the definitions given by Falconer et al. (Ap. J., 569, 1016, 2002), the
nonpotentiality parameters of AR8100 is computed based on the model outputs.
These non-potentiality parameters are:

Φ:                 The total magnetic flux content,

Lss:               The length of strong magnetic shear (>45º) and the
strong transverse field (>150 G) of the main neutral line.

LSG:               Strong Gradient Length (LOS magnetogram)

IN:                The net electric current, and
(α = μ I n   Φ) :
The flux normalized measure of the field twist.

Table I shows these parameters as a function of time and these values are
compared with values given by Falconer et al. (2002) at two specific times using
MSFC vectormagnetograph data. Reasonable agreement is reached.
B >0                                                                    B ≥ 100

All Fields (B > 0)                                              Strong Fields (B ≥ 100 gauss)
Magnetic Flux vs Time                                                 Magnetic Flux vs Time
40                                                                    40
observation                                                           observation
simulation                                                            simulation
30                                                                    30
Total Flux (1021Mx)

Total Flux (1021Mx)
20                                                                    20

10                                                                    10

0
0
10/29 10/30 10/31 11/1 11/2       11/3   11/4
Time (day)                                            10/29 10/30 10/31 11/1 11/2    11/3    11/4
(a)                                                                Time (day)
(b)
Δ

CMEs
Δ
Δ
Δ
Lss Evolution From 12:51UT to 17:39UT Oct 31, 1997 (AR 8100)

Time=12:51UT Oct 31                                                    Time=14:27UT Oct 31
1000

20                                                               20
Y(3.95arcsec)

Y(3.95arcsec)
15                                                               15
-400

-100

0
-40

-10
10                                 0                             10
20

-100
0                                                                                                        500

-100
200
5                                                                5
600

600
0                                                                0
0         10           20          30                            0           10         20        30
X(3.95arcsec)                                                    X(3.95arcsec)

0

Time=16:03UT Oct 31                                                    Time=17:39UT Oct 31

20                                                               20
Y(3.95arcsec)

Y(3.95arcsec)
15                                                               15                                              -500

0
-10
-100

10                                                               10
200                                                                   200

-100
-100

5                                                                5
600

600

0                                                                0                                             -1000
0         10           20          30                            0           10         20        30
X(3.95arcsec)                                                    X(3.95arcsec)
Non-Potential Parameters
AR8100
50
Lsg
Lss
flux
I
40                                              alfa
Normalized Units

30

20                                              Flux

10

0
0    200   400   600      800   1000   1200      1400
Time (min)
00:00 UT 11/3/97                                   24:00 UT 11/3/97

The solid vertical lines are represented 3 CMEs events. A weak CME occurred before this
computed period, and 8 CMEs occurred after this period of time.
AR10720
6
Lsg
Lss
flux
5
I
alfa

4

3

2

1

0
0                 20    40                60   80               100
Time (min)
22:09 UT 01/15/05                                      23:46 UT 01/15/05
AR 10720
250                                                                       250                                                                          250
22:09 UT 1/15/05
22:31 UT 1/15/05                                                           22:50 UT 1/15/05
200                                                                       200                                                                          200
Y(0.62arcsecs)

Y(0.62arcsecs)

Y(0.62arcsecs)
150                                                                       150                                                                          150

100                                                                       100                                                                          100

50                                                                        50                                                                           50

0                                                                         0                                                                            0
0   50   100   150        200   250   300    350                          0   50   100   150        200    250    300     350                          0   50   100   150        200   250    300    350
X(0.62arcsecs)                                                            X(0.62arcsecs)                                                               X(0.62arcsecs)

250                                                                       250                                                                          250

22:59 UT 1/15/05                                                          23:06 UT 1/15/05                                                             23:15 UT 1/15/05
200                                                                       200                                                                          200

Y(0.62arcsecs)
Y(0.62arcsecs)
Y(0.62arcsecs)

150                                                                       150                                                                          150

100                                                                       100                                                                          100

50                                                                        50                                                                           50

0                                                                            0
0
0   50   100    150        200    250     300     350                        0   50   100   150        200    250    300    350
0   50   100   150        200   250    300    350
X(0.62arcsecs)                                                              X(0.62arcsecs)
X(0.62arcsecs)
AR10646
100
Lsg
Lss
flux
80                                   I
alfa

60

40

20

0
0              500          1000          1500
Time (min)
16:42 UT 7/12/04                                 21:51 UT 7/13/04
AR 10646
300                                                                      300

16:42 UT 7/12/04                                                        18:54 UT 7/12/04
250                                                                      250

200                                                                      200

Y(0.62arcsecs)

Y(0.62arcsecs)
150                                                                      150

100                                                                      100

50                                                                       50

0                                                                       0
0   50   100        150         200   250   300                         0    50   100        150         200   250   300
X(0.62arcsecs)                                                           X(0.62arcsecs)

300                                                                      300

00:32 UT 7/13/04                                                          15:36 UT 7/13/04
250                                                                      250

200                                                                      200
Y(0.62arcsecs)

Y(0.62arcsecs)
150                                                                      150

100                                                                      100

50                                                                        50

0                                                                        0
0   50   100        150         200   250   300                          0   50   100        150         200   250   300
X(0.62arcsecs)                                                           X(0.62arcsecs)
AR10646
300                                                                     300

16:46 UT 7/13/04
250                                                                     250                 19:33 UT 7/13/04

200                                                                     200
Y(0.62arcsecs)

Y(0.62arcsecs)
150                                                                     150

100                                                                     100

50                                                                       50

0                                                                        0
0   50   100        150         200   250   300                          0   50   100         150          200    250    300
X(0.62arcsecs)                                                            X(0.62arcsecs)

300
300

21:27 UT 7/13/04
250                                                                                              21:51 UT 7/13/04
250

200
200
Y(0.62arcsecs)

Y(0.62arcsecs)
150
150

100                                                                     100

50                                                                       50

0                                                                        0
0   50   100        150         200   250   300                          0   50   100        150          200    250    300
X(0.62arcsecs)                                                           X(0.62arcsecs)
AR8210
35

30

25

20

15

Lsg
10                                          Lss
flux
I
alfa
5
0             50   100            150   200             250
Time (min)
17:00 UT 05/01/98                               21:15 UT 05/01/98
AR 8210, May 1, 1998

17:00 UT
Time=1                                                      17:19 UT
Time=2
40                                                                40                                                     1500
50                                                               50
50                                                                     50
30    50                                                          30 50
Y(1280km)

Y(1280km)
20                                                                20

50
0              50
-100                         450
10                                                                10
-100                                                            00
-1
0                                                                 0
0          10   20    30   40        50        60                 0     10   20    30   40                  50        60
X(1280km)                                                    X(1280km)
-600

17:48 UT
Time=3                                                       18:08 UT
Time=4
40                                                                40
50

50
50                                                                      50
30                                                                30 50
50
Y(1280km)

Y(1280km)
-1650
20                                                                20
500

50

-1000              50
10                                                                10
-100                                                                00
-1
0                                                                 0                                                     -2700
0          10   20    30   40        50        60                 0     10   20    30   40                  50        60
X(1280km)                                                    X(1280km)
AR 8210, May 1, 1998
18:24 UT
Time=5                                                                  17:19 UT
Time=6
40                                                                  40                                                     50
50                                                                       50
50                                                                          50
30                                                                  30    50
50

Y(1280km)

Y(1280km)
20                                                                  20
00                                                                        00

50
-10                                                                       -10                  50
10                    -100
10                         -100

0                                                                   0
0        10    20    30   40           50        60                 0              10   20    30   40               50         60
X(1280km)                                                                X(1280km)

18:55 UT
Time=7                                                                  19:27 UT
Time=8
40                                        50                        40                                                     50
50                                                                       50
50
30    50              -1000                                         30    50
Y(1280km)

Y(1280km)
20                                                                  20

500
-100
0               50
50

10                   -100                                           10                        -100

0                                                                   0
0        10    20    30   40           50        60                 0              10   20    30   40               50         60
X(1280km)                                                                X(1280km)
AR 8210 May 1, 1998
19:43 UT
Time=9                                                                    19:58 UT
Time=10
40                                                   50                    40                                                          1500
50                                                                         50
50                                                                         50
30 50                                                                      30    50
Y(1280km)

Y(1280km)
20                                                                         20

500
500                       -100
0             50                                                     -100
0              50        450
10                         -100
10                        -100

0                                                                          0
0              10    20    30   40            50        60                 0              10   20    30   40              50        60
X(1280km)                                                                 X(1280km)
-600

20:10 UT
Time=11                                                                  20:29 UT
Time=12
40                                                                         40
50                                                                             50
50                                                                          50
30    50                                                                   30 50
Y(1280km)

Y(1280km)
-1650
20                                                                         20
500

00               50                                                         00               50
-10                                                                        -10
10                         -100
10                       -100

0                                                                          0                                                          -2700
0              10    20    30   40            50        60                 0              10   20    30   40              50        60
X(1280km)                                                                 X(1280km)
VI. Concluding Remarks

• A data driven, 3D, MHD model for the initiation and propagation of solar disturbances is
presented. This model consist of two parts:
• Local (Rectangular) Coordinates/AR Study
• Global Coordinates/Solar Wind from the Sun-Heliosphere

•Time-dependent boundary conditions are fully implemented. Thus, the model is capable of
taking the observables as the inputs to drive the model such that the subphotosphere
effects are included.

Solar Interior              Corona                  Heliosphere
THANK YOU
AR 8100, October 31, 1997
Time=14:27UT Oct 31                                                     Time=16:03UT Oct 31
10

20                                                                20
Y(3.95arcsec)

Y(3.95arcsec)
15                                                                15
-4.0

0
-1.0

-1.
-4.0
10    2.0                                                         10

-1.0
5

-1.0
2.0
5                                                                 5
6.0

6.0
0                                                                 0
0         10           20           30                            0           10            20          30
X(3.95arcsec)                                                        X(3.95arcsec)

0

Time=17:39UT Oct 31                                                     Time=19:12UT Oct 31

20                                                                20
Y(3.95arcsec)

Y(3.95arcsec)
15                                                                15                                                   -5
-1.0
-1.0

10                                                                10
2.0                                                                   2.0

-1.0
-1.0

5                                                                 5
6.0

6.0

0                                                                 0                                                  -10
0         10           20           30                            0           10            20          30
X(3.95arcsec)                                                        X(3.95arcsec)
AR8100
50
Lsg
Lss
flux
I
40                                                         alfa

30

20

10

0
0            200   400   600      800   1000    1200        1400
Time (min)
00:00 UT 11/03/97                                     00:00 UT 11/04/97
AR 10646

08:00 UT 07/13/04
AR 10646

09:00 UT 07/13/04
AR 10646

10:00 UT 07/13/04
The left column of panels shows
the shock characteristics in the
corona (1-30 Rs), and the right
column indicates the evolution of
these characteristics into the
heliosphere (30-212 Rs).

(a) Position of MHD fast wave
front
(b) shock normal angle vs time
along four meridional angles
measured from the north
pole;
(c) plasma beta; and
(d) MHD fast wave Mach
number.
The solid lines = 0° (pole), dotted
lines = 45°, dotted-dashed lines =
70°, and dashed lines = 90°
(equator)
The disturbed velocity, density,
temperature, total magnetic
fields, and latitudinal components
(Bθ) at 1 A.U. as a function of time
at the equator (left panel)
and 45° away from the equator
(right panel), respectively, where
S, Sh and MC indicate the shock,
sheath, and magnetic cloud.
To study the evolution of an active region, the model is driven by differential
rotation, meridional flow and the measured magnetic field at the photospheric
level. To input the measured magnetic field the following procedure is used;
We take two consecutive sets of MDI Magnetograms subtract one
from the other, then incrementally increase at six second intervals
to cover these two sets of MDI data. Specifically, the expression
used is:

Bz ( x, y, t ) = Bz ( x, y, o ) + ∑∑ (Δ Bz )n ,
k    n

Since we have chosen our time step to be six seconds and the MDI
measurements are every 96 minutes, thus it takes 960 steps to fill the
period, hence
(ΔBz )n   =
(        )
Bz ( x, y, t k +1 ) − Bz x, y, t k
, n = 1, 2, ..., 960
960

Where Bz(x,y,tk+1) and Bz(x,y,tk) are obtained from MDI magnetograms.

This process is repeatedly carried out through the total simulation time.

```
To top