Constraints on near-source saturation models for avoiding over-saturation of response spectral ordinates in RVT-based stochastic ground-motion simulations

Inversions of empirical data and ground-motion models to find Fourier spectral parameters can result in parameter combinations that produce over-saturation of short-period response spectral ordinates. While some evidence for over-saturation in empirical data exists, most ground-motion modellers do not permit this scaling within their models. Host-to-target adjustments that are made to published ground-motion models for use in site-specific seismic hazard analyses frequently require the identification of an equivalent set of Fourier spectral parameters. In this context, when inverting response spectral models that do not exhibit over-saturation effects, it is desirable to impose constraints upon the Fourier parameters to match the scaling of the host-region model. The key parameters that determine whether over-saturation arises are the geometric spreading rate (γ) and the exponential rate within near-source saturation models (hβ). The article presents the derivation of simple nonlinear constraints that can be imposed to prevent over-saturation when undertaking Fourier spectral inversions.


Introduction
Non-ergodic hazard analyses typically require adjustments to be made to published ergodic ground-motion models in order to make them more appropriate for the target region and site location (Bommer and Stafford 2020). The necessary adjustments are typically performed using a parametric approach within the general framework of the Hybrid Empirical Method (Campbell 2003), although empirical approaches also exist (Atkinson 2008). The parametric approach involves the identification of two sets of parameters that describe the far-field Fourier amplitude spectra (FAS) of ground motions in the host and target regions. A comprehensive overview of these FAS parameters and the far-field model is provided by Boore (2003). For the target region, Fourier spectral parameters are inferred through inversions of empirical data. This empirical data is almost always dominated by relatively small-magnitude earthquake scenarios, as site-specific or region-specific instrumentation typically covers a relatively short observation period. In contrast, for the host region, the Fourier spectral parameters are either adopted from similar empirical inversions for a region assumed equivalent to that represented by the ergodic ground-motion model(s), or the parameters are obtained from direct inversion of the ergodic ground-motion model predictions. As discussed by Stafford et al. (2021), the internal consistency of obtaining FAS parameters that best represent the ergodic ground-motion model provides distinct advantages over the empirical route where one can only assume consistency. This is particularly the case where there is no unambiguous geographical "host region" associated with the ground-motion model.
Host-region models are selected using a number of criteria (Campbell 2003;Bommer and Stafford 2020), but a fundamental assumption is that these models implicitly capture the appropriate scaling of largemagnitude, near-source scenarios that are very rarely constrained by empirical data. The inversion of such models therefore requires the consideration of largemagnitude scenarios where the point-source model of the FAS does not strictly apply. However, Boore (2009) has demonstrated that the stochastic method implemented using point-source models can lead to very similar results to those from extended or finite source models-provided that an appropriate distance metric is adopted (this distance metric is referred to as the "equivalent point-source distance," and is defined shortly). The common basis for Fourier spectral inversions is to therefore adopt a model for the FAS based upon far-field point-source considerations, coupled with an equivalent point-source distance metric that approximates the finite-fault effects that are increasingly relevant as the magnitude increases. A key effect in this context is the saturation of groundmotion amplitudes with increasing magnitude. Saturation is the term used to describe the reduced slope in magnitude scaling of ground-motion models at large magnitudes. Over-saturation, of concern in the present article, is the condition where the derivative of ln Sa with respect to magnitude becomes negative at larger magnitudes.
The equivalent point-source distance, R P S (M), can be defined as a function of the rupture distance, R RU P , and the magnitude-dependent near-source saturation length, h(M): where the exponent n is most commonly set to n = 2 (Boore and Thompson 2015). However, the distance scaling in a model such as Chiou and Youngs (2014) would correspond to a case where n = 1 . The near-source saturation length is used to reflect the fact that the strongest contributions to the observed ground-motion intensity measure will not necessarily come from the point on the rupture that is closest to the site.
In recent years, considerable attention has been paid to accurately constrain models for the near-source saturation length, h(M), used within the computation of equivalent point-source distance metrics. In the literature, these saturation lengths are sometimes independent of magnitude, and are also referred to as finite-fault factors (Boore and Thompson 2015), effective depths (Atkinson et al. 2016), pseudo-depths (Yenier and Atkinson 2015), or fictitious depths (Scasserra et al. 2009). As h(M) has its most significant impact upon near-source scaling where empirical data is very limited, it is very challenging to appropriately constrain these terms. However, they can have a significant influence upon other FAS parameters. In particular, Boore (2012) has shown that estimates of the stress parameter, Δσ , depend very strongly upon the modelling of near-source saturation effects. Furthermore, as will be demonstrated in the present article, the scaling of h(M) plays a strong role in determining whether or not over-saturation of short-period response spectral ordinates arises.
Ground-motion models have included near-source saturation lengths for decades, with the values being inferred as part of standard regression analyses on empirical data. However, recent years have seen greater attention being given to the scaling of ground motions over a broader magnitude range, particularly extrapolating models down to smaller events. If inappropriate saturation lengths are adopted for these small events, the implied stress parameters associated with the events can be heavily influenced. Yenier and Atkinson (2014) introduced a model for h(M) that worked well for relatively large events, while Yenier and Atkinson (2015) proposed an alternative model that works better for smaller magnitude events. Boore and Thompson (2015) then developed a new model that merged the models from Yenier and Atkinson (2014) and Yenier and Atkinson (2015) to ensure strong performance over an extended magnitude range. Atkinson et al. (2016) focussed upon these finite-fault metrics further, with a focus on their applicability for relatively shallow, small magnitude, induced events (as source depth can play a role similar to the near-source saturation length). They concluded that the model from Yenier and Atkinson (2015) was suitable for application to induced events.
In most cases, near-source saturation models are developed using empirical observations over a magnitude range that is smaller than that relevant in common practical applications. As noted above, near-source saturation models have recently received heightened attention in the context of small-magnitude events (Atkinson et al. 2016). However, when these models are extended to very large magnitudes, their interaction with other components of the stochastic method can lead to over-saturation of earthquake ground motions, where predicted amplitudes for very large events are lower than those from smaller magnitude events.
The purpose of this short article is to define the conditions for which over-saturation occurs so that model developers can impose these constraints if they deem over-saturation to be an undesirable feature. This is primarily an issue for short-period response spectral ordinates. The derived conditions are useful for analysts inverting either empirical datasets, or ground-motion models, to identify underlying Fourier spectral parameters. The primary practical application where this work is relevant is when making hostto-target adjustments in the context of site-specific ground-motion modelling. Specifically, when using the Hybrid Empirical Method of Campbell (2003), one needs to invert empirical data in the target region and a ground-motion model from the host region in order to derive the final ground-motion model for the target region.

