Numerical simulations of solar spicules: Adiabatic and nonadiabatic studies
^{1} Group of Astrophysics, University of Maria CurieSkłodowska, ul. Radziszewskiego 10, 20031 Lublin, Poland
email: blazejkuzma1@o2.pl
^{2} IGAM, Institute of Physics, University of Graz, Universitätsplatz 5, 8010 Graz, Austria
^{3} Abastumani Astrophysical Observatory at Ilia State University, 0160 Tbilisi, Georgia
^{4} Space Research Institute, Austrian Academy of Sciences, Schmiedlstrasse 6, 8042 Graz, Austria
^{5} Dipartimento di Fisica Generale, Univesità di Torino, via Pietro Giuria 1, 10125 Torino, Italy
Received: 20 April 2016
Accepted: 30 September 2016
Aims. We aim to study the formation and evolution of solar spicules using numerical simulations of a vertical velocity pulse that is launched from the upper chromosphere.
Methods. With the use of the PLUTO code, we numerically solved adiabatic and nonadiabatic magnetohydrodynamic (MHD) equations in 2D cylindrical geometry. We followed the evolution of spicules triggered by pulses that are launched in a vertical velocity component from the upper chromosphere. Then we compared the results obtained with and without nonadiabatic terms in the MHD equations.
Results. Our numerical results reveal that the velocity pulse is steepened into a shock that propagates upward into the corona. The chromospheric cold and dense plasma follows the shock and rises into the corona with the mean speed of 20–25 km s^{1}. The nonlinear wake behind the pulse in the stratified atmosphere leads to quasiperiodic rebound shocks, which lead to quasiperiodic rising of chromospheric plasma into the corona with a period close to the acoustic cutoff period of the chromosphere. We found that the effect of nonadiabatic terms on spicule evolution is minor; the general properties of spicules such as their heights and risingtime remain slightly affected by these terms.
Conclusions. In the framework of the axisymmetric model we devised, we show that the solar spicules can be triggered by the vertical velocity pulses, and thermal conduction and radiative cooling terms do not exert any significant influence on the dynamics of these spicules.
Key words: Sun: activity / magnetohydrodynamics (MHD) / methods: numerical / Sun: corona / Sun: transition region
© ESO, 2017
1. Introduction
Spicules are thin, cool, and dense structures that are observed in the solar limb (Beckers 1968, 1972; Suematsu 1998; Sterling 2000; Zaqarashvili & Erdélyi 2009). They are seen to emerge from the chromospheric background at an altitude of about 2000 km above the solar surface where they reveal a speed of 25 km s^{1}, reach a maximum level, and then either disappear or sink down to the chromosphere. A typical lifetime of spicules is within the range of 5–15 min with an average value of ~7 min (Pasachoff et al. 2009). Spicules seem to consist of doublethread structures (Tanaka 1974; Dara et al. 1998; Suematsu et al. 2008), and they reveal bidirectional flows (Tsiropoula et al. 1994; Tziotziou et al. 2003, 2004; Pasachoff et al. 2009). Typical electron temperature and electron density in spicules are (15–17) × 10^{3} K and 2 × 10^{11}–3.5 × 10^{10} cm ^{3} at altitudes of 4–10 Mm above the solar surface (Beckers 1968) with a diameter estimated as 660 ± 200 km (Pasachoff et al. 2009). As a result, spicules are much cooler and denser than ambient coronal plasma. Highresolution observations by the Solar Optical Telescope onboard Hinode have revealed another type of spicules with many features different from those of classical limb spicules, and they are referred to as type II spicules (De Pontieu et al. 2007). The type II spicules are distinguished by (a) smaller diameters (≤200 km) in the Ca II H line and a significantly shorter height of 4 Mm; (b) a lifetime of 10–150 s; (c) the evolution, which shows an upflow and then disappears; and finally by (d) much higher speeds of 50–100 km s^{1} (e.g., De Pontieu et al. 2007, 2009; Rouppe van der Voort et al. 2009; Kuridze et al. 2015). Macrospicules were recently discussed by NóbregaSiverio et al. (2016), who considered whether the relevance of the entropy sources in the surges, such as the optically thin losses, can be applied to similar phenomena as macrospicules.
In spite of various theoretical models that have been brought forward to explain the spicule ejection in the lower solar atmosphere, many recent numerical methods have been developed to simulate the solar spicules or macrospicules with an energy input at their base in the photosphere such as a gas pressure pulse or an Alfvén wave that steepens into a shock wave (Sterling 2000, and references therein). Hansteen et al. (2006) and De Pontieu et al. (2007) simulated the formation of dynamic fibrils that are due to slow magnetoacoustic shocks through twodimensional (2D) numerical simulations. They suggested that these shocks are formed when acoustic waves generated by convective flows and global pmodes in the lower lying photosphere leak upward into the magnetized chromosphere. Heggland et al. (2007) used the initial periodic piston to drive the upward propagating shocks in 1D simulation, and MartinezSykora et al. (2009) considered the emergence of new magnetic flux, but the drivers of spicules originate from collapsing granules, energy release in the photosphere, or in the lower chromosphere. However, these simulations were unable to mimic the double structures and bidirectional flows in spicules. On the other hand, Murawski & Zaqarashvili (2010) performed 2D numerical simulations of magnetohydrodynamic (MHD) equations and showed that the 2D rebound shock model of Hollweg (1982) may explain both the double structures and bidirectional flows. They used a single initial velocity pulse, which led to the formation of consecutive shocks as a result of the nonlinear wake in the stratified atmosphere. However, they considered a simple model of atmospheric temperature that was approximated by a smoothed step function for the temperature profile.
Understanding exact drivers of spicules requires further investigation, and more than one mechanism may trigger their evolution depending on the local plasma and magnetic field conditions. The goal of this paper is to contribute to the abovementioned studies by performing simulations of the generation and evolution of the spicules in the solar atmosphere and compare our results for adiabatic and nonadiabatic MHD equations. The method we chose to trigger a spicule is a localized vertical velocity pulse launched from the upper chromosphere. This method is similar to the calculations performed by Shibata (1982), Sterling et al. (1993), Murawski & Zaqarashvili (2010), and Guerreiro et al. (2013). This approach differs from the models that attempt to model spicules with a disturbance in the photosphere (e.g., Suematsu et al. 1982; Hollweg 1982). Guerreiro et al. (2013) studied the midchromospheric energy inputs of earlier simulations by adding additional physics to the radiative loss term and including hydrogen ionization and recombination. They concluded that it would be difficult to produce spicules through those previously suggested mechanisms (specifically that of Sterling et al. 1993). Our simulations do not include the detailed energy losses of Guerrero et al. However, the energy input is different from that assumed by Guerrero et al. and Sterling et al., who adopted localized increase in the heating rate. Since the form of input energy is different, it is not clear whether the losses Guerrero et al. took into account would have any effect on the velocity of pulsedriven spicules. This will be tested in the future studies.
This paper is organized as follows. A numerical model is presented in Sect. 2, and the corresponding numerical results are shown in Sect. 3. Our paper is concluded by a summary of the numerical results in Sect. 4.
2. Physical model of the solar atmosphere
2.1. MHD equations
We consider a gravitationally stratified and magnetically confined solar plasma that is governed by the following set of nonadiabatic MHD equations: where ϱ is the mass density, p the gas pressure, V represents the plasma velocity, B is the magnetic field, T the temperature, q the anisotropic thermal conduction flux, L(ϱ,T) radiatively thin cooling terms (Mignone et al. 2007), and H(ϱ_{e},T_{e}) denotes the external heating therm that balances L and ∇·q at the equilibrium that is specified in Sect. 2.2. This term depends only on the equilibrium plasma quantities and therefore it does not vary in time. The symbol k_{B} denotes the Boltzmann constant, γ = 5/3 is the adiabatic index, m is the particle mass that is specified by a mean molecular weight of 0.6, and g = (0,−g,0) is the gravitational acceleration. The value of g is equal to 274 m s^{2}.
Fig. 1 Hydrostatic solar atmospheric temperature vs. height y. 

