VIEWS: 7 PAGES: 65 POSTED ON: 2/21/2012
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.