Research Article  Open Access
A Designer’s Approach for Estimation of NuclearAirBlastInduced Ground Motion
Abstract
A reliable estimate of freefield ground displacement induced by nuclearairblast is required for design of underground strategic structures. A generalized pseudostatic formulation is proposed to estimate nuclearairblastinduced ground displacement that takes into account nonlinear stressstrain behaviour of geomaterials, stressdependent wave propagation velocity, and stress wave attenuation. This proposed formulation is utilized to develop a closedform solution for linearly decaying blast load applied on a layered ground medium with bilinear hysteretic behaviour. Parametric studies of closedform solution indicated that selection of appropriate constrained modulus consistent with the overpressure is necessary for an accurate estimation of peak ground displacement. Stress wave attenuation affects the computed displacement under low overpressure, and stressdependent wave velocity affects mainly the occurrence time of peak displacement and not its magnitude. Further, peak displacements are estimated using the proposed model as well as the UFC manual and compared against the field data of atmospheric nuclear test carried out at Nevada test site. It is found that the proposed model is in good agreement with field data, whereas the UFC manual significantly underestimates the peak ground displacements under higher overpressures.
1. Introduction
Underground siting of strategic structures [1] is an option to enhance the safety against nuclearairblast [2]. Nuclear explosions generate rapidly moving airoverpressures capable of producing significant ground displacements [3]. Most severe loading occurs within a close vicinity around the ground zero (GZ) known as superseismic zone. This zone is subjected to significantly high overpressures. The velocity of moving airoverpressure fronts is also more than the Pwave velocity of the ground [3–6]. Hence, a reliable estimate of nuclearairblastinduced ground displacement is required for design in superseismic zone [5, 7, 8]. Several studies using numerical approaches that account for realistic stressstrain behaviour and boundary conditions are reported for calculating the freefield response [9–16]. Such studies require specialized expertise in numerical tools and are unattractive to practicing engineers. Before the advent of computational tools, Wilson and Sibley [17] conceptualized a onedimensional pseudostatic approach to estimate airblastinduced ground displacement. Whitman [4] and Baron et al. [18] have also shown that airblastinduced ground motion is predominantly onedimensional and vertical in superseismic zone. In the present work, a generalized onedimensional pseudostatic formulation is developed to estimate nuclearairblastinduced vertical ground displacement that accounts for (i) nonlinear stressstrain behaviour [e.g., 19, 20], (ii) stressdependent wave propagation velocity, and (iii) stress wave attenuation. Using the proposed formulation, a closedform solution is developed for linearly decaying blast load with negligible risetime applied on a layered ground medium with bilinear hysteretic stressstrain behaviour. The proposed model is validated against the field data of an atmospheric nuclear test conducted at Nevada test site [21], and various parametric studies are carried out. In addition, performance of the closedform solution is also compared with the model given in UFC [22].
2. Generalized Formulation
The developed formulation is explained in detail below.
2.1. Determination of Overstress Distribution (due to AirBlast) in the Ground
Airoverpressure generated due to nuclearairblast consists of an initial rising portion followed by a decaying portion (Figure 1). Overpressure timehistory can be visualized as consisting of multiple overpressure fronts. These overpressure fronts are transferred to the ground as Pwaves. On arrival of an airshock front at the point of interest, an initial pulse travels with seismic Pwave velocity and reaches a depth at time . Subsequently, other fronts arrive and penetrate into the ground. If the location of peak overpressure front at time is given as , then, the ground above the depth experiences compressive stresses lower than those caused by the peak overpressure front at that depth. Therefore, the depths above are referred to as “unloading zone.” However, the peak overpressure front is yet to reach at locations deeper than , and therefore, the zone lying beneath the depth is denoted as “loading zone” (Figure 1). It is noted that “unloading zone” corresponds to the “decay portion” of overpressure, whereas “loading zone” corresponds to the “rising portion.” Therefore, at times earlier than the risetime, there does not exist any unloading zone as overpressure fronts from the decay portion do not arrive before the end of risetime.
To determine the stress distribution in the ground at time , contribution from all stress fronts arrived before time is considered. If the kth overpressure front arrives at the point of interest on ground surface at time such that , then the kth overpressure front reaches at depth at time . Stress fronts propagating through ground are affected by attenuation (caused by hysteresis losses, viscosity of ground materials, and dispersion of energy in 3dimensional space) and interference with reflected wave fronts (generated due to impedance variation with depth). In the present study, stress wave attenuation is taken into account through a geometrical attenuation parameter, whereas interference of incident and reflected stress waves is ignored. If the attenuation coefficient at depth is given as , then the stress generated at depth can be written as , where is the magnitude of the kth overpressure front. Depth of the penetration depends upon wave propagation velocity of the stress front, which in turn depends upon the tangential modulus of the geomaterial at the stress level at the depth of interest. If wave velocity of the kth stress front is denoted as (a function of depth “z”), then (1) can be solved to obtain as a function of , , and (2),where is the function obtained by solving (1). It is noted that for a given , is effectively a function of only because is also a function of . As soon as the airblast slaps the ground surface, ground deforms elastically and an elastic wave is generated which propagates with seismic Pwave velocity. Upon subsequent arrival of higher overpressure fronts, inelastic waves propagate at velocities smaller than Pwave velocity. If the ratio of Pwave velocity to wave velocity corresponding to the kth overpressure front at depth “z” is denoted as , then at initial stress levels is equal to one. With increasing stress level, decreases [17, 23]; hence, it is logical to assume that increases with increasing stress level and the maximum value is attained at a peak stress level corresponding to the peak overpressure front. Therefore, a simplest choice is to assume a linearly increasing in direct proportion to overpressure (3). Overpressure fronts from the decay portion propagate through the media which is already stressed in a nonlinear range due to passage of the peak overpressure front. Therefore, it is assumed that wave velocity of overpressure fronts in the decay portion is the same as the peak overpressure front, and hence is the same as :where subscript “z” represents the depth coordinate. Using the average loading rate as , can be approximated as (3). As wave propagation velocity through ground media is proportional to the square root of the tangential modulus [24], is given by the following equation:where is the tangential modulus at initial stress level at depth “z” and is the tangential modulus at peak stress level at depth “z.”
It should be noted that the ground is in equilibrium under geostatic stresses and at rest before nuclear explosion and only the overstresses caused by nuclearairblast causes the ground displacement.
2.2. Determination of Strain Distribution in Ground
Strain distribution with depth can be determined using stress distribution and an appropriate stressstrain relationship [e.g., 25–27] as given in the following equation:where functions and denote the loading and unloading branches of the stressstrain curve, respectively (Figure 2).
2.3. Integration of Strains
To obtain the vertical ground displacement at time instant , the strain distribution is integrated (6) from the ground surface to the penetrated wavelength in the ground up to time ,
Substituting in (6) leads to the following equation:
Thus, the displacement timehistory can be estimated using (7).
3. ClosedForm Solution
Using (7), closedform solutions can be obtained for several simplified cases [28]. In this article, a closedform approximation is developed for the following simplifications:(a)Linearly decaying overpressure timehistory with zero risetime (Figure 3, (8)),where and are peak overpressure and equivalent positive phase duration, respectively. The reason behind choosing the special case with zero risetime is the popularity of approximating the actual decay of the incidental pressure by an equivalent triangular pressure pulse among practicing engineers [8, 22]. It is to mention that design charts and empirical relations between the weapon yield, peak overpressure, and the equivalent positive phase durations are also available for equivalent triangular pulses with zero risetime [5, 8, 29]. Thus, for practicing engineers, a linear decay model with zero risetime is more useful.(b)Bilinear hysteretic stressstrain model [8] with the loading secant modulus (M_{L}) and strain recovery ratio (r) as parameters (Figure 2).(c)Attenuation coefficient is given by the following equation [8, 30]:where W is the yield of the explosion in kiloton and V_{L} is the wave propagation velocity of the peak overpressure front.(d)Using (4), is determined as 1 for the bilinear stressstrain model. However, due to nonlinear stressstrain behaviour, is usually greater than 1. Wilson and Sibley [17] and Batdorf [23] recommended a range of 1.5 to 2. Hereafter, is denoted as and assumed to be constant with depth.(e)A constant but representative Pwave velocity of ground media is assumed.
Using assumptions (a)–(e), integral for loading fronts in (7) (with ) between general time instants and is written as , where can be shown to be given by the following equation:
It is to clarify that the overpressure timehistory defined by the piecewise function such that for and for converges to (8) in the limiting case when . Therefore, the expression in (10) is first obtained for a case of finite and then it is evaluated for the case when to determine the solution for overpressure timehistory of (8). Similarly, the integral for unloading fronts in (7) (with ) between general time instants and can be written as , where can be shown to be given by the following equation:
By setting the appropriate values of and in (10) and (11), the closedform integrals between desired time instants can easily be evaluated. Thus, the closedform solution is given by the following equation:
The closedform solution is further extended to accommodate multiple ground layers. This requires the determination of time instants and during which the wave fronts pass through a particular layer. For illustration, a layered ground medium with an interface at depth with the modulus of top layer as and modulus of bottom layer as is considered as shown in Figure 3. This leads to the three cases (Table 1) and corresponding displacement solutions (shown in Table 1).