Open with DEXTER 
Fig. 2 Vertical profile of plasma β (left) and sound speed c_{s} (right) at the plasma equilibrium. 

Open with DEXTER 
2.2. Equilibrium solar atmosphere
In a static solar atmosphere all plasma quantities are timeinvariant, which means that ∂f_{e}/∂t = 0, where f_{e} denotes a plasma quantity and the subscript e corresponds to the equilibrium. Then, from Eqs. (1)–(4) it follows that for a still (V_{e} = 0) medium the Lorentz force must be balanced by the gravity force and the gas pressure gradient, (5)and the heating term must compensate for the radiative losses and thermal conduction, (6)This model of the solar atmosphere corresponds to a quiet Sun.
Fig. 3 Magnetic field lines at the plasma equilibrium. 

Open with DEXTER 
2.2.1. Forcefree magnetic field of the hydrostatic atmosphere
A hydrostatic atmosphere corresponds to the forcefree ((∇ × B_{e}) × B_{e} = 0) magnetic field. We additionally assume a currentfree (∇ × B_{e} = 0) magnetic field whose radial B_{er}, azimuthal B_{eθ}, and vertical B_{ey} components are given as (7)where a and S are free parameters corresponding to the vertical location of the singularity in the magnetic field and the magnetic field strength, respectively. We set a = −1 Mm and S in such way that at the reference point (r = 0,y = y_{r} = 6) Mm the magnitude of magnetic field B_{e} = 9.5 Gauss. The corresponding magnetic field lines are displayed in Fig. 3. We note that the magnetic lines diverge with height and B_{e} is vertical along the symmetry axis, r = 0 Mm.
For a forcefree magnetic field it follows from Eq. (5) that the gas pressure gradient has to be balanced by the gravity force, (8)The subscript h corresponds to a hydrostatic quantity. With the use of the ideal gas law and the vertical ycomponent of Eq. (8), we express the hydrostatic gas pressure and mass density as (9)where (10)is the pressure scale height, and p_{0} denotes the gas pressure at the reference level, y = y_{r}.
For simplicity reasons we assume that T_{h} varies with height y only, and it specifies a hydrostatic atmosphere that is determined by the semiempirical model of Avrett & Loeser (2008) that is extrapolated into the solar corona (Fig. 1). In this model, the temperature attains a value of about 7 × 10^{3} K at the top of the chromosphere, y ≈ 2.0 Mm. At the transition region, which is located at y ≃ 2.1 Mm, T_{h} exhibits an abrupt jump (Fig. 1) and grows to about 1.0 × 10^{6} K in the solar corona at y = 10 Mm. Higher up in the solar corona, the temperature increases very slowly, tending to its asymptotic value of about 2 MK at y = 40 Mm. The temperature profile uniquely determines the equilibrium mass density and gas pressure profiles, which decrease with height (not shown).
We specify the plasma β as the ratio of gas to magnetic pressures, (11)The vertical profile of plasma β is illustrated in Fig. 2 (left panel). We note that for the coronal plasma, the value of plasma β is lower than 1 within the displayed region. The vertical profile of the sound speed (12)is displayed in the right panel of Fig. 2. Below the transition region, c_{s} ≈ 10 km s^{1}. Higher up c_{s} grows, first fast right above the transition region, and higher up slowly, reaching a value of about 100 km s^{1} at y ≈ 20 Mm.
Fig. 4 Temporal evolution of log (ϱ(r,y)) at t = 20 s, t = 50 s, t = 110 s, t = 175 s, t = 220 s, t = 300 s (from top left to bottom right), for adiabatic MHD equations and for the case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. Arrows represent velocity vectors in the r−y plane, [V_{r},V_{y}]. 

