Interactive comment on “ Verification of SpacePy ’ s radial diffusion radiation belt model ”

Abstract. Model verification, or the process of ensuring that the prescribed equations are properly solved, is a necessary step in code development. Careful, quantitative verification guides users when selecting grid resolution and time step and gives confidence to code developers that existing code is properly instituted. This work introduces the RadBelt radiation belt model, a new, open-source version of the Dynamic Radiation Environment Assimilation Model (DREAM) and uses the Method of Manufactured Solutions (MMS) to quantitatively verify it. Order of convergence is investigated for a plethora of code configurations and source terms. The ability to apply many different diffusion coefficients, including time constant and time varying, is thoroughly investigated. The model passes all of the tests, demonstrating correct implementation of the numerical solver. The importance of DLL and source term dynamics on the selection of time step and grid size is also explored. Finally, an alternative method to apply the source term is examined to illustrate additional considerations required when non-linear source terms are used.


Introduction
The terrestrial radiation belts are regions of near-Earth outer space where relativistic electrons and ions are electromagnetically orbiting the planet.The belts naturally organize into two tori: a stable inner belt lying within ≈2 R E and a far more dynamic outer belt located outside of ≈3 R E (van Allen and Frank, 1959).Since their discovery, the belts have been the focus of intense research due to the innumerable unknowns concerning their behavior and their damaging effects on spacecraft, both transient (Baker, 2000;Feynman and Gabriel, 2000;Koons et al., 1999;Pirjola et al., 2005, etc.) and accumulated over a satellite's lifetime (Gubby and Evans, 2002;Welling, 2010).
During slowly changing conditions, radiation belt particles undergo three types of periodic motion, each with its own corresponding adiabatic invariant: gyration about field lines, bounce along field lines, and drift about the Earth.During active times, when conditions change on time scales shorter than the periods of motion, adiabaticity can be broken and particle motion can no longer be described by a simple sum of these three.Casting the particle evolution in terms of the three invariants allows non-adiabatic motion to be modeled diffusively via the Fokker-Planck equation, yielding a powerful description of belt dynamics (Schulz and Lanzerotti, 1974).
The third invariant represents the total magnetic flux enclosed within a full particle orbit.It is common to use a normalized form of this, L * (referred to as simply L herein), described by Roederer (1970).It is analogous to radial distance from the center of the Earth (in Earth radii) to the equatorial crossing point of the bouncing particle, exactly so if the terrestrial magnetic field is adiabatically relaxed to a simple dipole geometry.Diffusion in L alone (the other invariants shall be considered conserved) accounts for the capture and inward radial transport of radiation belt particles (Roederer, 1970;Fälthammar, 1965).The long period of a particle orbit makes this invariant the most readily broken.For these reasons, modeling of phase space density in L space is both fruitful (a great deal of dynamics are captured) and easily achieved (the long particle orbit time allows for larger time steps).
The evolution of the phase space density, f , of a radiation belt population of constant µ and K (values corresponding to the first and second adiabatic invariants, respectively) for a given diffusion rate, D LL , is expressed as (Lyons and Schulz, 1989), Published by Copernicus Publications on behalf of the European Geosciences Union.
where Q represents combined sources and losses.Solving this equation for f has been the focal point of many important radiation belt studies.Brautigam and Albert (2000) and Miyoshi et al. (2003), using different boundary conditions at high L and different loss processes, demonstrated the importance of radial diffusion in radiation belt dynamics.Shprits et al. (2006) showed that radial diffusion and magnetopause shadowing (particle loss as the magnetopause intersects drift paths) can account for observations of rapid flux dropouts at the onset of particular storms.Lam et al. (2007) used a radial diffusion model to explore the importance of plasmaspheric hiss driven losses.The results compared favorably to measurements of the radiation belts around L ≈ 4.These are but a few examples; a more complete review can be found in Shprits et al. (2008).Clearly, radial diffusion representation of the belts is a powerful tool for investigation of their dynamics.
Many models and studies have moved beyond simple onedimensional diffusion.Three-dimensional diffusion models have arose, which allow for not only radial but also pitch angle and energy space diffusion.The Salammbo model (Beutier and Boscher, 1995) is an early example of such.More recent 3-D diffusion models now include cross-diffusion terms and have shown that such terms can be important to the evolution of the radiation belts (Xiao et al., 2010;Subbotin et al., 2010).The Dynamic Radiation Environment Assimilation Model (DREAM, Reeves et al., 2005) incorporates data assimilation to drive radial diffusion results towards more realistic values (Koller et al., 2007).Other models, such as the Radiation Belt Environment (RBE) model (Fok et al., 2008), use bounce-averaged kinetic representations of the belts along with time-accurate magnetic and electric field specifications instead of simpler diffusion equations.Despite these important modeling advances, radial diffusion still remains a core part of radiation belt investigations.
This work introduces a new, open source version of the DREAM radial diffusion radiation belt code and, as part of the code's development, verifies the model.Code verification is the process of testing for proper implementation of a code's numerics and other processes.Verification asks, "Is the model solving the prescribed equations correctly?"and is separate from model validation, which asks, "Are the equations solved representative of the real world?";examples of quantitative model validation can be found in Wang et al. (2008) and Welling and Ridley (2010).Verification studies are key steps in code development and guide user decisions in terms of code configuration.Proper verification also expedites future development; as additional features and expanded physics are added to a model, it is critical to ensure that the existing features are instituted properly.Because RadBelt is currently being expanded to include diffusion in pitch angle and energy space, verification at this early stage is crucial.
This study leverages the Method of Manufactured Solutions, described below, to quantitatively assess the code's performance and ensure that the results are accurate to typical machine limitations.Code exercises are repeated for a plethora of configurations to demonstrate proper implementation and convergence.Alternative source functions and numerical source applications are applied to test all facets of the model.