4. Validation
The proposed closedform solution is validated against the field data from an atmospheric nuclear test conducted at Frenchman Flat (Nevada). Perret [21] presented the measured overpressure timehistories (Figure 4) along with corresponding peak ground displacements at different distances from GZ (Table 2) for a 37 kt nuclear explosion at a height of 214 m. Model input parameters for above mentioned nuclear test are determined as follows.
(a)
(b)
(c)
(d)

The recorded overpressure timehistories are converted to equivalent linearly decaying overpressure timehistories (Figure 4). Equivalent overpressure timehistories are given by (8) such that the peak overpressure is equal to the recorded peak overpressure, and equivalent positive phase duration is given by the following equation:where is the positive phase impulse (area under the recorded overpressure timehistory). It is worth mentioning, though the risetime is taken to be zero in the equivalent overpressure timehistory in (8), that the effect of nonzero rising time has been accounted for by considering the total positive phase impulse (in (13)) which includes the area under the rising and decaying portions of actual overpressure timehistory. Therefore, it is expected that setting the risetime as zero in (8) would have a negligible impact on the magnitude of peak displacement. Furthermore, the ratios of length of risetime to total positive phase duration for the four cases (P1, P2, P3, and P4) are 0.10, 0.21, 0.23, and 0.27, respectively. This indicates that risetime is reasonably small near ground zero; however, the relative length of risetime increases with increasing distance from ground zero. Thus, it is expected that the errors (if any) due to setting risetime to zero would be more pronounced only at distances far away from ground zero. Even in such cases, a closedform solution can be developed using (7) by setting a finite risetime. However, it is worth noting that at far away distances (from ground zero), other important factors related to outrunning ground motion would start governing [31].
The initial pulse moving at seismic Pwave velocity penetrates to a depth of for equivalent case, whereas for actual case, it penetrates to a depth of . It is noted that the depth of the stressed zone in the equivalent case is smaller compared to the actual case as equivalent duration is smaller than the actual positive phase duration (Figure 4). Thus, all depthdependent parameters are scaled by a factor , such that the total depth of the stressed zone becomes with the depth of the top layer being , and with a characteristic attenuation length of . Based on the variation of Pwave velocity with depth at Frenchman Flat (Figure 5), the average Pwave velocity is estimated as 658.69 m/s. Whitman [4] presented the representative constrained modulus for Frenchman Flat interpreted from different experimental techniques (Table 3). Whitman [4] also emphasized that all aspects of stressstrain behaviour of geomaterials under blast loading are not captured by a single test and advocated to choose the modulus judiciously. The stress attenuation coefficient is calculated using (9) with scaled (Table 4). Using the closedform solution (shown in Table 1), parametric variations are studied (Table 5) and peak displacements (Table 6) are estimated.
(a)
(b)
(a)
(b)
(c)
(d)