Should over-saturation be prevented?
The purpose of this article is to define constraints explicitly to avoid the occurrence of over-saturation for large-magnitude, short-period, response spectral predictions. However, some empirical data provides evidence of over-saturation , and a valid question is whether one should actively attempt to prevent over-saturation. This issue is primarily of concern for ground-motion model developers, who must decide whether to let sparse data drive the scaling of their models, or impose some constraints to control this scaling (which may allow or prevent oversaturation). For the present article, this issue is of tangential interest, as the objective when performing FAS inversions of a ground-motion model is to identify the parameter combinations that collectively describe the scaling embedded within the model. While ensuring that these parameters are consistent with physical expectations is desirable, it is not the primary objective. Therefore, if a ground-motion model does not allow over-saturation to occur, then the FAS parameter sets should be obtained to respect this condition.
Apparent over-saturation can arise for a number of reasons, such as imperfect calibration of nonlinear site response functions, or even legitimate effects of nonlinear site response. However, the primary reason why over-saturation can arise is due to the way in which source-to-site distances are defined in groundmotion models. It is now standard to use finite-fault distance metrics such as the Joyner-Boore distance, R J B , or the rupture distance, R RU P , within groundmotion models (Abrahamson and Shedlock 1997). The rupture distance is also recommended for use within the stochastic ground-motion method (Boore 2003), when combined with an appropriate finite-fault factor (Boore and Thompson 2015), as seen earlier in Eq. 1.
However, for the very large ruptures associated with large-magnitude events, ground motions may be driven by asperities that are located some distance from the point implied by these distance metrics. Finite-fault factors are defined precisely for this purpose, to increase the effective rupture distance as the magnitude increases, but the effectiveness of these factors when applied to empirical data depends strongly upon the particular spatial distribution of recording stations around the rupture. Alternative distance metrics, such as the mean rupture distance proposed by (Thompson and Baltay 2018), aim to remove the need for these finite-fault factors, but come with additional issues, and also lead to over-saturation (depending upon the adopted power p used in the mean distance definition).