Model description
The code being investigated here is the RadBelt module of the SpacePy software library (Morley et al., 2011).This model solves Eq. ( 1) for a single µ and K.The bulk of the code is written in Python, making the model unique among others in terms of flexibility and capability; the model is initialized, configured, executed, and visualized all through a Python interface.Code domain (typically L = [1 : 10]), sources and losses, D LL , time step and grid size are all specified in an object-oriented manner, e.g., setting object attributes.Required input data, such as Kp indices, are obtained automatically by SpacePy from the Qin et al. (2007) input set housed at the Virtual Radiation Belt Observatory (ViRBO).Results can be quickly visualized through built-in object methods.The RadBelt module is currently being expanded to include full three-dimensional diffusion and data assimilation capabilities.
Several built-in empirical relationships act as default sources and losses in RadBelt.The default electron source is an empirical acceleration term given by, where γ is the magnitude of the source (typically 0.1 days −1 ), L center is the center of the source curve (typically L = 5.6), and L width is the width of the source curve (typically 0.3).This function represents electrons being accelerated into the current µ − K slice as a function of magnetospheric activity.This simple source model was instituted to provide users with a quick, built-in source of phase space density; more accurate, physics-based source models will be developed in the future.five minute time resolution.L max was calculated using the Tsyganenko 2001 storm-time empirical magnetic field model (Tsyganenko, 2002;Tsyganenko et al., 2003).This magnetic field model was selected to balance accuracy with speed of execution, making the six-month calculation time feasible while returning reasonable results.The root-mean-squared error of this relationship is 0.64 R E .Loss due to plasmaspheric hiss (Millan and Thorne, 2007), important for developing the slot region, is included by using a loss lifetime of 10 days.The radial extent of the plasmapause is obtained via the Kp-dependent relationship defined in Carpenter and Anderson (1992).Though these relationships are the defaults, RadBelt is flexible and any arbitrary functions can be used.
An example of RadBelt output is shown in Fig. 1, which shows the results from a simulation of the well known "Halloween Storms" occurring from 20 October 2003 to 5 December 2003.The simulation was run using the Brautigam and Albert (2000) specification of D LL , the default source and losses described above, no density at the boundaries, and initially empty radiation belts.The top frame shows electron phase space density for µ = 2083 MeV G −1 and K = 0.03

