Shock-induced chemistry and high strain-rate viscoelastic behavior of a phenolic polymer

The high strain rate behavior of phenolic is relevant to a wide range of reinforced plastics,^1–6 including phenolic foams and aerogels,⁷ and as binders in numerous laminates, composites, abrasives, and refractories.^8–11 Accelerated composite development has been partially driven by the growing commercial space market concerned with advanced processing and reliability under extreme conditions encountered during rocket load, planetary reentry, or orbital debris impact,^12–18 as well as aviation.^2,8 Extreme impact scenarios involving phenolic composites can include encounters of spacecraft by micrometeorites with velocities up to ∼70 km/s.¹⁸ 

Despite the fact that composite performance is intimately linked to interphase properties, studies of the high-rate compression of phenolic polymers are comparatively scarce.^19–21 Numerous studies of shock compression have been conducted for other cross-linked aromatics, including epoxies.^2,20,22,23 Investigations of the quasi-static behavior of phenolic have indicated sensitivity to the choice of resin and curing method.^24–28 For example, Hartmann²⁶ suggested a strong relationship between the acoustic wave speed and density and little relationship between wave speed and the extent of cross-linking, although the corpus of work indicates variation with respect to all three properties and is insufficient to define a correlation.^19–21 

Typical cure reactions can include the saturation of phenol with methylol or hydroxymethyl groups, dibenzyl ether linker formation via dehydration reactions, and breakdown of dibenzyl ether linkages to form methylene linkages.^29–31 The number of cross-links and ratio of linker types is, thus, determined by the starting resin and cure process, leading to a potential source of significant variability in phenolic composite properties. Dynamic relaxation during compression involves both fast relaxations (e.g., dilation, side-group rotation, scission) and slow relaxations (e.g., chain reptation) and may be sensitive to these structural differences. Arman et al.³² used molecular dynamics (MD) simulations to explore viscoelastic/viscoplastic mechanisms during shock compression of phenolic; however, the disparate timescales necessitate experiments to unify the quasi-static and shock regimes. Beyond mechanical strain, several aromatic polymers have also been shown to undergo reaction under sufficient shock compression,^33–37 while mechanistic details are generally lacking.²² Indeed, the reaction of the phenolic binder has been shown to influence the viscoelastic behavior of some composites.³⁸ 

We report here a comprehensive investigation of phenolic behavior under high strain rate compression (∼10⁵ s⁻¹) by exploring both unreactive and reactive regimes. First, experiments and finite element simulations are used to show that viscoelastic behavior observed at ordinary laboratory rates (≲1 s⁻¹) extends to the high-rate regime for stresses up to 1.2 GPa and to determine the spall strength. These findings are used to interpret shock Hugoniot experiments, which are compared to MD simulations, from which the Hugoniot elastic limit is also determined. Finally, the MD technique is used to confirm a chemical transition on sufficient shock compression and to provide insight into the densification mechanism for stresses up to 60 GPa. This study provides an essential basis for interpreting the high strain rate response of phenolic phases, needed for designing composite materials, with implications for understanding the shock compression of aromatic polymers.

All samples were derived from the same batch of resole-type formaldehyde phenolic resin (PSR133, Barrday, Millbury, MA; Lot No. J47740) typical for high-temperature applications. Large (∼50 mm diameter), homogenous pucks free of visible vapor inclusions were produced by pouring resin into aluminum pans, allowing excess solvent to evaporate at room temperature for several hours, and then curing for 6 days with temperature stepped over 80–95 °C, followed by a 24 h post-cure to achieve 95%–98% of the fully cured state as assessed by the residual cure enthalpy (differential scanning calorimetry) and with residual volatiles <3 wt. % (thermal gravimetric analysis). Specimens were cut and diamond-turned on both faces using a room temperature, solvent-free process to maintain the original material structure and remove oxidized exterior layers. Samples ≤∼1 mm thickness were selected to ensure that only targets with full density (absent inclusions) were used in the study.