The parametric studies (Table 5) highlight that the proposed model is most sensitive to the constrained modulus compared to other parameters and choice of the modulus is also subjective for engineers. However, based on the parametric studies, the last column “Recommended Value” of Table 5 provides some guidelines to select the constrained modulus. With increasing availability of similar case studies, these guidelines can further be refined.
5. Comparison with UFC Model
The UFC manual [22] provides an expression for airblastinduced peak vertical displacement based on onedimensional elastic wave propagation: where is the bulk density of geomaterials. To estimate peak displacements using (14), positive phase impulse is taken from Table 2. Representative Pwave velocity and density are taken as 658.69 m/s and 1331 kg/m^{3}, respectively, taking into account the variation with depth (Figure 5). Computed peak displacements using the proposed model and UFC model are shown in Figure 9 along with the measured field values. It can be clearly seen that the predicted values of the proposed model are in good agreement with the measured values and the UFC model significantly underestimates the peak displacements under high overpressures (with increasing war head capacity).
6. Conclusions
A closedform expression is developed to estimate nuclearairblastinduced freefield ground displacement that takes into account peak overpressure, positive phase impulse, depth of layer interface, representative constrained modulus of each layer, strain recovery, stress attenuation, Pwave velocity, and velocity ratio. The solution is validated against a nuclear test conducted at Frenchman Flat (Nevada) and the following conclusions are arrived at:(i)Peak ground displacement estimates are quite sensitive to the constrained modulus, and a judicious selection of appropriate constrained modulus based on the magnitude of applied overpressure is recommended.(ii)Based on the presented case study, two guidelines are recommended to select the appropriate modulus: (a) shallow ground layers are likely to have modulus values determined by unconfined or triaxial compression test and deep ground layers are likely to have modulus values computed from seismic velocity test, and (b) constrained modulus increases with decreasing ratio of applied stress to overburden.(iii)Effect of attenuation should be accounted for low overpressures and may be neglected at higher overpressures (or high war head capacities).(iv)The velocity ratio affects mainly the time of occurrence of peak ground displacement and not the magnitude under higher overpressures. Therefore, the velocity ratio becomes important when design calculations utilize the complete displacement timehistory such as in case of shock spectra.(v)A complete strain recovery is a better representation of actual conditions under low overpressure zones, and under higher overpressures, a partial strain recovery is recommended.(vi)The proposed model closely estimates the experimental values at all overpressures, whereas the UFC model significantly underestimates the peak ground displacements at high overpressures.
Notations
:  Attenuation coefficient at depth 
:  Strain at time at depth “z” 
:  Peak strain level in geomaterial at depth “z” 
:  Residual strain level in geomaterial at depth “z” 
:  Representative velocity ratio (between Pwave velocity and peak stress velocity) of ground media 
:  Ratio of Pwave velocity to wave velocity corresponding to the kth overpressure front at depth “z” 
:  Functional form for the loading branch of the stressstrain curve 
:  Ratio of Pwave velocity to wave velocity corresponding to the peak overpressure front at depth “z” 
:  Functional form for the unloading branch of the stressstrain curve 
:  Depth of the top ground layer 
:  Positive phase impulse of overpressure timehistory 
:  Characteristic attenuation length 
:  Loading secant modulus 
:  Modulus of bottom layer 
:  Modulus of top layer 
:  Peak overpressure 
:  Pressure magnitude of the kth overpressure front 
:  Strain recovery ratio 
:  Bulkdensity of geomaterials 
:  Scale factor 
:  Peak stress level in geomaterial at depth “z” 
:  Equivalent positive phase duration of linearly decaying overpressure timehistory 
:  Time when the kth overpressure front arrives at the point of interest on ground surface 
:  Time of arrival of the last loading front to reach at depth H (by the time ) 
:  Positive phase duration of overpressure timehistory 
:  Risetime to peak overpressure in overpressure timehistory 
:  Vertical ground displacement at time 
:  Representative Pwave velocity of ground media 
:  Representative wave propagation velocity of the peak overpressure front in ground 
:  Wave velocity of the kth stress front at depth “z” 
:  Yield of the explosion 
:  Depth penetrated by the kth overpressure front at time 
:  Stress at time at depth “z” 
:  Depth penetrated by initial pulse at time 
:  Depth penetrated by the peak overpressure front at time 
:  Tangential modulus at initial stress level at depth “z” 
:  Tangential modulus at peak stress level at depth “z.” 
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
References
 S. Pinto, “Underground construction of nuclear power reactors,” Nuclear Engineering and Design, vol. 61, no. 3, pp. 441–458, 1980. View at: Publisher Site  Google Scholar
 C. V. Chester and G. P. Zimmerman, “Civil defense shelters: a stateoftheart assessment,” Tunnelling and Underground Space Technology, vol. 2, no. 4, pp. 401–428, 1987. View at: Publisher Site  Google Scholar
 S. Glasstone and P. J. Dolan, The Effects of Nuclear Weapons, Castle House Publications, London, UK, 1980.
 R. V. Whitman, The Response of Soils to Dynamic Loadings; Report 26, Final Report, Massachusetts Institute of Technology, Cambridge, MA, USA, 1970.
 R. E. Crawford, C. J. Higgins, and E. H. Bultmann, The Air Force Manual for Design and Analysis of Hardened Structures, Final Report Jun 71May 74 (No. ADB004152/5), Civil Nuclear Systems Corporation, Albuquerque, NM, USA, 1974.
 B. M. Das, Fundamentals of SoilDynamics, Elsevier, Amsterdam, Netherlands, 1983, Chapter 4.
 H. L. Brode, A Review of Nuclear Explosion Phenomena Pertinent to Protective Construction, R 425 PR, The Rand Corporation, Santa Monica, CA, USA, 1964.
 M. S. Agbabian, Design of Structures to Resist Nuclear Weapon Effects, ASCE Manual on Engineering Practice, vol. 42, ASCE, Reston, VA, USA, 1985.
 J. H. Prevost, DYNAFLOW: A Nonlinear Transient Finite Element Analysis Program, Princeton University, Princeton, NJ, USA, 1981.
 F. Rischbieter and K. Michelmann, Bodenwelle bei Gasexplosionen. Beitrag der induzierten Bodenwelle zur Belastung von KKWStrukturen, BatelleInstitut, Frankfurt, Germany, 1984.
 K. J. Kim and S. E. Blouin, Response of Saturated Porous Nonlinear Materials to Dynamic Loadings, Applied Research Associates Inc., South Royalton, VT, USA, 1984.
 H. Werkle and G. Waas, “Analysis of ground motion caused by propagating air pressure waves,” Soil Dynamics and Earthquake Engineering, vol. 6, no. 4, pp. 194–202, 1987. View at: Publisher Site  Google Scholar
 Z. Yang, “Finite element simulation of response of buried shelters to blast loadings,” Finite Elements in Analysis and Design, vol. 24, no. 3, pp. 113–132, 1997. View at: Publisher Site  Google Scholar
 Z. Wang, H. Hao, and Y. Lu, “A threephase soil model for simulating stress wave propagation due to blast loading,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 28, no. 1, pp. 33–56, 2004. View at: Publisher Site  Google Scholar
 Y. Lu and Z. Wang, “Characterization of structural effects from aboveground explosion using coupled numerical simulation,” Computers & Structures, vol. 84, no. 28, pp. 1729–1742, 2006. View at: Publisher Site  Google Scholar
 T. Bierer and C. Bode, “A semianalytical model in time domain for moving loads,” Soil Dynamics and Earthquake Engineering, vol. 27, no. 12, pp. 1073–1081, 2007. View at: Publisher Site  Google Scholar
 S. D. Wilson and E. A. Sibley, “Ground displacements from airblast loading,” Journal of the Soil Mechanics and Foundations Division, vol. 88, no. 6, pp. 1–32, 1962. View at: Google Scholar
 M. L. Baron, I. Nelson, and I. Sandler, Investigation of Air Induced Ground Shock Effect Resulting from Various Explosive Sources. Report 2, Influence of Constitutive Models on Ground Motion Predictions, Weidlinger Associates, New York, NY, USA, 1971.
 Z. Wang and Y. Lu, “Numerical analysis on dynamic deformation mechanism of soils under blast loading,” Soil Dynamics and Earthquake Engineering, vol. 23, no. 8, pp. 705–714, 2003. View at: Publisher Site  Google Scholar
 S. Wang, L. Shen, F. Maggi, A. ElZein, and G. D. Nguyen, “Uniaxial compressive behavior of partially saturated granular media under high strain rates,” International Journal of Impact Engineering, vol. 102, pp. 156–168, 2017. View at: Publisher Site  Google Scholar
 W. R. Perret, Ground Motion Studies at High Incident Overpressure, Sandia Corporation, Albuquerque, NM, USA, 1960.
 UFC Manual, Structures to Resist the Effects of Accidental Explosions, US DoD, Washington, DC, USA, 2008, UFC 334002.
 S. B. Batdorf, “Airinduced ground shock: a onedimensional theory,” in Proceedings of International Symposium on Wave Propagation and Dynamic Properties of Earth Materials, Albuquerque, NM, USA, August 1967, G. Triandafilidis, Ed., pp. 423–432, Aerospace Corp., San Bernardino, CA, USA, 1969. View at: Google Scholar
 S. L. Kramer, Geotechnical Earthquake Engineering, Pearson Education, Delhi, India, 1996.
 X. Tong and C. Tuan, “Viscoplastic cap model for soils under high strain rate loading,” Journal of Geotechnical and Geoenvironmental Engineering, vol. 2, no. 206, pp. 206–214, 2007. View at: Publisher Site  Google Scholar
 J. An, C. Y. Tuan, B. A. Cheeseman, and G. A. Gazonas, “Simulation of soil behavior under blast loading,” International Journal of Geomechanics, vol. 11, no. 4, pp. 323–334, 2011. View at: Publisher Site  Google Scholar
 P. Y. Hicher and J. F. Shao, Constitutive Modeling of Soils and Rocks, John Wiley & Sons, Hoboken, NJ, USA, 2013.
 S. Pathak and G. V. Ramana, “Airblast induced ground displacement,” Procedia Engineering, vol. 173, pp. 555–562, 2017. View at: Publisher Site  Google Scholar
 N. M. Newmark and J. D. Haltiwanger, Air Force Design Manual, Principles and Practices for Design of Hardened Structures, University of Illinois at Urbana–Champaign, Champaign, IL, USA, 1962.
 A. J. Hendron Jr. and H. E. Auld, “Effect of soil properties on the attenuation of airblastinduced ground motions,” in Proceedings of International Symposium on Wave Propagation and Dynamic Properties of Earth Materials, Albuquerque, NM, USA, August 1967, G. Triandafilidis, Ed., pp. 29–47, University of Illinois, Urbana, IL, USA, 1969. View at: Google Scholar
 I. Nelson and G. Y. Baladi, “Outrunning ground shock computed with different models,” Journal of the Engineering Mechanics Division, vol. 103, no. 3, pp. 377–393, 1977. View at: Google Scholar
Copyright
Copyright © 2018 Shashank Pathak and G. V. Ramana. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.