√
G R E over the L domain.These values will be used throughout this study.The over-plotted white line shows the location of the last closed drift shell for this particular µ and K combination.The bottom frame shows the Dst and Kp indices during the simulated period.These values drive both the D LL and the empirical functions described above.
Throughout the simulated period, the last closed drift shell confines the outer radiation belt as phase space density outside of L max decays quickly (Fig. 1, top frame).During quieter periods, L max grows and the belts diffuse outward to fill the newly closed drift shells.The limitations of the simple source model are evident, as there is little variability  Formulation α β S1997 1.9 × 10 −10 11.7 BA2000 10 0.506Kp−9.32510.0 FC2006 1.5 × 10 −6 8.5 U2008 7.7 × 10 −6 6.0 within the belts and phase space density values remain elevated throughout the simulation.
The driving physics of any radial diffusion model is contained in the diffusion coefficient, D LL .Recognizing this, RadBelt presently includes several different D LL formulations.Adding a new formulation is as simple as defining a new function and assigning it as a RadBelt object attribute.The currently included D LL models are taken from Selesnick et al. (1997), Brautigam and Albert (2000), Fei et al. (2006), and Ukhorskiy and Sitnov (2008), denoted as S1997, BA2000, FC2006, and U2008 herein for convenience.Each of these has the generic form, The values of α and β are summarized in Table 1; each is compared in Fig. 2 over the typical spatial domain of Rad-Belt.
RadBelt uses a modified Crank-Nicolson implicit finite difference solver which is second-order accurate in space and time (Crank and Nicolson, 1947).Equation (1) has several characteristics that make the standard Crank-Nicolson approach inappropriate, namely a space (and potentially time) dependent diffusion coefficient as well as a right hand side that contains factors of L 2 in both the numerator and  8), for a smaller α factor, the source closely resembles the analytic solution, but scaled and phase-shifted (center frame).As α becomes stronger, the source becomes more complicated to counteract the strong diffusion at high L (bottom frame).
denominator.The traditional solver is modified by combining the factor of L −2 with D LL and discretizing the right hand side via a second order difference scheme that takes into account the spatial dependence of D LL , outlined in Press et al. (2007).The possibility of a time-dependent D LL is accounted for by computing it at the midpoint between t and t + t.This approach maintains unconditional stability of the scheme as demonstrated by Tadjeran (2007).The resulting discretized form of Eq. ( 1) is, where indexes n and j represent discretization in time and space, respectively.The derivation then follows the usual path of separating f n+1 terms from f n terms, resulting in the familiar arrangement, where A and B are tridiagonal matrices of coefficients that include factors of L, t, L, and D LL .In RadBelt, routines to invert A via standard LU tridiagonal matrix decomposition (Press et al., 2007) to advance the solution forward in time are written in the C programming language to obtain the fastest possible execution speed.The portions of SpacePy written in C are compiled to a shared object library which Python interfaces using its built-in C-Types module.

Verification technique
The prototypical code verification approach begins by obtaining an analytical solution to the equation or system of equations.The known solution is compared to the numeric solution in a systematic way covering a range of grid spacings and time steps, demonstrating the code's ability to properly converge to the correct solution as expected given the implemented solver.For complicated systems, however, obtaining an analytical solution becomes either prohibitively difficult or outright impossible.
One way to overcome this difficulty is to employ the Method of Manufactured Solutions (MMS) (Roache, 2002).This method begins by selecting a solution independent of the governing equations.The manufactured solution is substituted into the governing equation to produce a source term that, when added to the original governing equation, yields a new equation for which the manufactured solution is an exact solution.The manufactured solution can now be applied to convergence studies of the code, which, if passes, can be considered verified -it correctly solves the governing equations and properly applies the given numerical scheme (Roache, 1998) This process is illustrated clearly for the situation at hand.The arbitrary solution is selected with no regard for Eq. ( 1) outside of the variables we wish to exercise (t and L): where a and b will be set to satisfy prescribed boundary conditions.This "solution" is >2 times differentiable in both t and L, which is required once substituted into Eq.( 1).Substituting Eqs. ( 7) and (4) into Eq.( 1) yields the source term, For this exercise, Dirichlet boundary conditions will be used in L (set to zero for convenience using a typical domain of L = [1 : 10]) while the initial conditions will be no phase space density.This restricts the choice of a to 2π 9 .The constant b is then chosen to ensure a source that has a short period so that the verification simulations can be performed quickly.Given a simulation run to t = 600 s, b is set to 2π 300 such that two full periods are completed in one simulation.
The behavior of both the manufactured solution and the associated source term for the chosen values of a and b are shown in Fig. 3.As desired, the solution shows plenty of variation in just a short time period (Fig. 3, top frame).For a moderate D LL (S1997, whose behavior is shown as a black dashed-dotted line in Fig. 2), the first term of Eq. ( 8) dominates and the artificial source is similar to the manufactured solution, but phase-shifted and scaled (Fig. 3, center frame).However, for a strong D LL , the pattern in the source term changes drastically (bottom frame).This occurs for the BA2000 D LL when Kp is high.As α increases by several orders of magnitude, the L-dependent second and third terms of Q MMS (Eq.8) become dominant.The strong source term values at high L overcomes the strong, L-dependent storm time diffusion of electrons and maintains the shape of the manufactured solution and enforces the prescribed conditions.
With this manufactured solution and source term, the Rad-Belt model was run thousands of times with different combinations of t and L. Each run covered ten minutes of simulation time, or two full periods of the manufactured solution, over the typical radial domain of L = [1 : 10].This was repeated for each D LL model currently included in the code.plot.A quick comparison between the time steps used and the corresponding maximum error suggests that the code is converging at the expected order (second) given the selection of solver: when the time step is decreased by a factor of x, the error decreases by a factor of x 2 .Figure 5 illustrates this clearly and quantitatively.Time steps and associated error values are plotted on a log-log scale for 100 separate simulations (top frame).These results are naturally sorted into three distinct regimes.On the far right is the regime where the time step is too large to properly capture the dynamics.This is illustrated in the top row of Fig. 4, where t = 300 s misses most of the dynamics of the analytical solution and the error is greater than the amplitude of the forcing source term.Because the domain is resolved so poorly, the error varies wildly with spurious peaks and valleys developing.Such oscillations are caused by a variety of factors, such as the time step or grid resolution harmonizing with the spatial and temporal frequencies of Q MMS , but are ultimately the result of very poor resolution and should be disregarded.On the far left is the regime where the error is dominated by floating point error.Though the time step is refined throughout this regime, the maximum error is not reduced as it is either of machine origin or requires further refinement in L. This is observed in the bottom two