Phenolic density of ρ₀ = 1252 ± 4 kg/m³ (avg. ± s.d.) indicated sample-to-sample reproducibility, along with pulsed ultrasound (5–10 MHz) measurements of longitudinal c_L (3.05 ± 0.02 km/s), transverse (1.41 ± 0.01 km/s), and bulk (2.58 ± 0.02 km/s) sound speeds, Poisson's ratio (0.364 ± 0.003), shear modulus (2.50 ± 0.03 GPa), bulk modulus (8.34 ± 0.15 GPa), and Young's modulus (6.81 ± 0.06 GPa). Differential scanning calorimetry was used to confirm a typical glass transition temperature of T_g ≈ 150–170 °C and a linear expansion coefficient of 57 ± 4 ppm/K, and calculated Grüneisen coefficient of 0.88 ± 0.08.²⁶ Activation energies of 75 ± 4 kJ/mol for residual cure and 146 ± 26 kJ/mol for thermal degradation (TGA) were also as expected.^39–41 

Impact tests were performed on the single-stage gas gun at the Dynamic Integrated Compression Experimental (DICE) Facility (Sandia National Laboratories, Albuquerque, NM) by mounting impactors to ∼1.5 cm thick Al6061 flyers. Four shots were performed in the low-stress (≤1.2 GPa) regime to investigate the shock Hugoniot, stress-wave transmission, and spallation of phenolic. Each shot involved simultaneous shock compression and measurement at three specimen stations to maximize the number of interrogations per shot (Fig. 1 and Table I). Impactor velocity and tilt relative to the target assembly were determined from arrival-time signals from an array of electrically charged pins located around the target periphery. In the forward-ballistic (F) configuration, phenolic targets were backed by poly(methyl methacrylate) (PMMA). One specimen was unbacked to permit the measurement of spall strength. In the reverse-ballistic (R) configuration, a phenolic impactor simultaneously struck three observation window materials to yield six measurements of shock Hugoniot states in two shots.

The front face (i.e., closest to the impact plane) of each window was coated with a thin aluminum film that reflected laser light used for time-resolved velocity interferometry (VISAR)^42–44 measurements of the motion at the target/window or impactor/window interface. For the non-windowed specimen (shot F1, station B), the reflective aluminum coating was applied directly to the sample and the sample's free-surface velocity was tracked by VISAR. Velocity amplitudes were adjusted on the basis of optical corrections that have previously been published for the individual window materials, namely, PMMA,⁴² lithium fluoride (LiF),⁴⁵ and sapphire.⁴⁶ The initial time, $t 0$⁠, for each velocity history was set at the first detectable motion of each reflective surface. A dual delay VISAR configuration featuring two different velocity-per-fringe (VPF) constants of 0.090 and 0.050 km/s was used for these shots, allowing two phenolic Hugoniot points to be reported for each impact condition. Slight variations are indicative of uncertainties in the interferometric measurements.

Hugoniot states were calculated using a standard impedance-matching analysis,⁴⁷ based on (i) the measured impact velocity, (ii) the magnitude of the initial interface-velocity plateau, and (iii) published Hugoniot data as follows. For sapphire, previously published linear fits with uncertainties were used as-is,⁴⁶ i.e., $D$ $=0.89( ± 0.07)u+11.208( ± 0.013)$ for u in km/s. Although sapphire exhibits nonlinear elastic behavior in this regime, its Hugoniot is well-defined.^42,46 For LiF, previously published linear fit parameters were used;⁴⁸ however, it was necessary to refit the data to estimate parameter uncertainties, yielding $D=1.353( ± 0.012)u+5.148$ $( ± 0.030)$⁠. For PMMA, data from⁴² were slightly nonlinear and can be approximated by the second-order polynomial fit, $D=−5.063$ $( ± 0.577) u 2+3.099( ± 0.205)u+2.758( ± 0.013)$ for $u≤0.32 km / s$⁠. For (ii), velocity histories along the shock plateau, i.e., after equilibration but before release, were fit to a line to determine an average particle velocity. The impedance matching followed the analytical procedure of Ref. 49 with the exception that it was necessary to solve for the reflected wave intercepts numerically for the nonlinear PMMA Hugoniot. Uncertainties in window Hugoniot parameters, impactor density (±0.3%), window density (∼±1%), flyer velocity (∼±1 m/s), and measured plateau velocity (∼0.1–1.0 m/s), were propagated through the impedance matching analysis by evaluating the sensitivity of the solution with respect to each parameter,⁵⁰ yielding uncertainties for each phenolic Hugoniot state parameter.