Mathematical framework
Within the RVT framework (Boore 2003), response spectral ordinates for a given natural period, T n , are computed as: where Sa(f n , ζ n ) is the spectral acceleration of an oscillator with natural frequency f n = 1/T n and damping ratio ζ n , ψ is the peak factor (the ratio of the peak response to the root-mean-square, RMS, response), m 0 is the zeroth spectral moment of the oscillator response, and D rms is the root-mean-square duration. The zeroth moment is computed from using k = 0 in the general expression: (3) When deriving ground-motion models, it is most common to work with the logarithmic spectral ordinates, and we can write: ln The RMS duration, D rms , is computed from the excitation duration, D ex , as: where Γ is a function of the excitation duration, D ex , the natural period of the oscillator, T n , and the damping ratio, ζ n . The term Γ is defined and computed as the adjustment to the excitation duration that is required in order that RVT-based estimates of peak oscillator response match those obtained from timedomain simulations. Boore and Thompson (2012) provide a brief history of the development of these adjustment factors. For the current article, it is important to note that these factors are magnitude and distance dependent, but that the magnitude dependence is relatively mild. There are various methods for computing the peak factor, ψ. The most recent recommendation of Boore and Thompson (2015) is to use the peak factor approach of Der Kiureghian (1980) that builds upon the work of Vanmarcke (1975). This particular formulation is a function of the first three spectral moments, i.e., m 0 , m 1 , and m 2 , as well as the excitation duration. As such, the peak factor is also magnitude and distance dependent. However, like the Γ factor, the dependence upon magnitude is relatively mild.
Dropping the dependence upon f n and ζ n for brevity, but introducing explicit dependence upon magnitude M, we can write: The validity of the approximation in Eq. 7 can be appreciated from Fig. 1, where the partial derivatives Fig. 1 Derivatives of the logarithmic spectral ordinate, Sa(T n = 0.01, ζ n = 0.05), contributions in Eq. 6 with respect to magnitude. All derivatives are computed for R RU P = 1 km, Δσ = 10 MPa, perfect spherical geometric spreading (the geometric spreading rate is γ = 1), the anelastic attenuation of Q(f ) = 200f 0.5 , and κ 0 = 0.035 s. The recommendations of Boore and Thompson (2015) for the peak factor, finite-fault factor, and RMS and excitation durations are also adopted of each term in Eq. 6 are shown. The partial derivatives associated with ψ(M) and Γ (M) are far weaker than those for m 0 (M) and D ex (M), particularly at large magnitudes where over-saturation may arise.
If we wish to prevent over-saturation, we require that: and this condition is most relevant for large magnitudes and short periods, where over-saturation can occur.
Assuming that the magnitude dependence of ln Sa(M) is driven by the scaling within m 0 (M) and D ex (M), we can write: which is equivalent to: For the first term, we can write: The frequency response function of the oscillator is not a function of magnitude, so the derivative of the above expression is: We can also recognize that for short periods, f n is very high, and |H (f ; f n , ζ n )| ≈ 1 over the range of frequencies that dominate the summation (or integral-referring back to Eq. 3) (Bora et al. 2016). Therefore: Assuming a single corner-frequency ω 2 source spectrum, with corner-frequency dependent upon magnitude (Aki 1967;Brune 1970;1971), the Fourier amplitude spectrum of acceleration can be written as: where f c (M) is the magnitude-dependent corner frequency, M 0 (M) is the seismic moment, R 0 is the reference distance at which the source amplitude is defined, g(R; M) is the geometric spreading function, q(f, R; M) is the anelastic attenuation filter, and s(f ) is the combination of impedance and damping effects at the site. The term C is defined as: with R θφ = 0.55 being the average radiation pattern, V = 1/ √ 2 partitions energy into horizontal components, F = 2 is a free-surface factor, ρ s and β s are the density and shear-wave velocity at the source, and R 0 = 1 km is a reference distance (and is interpreted as a rupture distance, i.e., R 0 ≡ R RU P ,0 = 1 km) (Boore 2003).
In the present article, we are most concerned with the situation in which we are at very short distances, ideally at the reference distance where the source motions are defined. The corresponding point-source distance at the source is therefore: At these short distances, geometric spreading functions are also typically characterized as: where γ is, again, the geometric spreading rate. The anelastic attenuation filter is normally written as: with Q 0 being the quality factor, η being the quality exponent (which controls the frequency dependence of the damping along the propagation path), and c Q being the crustal velocity used to define Q 0 . Using the equivalent point-source distance, this is: The site effects are expressed as: with S I (f ) representing the impedance effects (generally site amplification, increasing with frequency), from a model such as Boore (2016), and S K (f ) representing damping effects. The overall expression for the FAS at the reference distance R 0 can then be expressed as: where it can be appreciated that magnitude dependence arises from the source model, and path components given that the equivalent point-source distance used within the geometric spreading and anelastic attenuation filters is magnitude dependent.

