Automated Shock Analysis


Automated Shock Detection and Analysis Algorithm

-- Details published in Space Weather Journal --


Interplanetary shocks appear to have "...[density] compresion ratio between 1 and 8 with an average close to 2... magnetic compression [ratio] with an average at 1.9... Alfvén Mach number is between 1 and 13 with an average at 1.7." [1] The averages tend to be low, thus showing that weak shocks dominate the interplanetary space. Wind spacecraft from 1994 to 1997 confirms that weak shocks with low Mach number constitute the most common observation. [2] Unfortunately, there exist other events in the SWEPAM data, including 'waves', that have comparable observed parameters. Therefore, the foremost complication of detecting shocks in an automated fashion is rejecting false positives while keeping as many real shocks as possible. Automated shock analysis has proved to reject false positives at the price of also rejecting weak shocks that have low impact on space weather. This trade off is acceptable so long as space weather interests continue to focus exclusively on the strongest shocks.

The automated detection and analysis is a 7-step process. The first 3 steps contribute to finding likely shock candidates and selecting upstream and downstream data for the analysis. The next 2 steps involve the fitting of the data to the Rankine-Hugoniot jump conditions to obtain a solution that includes propagation speed, its direction, density, and magnetic compression. We use the methods of Viñas and Scudder [3] and Szabo [4] to obtain solutions to the shock equations, but the technique described here is applicable to any method that yields valid shock parameters. The remaining 2 steps are what is normally performed in some comparable manner by the operator. They involve passing judgement on the quality of the solution. When the solution obtained in this manner is judged to be unacceptable the investigator will either judge the interval to be something other than a shock or select other upstream and downstream points to represent the plasma conditions. The automated analysis must perform these steps in a reasonable manner so as to produce results in good agreement with the best interactive solutions while reliably determining whether the candidate event is in fact a shock or is possibly something else to be discarded or noted diferently.

The code was tested on science-quality Level-2 data from the Advanced Composition Explorer (ACE) mission [5] that merges thermal ion data from the SWEPAM instrument [6] with magnetic field data [7]. These data have already been analyzed [8], allowing for comparison of solutions: automated vs. interactive. The data files contain date, time, proton density, temperature, bulk speed, three components of the solar wind velocity in RTN coordinates, three components of average interplanetary magnetic field (IMF) also in RTN coordinates, and the average IMF scalar magnitude. Data are recovered every 64 seconds and are handled in a manner that closely resembles the 1-min environment of the NOAA real-time data stream.

Step 1: Identifying Shock Candidates

The first step is to find a suitable candidate. A shock finding subroutine compares each of the two consecutive data points or with a single gap of one point in case of missing data. The subroutine looks for jumps in velocity, temperature, and proton density. Shock identification implementation does not look across data gaps larger than 1 point. The criteria for a set of data to pass step 1 are: 1.5% jump in velocity, 15% jump in temperature, 20% jump in proton density, and at least 68 for the sum of the percentage values. These criteria were chosen based on the weakest shocks in 1999 ACE data.

Step 2: Refining Shock Candidates

In the second step we use averages of velocity, temperature, and density. Ten data points immediately before and ten data points after the shock candidate time are selected (roughly 11 minutes on either side of the shock), and bad data are disregarded. Examining wind speed, density, and temperature separately, we select three points from before and after the shock that are not necessarily consecutive, but are closest in value, and average the values on either side of the discontinuity. The points chosen for one parameter may not be the same points chosen for another, and the averages are obtained. We then check whether the averages are within 50% of the single-point jump parameters, and apply a similar test on averages as we applied in step 1.

Step 3: Selecting Upstream and Downstream Data Points

Because the final shock solutions must be based on a common selection of time tags for all variables and because the wind speed can vary greatly and independetly from the shock speeds, we use the data points from the density selection in Step 2, together with all variable from those data times, as input for our obtaining of shock solutions.

Step 4: Obtaining Shock Normal and Speed

In this step, we perform a solution to Rankine-Hugoniot equations. [3] [4] The algorithm involves an iteration procedure starting from an initial guess of (θ, φ), the angular direction of the normal of a planar shock, and often the initial guess is crucial to ensuring the convergence of the solution. A well-educated guess comes from an extra step in which iterations are carried out for selected (θ, φ) parameters over a finite 2D grid spanning the whole parameter space. Each yields a χ², the least-squares minimization of deviations between the theoretically predicted values and actual in situ spacecraft measurements. Then a pair of (θ, φ) that is close to the one of minimum χ² is used as the initial guess to start teh iteration to converge to the best-fit solution. While in case of an infinite loop, the initial guess of (θ, φ) has to be adjusted or different points need to be selected, the automated code marks the solution as poor. A few candidates in the 1999-2005 data set ran into this problem, none of which appeared to be real shocks. Once the shock normal is obtained, the shock speed is provided using the observed upstream and downstream plasma velocity and density data.

Step 5: Obtaining Asymptotic Magnetofluid State

In this step, the Rankine-Hugoniot conservation constants and the self-consistent asymptotic magnetofluid states upstream and downstream of the shock are obtained. These involve a one-dimensional nonlinear least-squares procedure that minimizes a system of functions of a single variable, ρ, the mass density, by applying the same algorithm as in Step 4.