Open with DEXTER 
2.2.2. Perturbation
Initially, at t = 0 s we perturb the model equilibrium by the initial pulse in the ycomponent of velocity, which is expressed as follows: (13)where y_{0} is the vertical position of the initial pulse, w is its width, and A_{V} its amplitude. We set and hold fixed w = 0.25 Mm, while allowing other parameters to vary. For our studies, the initial position of y_{0} varies between 1.5 Mm and 1.75 Mm, and the amplitude A_{V} varies between 30 km s^{1} and 50 km s^{1}. The detailed studies were performed for the case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. The value of A_{V} may be associated with reconnection of magnetic field lines.
3. Numerical simulations of MHD equations
To solve Eqs. (1)–(3) numerically, we used the PLUTO code (Mignone et al. 2007, 2012). In our problem, we set the CourantFriedrichsLevy number equal to 0.3 and chose piecewise TVD linear interpolation in a secondorder RungeKutta method, which leads to secondorder accuracy in space and time. Additionally, we adopted the HartenLaxvan Leer discontinuities (HLLD) approximate Riemann solver (Miyoshi & Kusano 2005).
Our simulation box in (r,y) was set as (0.0,5.12) Mm ×(1.0,40.0) Mm, where y = 0 denotes the bottom of the photosphere. For our study we used the uniform grid within the region (0.0 ≤ r ≤ 5.12) Mm ×(1.0 ≤ y ≤ 11.24) Mm, which is covered by 1024 × 2048 grid points. This leads to a resolution of 5 km in the lower regions of the simulation box. Above this region, namely within the box (0.0 ≤ r ≤ 5.12) Mm ×(11.24 ≤ y ≤ 40.0) Mm, we implemented a stretched grid along the ydirection; this box was divided into 648 cells whose size grows with y. Such a stretched grid plays the role of a sponge as it absorbs incoming signal and allows us to avoid significant reflections from upper boundary. We imposed open boundary conditions at r = 5 Mm, but at the bottom and top, we fixed all plasma quantities to their equilibrium values. The left boundary (r = 0) was set as axisymmetric.
The heating source term H in Eq. (3) was implemented as follows. The code computed the residuum R^{0} using the initial condition and then it subtracted R^{0} also at later times during the update over time step Δt as (14)where U^{n} is a plasma vector state at time t = nΔt, n = 1,2,... This automatically mimics the heating term without any need to explicitly write it.
Fig. 5 Temporal evolution of ϱ(r = 0,y) (top) and V_{y}(r = 0,y) (bottom) for the case of adiabatic (left) and nonadiabatic (right) MHD equations and for A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. 