We use classical molecular dynamics (MD) simulations to investigate high-pressure material response in phenolic resins, looking to differentiate signatures of pressure-driven chain compaction and shock-induced chemistry. Simulations were performed using Sandia's LAMMPS code,⁵¹ with the ReaxFF force field,^52,53 over a range of shock pressures ranging from 300 MPa to 60 GPa. This method has been used widely to explore the thermal decomposition of polymers.^54,55 For this work, two sets of ReaxFF parameterizations were used. “MD1” simulations used parameters recently modified to improve the thermal expansion behavior of phenolic.⁵⁶ “MD2” simulations used the CHO parameters, which are well known to capture shock-driven chemistry.⁵² 

Cross-linked phenolic structures can be complex. Experimentally, they are produced by the reaction of phenol monomers with formaldehyde to first saturate the ortho and para positions with hydroxymethyl or methylol groups. Cross-linking extent depends primarily on the formaldehyde to phenol ratio but also on the secondary reaction of methylol groups reacting to form a dibenzyl ether or a methylene linkage. Longer reaction times lead to further converted linkages to methylene.^29–31 This experimental end-state is our target simulated structure [Fig. 2(a)].

Cross-linked structures are produced using the bond/react procedure,⁵⁷ non-reactive PCFF,^32,58 and Lennard–Jones potentials.⁵⁹ The cross-linking is facilitated by temperature and pressure in three stages: (1) 2000 K and 0 GPa for 5 ps, (2) 1000 K and 2 GPa for 8 ps, and (3) 1300 K and 2 GPa until the desired cross-link extent is reached. The non-reactive potentials allow cross-link sites to be brought into proximity with pressure and temperatures that might otherwise cause phenolic decomposition. After the cross-link bonds are formed, the structures are thoroughly equilibrated with the ReaxFF potential to relax any residual stresses. After a brief push-off using LAMMPS fix nve/limit, these equilibrations are completed in stages: (1) 300 K and 2 GPa for 5 ps and then 700 K for 5 ps, (2) then temperature is ramped down to 300 K over 2.5 ps and held for another 2.5 ps, (3) then pressure is ramped down to ambient using an NPT ensemble at 300 K for 2.5 ps, and (4) then a system is equilibrated for 5 ps under ambient conditions. Any volatiles that were produced are removed and the final 5 ps equilibration repeated, as necessary.

The final cross-linked structure contained ∼4 × 10⁴ atoms with density 1255 kg/m³, importantly near the experimental density. A second simpler linear structure consisting of 432 linear phenolic chains with eight monomeric subunits, ∼5 × 10⁴ atoms, was produced with a density of 1400 kg/m³ for comparison [Fig. 2(b)]. These system sizes were chosen such that linear chains were no more than a third of the total box size, to avoid self-interaction and finite size effects. For the cross-linked system, a balance was made between accuracy and efficiency consistent with recent research on system sizes in epoxy systems.⁶⁰ The linear systems were annealed near ambient pressure at 700 K and then equilibrated for 150 ps, followed by a thermal ramp down by 50 K at a rate of 1 K/ps—repeated until systems reached 300 K.

To obtain shock states, we employ the equilibrium, non-propagating constant stress Hugoniostat method of Ravelo et al.,⁶¹ wherein the system is compressed uniaxially until the target pressure is reached and the system is thermostatted dynamically to satisfy the Rankine–Hugoniot relations. This method, in combination with ReaxFF, has been used successfully in the shock simulation of other polymers, e.g., Ref. 53. The systems were shocked to various pressures using temperature and pressure damping coefficients of 5 ps and 0.05 fs timesteps. These simulations are expected to reasonably model experiments, which are reasonably homogenous over the scale of the experiments. No parameters of the MD model were tuned other than to ensure that all starting structures had initial density comparable to the experiments.

Shock waves generated from sapphire impact transmitted through the phenolic of different thicknesses were seen to decrease in the amplitude with increased sample thickness (Fig. 3). In addition, the terminal region of the velocity plateau shows dynamic behavior indicative of shock state evolution. This “rounding” of shock wave profiles is a characteristic of nonlinear waves with dissipation,⁶² such as from viscoelasticity,^63,64 influencing interpretation of shock states.⁶⁵ While additional data would be required to differentiate the response of each thickness, the three waveforms are nonetheless consistent in showing this behavior.