Max Error
Fig. 6.Error map for 8000 simulations as a function of time step (vertical axis) and grid size (horizontal axis).White lines denote contours of constant error at each order-of-magnitude.These results were obtained using the S1997 D LL model.A cut of this map along a line of constant L ( t) will yield a curve similar to what is seen in the top (bottom) frame of Fig. 5 with the slope indicating the order of convergence.The proper behavior is observed here over a large range of time steps and grid sizes.rows in Fig. 4 where the error reduction is less than expected given the second order solver.The center region of Fig. 5, delimited by the vertical dashed red lines, is known as the asymptotic regime.Here, the time step and grid size are both appropriate for the conditions being simulated and the error drops off as expected given the selected numerical scheme.The slope of this region is nearly equal to the order of the solver, indicating that the code is correctly solving the system as intended.The lower frame of Fig. 5 shows that the same is true when the grid resolution is refined.
Figure 6 takes these results one step further.The number of simulations is now increased to 8000, and the error for each t-L combination is sorted into a two-dimensional color map.The white contour lines show order-of-magnitude boundaries clearly.Figure 6 not only confirms the results from above, but also reveals more about the code's behavior for a particular D LL and Q MMS .For t > 1 s, refining the grid has almost no effect on the overall error.In this region, the contours are parallel to the L axis.This is because for a moderate value of D LL and a quickly varying source function, it is much more important to have a properly refined time grid as opposed to a fine L grid.Only when the periodic source term is properly resolved does reducing the grid size improve results.The curve that passes through the intersection of the vertical and horizontal contours (not shown) traces the most efficient combinations of t and L. This information is important to code users who want to maximize code performance in terms of controlling error and maintaining fast run speeds; understanding the balance between the rate of variation in time and space is key.Results for DLL=BA2000, Kp=9.0 10 -8 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 Max Error This process was repeated for all other D LL models.For BA2000, the simulations were performed for Kp = [0 : 9] for a total of fourteen separate verification sets.The results are summarized in Fig. 7. Returning to Fig. 2, it would be expected that the error maps for FC2006, U2008, and BA2000 for low Kp should appear similar to the error map for S1997 (Fig. 6) as the D LL models are similar in shape and have weak to moderate values.As seen in the top four frames of Fig. 7, this is indeed the case.As Kp increases, however, BA2000 D LL becomes far stronger than any other values investigated to this point.This is evident in the bottom two frames of Fig. 7, as the region where L refinement is most important (contour lines are dominantly parallel to the t axis) grows as Kp increases from 6 (lower left frame) to 9 (lower right frame).Again, it is important for code users to understand the balance between time dynamics and diffusive L dynamics when selecting time step and grid resolution.