Open with DEXTER 
To fully understand the influence of nonadiabatic effects on simulated spicules, we performed first simulations of the adiabatic case, that is, without the thermal conduction and cooling terms. Figure 4 shows the spatial profiles of logarithm of ϱ(r,y) at six instants of time. We note that the system is axisymmetric along r = 0 Mm. The initial pulse in V_{y} is launched from r = 0 Mm, y = 1.75 Mm, which is located about 0.35 Mm below the transition region. The amplitude of the initial pulse is A_{V} = 40 km s^{1}. The shock front that results from the initially launched pulse arrives at the transition region and triggers the plasma jet (top panels), which reaches its maximum height of ≈4.3 Mm (top right panel). At later moments in time, the chromospheric plasma injected into the corona begins to fall toward the transition region (bottom left panel). The falling plasma then triggers the second pulse (bottom middle panel), which results in the second spicule (bottom right panel).
Fig. 6 Temporal evolution of ϱ(r = 0,y = 5) Mm for the case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. 

Open with DEXTER 
Figure 5 shows the temporal evolution of ϱ(r = 0,y) (top panels) and V_{y}(r = 0,y) (bottom panels) for adiabatic (left panels) and nonadiabatic (right panels) cases. The rise time of the spicule to its maximum height in the adiabatic case is ≈110 s, which is shorter by about 15 s than in the nonadiabatic case (≈125 s), and the chromospheric plasma reaches a height of ≈4.3 Mm. This value is lower than in nonadiabatic case (≈4.7 Mm). In the adiabatic case we spot three waves (Fig. 5, top left panel): the leading wave is a shock wave, which travels with velocity ≈170 km s^{1}; this shock is followed by the contact wave, and the rarefaction wave. We note that all plasma quantities are continuous across a contact wave; the only exception is the mass density, which is discontinuous. As a result, all three waves are seen on the mass density profiles, while V_{y}(r = 0,y) profiles reveals only two of these waves. As in the nonadiabatic case, in later moments, the system becomes thermally unstable. The right panels show the initial phase of the system evolution. The crosssection of Fig. 5 for the adiabatic case for a fixed value of y = 3.5 Mm is presented in Fig. 6. The corresponding three waves (shock, contact, and rarefaction) are represented by steep gradients at t ≈ 30 s, t ≈ 55 s, and t ≈ 75 s. As the formation mechanism of the spicule is a shock wave, we assume that it is the crestshocktype jet (Shibata 1982). The mass density at the top of the spicule is ~10^{2} higher than in the ambient coronal plasma (Fig. 6, t ≈ 75 s), and it matches the predicted value for the crestshocktype jet.
Fig. 7 Radial profile of ϱ(r,y = 3.1 Mm) at t = 35 s (solid line), t = 80 s (dashed line), and t = 155 s (dotted line) in the case of adiabatic MHD equations. The case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. 