To show that the dissipation at these strain rates of ∼10⁵ s⁻¹ can be explained by viscoelasticity—a known characteristic of phenolic at strain rates ≲1 s⁻¹ (Refs. 27 and 28)—simulations of stress-wave propagation were performed using CTH v.10.3,⁶⁶ with 1-D rectilinear meshes representing the flyer, impactor, target, and window. CTH is a three-dimensional, multi-material Eulerian hydrocode capable of modeling high strain rates characterized by high velocity impact. All materials in this analysis were modeled with a pressure-dependent Mie–Grüneisen equation of state formulated from Hugoniot curve parameters including the measured Hugoniot for phenolic (see Sec. IV) and reference values for sapphire,⁴⁶ and for the Al6061 in the bounding fixture.⁶⁶ 

To parameterize the M–B portion of the model, we measured the storage and loss moduli for the phenolic over a temperature range of 140–280 K and frequency range of 0.2–20 Hz using Dynamic Mechanical Analysis (DMA). Linear regression of the DMA data was used to determine the shift of the β-transition as 9.2 ± 3.4 °C/decade in the strain rate, which was somewhat lower than the values of 15.3 and 25.2 °C/decade for PMMA and polycarbonate, respectively.⁶⁷ With this value for phenolic fixed, we then varied the shear modulus vs effective temperature relationship, using piecewise linear splines to obtain a set of ${ C β , 1 , C β , 2}$ parameters for each spline, as a numerical representation for the temperature and strain rate dependence of the shear modulus. All other parameters for the viscoelastic model were taken from Refs. 67 and 68 for PMMA, noting that Poisson's ratio measured for the phenolic (0.364 ± 0.003) was comparable to expected values for PMMA (0.32–0.36) over microsecond relaxation times and temperatures up to ∼125 °C.⁶⁹ 

The shear modulus representation, shown in the inset in Fig. 4 in comparison with PMMA, was varied to obtain the best fit to the experiment (Fig. 4). We found that the location of the “knee” in the shear modulus vs effective temperature relationship influences the strain rate required to see a significant change in the yield behavior under shock compression of phenolic, as seen also for PMMA.⁶⁷ To further refine the model, the dynamic spall strength was varied within the CTH model over the range 82–112 MPa, with the best overall fit to the plate impact data given by a value of 102 MPa. This value derived from the model incorporating viscoelastic effects is 11% larger than that determined directly from our data using the pull-back signal method ignoring viscoelasticity and nonlinear compressibility.⁷⁰ While the direct determination method is known to be approximate,⁷⁰ the agreement nonetheless indicates overall consistency. The sharpness of the pullback signals may arise from the use of comparatively thin samples.^70–72 From the CTH model, we infer that the spall plane was located ∼28 μm from the front surface of the target such that the spall scab had an estimated thickness of 0.764 mm or ∼97% of the original target thickness. The model successfully captures the inelastic deformations on both initial shock compression and spall fracture of a non-windowed phenolic sample, as well as the associated stress-wave reverberations within the flyer, target, and spalled sections as exampled in Fig. 4. We note that while the shock Hugoniot, shear modulus, and spall strength used in the CTH model either derive from, or are related to, the plate impact data, all other parameters, including the strain rate dependence of the β-transition, were determined through independent measurements. These results show the extensibility of the viscoelastic behavior of phenolic across strain rates of ∼1–10⁵ s⁻¹ and justify extracting the Hugoniot state following velocity equalization, which we show in Sec. III B.

Having demonstrated viscoelasticity as the primary relaxation mechanism in shock compressed phenolic, we proceed to identify the equilibrium shock state from the terminal velocity reached along the shock plateau and by applying the Rankine–Hugoniot jump conditions, with shock velocity D, particle velocity u, shocked density $ρ H$⁠, and Hugoniot stress $σ H$ (Table II). Reasonable consistency is seen between each pair of Hugoniot states determined with different VPF constants. The spread in values provides an estimate of measurement uncertainty.