Time dependent D LL
Real world applications of RadBelt will involve a timevarying Kp.When using the BA2000 model, which is Kp dependent, D LL will vary throughout the simulation time domain.Verification of RadBelt's ability to properly capture this new complexity must be performed.The new test case is summarized in the top frame of Fig. 8. Kp is allowed to change from 9 to 1 in a smooth, linear fashion over the simulated temporal domain (shown in green with the scale on the left).This results in an exponentially decreasing D LL (shown in blue with the scale to the right).Real Kp values vary as a step function with a three-hour window; it is often linearly interpolated to create a smoothly varying function when used as a model input.As such, the situation selected in the top frame of Fig. 8 represents conditions found in science applications.
The bottom frame of Fig. 8 illustrates that the code indeed converges properly as t decreases (constant L = 0.01) when a time-varying D LL is used.Because the maximum error is bounded by what can accumulate when D LL is at its maximum, the curve closely matches what is observed for a constant BA2000 D LL (K P = 9) (Fig. 7, lower right frame, values along L = 0.01).Similar tests for an impulsive Kp that jumps from 1 to 9 at t = 300 s were performed with similar, positive results (not shown).These results demonstrate that RadBelt properly handles a time-varying D LL model.

Alternative source function
Though unlikely (Roache, 1998), false positives (results that indicate the code converges properly even though it is not) may be possible when using the MMS method.To test the veracity of the results to this point, portions of the above exercises are repeated using a new manufactured solution.This Figures 9 and 10 summarize the results when using the second manufactured solution and S1997 D LL .Figure 9, similar to Fig. 4, demonstrates that the code does indeed converge on the analytic solution as t is progressively refined.The error values are much lower than before, likely a function of Q 2nd MMS (Eq.12) varying far more slowly in time than the initial Q MMS (Eq.8).With the overall error drastically reduced, residual error at high L (lower right frame) becomes evident.This occurs in the region where D LL is strongest.Because the new source term varies more slowly in time compared to the old, these L-dependent features can be seen at a higher t than before.In other words, the results for this new source term are more sensitive to L refinement than t.This feature is illustrated clearly in Fig. 10, where the region of contours that follow constant T values is reduced compared to what was seen in Fig. 6.The results demonstrate proper convergence, but emphasize that the source term must

Max Error
Fig. 10.Similar to Fig. 6 but showing results when the second, non-periodic manufactured solution is used (Eq.9).be taken into consideration when selecting time and space discretizations.

Code performance
Figure 11 benchmarks the code performance using a single core on an Intel Xeon 2.6 GHz CPU.The color map shows the time to complete each of the 8000 simulations for the BA2000 D LL (Kp = 9) set.The black contour lines show the same results in terms of simulation time (always ten minutes) divided by total CPU time.RadBelt runs very quickly on a modern machine, with most simulations finishing at >100times real time.These results demonstrate that it is possible to combine the benefits of a modern interpreted programming language with the light weight, fast number-crunching abilities of a compiled language.

Source term application considerations
The source term used thus far has been simple enough to include in the derivation of the Crank-Nicolson solver.This leaves nothing to the imagination when the source term is applied.In real-world, non-verification situations, more complex source terms may arise that require a more sophisticated application, for example, a non-linear source term that depends on the phase space density.A popular approach is to split out the source term as one would when applying Strang splitting (Strang, 1968) to a multi-dimensional problem.This approach splits the differential equation into two parts: a partial differential equation and an ordinary differential equation.The two subproblems, the first representing the diffusion of f and the second representing the source of f , can be solved in an alternating fashion that is mathematically equivalent to solving the whole system (Toro, 1999).A second order accurate method for this approach can be written as, (seconds) (seconds) Fig. 11.Run speed benchmark for all 8000 simulations performed using BA2000 D LL (Kp = 9).The color map shows raw simulation completion time (seconds) while black contour lines show the ratio of simulation time to CPU time.For example, for a total simulation time of 600 s and a CPU run time of 6 s (near-white on the color scale), the code simulates the radiation belts 100 times faster than real time and has a simulation-to-CPU time ratio (black contour lines) of 100X.
where f is the phase space density, S t 2 is the operator that advances the source subproblem forward in time by onehalf time step, and C t is the Crank-Nicolson operator that advances the diffusion subproblem forward in time a full time step.Such splitting yields a powerful, generalized algorithm for tackling complicated source terms but requires additional considerations not previously taken into account, namely careful selection of the source term solver.
Figure 12 illustrates the impact of the selection of source solver if splitting were applied to the artificial source terms used in this work.Source splitting is not necessary for Q MMS as it is only a function of t and L, but is applied to illustrate the complications that arise when it is used.The blue line labeled "Case 1" shows that employing a simple, second order accurate trapezoidal integration to solve the source subproblem yields the desired results: a model that converges cleanly with the prescribed accuracy.But what if a higher accuracy quadrature method is used?Case 2 uses the fourth-order accurate Simpson's Rule to integrate the source term; the results are displayed as the green line in Fig. 12.The error is reduced from Case 1, however the code still converges in a second order manner.This is because less error is introduced from the integration of the source term but the total error is still dominated by the Crank-Nicholson solver and converges along that lower-order rate.Finally, Case 3 (red line) is identical to Case 2 but the magnitude of D LL is reduced by an order of magnitude.Because the rate of diffusion is smaller, a larger time step is better suited to the problem and the overall error should drop.This proves true for t < 30.0 s in Fig. 12, but for larger time steps, the error quickly rises to match Case 2. What is observed is the two solution operators trading the bulk of the error.When the source term solver accounts for the greater portion of the error, the code converges with fourth order accuracy.Because this drops below the error inherent in the diffusion operator quickly, the code then begins to converge more slowly, and all three lines become parallel.All three of these curves result from slight variations of a Crank-Nicolson, split-source term implementation.
This example illustrates the additional complexities that must be taken under consideration when more complicated source and loss terms are included.The selected source integration scheme has dramatic implications on the overall error.Changing diffusion rates, which happens commonly when the BA2000 model is used, can push the model into different error regimes and change the behavior of the model.These effects must be kept in mind as the model is developed in the future.