Open with DEXTER 
Fig. 8 Maximum height of spicule y_{max} vs. amplitude of the initial pulse A_{V} for y_{0} = 1.75 Mm. The case of adiabatic (+) and nonadiabatic (+?) MHD equations. 

Open with DEXTER 
Figure 7 shows the effect caused by the rarefaction wave. The horizontal plane is located at y = 3.5 Mm. The plasma inside the spiculelike structure becomes rarefied. The time at which this rarefaction occurs is clearly shown: the dense plasma at t = 35 s (solid line) experiences rarefaction within the central part of the spicule at t = 80 s (dashed line), which keeps its trend up to t = 155 s (dotted line). At this time, plasma is so rarefied that it may be missed by the detectors. This may mimic the disappearance of II type spicules.
Figure 8 illustrates spicule height y_{max} vs. amplitude of the initial pulse A_{V} with a fixed position of the initial pulse y_{0} = 1.75 Mm. Hence we infer that a larger amplitude pulse results in chromospheric plasma reaching a higher maximum altitude, which is intuitively expected as a larger amplitude corresponds to more energy being launched initially. This more energetic pulse leads to higher jets. Figure 9 shows a similar trend in the arise time of the spicule for a growing amplitude of the pulse. The comparison of results for adiabatic (+) and nonadibatic (*) MHD equations reveals that the nonadiabatic effects have a minor influence on both spicule height and rise time of the spicule. The nonadiabatic terms result in a growth of the maximum height by about ~10–15%, while the rise time of the spicule grows by about ~15–20%.
Figure 10 illustrates the influence of the vertical position of the initial pulse on the maximum height of the spicule. For a pulse launched closer to the transition region, the mass density of the upwardpushed plasma is lower. Thus such a pulse exhibits a lower momentum, and we can expect it to penetrate lower coronal regions. We also expect that the rise time of the spicule becomes shorter with a higher value of y_{0}. Indeed, Fig. 11 illustrates this declining trend. The obtained values of spicule height lie within the range of the observed height of the spicules.
Fig. 9 Rise time of spicule t_{r} vs. amplitude of the initial pulse A_{V} for y_{0} = 1.75 Mm. The case of adiabatic (+) and nonadiabatic (+?) MHD equations. 