Partial derivative of the zeroth moment
Equation 13 showed how the derivative of the zeroth moment is directly related to the derivative of the FAS, |A(f, M; R 0 )|. We therefore require the specification of the partial derivative of |A(f, M; R 0 )| with respect to magnitude. As |A(f, M; R 0 )| is a product of four distinct components that have magnitude dependence, we make use of the chain rule to compute the overall derivative.
The seismic moment can be written as: where α = ln(10) 3 2 and β = ln(10) × 16.05. These constants are simply conversions of the original constants presented by Hanks and Kanamori (1979) to account for the change from the base-10 to natural logarithm (hence the ln(10) multiplicative factors). The derivative of the seismic moment can then be defined as: Following Brune (1970), the source cornerfrequency can be defined in terms of the seismic moment and stress parameter, Δσ , as: where F is a constant that embeds the effects of the rupture velocity, source geometry, and slip distribution. The corresponding partial derivative of the corner frequency is: For the equivalent point-source distance, we have the general expression: For the most common case, where n = 2, we have: and for large magnitude events where over-saturation is potentially of concern, we will have h(M) 1, so that the above expression reduces to: When n = 1, we will always have the equivalence shown in Eq. 28 regardless of how large h(M) is. So, for typical cases, we have the derivative of the distance being equivalent to the derivative of the saturation distance. As the saturation distance scales according to: over the large magnitude range (e.g., Chiou and Youngs (2014), Yenier and Atkinson (2014), Yenier and Atkinson (2015), and Boore and Thompson (2015)), we have: where the final approximation stems from the condition that h(M) 1 for the scenarios of interest. Finally, for the anelastic attenuation filter, we have: from which we can write: Due to the underlying exponential nature of all of the components of Eq. 21, we can express their derivatives in terms of a factor multiplied by the original component. This makes the application of the chain rule lead to a relatively compact expression for the derivative of a FAS ordinate, as shown in Eq. 33.
Now, recalling that: and that the moment expression (for very high f n , where |H (f ; f n , ζ n )| ≈ 1 over the frequencies that are of most importance) is: we can write the contribution to the derivative of ln Sa as: The final term involving the summations (integrals) provides a weighted average of the terms in square brackets, with the weighting proportional to the squared FAS. Within the square brackets, the first term involving the corner frequency tends to low values below f c , and tends to unity as f i f c . Therefore, for the frequencies where |A(f i , M)| 2 is strong, the first term is essentially a constant equal to 2 3 α. Equation 36 can therefore be simplified to: In this light, we can see that the second term provides a reduction based upon a weighting of the FAS. Replacing the elaborate summation terms by a function Ψ Q , for brevity, we have: The nature of the Ψ Q function will be discussed subsequently in Section 3.3. For now, however, it is worth stating that it provides a modest reduction to the main α/3 − γ h β term, and that the magnitude of Ψ Q depends upon a number of FAS parameters-as noted by the function arguments.

Partial derivative of the excitation duration
The other key contribution to the derivative of ln Sa is the partial derivative of the excitation duration. The excitation duration is composed of source and path contributions. A number of different models for the excitation duration have been proposed, but they generally take the form of a source contribution that is related to the source rise time (the reciprocal of the corner frequency), and piecewise linear path scaling components. While the particular linear segments of the path scaling can change over different distances depending upon the region (e.g., Boore and Thompson (2014)), in the near-source region a single segment with slope d β can be defined. Such a model is shown in Eq. 39.
Equation 39 includes a constant term d α that may be zero, but otherwise represents any linear segment of the path duration model defined prior to that relevant for distance R P S (M).
The partial derivative of the excitation duration can then be defined as: where results from the previous section have been used.
The contribution to the derivative of ln Sa can then be defined as: or, after factoring out the α/3 term, and substituting D src (M) = 1/f c (M). where The behavior of Ψ D is explored in Section 3.3. It is primarily dependent upon the properties of the excitation duration model, but it is also coupled to the FAS parameter Δσ given that source duration is directly linked to this parameter.