Step 6: Judging the Quality of the Solutions

This step is needed to discard false positive events based on the candidate's asymptotic solution to the R-H equations. As presently configured (and subject to modification), the test accepts shock candidates if MA > 0.7 (most real shocks have a higher Mach number), Vs/uncertainty > 2.9, ρVp/uncertainty > 3.0, and the tangential momentum divided by its uncertainty > 0.5. False positive that pass this test tend to be weak. Real shocks that fail this test tend to be similarly weak and are not of great importance to space weather applications.

Step 7: Raking by Shock Solution Quality

While Step 6 performs reasonably well, it still tends to admit many false positive solutions and many poor solutions for reliable space weather applications. Therefore, this final step is designed to remove false psitives previously accepted while keeping the strongest of the shocks. We make sure that the shocks accepted are indeed strong shocks with self-consistent solutions by assigning points to each aspect of the shock solution. This point system for evaluating the quality of the solution and the strength of the shock are used in this final step. Most shock candidates with point values of 40 or above appear to be real shocks, most of which are strong.

Points Criteria
5 Mass Flux/uncertainty > 4.5
5 Tang Mom/uncertainty > 1.5
10 MA - uncertainty > 4.2
15 Vs/uncertainty > 14
10 Vs uncertainty < 35
10 Th_Bn uncertainty < 1.7
20 r_n - uncertainty > 2.9
25 r_b - uncertainty > 1.9


Funding provided by NASA grants NNG04GMO5G, NAG5-12492, and Caltech subcontract 44A-1062037 in support of the ACE/MAG experiment. Support at LANL was provided under the auspices of the U.S. Department of Energy, with financial support from the NASA ACE program.


1. Kallenrode, M.-B. (2001), Space Physics, Springer-Verlag, New York.
2. Berdichevsky, D.B., A. Szabo, R.P. Lepping, A.F. Vinas, and F. Mariani (2000), Interplanetary fast shocks and associated drivers observed through the 23rd solar minimum by Wind over its first 2.5 years, J. Geophys. Res., 105, 27,289-27,314.
3. Viñas A.F., and J.D.Scudder (1986), Fast and optimal solution to the "Rankine-Hugoniot problem", J. Geophys. Res., 91, 39-58.
4. Szabo, A. (1999), An improved solution to the "Rankine-Hugoniot" problem, J. Geophys. Res., 99, 14,737-14,746.
5. Stone, E.C., A.M.Frandsen, R.A.Mewaldt, E.R.Christian, D.Margolies, J.F.Ormes, and F.Snow (1998), The Advanced Composition Explorer, Space Sci. Rev., 86 (1-4), 1-22.
6. McComas, D.J., et al. (1998), Solar wind electron proton alpha monitor (SWEPAM) for the Advanced Composition Explorer, Space Science Rev., 86 (1-4), 563-612.
7. Smith, C.W., M.H.Acuña, L.F.Burlaga, J.L'Heureux, N.F.Ness, and J.Scheifele (1998), The ACE magnetic field experiment, Space Sci. Rev., 86 (1-4), 613-632.
8. Hu, Qiang. Shock solutions 1998-2007

Note: Velocity data is not provided in real-time fashion. Only bulk speed is provided to public online. For the purpose of this analysis, an assumption was made that the x-component of the velocity is equal to the bulk speed.

The following figures show the comparison of shock solutions using the old automated code versus the new automated code. The old code uses three components of the shock velocity while the new code assumes that the radial velocity component is equal to the bulk speed, holding the other two components zero. While the assumption is not true, it allows for the same old code to analyze raw data, where only the bulk speed is given.

2000 data set was used for this comparison. It was observed that many solutions did not agree. However, the strongest shock candidates, i.e. with point values greater than 40 points appeared to have similar solutions. The figures below show the comparison (old vs. new automated code) of shock candidates with point values of 40 or above. The horizontal axis represents values obtained using the old method of using all three velocity components. The vertical axis represents values obtained using the new method of using bulk speed as a single velocity component. The red line (y=x) represents where all points should be.

Speed in Spacecraft Frame (Vs)
Speed in Plasma Frame (Vp)
Density Compression Ratio (Rn)
Magnetic Field Strength Compression Ratio (Rb)
Alfven Mach Number (MA)

Automated Shock Analysis Authors

Vassili S. Vorotnikov
Department of Chemical Engineering
University of New Hampshire

Charles W. Smith
Space Science Center, Morse Hall
University of New Hampshire

Qiang Hu
Institute of Geophysics and Planetary Physics
University of California at Riverside

Adam Szabo
Goddard Space Flight Center
Greenbelt, Maryland

Ruth M. Skoug
Los Alamos National Laboratory
Los Alamos, New Mexico

Christina M.S. Cohen
Space Radiation Laboratory
Caltech, Pasadena, CA

Page Designed by: Vassili S. Vorotnikov
Department of Chemical Engineering, UNH

Supervised by: Charles W. Smith
Department of Physics and Space Science Center, UNH

Page last modified: 09/30/2008