Open with DEXTER 
Fig. 10 Maximum height of spicule y_{max} vs. vertical position of the initial pulse y_{0} for A_{V} = 40 km s^{1} and in the case of adiabatic MHD equations. 

Open with DEXTER 
Fig. 11 Rise time of spicule t_{r} vs. vertical position of the initial pulse y_{0} for A_{V} = 40 km s^{1} and in the case of adiabatic MHD equations. 

Open with DEXTER 
Fig. 12 Spicule velocity V_{y} vs. time for y_{0} = 1.75 Mm, A_{V} = 40 km s^{1} and adiabatic MHD equations. 

Open with DEXTER 
Figure 12 displays the dependence of spicule velocity on time. The location of the top of the spicule is determined by the vertical position of the contact wave. The upward movement of the transition region plasma starts at t ≈ 10 s, and at this time the spicule reaches its speed of about 45 km s^{1}. At later moments in time, this spicule decelerates until at t ≈ 100 s the spicule speed attains zero, which occurs at the highest location of the spicule. At the lowest location of the spicule, which is reached at t = 200 s, the vertical velocity of the decreasing plasma is about 30 km s^{1}. After the first rise and fall, the transition region subsequently experiences oscillations.
Fig. 13 Timedistance plot of ϱ(r = 0,y,t) vs. time for y_{0} = 1.75 Mm, A_{V} = 30 km s^{1} and adiabatic MHD equations. 

Open with DEXTER 
Fig. 14 Acoustic cutoff wave period P_{ac} vs. altitude y. 