Partial derivative of the log spectral acceleration
Returning to Eq. 9, we can now state that the partial derivative of ln Sa with respect to magnitude is defined by the approximate expression: As neither Ψ D nor Ψ Q depends upon the geometric spreading rate γ , we can specify a condition upon this rate such that ∂ ln Sa(M)/∂M ≥ 0. This condition is: That is, provided that the geometric spreading rate γ is less than or equal to the expression on the righthand side of Eq. 45, over-saturation of short-period response spectral ordinates should not arise. Here, "short-period" refers to cases where the transfer function of the SDOF oscillator is essentially flat and equal to unity over the frequency range that dominates the computation of m 0 . For the purposes of performing an inversion of a ground-motion model, or of empirical data, the condition in Eq. 45 is not particularly convenient as it involves all of the FAS parameters. To determine whether simplifications can be made to this expression, Figs. 2 and 3 show the typical magnitudes of the Ψ D and Ψ Q terms, respectively, as well as how these terms vary with the FAS parameters.
In Fig. 2, we observe some sensitivity to both the stress parameter and the coefficients of the excitation duration model (as evidenced by the different values of Ψ D for active or stable crustal regions). In theory, as M → ∞, Ψ D → 1 2 , regardless of the path duration model. However, over the magnitudes of practical interest shown in Fig. 2, we can see that Ψ Q ≈ 0.4 − 0.45 for active crustal regions, and Ψ Q ≈ 0.3 − 0.4 for stable crustal regions. Making the simplification that Ψ D = 1 2 is therefore a conservative approximation that limits the value of γ below the levels strictly implied by Eq. 45. The performance of this approximation is assessed shortly.
For the scaling of Ψ Q , Fig. 3 shows that magnitude dependence is relatively weak, and that large variations in the FAS parameters generally do not lead to major variations from a representative value of about Ψ Q ≈ 0.1. Whereas making an assumption that Ψ D = 1 2 leads to conservative constraints upon γ , neglecting the Ψ Q contribution by assuming Ψ Q ≡ 0 is unconservative. However, it is also important to reconsider Fig. 1 in which the various contributions to the overall derivative of ln Sa were shown. From consideration of that figure, it was decided to ignore the influence of the peak factor, ∂ ln ψ(M)/∂M, given  that its magnitude was significantly lower than the contributions from ln m 0 (M) and ln D ex (M). Closer inspection of Fig. 1 shows that the peak factor contribution is ∂ ln ψ(M)/∂M ≈ 0.1 over the range of magnitudes shown. The reason why this is important is that ∂ ln ψ(M)/∂M and Ψ Q have similar magnitudes and opposite signs. As such, neglecting both allows the overall derivative of ln Sa to still be reasonably well approximated, but also permits Eq. 45 to be significantly simplified.
Making the assumption that Ψ D ≡ 1 2 and that Ψ Q ≡ 0 in Eq. 45 therefore allows the following condition to be established.
This very simple condition has the significant advantage of only involving two key FAS parameters.
Furthermore, despite its simplicity, it has proven effective in ensuring that over-saturation does not arise when performing FAS inversions of ground-motion models .