The MD simulations in the low-pressure regime demonstrate reasonable agreement with the experimental data (Fig. 5). Simulations using several other initial densities were also performed, resulting in poorer agreement (omitted for brevity). From simulations, we can also infer the maximum temperature achieved in the impact experiments was ∼338 K, low enough to preclude any appreciable chemical reaction.⁷³ By observing the presence or absence of shear stress relaxation in the temporal profiles from the simulated shock compression, we estimate the Hugoniot elastic limit (HEL) as between 1 and 2 GPa, consistent with the ∼1.5 GPa plasticity threshold determined by Arman et al.,³² for SC1008 phenolic by evaluating shear stress relaxations behind the shock front within MD simulations, suggesting that only viscoelastic deformation occurs in the experiments. We also note that these HEL values are similar to earlier experimental observations, e.g., for epoxy resin (∼0.6 GPa),⁷⁴ PMMA (0.7–0.9 GPa),⁴² and polyether ether ketone (∼1.0 GPa).⁷⁵ 

In Fig. 6, the data are plotted with respect to relative compression $η=1− ρ 0/ ρ H$ to facilitate comparison to phenolics of different initial densities.^19–21 Because the measurements of Carter and Marsh²⁰ were conducted over 3.8–66.5 GPa, their data at higher pressure were fit to a polynomial and extrapolated to lower pressure for the purpose of comparison (dashed line in Fig. 6). Data from the MD simulations are tabulated in  Appendix, Table III. On balance, these results show that the choice of starting resin and cure process, which affect polymer structure and density, can significantly influence the shock compression behavior of phenolic. The two sets of MD simulations, which also reflect different cross-link structures, similarly show significant differences.

The relative residuals had a mean of ±12% (bidirectional), which restricts prediction to a span of ∼24% of the local mean absent knowledge of the phenolic structure. Within this variation, our experimental data for PSR133 phenolic is in reasonable agreement with prior work including measurements for Durite HR 300-Borden phenolic examined by Carter and Marsh,²⁰ FF-17 phenolic examined by Dandekar and Lamothe,²¹ and SC1008 phenolic examined by Wood et al.¹⁹ 

When our Hugoniot data, which extends to lower shock stress, are reparametrized using the jump conditions, a first-order transition in the shock velocity—particle velocity relationship in the limit of weak compression becomes apparent that has not been previously observed (Fig. 6, inset), presumably due to an ordering process. It can be shown that this discontinuity can manifest as a deviation from a simple quadratic in the $σ H$— $η$ relationship (c.f. fit to Carter and Marsh²⁰ in Fig. 6), however, the discontinuity is accentuated in the D-u plane (Fig. 6 inset). Excellent fits to the linear form $D= c 0+bu$ were obtained above this transition with b = 1.24 ± 0.15 and c₀ = 3.16 ± 0.03 km/s. A fit below the transition gave b ≈ 0 within uncertainty and c₀ = 3.09 ± 0.03 km/s, similar to the bulk sound speed as seen for many materials at small compression.⁷⁶ Similarly, Carter and Marsh²⁰ report fit parameters of b = 1.39 and c₀ = 2.98 km/s for 0.6 ≤ u ≤ 2.6 km/s (0.16 ≤ η ≤ 0.39) and ρ₀ = 1385 kg/m³. Wood et al.,¹⁹ fit Hugoniot data for SC-1008 phenolic to the form $D$ $= c 0+bu+c u 2$ for u ≤ 0.8 km/s, which was approximately linear with c₀ = 2.16 km/s, b = 3.29, and c = 0 for u ≤ 0.3 km/s and ρ₀ = 1180 kg/m³. A similar fit by Dandekar and Lamothe²¹ gave c₀ = 2.970 km/s, b = 3.225 and c = −2.726 for u ≤ 0.5 km/s and ρ₀ = 1290 kg/m³. These results involving different phenolics are in coarse agreement, with differences in shock stress ≲30% arising from the choice of resin and cure process.

To consider the upper limit of phenolic behavior, e.g., for micrometeorite impacts, the shock Hugoniot of phenolic was predicted using MD for shock stresses up to ∼60 GPa. Example results are shown in Fig. 7(a) at multiple times following the initial shock to various pressures. Simulations of several cross-linked phenolic structures with methylene linkages with or without methylol side-groups gave similar results. Importantly, the shocked density was found to evolve over 10–100 s of ps at stresses ≳20 GPa (≳1600 K). This transition, which occurs at stresses above the instantaneous plastic limit, implies a rate-dependent densification mechanism. The simulations are in excellent agreement with measurements by Carter and Marsh²⁰ [Fig. 7(a), inset] up to this pressure, also implying a departure from steric compression.