Open with DEXTER 
The transition region exhibits oscillations that are due to the rebound shocks (Hollweg 1982). It is interesting that these oscillations occur with a wave period of about 250 s (Fig. 13). This wave period is close to the acoustic cutoff period, (15)which for y = y_{0} = 1.75 Mm attains a value of about 259 s (Fig. 14).
4. Summary and conclusions
We performed numerical simulations of a spicule by launching initially (at t = 0 s) a localized pulse in the ycomponent of a plasma velocity in the upper chromosphere. The initial magnetic field configuration was forcefree and the initial stratification remained in static equilibrium. We simulated both ideal and nonideal MHD equations to reveal nonadiabatic effects such as thermal conduction and radiation. Numerical simulations showed that an upwardpropagating signal quickly steepened into a shock that propagated into the corona along the magnetic field lines. This shock was followed by the cold and dense chromospheric plasma jet, which exhibited properties of a contact wave and reached a certain height (typically 4–5 Mm) and then returned to the chromosphere. The mean upflow speed was 20–25 km s^{1}. The obtained values match those given for the spicules observed by Beckers (1968, 1972). However, subsequent shocks come from the chromosphere with a periodicity of almost the chromospheric acoustic cutoff period, which is the result of the nonlinear wake behind the pulse that propagates in the stratified atmosphere (Kuridze et al. 2009). These shocks again lift up the chromospheric plasma into the corona and cause the quasiperiodic appearance of plasma jets. This is consistent with the rebound shock model suggested by Hollweg (1982). Numerical solutions show that the rarefaction wave that follows the shockwave results in a reduction of the mass density inside the structure. Both ideal and nonideal simulations give similar results for the maximum height, upward speed, rise time, and periodicity. Therefore, nonadiabatic effects do not significantly affect the dynamics of jets, but slightly increase the maximum heights and rise time.
In conclusion, our numerical simulations of the solar spicules approximately mimicked the observed properties of the spicules. Since we did not attempt to synthesize observables, an exact comparison with observations cannot be achieved here.
Acknowledgments
This work (B.K. and K.M.) was financially supported by the project from the Polish National Foundation (NCN) Grant No. 2014/15/B/ST9/00106. The work of T.V.Z. was supported by the Austrian “Fonds zur Förderung der Wissenschaftlichen Forschung” (FWF) under projects P26181N27 and P25640N27, by Austrian ScientificTechnology Collaboration (WTZ) Grant PL15/2013 and by FP7PEOPLE2010IRSES269299 project SOLSPANET. The PLUTO code used in this work was developed by Andrea Mignone. Numerical simulations were performed on the Lunar/Solaris cluster at Institute of Mathematics of M. CurieSkłodowska University. The visualizations of the simulation variables have been carried out using the IDL (Interactive Data Language) and VisIt software packages.
References
 Avrett, E. H., & Loeser, R. 2008, ApJS, 175, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Beckers, J. M. 1968, Sol. Phys., 3, 367 [NASA ADS] [Google Scholar]
 Beckers, J. M. 1972, ARA&A, 10, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Dara, H. C., Koutchmy, S., & Sueamatsu, Y. 1998, in Solar Jets and Coronal Plumes, ESA SP, 421, 255 [NASA ADS] [Google Scholar]
 De Pontieu, B., McIntosh, S. W., & Hansteen, V. H., et al. 2007, PASJ, 59, 655 [NASA ADS] [CrossRef] [Google Scholar]
 De Pontieu, B., McIntosh, S. W., Hansteen, V. H., & Schrijver, C. J. 2009, ApJ, 701, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Guerreiro, N., Carlsson, M., & Hansteen, V. 2013, APJ, 766, 128 [NASA ADS] [CrossRef] [Google Scholar]
 Hansteen, V. H., De Pontieu, B., Rouppe van der Voort, L. H. M., van Noort, M., & Carlsson, M. 2006, ApJ, 647, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Heggland, L., De Pontieu, B., & Hansteen, V. H. 2007, ApJ, 666, 1267 [NASA ADS] [CrossRef] [Google Scholar]
 Hollweg, J. V. 1982, ApJ, 257, 345 [NASA ADS] [CrossRef] [Google Scholar]
 Kuridze, D., Zaqarashvili, T. V., Shergelashvili, B. M., & Poedts, S. 2009, A&A, 505, 763 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kuridze, D., Henriques, V., Mathioudakis, M., et al. 2015, ApJ, 802, 26 [NASA ADS] [CrossRef] [Google Scholar]
 MartinezSykora, J., Hansteen, V., De Pontieu, B., & Carlsson, M. 2009, ApJ, 701, 1569 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228 [NASA ADS] [CrossRef] [Google Scholar]
 Mignone, A., Zanni, C., Tzeferacos, P., et al. 2012, ApJS, 198, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Miyoshi, T., & Kusano, K. 2005, JCP, 208, 315 [NASA ADS] [CrossRef] [Google Scholar]
 Murawski, K., & Zaqarashvili, T. V. 2010, A&A, 519, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 NóbregaSiverio, D., MorenoInsertis, F., & MartínezSykora, J. 2016, ApJ, 822, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Pasachoff, J. M., Jacobson, W. A., & Sterling, A. C. 2009, Sol. Phys., 260, 59 [NASA ADS] [CrossRef] [Google Scholar]
 Rouppe van der Voort, L. H. M., Leenaarts, J., De Pontieu, B., Carlsson, M., & Vissers, G. 2009, ApJ, 705, 272 [NASA ADS] [CrossRef] [Google Scholar]
 Shibata, K. 1982, Sol. Phys., 81, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Sterling, A. C. 2000, Sol. Phys., 196, 79 [NASA ADS] [CrossRef] [Google Scholar]
 Sterling, A. C., Shibata, K., & Mariska, J. T. 1993, ApJ, 407, 778 [NASA ADS] [CrossRef] [Google Scholar]
 Suematsu, Y. 1998, in Solar Jets and Coronal Plumes, ESA SP, 421, 19 [NASA ADS] [Google Scholar]
 Suematsu, Y., Shibata, K., Neshikawa, T., & Kitai, R. 1982, Sol. Phys., 75, 99 [NASA ADS] [CrossRef] [Google Scholar]
 Suematsu, Y., Ichimoto, K., Katsukawa, Y., et al. 2008, ASP Conf. Ser., 397, 27 [NASA ADS] [Google Scholar]
 Tanaka, K. 1974, in Chromospheric Fine Structure, ed. G. Athay, IAU Symp., 56, 239 [Google Scholar]
 Tsiropoula, G., Alissandrakis, C. E., & Schmieder, B. 1994, A&A, 290, 294 [NASA ADS] [Google Scholar]
 Tziotziou, K., Tsiropoula, G., & Mein, P. 2003, A&A, 402, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tziotziou, K., Tsiropoula, G., & Mein, P. 2004, A&A, 423, 1133 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Zaqarashvili, T. V., & Erdélyi, R. 2009, Space Sci. Rev., 149, 355 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1 Hydrostatic solar atmospheric temperature vs. height y. 