Conclusions
This work introduced and verified the RadBelt radial diffusion model that is part of the SpacePy software library.The method of manufactured solutions was employed to provide an analytical solution where one could not be trivially obtained.The model's flexibility, especially in terms of choice of D LL and sources, was put to many rigorous tests.It was found that the solver and many code features have been implemented correctly and robustly.
The results here are of special interest to code users.Repeatedly, it was demonstrated that code convergence is beholden to the choice of D LL and the rate of change of the source function.While selecting t and L, it is easy to find oneself in a regime where refining one has little impact on the final error.How the source and diffusion coefficient affect these semi-stagnant convergence regimes is vital knowledge for code users who want to maximize code performance while minimizing error.
Finally, it must be emphasized that performance in terms of convergence regimes and overall error is situation specific.The time and grid discretizations used in this study are far smaller than the typical scientific user will require.It is suggested that, when employing new source, loss or diffusion coefficient models, the user perform a quick study to ensure proper refinement for the most extreme situation they choose to encounter.

Fig. 1 .
Fig. 1.Example output from a RadBelt simulation of the well known "Halloween Storms".The top frame shows phase space density with L max location over plotted in white.The bottom frame shows the geomagnetic indices used to drive D LL and the empirical source and loss functions.

Fig. 2 .
Fig. 2. Comparison of the four formulations of D LL used in this study.Because the formulation of Brautigam and Albert (2000) is Kp dependent, it is shown for several different Kp values.

Fig. 3 .
Fig. 3. Behavior of the manufactured solution (top frame), Q MMS when S1997 D LL is employed (center frame), and Q MMS when BA2000 D LL (Kp = 9) is used (bottom frame).Following Eq. (8), for a smaller α factor, the source closely resembles the analytic solution, but scaled and phase-shifted (center frame).As α becomes stronger, the source becomes more complicated to counteract the strong diffusion at high L (bottom frame).

Fig. 4 .
Fig. 4. A demonstration of the code converging towards the analytical solution as the time step is refined.The left column shows the numerical solution while the right column shows the corresponding error.Because the maximum error, listed in the title of each plot, reduces as the square of the time step, this figure hints at the expected convergence rate given the numerical scheme used.

Figure 4 Fig. 5 .
Figure4illustrates the model converging towards the analytical solution as t is refined for a constant L = 0.1 with S1997 D LL .The left column displays the solution for five different values of t, decreasing from top to bottom, while the right column shows the corresponding error (defined as the magnitude of the difference between the analytical solution and the numerical solution) throughout the domain.The maximum error for each is listed at the top of each error

Fig. 8 .
Fig. 8. Results for time verification when using a time-varying D LL .L = 0.01 for this set.The top frame shows the values of Kp (green line) and BA2000 D LL (blue line).The bottom frame shows the expected convergence of the model in the same manner as the top frame of Fig. 5.

Fig. 12 .
Fig. 12.Convergence curves for three separate simulations that split out the source term from the Crank-Nicolson solver.Case 1 uses a simple second-order accurate source solver, while Cases 2 and 3 use a fourth order source solver.The variety in the curves demonstrate the importance of properly selecting a solver along with an appropriate grid size and time step.

Table 1 .
A comparison of the different D LL formulations used in this study in terms of α and β.