The phenolic evolution above approximately 20 GPa can be related to chemical reaction kinetics as witnessed by the population of carbon coordination number C_i [Fig. 7(b)]. Whereas most carbons are initially bound to two other carbons (C₂), additional cross-linking on shock compression results in an increase of C₃ and C₄ populations. The average carbon coordination number is strongly correlated with the increase in the density through the evolution of the shock state [Fig. 7(c)].

Consideration of the C–C–C angular distribution functions (ADFs) in Fig. 8(a) indicates a transition from the initial planar, aromatic, sp2 structure to tetrahedral sp3 configurations, as also seen during shock compression in benzene.³⁵ The primary chemical mechanism facilitating the densification appears to be dehydrogenation followed by the recombination of carbon radicals, as shown in Figs. 8(b) and 8(c). For example, at 40.5 GPa, this occurred within 0.1 ps for 75% of dehydrogenation events. Other reactions also occur, notably the formation of H₂O from dehydration reactions between neighboring H and OH groups, and the production of other light volatiles, including CO, CO₂, H₂, and CH₄. The simulations, thus, indicate that the local bonding environment shifts from cross-linked phenolic rings to a denser, non-aromatic structure on sufficient shock compression.

These findings explain the discontinuous slope in the shock Hugoniot of phenolic at 23.2 GPa observed experimentally by Carter and Marsh²⁰ [Fig. 7(a), inset], who reasoned by comparison with other aromatic polymers that a chemically driven phase transition should occur. A discontinuity, or “cusp,” has been observed around ∼20 − 30 GPa, and attributed to chemistry, in several shock-compressed resins and carbon-resin composites,^20,77–80 and other polymers.^{20,77,81–83} We note, however, that at the highest shock stress in our simulations, equilibration required more than 250 ps, approaching a practical computational limit, whereas reaction times up to 100 s of ns following shock compression of other polymers have been observed experimentally.^22,81,83,84 Consequently, the existence of other reactions and/or bulk relaxation processes remains plausible and may explain the discrepancy with their experiments at the highest pressures.

Our study finds that the shock compression of phenolic at high strain rates (∼10⁵ s⁻¹) and low pressures can be characterized by steric compaction, viscoelastic wave dampening, spall fracture (0.10 GPa), first-order Hugoniot transition (∼0.36 GPa), and viscoplastic yield (1–2 GPa), all independent of chemistry. Through comparison to the literature and our own simulations, we have shown the shock Hugoniot to be highly sensitive to the choice of starting resin and cure process. This sensitivity could propagate in determining properties of phenolic-based composites (e.g., density, cross-link structure, and shock compressibility). In contrast, chemical processes become the primary driver at pressures above approximately 20 GPa of the rate-dependent relaxation process during shock compression. The densification is facilitated by dehydrogenation and dehydration reactions, which allow carbon atoms to recombine into denser, non-aromatic structures with interstitial volatiles. These findings provide an essential basis for interpreting the high strain rate response of phenolic phases to support composite design and assessment for the growing number of commercial space applications relying on phenolic composites, while highlighting compression mechanisms relevant to other polymer systems.

We thank Benjamin Anderson, Mark Stavig, Gary Chantler, Dennis Lobley, Kevin Rolfe, Randy Hickman, and the DICE Facility team for experimental support. Partial support was provided by the Laboratory Directed Research & Development program at Sandia National Laboratories. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy's National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC, under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns all right, title, and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (https://www.energy.gov/downloads/doe-public-access-plan).

The authors have no conflicts to disclose.

Nathan W. Moore: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Keith A. Jones: Formal analysis (equal); Investigation (equal); Software (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Jack L. Wise: Formal analysis (equal); Investigation (equal); Validation (equal). Darren G. Talley: Formal analysis (equal); Software (equal); Validation (equal). J. Matthew D. Lane: Formal analysis (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author on reasonable request.

Table III