Open with DEXTER  
In the text 
Fig. 2 Vertical profile of plasma β (left) and sound speed c_{s} (right) at the plasma equilibrium. 

Open with DEXTER  
In the text 
Fig. 3 Magnetic field lines at the plasma equilibrium. 

Open with DEXTER  
In the text 
Fig. 4 Temporal evolution of log (ϱ(r,y)) at t = 20 s, t = 50 s, t = 110 s, t = 175 s, t = 220 s, t = 300 s (from top left to bottom right), for adiabatic MHD equations and for the case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. Arrows represent velocity vectors in the r−y plane, [V_{r},V_{y}]. 

Open with DEXTER  
In the text 
Fig. 5 Temporal evolution of ϱ(r = 0,y) (top) and V_{y}(r = 0,y) (bottom) for the case of adiabatic (left) and nonadiabatic (right) MHD equations and for A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. 

Open with DEXTER  
In the text 
Fig. 6 Temporal evolution of ϱ(r = 0,y = 5) Mm for the case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. 

Open with DEXTER  
In the text 
Fig. 7 Radial profile of ϱ(r,y = 3.1 Mm) at t = 35 s (solid line), t = 80 s (dashed line), and t = 155 s (dotted line) in the case of adiabatic MHD equations. The case of A_{V} = 40 km s^{1} and y_{0} = 1.75 Mm. 

Open with DEXTER  
In the text 
Fig. 8 Maximum height of spicule y_{max} vs. amplitude of the initial pulse A_{V} for y_{0} = 1.75 Mm. The case of adiabatic (+) and nonadiabatic (+?) MHD equations. 

Open with DEXTER  
In the text 
Fig. 9 Rise time of spicule t_{r} vs. amplitude of the initial pulse A_{V} for y_{0} = 1.75 Mm. The case of adiabatic (+) and nonadiabatic (+?) MHD equations. 

Open with DEXTER  
In the text 
Fig. 10 Maximum height of spicule y_{max} vs. vertical position of the initial pulse y_{0} for A_{V} = 40 km s^{1} and in the case of adiabatic MHD equations. 

Open with DEXTER  
In the text 
Fig. 11 Rise time of spicule t_{r} vs. vertical position of the initial pulse y_{0} for A_{V} = 40 km s^{1} and in the case of adiabatic MHD equations. 

Open with DEXTER  
In the text 
Fig. 12 Spicule velocity V_{y} vs. time for y_{0} = 1.75 Mm, A_{V} = 40 km s^{1} and adiabatic MHD equations. 

Open with DEXTER  
In the text 
Fig. 13 Timedistance plot of ϱ(r = 0,y,t) vs. time for y_{0} = 1.75 Mm, A_{V} = 30 km s^{1} and adiabatic MHD equations. 

Open with DEXTER  
In the text 
Fig. 14 Acoustic cutoff wave period P_{ac} vs. altitude y. 

Open with DEXTER  
In the text 