Performance of the derived constraints
In order to demonstrate the performance of the elaborate and simplified constraints from Eqs. 45 and 46, it is necessary to work with a near-source saturation model that allows h β to vary in a sensible manner. Stafford et al. (2021) presented a near-source saturation model with the functional form shown in Eq. 47.
As described by Stafford et al. (2021), the functional form of this expression is adopted from the magnitude-scaling parameterization of the Chiou and Youngs (2014) ground-motion model. The function is effectively composed of two linear portions (in ln h(M)-M space), connected by a smooth transition. It has many similarities with the model of Boore and Thompson (2015) in this respect. At small magnitudes, the slope is equal to h γ , while at large magnitudes the slope is equal to h β . As these linear portions relate to the scaling of ln h(M), the slopes can be regarded as exponential saturation rates. The transition between the two linear portions occurs around a magnitude centered at h ε , and h δ controls how quickly this transition takes place. When varying the value of h β , the overall function in Eq. 47 will shift. Therefore, to preserve the scaling implied for smaller magnitudes an adjustment to h α needs to be made. If we define Δh β as a shift from some nominal value h β,0 , such that h β = h β,0 + Δh β , Eq. 47 can be expressed as: where h α = h α − Δh β h ε . This approach has been adopted to generate the alternative models shown in Fig. 4. The value of h β = 0.5 has been selected to be similar to both the Yenier and Atkinson (2015) and Fig. 4 Comparison of near-source saturation models. The generic saturation model of Eq. 47 is implemented here with h α = −0.9, h γ = 1.15, h δ = 2.5, and h ε = 6.5, to approximately mimic the scaling of Boore and Thompson (2015) Boore and Thompson (2015) models at large magnitudes (to match exactly, h β = 0.541), while the other parameters noted in the figure caption have been selected (from visual comparison) to provide similar scaling to Boore and Thompson (2015). A large number of FAS parameter combinations are then considered, and the exact partial derivative of ln Sa is computed for different large-magnitude values. Specifically, representative values of Δσ = 10 MPa, Q 0 = 200, η = 0.5, and κ 0 = 0.035 s are chosen, while γ and h β are allowed to vary. The nominal FAS parameters of {Δσ, Q 0 , η, κ 0 } are comparable (but not identical) to values obtained from inversion of the Chiou and Youngs (2014) ground-motion model .
The exact derivatives are computed using automatic differentiation (Molkenthin et al. 2014). This is achieved through the Julia (Bezanson et al. 2017) programming language, and specifically using the package ForwardDiff (Revels et al. 2016) in conjunction with the bespoke RVT-based package StochasticGroundMotionSimulation (Stafford 2021). Figure 5 shows the results for four magnitude values and broad ranges of h β and γ . The exact boundary where ∂ ln Sa(M)/∂M = 0 is shown with a labelled contour. The region to the upper right of this contour corresponds to parameter combinations where oversaturation will occur. In addition to Eqs. 45 and 46, two additional cases are shown in Fig. 5. The dashed line shows the case where Ψ D is considered, but Ψ Q is neglected. The dotted line shows the case where Ψ D = 1 2 and Ψ Q is considered. These two options for simplifying the complete expression in Eq. 45 bound the more appropriate constraints given in Eqs. 45 and 46, with the dashed lines being unconservative, and the dotted lines being conservative.
In principle, the most accurate constraint should be obtained using Eq. 45. However, Fig. 5 shows that this is not always the case, particularly at the lower magnitudes considered. Recall that a number of approximations were made throughout the derivation and that these approximations should improve as the magnitude increases. This is seen in Fig. 5 where for the largest magnitudes in the lower panels the expression of Eq. 45 is essentially parallel to the exact zero contour, and is slightly on the conservative side. However, what is more pleasing is that the very simple Comparison of the exact derivatives for ranges of h β and γ with the approximate constraints presented in the text. The nominal FAS parameters that are not varied are Δσ = 10 MPa, Q 0 = 200, η = 0.5, and κ 0 = 0.035 s constraint of Eq. 46 does a very good (arguably the best) job at M7.5 and also performs reasonably well for most γ and h β combinations at larger magnitudes. For low geometric spreading rates of γ ≤ 1.0, the simple constraint becomes slightly unconservative, but otherwise it provides an excellent approximation to the true boundary.
Because the constraint of Eq. 46 only involves γ and h β , it is very simple to specify this condition as a nonlinear constraint within an inversion. For this reason, it is recommended that this simple constraint be favored over the more elaborate expression of Eq. 45. Note that in the recent inversion of the Chiou and Youngs (2014) ground-motion model by Stafford et al. (2021), this simple constraint was adopted, and served its purpose well.

Conclusions
When inverting empirical ground-motion data, or ground-motion models, it is often desirable to prevent the over-saturation of short-period response spectral ordinates. Depending upon how the rupture scenarios within the empirical data are distributed, or how the parameter space is sampled within an inversion of a ground-motion model, over-saturation can arise even when the underlying data or model does not support this. To prevent over-saturation, the optimization process used to identify the FAS parameters can be conducted with the use of nonlinear constraints among the parameters. This article has presented two such constraints. The elaborate expression of Eq. 45 provides a good, slightly conservative, approximation to the limiting condition of ∂ ln Sa(M)/∂M = 0 for the largest magnitudes of typical interest. However, this expression requires the specification of very complex nonlinear constraints within any inversion procedure. On the other hand, the extremely simple result of γ h β ≤ α/6 (see Eq. 46) performs extremely well for the vast majority of cases of practical interest. This simplified constraint places a direct limiting relationship upon the near-field geometric spreading rate, γ , and the large-magnitude saturation rate, h β , of near-source saturation models. Therefore, for analysts making general forward simulations under the stochastic method, respecting the constraints derived in the present article will ensure that over-saturation of short-period response spectral amplitudes will be avoided.