The empirical probability of submarine mass failure is quantified from a sequence of dated mass-transport deposits. Several different techniques are described to estimate the parameters for a suite of candidate probability models. The techniques, previously developed for analyzing paleoseismic data, include maximum likelihood and Type II (Bayesian) maximum likelihood methods derived from renewal process theory and Monte Carlo methods. The estimated mean return time from these methods, unlike estimates from a simple arithmetic mean of the center age dates and standard likelihood methods, includes the effects of age-dating uncertainty and of open time intervals before the first and after the last event. The likelihood techniques are evaluated using Akaike’s Information Criterion (AIC) and Akaike’s Bayesian Information Criterion (ABIC) to select the optimal model. The techniques are applied to mass transport deposits recorded in two Integrated Ocean Drilling Program (IODP) drill sites located in the Ursa Basin, northern Gulf of Mexico. Dates of the deposits were constrained by regional bio- and magnetostratigraphy from a previous study. Results of the analysis indicate that submarine mass failures in this location occur primarily according to a Poisson process in which failures are independent and return times follow an exponential distribution. However, some of the model results suggest that submarine mass failures may occur quasiperiodically at one of the sites (U1324). The suite of techniques described in this study provides quantitative probability estimates of submarine mass failure occurrence, for any number of deposits and age uncertainty distributions.

Submarine mass failures present a significant hazard for submarine infrastructure, for offshore energy development, and to coastal regions via the generation of tsunamis (Locat and Lee, 2002; Masson et al., 2006; Sawyer et al., 2009; ten Brink et al., 2009b; Jackson, 2012). One of the critical components of assessing natural hazards is determining the probability of occurrence. Empirically estimating the probability of submarine mass failures is particularly difficult, owing to the paucity of age dates for individual events. Typical problems in dating submarine mass failures and their deposits include lack of overlying strata, the inability to core underlying strata, and lack of datable material. William R. Normark, to whom the “Exploring the Deep Sea and Beyond” themed issue is dedicated, was heavily involved in several studies dating submarine mass failures (Lipman et al., 1988; Normark, 1990; Laursen and Normark, 2002; Lee et al., 2004; Fisher et al., 2005; Normark and Gutmacher, 1988). In particular, for the Palos Verdes debris avalanche in southern California, Normark et al. (2004) presented an ideal case where radiocarbon ages from piston cores of distal strata above and below an acoustically transparent layer continuous with the debris avalanche were used to constrain the age of the failure. This and similar studies (e.g., Haflidason et al., 2005) provide the necessary data with which to determine the probability of submarine mass failures.

For an empirical determination of failure probability, there are very few places in the world where repeated submarine mass failures have been individually dated. Cores from two Integrated Ocean Drilling Program (IODP) sites in the Ursa Basin, northern Gulf of Mexico (Fig. 1), penetrated a sequence of mass transport deposits (MTDs) (Expedition 308 Scientists, 2006a, 2006b). These deposits are interpreted to have been emplaced during episodes of retrogressive submarine mass failures upslope from the drill sites (Sawyer et al., 2009). As many as 14 MTDs have been identified in both drill holes and have been assigned ages based on microfossils and magnetostratigraphy (Urgeles et al., 2007).

The objective of this study is to develop a methodology to estimate the probability of submarine mass failures from a sequence of dated MTDs, such as found in the Ursa Basin. A previous study (Geist and Parsons, 2010) focused on the more common situation in which a sequence of MTDs are identified in seismic reflection records, but only a basal horizon beneath the oldest event is dated (e.g., Fisher et al., 2005). This previous method assumes that submarine mass failures occur as a stationary Poisson process. A gamma distribution is used as a conjugate prior to the Poisson distribution in a Bayesian framework, in order to estimate the Poisson intensity or rate parameter and its uncertainty. In contrast, it is not assumed from the outset of this study that submarine mass failures occur according to a Poisson process; several other time-dependent probability models that include clustering and quasiperiodic behavior are also examined. Other Bayesian techniques that have been adapted for small data sets have been proposed by, for example, Nomura et al. (2011).

Paleoseismic records represent an analogous situation where a sequence of earthquake events has been dated. While similar in many respects, the stratigraphic record of paleoseismic and mass transport deposits differ in terms of the scale of their stratigraphic footprint and the effects of erosion and deposition processes on the retention of the event record. Individual paleoseismic events are generally marked by the presence of distinct colluvial wedges, angular unconformities, fissures, sand blows, upward termination of fault displacements, and abrupt lateral changes in bed thickness (McCalpin, 1996). Often, these stratigraphic indicators are present within close proximity to one another and overlap with similar features from previous seismic events, leading to a tightly constrained, aggraded event record. Mass transport deposits, on the other hand, are often transported significant distances from their source regions, and, depending on the properties of the sediment, water content, and flow energy, can undergo significant transformation prior to deposition (e.g., Normark and Gutmacher, 1988). During transport, these sediment flows can be extremely erosive, scouring underlying MTDs and intervening hemipelagic sediment. Limited stratigraphic data, the difficulty in quantifying the effect of this erosion, and identifying the subdued change from MTD to post-event sediments can significantly impact the fidelity of the MTD age record. Although the absolute uncertainty of MTD age dates is often greater than that of paleoseismic horizons, the uncertainty relative to the mean return time may be less for MTDs compared to paleoseismic ages along active margins.

Probabilistic methods applied to paleoseismic data have included estimating distribution parameters from a sequence of event times (Utsu, 1984; Nishenko and Buland, 1987) and more recent methods that take into account the uncertainty of age dating and open time intervals before the oldest event and after the last event up to the present (Rhoades et al., 1994; Ogata, 1999b; Parsons, 2008). We adapt several methods from the latter studies to estimate the probability of submarine mass failures. The size of the failure is not explicitly considered and it is assumed that the probability estimates apply to all sizes as expressed by the observed MTDs.

In many studies, the arithmetic (sample) mean using the center age of dated events is commonly cited as the mean return time for a natural hazard. This statistic does not take into account the open time intervals and the uncertainty in age dating. Moreover, uncertainty is commonly reported in terms of the standard deviation of the center ages, which is not appropriate for events that do not occur according to a Gaussian distribution (Parsons, 2008). Using the methods described in this paper, confidence intervals can be established for each of the parameters in the probability models to better characterize uncertainty in the hazard estimate.

For natural hazards where there has been a lack of probability models developed, such as for submarine mass failures, the least astonishing and simplest hypothesis is to assume that events occur as a Poisson process. For a Poisson process, the time between individual events (i.e., inter-event or return time) is independent of the past history of events and the return times are distributed according to an exponential distribution parameterized by a rate or intensity parameter. Moreover, a Poisson process is in general the one with the maximum entropy and an ideal reference model against which the information content in other probability models can be gauged (Kagan and Knopoff, 1987). Other models that are developed for natural hazards are most often evaluated against the exponential model using various statistical tests and measures of goodness of fit.

One non-Poissonian probability model that has been developed for subaerial mass failures is the Weibull model. The Weibull model is one of several models that includes an additional parameter quantifying whether a sequence of events is clustered or quasiperiodic in time. Other probability models that are commonly invoked are the lognormal and gamma models. In Griffiths’s (1993) study, the author used the Weibull model and parameter estimation using a method of moments to indicate that large rock avalanches are quasiperiodic. In contrast, Witt et al. (2010) provided detailed evidence for the clustering of small subaerial mass failures from historical records. One other model that has been recently employed for earthquakes is the Brownian Passage Time model (BPT) (Matthews et al., 2002), which has the characteristic of an asymptotically constant hazard function, as opposed to the lognormal (asymptotically decreasing hazard) and the Weibull model (asymptotically increasing hazard).

One may intuitively expect that submarine mass failures do not occur independently, owing to the fact that one event may destabilize the slope surrounding the failure, making another failure more likely to occur (e.g., Pratson and Coakley, 1996; Sawyer et al., 2009). This suggests, on the one hand, that submarine mass failures cluster in time. On the other hand, the long-term occurrence of large submarine landslides is likely dependent on sediment supply, particularly in seismically active regions, where the occurrence is dependent on the buildup of sufficient sediment. This line of reasoning would suggest that submarine mass failures are quasiperiodic in time (Griffiths, 1993). The relationship between variable sedimentation rate, pore pressures, and slope failures, however, is complex (Stigall and Dugan, 2010; You et al., 2012).

We test the aforementioned probability models against the data presented in the Ursa Basin IODP holes. Assignment of the range of dates associated with each failure is described in the Age Dates of Mass Transport Deposits section. Expressions for the probability models under consideration are described in the Probability Models section. In the Maximum Likelihood Methods and the Monte Carlo Method sections, we describe two general types of statistical inference methods that estimate the parameters of the probability models. In the Model Selection section, we compare each of the probability models, and we estimate the hazard function for each model in the Hazard Estimation section. Overall, this study provides a framework of inferring probability models for other areas in which a sequence of MTDs has been dated.

The approximate ages of nine MTDs from site U1322 (Table 1) and five MTDs from site U1324 (Table 2) are used to test the probability models. While the morphology of MTDs can change rapidly over very small distances leading to incomplete records in recovered stratigraphic sections, the MTD stratigraphic record at both sites is assumed complete for the areas local to the holes based on the interpretation of seismic reflection data from the region (Winker and Shipp, 2002; Expedition 308 Scientists, 2006b; Sawyer et al., 2007). Age ranges for each MTD were estimated from the age model of Urgeles et al. (2007), which was established from microfossil zonations and magnetostratigraphic datums tied to global δ18O records (Expedition 308 Scientists, 2006b). Events 1b and 1d in both sites occur in lithostratigraphic unit 1 (subunits 1b and 1d), while events 2a–2g in U1322 (Table 1) are thought to be time correlative with subunits IE through IIB at site U1324 (Expedition 308 Scientists, 2006a).

The MTD sequences at both sites overlie and erode into the regionally identifiable Blue Unit, a thick layer of interbedded sand and mud, likely emplaced by turbidity currents in a basin-floor fan setting during the Marine Isotope Stage (MIS) 4 sea-level lowstand (Sawyer et al., 2009). As the top of the Blue Unit has been substantially eroded, for the purposes of this study the base of the Blue Unit is assigned an age of 74 ka, the maximum age of the MIS 4 interval (Martinson et al., 1987). This date is used as the earliest time when no submarine mass failures had occurred prior to the first MTD (i.e., the “S” date in Ogata, 1999b and Maximum Likelihoods Section). Given the lack of absolute age information and therefore low fidelity/resolution of the age model, we apply no formal error estimates to the upper and lower bounds of MTD ages in Tables 1 and 2.

Models of submarine mass failure occurrence assume that the return times are distributed according to a particular probability distribution. Consider a point process consisting of a sequence of n events, in this case submarine mass failures, that are ordered in time: t1 < t2 <…< ti<…< tn. The time between events is termed the inter-event or return time τ and is given by τi = ti+1ti. The process that gives rise to these events is considered a renewal process if the return times are independent and identically distributed according to a particular probability distribution (Daley and Vere-Jones, 2003). A renewal process is considered “memoryless” in that the probability of a future event does not depend on the past history of events, only the elapse time since the last event. Moreover, a stationary Poisson process is one in which the probability of a future event is constant for a given time interval and depends on neither the past history of events nor the elapsed time since the last event. For a Poisson process, the number of events in a finite interval follow a Poisson distribution and the probability density function (pdf) of return times follows an exponential distribution:
where λ is the rate or intensity parameter. The mean return time for this distribution is given by 1/λ. A stationary or temporally homogeneous Poisson process is one in which λ is constant with time. All probability models considered in this study are assumed to be stationary.
Several other renewal-process probability distributions that account for increased clustering or quasiperiodic behavior can be compared with the exponential distribution. The pdf of return times described by a gamma distribution is given by
where Γ is the gamma function, β is a scale parameter, and α is a shape parameter, such that α < 1 is associated with a cluster process and α > 1 is associated with a quasiperiodic process (α = 1 recovers the exponential distribution). The mean return time for the gamma distribution is given by αβ.
The pdf for the Weibull distribution, used by Griffiths (1993) and Witt et al. (2010) as a model for subaerial landslides, is given by
where β is a scale parameter and α is a shape parameter (α < 1 and α > 1 represent cluster and quasiperiodic processes, respectively, as with the gamma distribution). The mean for this distribution is given by βΓ(1 + 1/α).
The pdf for the lognormal distribution that commonly characterizes quasiperiodic processes is given by
where μ is a location parameter and σ is a shape parameter. The mean return time for the lognormal distribution is given by exp(μ + σ2/2).
Finally, the BPT distribution, developed by Matthews et al. (2002) for earthquake sequences, has a pdf given by
where μ is a location parameter and α is the coefficient of variation or aperiodicity parameter. The mean is simply given by μ. This distribution is equivalent to the inverse Gaussian distribution or Wald distribution given by
where forumla. The BPT distribution is nearly indistinguishable from the lognormal distribution described above. Both the lognormal and BPT models have the property of zero probability of repeated failure immediately following an event. In addition, at long times since the last event relative to the mean return time, the hazard rate for the BPT distribution is quasi-stationary, as is the hazard rate for the gamma distribution.
Aside from the density distribution of return times, the parameter estimation methods below are based on functions that can be derived from the pdf. These functions include the survivor function ϕ(τ) = 1 – F(τ), where F(τ) is the cumulative distribution function, and the hazard-rate function υ(τ) given below:

Note that according to this definition and Equation 1, υ(τ) for a Poisson process simply equals the intensity parameter λ and is time independent. The hazard functions associated with the other probability models are all time dependent. A key aspect in the development of likelihood functions for a point process in general is the idea of conditional intensity of the form λ(t|t1tn) (t > tn) (Daley and Vere-Jones, 2003) (see Appendix). For a renewal process, the conditional intensity depends only on the time since the last event and equals υ(ttn) (Ogata, 1999b; Daley and Vere-Jones, 2003).

The likelihood function represents the probability that a fixed set of observations would be observed as a function of a given set of parameters (θ) for a particular distribution (Aitkin, 2010). The parameters that maximize the log-likelihood function (forumla) are termed the maximum likelihood estimate (MLE). Standard error estimates of inferred parameters based on the likelihood function depend on how well the function can be approximated by a quadratic function (Pawitan, 2001).

MLE Considering Open Intervals and Exact Event Times

Standard expressions for the log-likelihood function and MLE are known for each of the probability models described above, assuming that the occurrence times are known exactly and the open intervals are ignored. However, these expressions may be significantly biased for small sample numbers. Daley and Vere-Jones (2003) and Ogata (1999a, 1999b) describe expressions of the likelihood function, applicable to an arbitrary probability model, from point-process theory and the conditional intensity function. Considering the open intervals and events ti over S < t1 < t2 <…< tn < T, the log-likelihood function for a stationary renewal process (Ogata, 1999b; Daley and Vere-Jones, 2003) is
where νθ is the mean return time over the entire distribution (not the arithmetic mean of the samples) and is given by:

See Appendix for a summary of how the log-likelihood function shown in Equation 8 is derived. The log-likelihood function expressed by Equation 8 differs from standard log-likelihood formulas associated with the distributions and from those of previous paleoseismic studies that do not include the open interval before the oldest dated event and the –ln υθ term (Ogata, 1999b).

Using the probability models described above, we determine the MLE of the model parameters, assuming an exact occurrence age at the center of each interval indicated in Tables 1 and 2. The log-likelihood function is formulated from Equations 8 and 9. The maximum in the function is found by the Nelder-Mead direct search method of optimization (Nelder and Mead, 1965; Press et al., 2007). Results of parameter estimation for each probability model are shown in Tables 3 and 4.

Confidence intervals for the distribution parameters can be calculated from the likelihood profile, holding the other parameter at its MLE. The 95% confidence interval is obtained by comparing likelihood ratios to the chi-squared distribution with one degree of freedom (χ295%[1]):

The gamma and Weibull models have shape parameters whose domains span clustering and quasiperiodic behavior. For site U1322, the 95% confidence interval for α spans across the boundary α = 1 that divides clustering and quasiperiodic behavior for both models (Fig. 2). For site U1324, the lower bound of the 95% confidence interval for α is 3.59 for the gamma model and slightly greater than one (1.04) for the Weibull model (Fig. 3), suggesting quasiperiodic behavior. Sample size, however, needs to be considered in overall model selection (discussed in the Model Selection section).

Potential bias in standard log-likelihood functions based only on the return times can be determined by comparing mean return times from parameter estimates shown in Tables 5 and 6 with those from standard MLE expressions. For example, for the exponential distribution, the MLE using the standard expression equals the arithmetic or sample mean from the data (5.81 k.y.), whereas the MLE using Ogata’s (1999b) derived log-likelihood function that considers open intervals is significantly longer (8.22 k.y.). In each case, the sample mean and the standard MLE underestimate the mean return time (Parsons, 2008).

Type II MLE Considering Open Intervals and Uncertain Event Times

A Bayesian method was developed by Ogata (1999b) to account for uncertain event times, in addition to the effect of open intervals. Unlike traditional Bayesian inference techniques in which the prior and posterior distribution are a function of the parameter space (θ), Ogata’s (1999b) method was formulated such that the prior and posterior distributions are the probability densities of each event age distribution [ψi(ti); i = 1,…, n]. Thus, both the event ages and models parameters are uncertain. The prior distributions for event ages used by Ogata (1999b) include Dirac delta functions, and uniform and triangular distributions. Bayes’ theorem is then written as
where φ(t1,…,tn|θ) is the posterior distribution and ℒ[S,T](θ) is termed the integrated or marginal likelihood (Ogata et al., 2000; Aitkin, 2010). ℒ[S,T](θ) integrates the likelihood function with the joint distribution of the priors:

Ogata (1999b) provided details on the formulation of ℒ[S,T](θ) in a manner similar to the development of Equation 8. Using center-age Dirac delta functions for the priors recovers the likelihood expression given in Equation 8 exactly.

The model parameters forumla that maximize the integrated likelihood (Equation 12) is termed, for this type of statistical inference, Type II maximum likelihood (Good, 1965; Robert, 2007). The Type II MLEs for the different probability models are indicated in Tables 7 and 8 for the two IODP sites. Estimates of the rate parameter for exponential distribution using the center-age (Tables 3 and 4) and Type II likelihood methods are identical (Ogata, 1999b). For the other distributions, there appears to be more variation between the two methods in estimates of the shape parameters compared to the scale parameters. The mean return times that correspond to the Type II MLEs are slightly longer for site U1322 and nearly the same for site U1324 compared to the center-age MLEs (Tables 5 and 6). As indicated by Ogata (1999b), the center-age MLEs can be considered Type II MLEs in which Dirac functions are used as priors. Type II MLEs using different priors, such as triangular functions, will likely result in estimates intermediate between the center-age MLEs and the Type II MLEs using a uniform prior as described in this study.

Monte Carlo methods are an important alternative to the likelihood methods, particularly where the likelihood functions and their maximum are difficult to compute. The Monte Carlo method described in detail by Parsons (2008) takes into account age-dating uncertainty and the open intervals, as in the Type II likelihood methods developed by Ogata (1999b). Two probability models representing time dependence and time independence are compared using a Monte Carlo method that tallies the relative success of different distributions against sets of observed intervals. Time independence is represented by exponential distributions (Poisson process), and Brownian Passage Time (inverse Gaussian) distributions are taken to represent time dependence.

To make the test, a series of distributions that covers all reasonable mean recurrence intervals is developed (100 yr to 20,000 yr for the Ursa Basin MTDs). Time-dependent distributions are characterized by two parameters (Equation 5), and are thus also constructed across coefficient of variation (shape parameter) values between 0.01 and 0.99 for each mean return time. Groups of MTD times are randomly drawn 5 million times from each possible distribution and assembled into sequences. With this method, sequence means are identified directly from the parameters of parent distributions rather than from taking arithmetic means of sequences.

The Monte Carlo calculations use a uniform distribution for the event time windows (Tables 1, 2), and an event that happens at any time within the window is considered a match. The Monte Carlo sequences begin with an event that is given freedom to happen at any time prior to the first observed MTD time window. In this regard, it is different from Ogata’s (1999a) method in which no event is assumed to occur between time S and the first deposit. The extra event in the Monte Carlo method contributes nothing other than a starting point for the sampling. This is needed because the first observed time window has some range within which the event might have happened, whereas Monte Carlo simulation must begin at a point in time. A second open interval extending from the end of the last observed event window to the present is utilized; if a simulated event falls into this interval, its sequence is not considered a match.

To summarize the technique, exponential distributions with mean return times from 100 to 20,000 yr are generated in 100 yr increments, and a group of Brownian Passage Time distributions is also generated with mean return times from 100 to 20,000 yr in 100 yr increments. The time-dependent distributions also range through coefficients of variation between 0.01 and 0.99 (0.1 increments) for every mean return time. Each individual distribution is given 5 million chances to reproduce each MTD record within the uncertainties (Tables 1, 2). The number of successful matches is tallied for every possible distribution. The most likely distributions to represent Ursa Basin submarine mass failure occurrence are taken to be the ones that most frequently match the records.

Results from Monte Carlo analysis are shown in Figure 4, where the number of matches for each potential return time distribution is displayed. Generally, the record from site U1324 is easier to match than the U1322 sequence because there are fewer total events, and they are observed within wider time intervals (Tables 1, 2). However, the results of comparing a Poisson process against time-dependent models are the same; in both cases, many more matches are attained from exponential distributions than from time-dependent distributions (factor of 2:1 at the U1322 site, and 20:1 at the U1324 site; Fig. 4), even when considering the difference in number of model parameters (see the Model Selection section). Best-fit mean return times are given in Tables 9 and 10.

In summarizing results from Monte Carlo modeling, we note that it is possible to fit the MTD sequences found at the two IODP sites from the Ursa Basin to time-dependent models, particularly for higher (more random) coefficients of variation (Fig. 4, Tables 9 and 10). However, a better overall match can be made with random-in-time exponential distributions.

Comparison of different probability models for submarine mass failure occurrence is made using Akaike’s Information Criterion (AIC). AIC is a relative measure of the information lost when a particular model is used as an approximation to the true model. AIC balances the goodness of fit with the number of parameters needed to describe the model (Burnham and Anderson, 2010). The preferred model, among competing models for the same data and constraints, is the one with the lower AIC. The formula for AIC is given by
where K = dim{θ} is the number of parameters in the model. For small samples, there is a bias correction given approximately by (Hurvich and Tsai, 1989; Burnham and Anderson, 2010)
where n is the number of samples, in this case the number of MTDs.

AIC and AICc are used to evaluate the different models (Tables 3 and 4) from the MLE using center-age event times. Using AIC without the bias correction indicates that some of the quasiperiodic models perform better than the exponential model (BPT model optimal for both sites U1322 and U1324). This is consistent with α = 1 being outside the 95% confidence interval for the gamma and Weibull models (site U1324) discussed in the Maximum Likelihood Methods section (Fig. 3). However, using AICc with the small-sample bias correction indicates that the exponential model associated with a stationary Poisson process is the optimal model for both sites.

For the Type II MLE results based on Ogata’s (1999b) Bayesian inference method, Akaike’s Bayesian Information Criterion (ABIC) is used to evaluate the different models (Tables 7 and 8) (Akaike, 1980). ABIC is of similar form to AIC:

However, in this case, ℒ is the integrated likelihood function in the Bayesian method (Equation 12). The optimal model for submarine mass failure occurrence using Type II MLEs is the BPT model for both sites. It is unclear whether ABIC is biased for small sample sizes.

For the Monte Carlo method, one can view the number of successful matches over a finite number of possibilities as an approximation to the likelihood function. In this case, discrete values of the parameter space are used as explained in the previous section. In this view, the many more successful matches associated with the exponential distribution than the BPT solution suggests that maximum log-likelihood is greater for the exponential distribution, even without considering the number of parameters in each model.

The hazard-rate function (Equation 7) is calculated using the center-age MLE for each probability model and each site (Fig. 5). The dashed lines indicate the range in elapsed time from the last event to the present, given the uncertainties in dating. Formal techniques for incorporating this uncertainty in the hazard-rate function are given by Rhoades et al. (1994) and Ogata (1999b). At site U1322, the gamma and Weibull models exhibit increasing hazard rate at the present time, whereas the lognormal and BPT models exhibit decreasing hazard rate (Fig. 5A). For all cases, the hazard rate for the exponential model (Poisson process) is constant. As discussed in the previous section, the present-day hazard rate for the exponential model is optimal according to the AICc measure and Monte Carlo results, with the associated hazard rate intermediate at the present time between the two classes of quasiperiodic models.

The difference between the hazard functions at the two sites relates to the mean return time relative to the time since the last event. At site U1324, all of the quasiperiodic models exhibit increasing hazard rate at the present time (Fig. 5B). The time since the last event at site U1322 is more than twice the mean return time, such that the hazard-rate functions are approaching their asymptotic values. In contrast, the time since the last event at site U1324 is comparable to the mean return time. As such, there is a significant difference between the present-day hazard rate estimated by the optimal exponential model and the quasiperiodic models at site U1324. If it were found that the quasiperiodic models are optimal in terms of estimating the occurrence of submarine mass failures at site U1324, as suggested by the Type II MLE results where the BPT model is optimal, then the present-day hazard rate would be considerably higher compared to the hazard rate associated with the exponential model.


Both the maximum likelihood techniques and the Monte Carlo technique for estimating the occurrence probability of submarine mass failures include the effect of the open time intervals on parameter estimation. The effect of the open intervals alone can be examined by comparing the standard MLE with the center-age MLE derived by Ogata (1999b). For both drill sites and all probability models, inclusion of the open intervals increases the estimated mean return time for submarine mass failures (Tables 5 and 6). The Monte Carlo technique and the Type II maximum likelihood technique also account for uncertainty in age dating events. Uncertainty is endemic in estimating event ages from geologic samples, such as the age of mass transport deposits sampled by drill holes as discussed in the age dates section above. The effect of this uncertainty partly depends on how the uncertainty is distributed for each event, incorporated in the Type II likelihood method as prior distributions in a Bayesian framework. End members examined in this study are uniform prior and the Dirac prior, the latter of which reduces to the center-age MLE. Age uncertainty using a uniform prior tends to slightly increase the mean return time for site U1322 compared to center-age MLEs, though there is little effect for site U1324. For each method and probability model, it is possible to define parameter confidence intervals as shown in Figures 2 and 3.

Differences between the Type II likelihood and Monte Carlo methods can be examined by comparing the results from the BPT model. For the Type II likelihood method, the mean return time is 8.62 and 13.9 k.y. for sites U1322 and U1324, respectively (Tables 5 and 6). In comparison, the mean return time estimated by the Monte Carlo method is 4.03–4.37 and 8.25 k.y. for sites U1322 and U1324, respectively (Tables 9 and 10). In addition, the Type II likelihood technique appears to be estimating a significantly lower value for the shape parameter (coefficient of variation, α) at site U1324, than the Monte Carlo technique (Tables 8 and 10). The difference in the results from the two methods may relate to how the open intervals are handled. In the case of the likelihood methods, no event is assumed to occur between age S (base of the Blue Unit, estimated to be 74 ka) and the first deposit, whereas in the case of the Monte Carlo method an event can occur at any time prior to the first deposit. Because of the small number of events, particularly at site U1324, the latter may result in more Poisson-like behavior and shorter estimates for the mean return time than the likelihood methods for this case study.

The computational requirement of each technique depends on the number of events and event-age uncertainty relative to the mean return time and the probability model used. For the Monte Carlo results shown in Figure 4, twenty times the number of runs needed to be made for site U1322 than for site U1324 (100 million runs versus 5 million runs for each pair of distribution parameters). In contrast, the computational requirement for the Type II likelihood method depends on the number of events, in terms of the number of integrations required to define the likelihood function (Equation 12) and the analytical form of the probability model.


Using these methods, we can assess whether submarine mass failures recorded at the two IODP sites occur randomly (i.e., according to a Poisson process, exhibiting an exponential distribution of inter-event times) or whether we can say with confidence that there is a clustering or quasiperiodic aspect to these events. For the center-age MLE, the exponential distribution of return times associated with a Poisson process is the optimal model for both sites according to the AICc model selection criterion. Because of the small number of samples, the correction factor to AIC is particularly large, in favor of the exponential model. If more events were identified, it is possible that the BPT model would have been selected, as indicated by the raw AIC values in Tables 3 and 4. The ABIC model selection criterion for the Type II MLE method indicates that the BPT model is optimal for both sites, however no small-sample bias correction is applied to ABIC as it is for AICc. The Monte Carlo method that considers both age-dating uncertainty and open intervals indicates that the exponential model is optimal. As indicated above, the difference between the likelihood and Monte Carlo results likely relates to how the open interval before the first event is treated. The difference in the resulting hazard function between the exponential model and quasiperiodic models can be significant, as shown in Figure 5.

The difference in the overall rates and mean return times between the two sites relates to the simple fact that fewer MTDs were present at site U1324 (Urgeles et al., 2007). There may also be differences in the physical processes relating to the long-term occurrence of submarine mass failures. The two sites are in slightly different geographic positions relative to the ancestral canyons and their levee deposits (Sawyer et al., 2009; Stigall and Dugan, 2010). Hence, the record of MTDs at each site may be a function of the basin architecture. In addition, physical properties such as pore pressure between the two sites are slightly different, indicative of regional variation of sediment accumulation and fluid flow focusing (Flemings et al., 2008; Stigall and Dugan, 2010).

Differences between the occurrence results of subaerial landslides that show clustering behavior (Witt et al., 2010) and the results of this study that primarily show Poisson/quasiperiodic behavior may relate to differences in the detection level between the two studies. Although it is difficult to assess the range of landslide sizes among the studies, it is certain that subaerial landslides examined by Witt et al. (2010) (with time scales measured in days) are smaller in size than the submarine mass failures considered here and in Griffiths (1993) (with time scales measured in thousands of years). Smaller landslides are likely conditional to a larger predecessor landslide (e.g., secondary failures along a scarp from a predecessor landslide) and therefore tend to cluster temporally. The MTDs record larger events, possibly from independent sources. There does not appear to be a strong time-dependent failure process underlying the timing of the MTDs in this setting. Although some events may be part of a larger retrogressive system (Sawyer et al., 2009), this does not necessarily imply that events are conditional on preceding events, in terms of temporal occurrence. It should be noted, however, that there may be fundamental differences in the occurrence of subaerial landslides compared to submarine mass failures, analogous to the differences in size distributions in the two environments (ten Brink et al., 2009a).

Finally, it is assumed throughout this study that the probability models are stationary. The driving force for basin formation is salt withdrawal, which is likely stationary relative to the geologic record of the MTDs. However, sediment input and accumulation may be dependent on periods of glaciation and changes in sea level (Kolla and Perlmutter, 1993; Lee, 2009; Piper and Normark, 2009; Korup et al., 2012) and thus may not be stationary. Testing of the stationary assumption is likely only possible for global data sets until more age dates of mass transport deposits are acquired.

Techniques have been developed to empirically determine the probability distribution of submarine mass failures from dated mass transport deposits. Adapted from paleoseismic methods, these techniques include the uncertainty in age dating the deposits and the open time intervals before the first and after the last deposit. The renewal-process likelihood techniques include those for center-age occurrence, accounting for open intervals, and Type II likelihood methods under a Bayesian framework to also accommodate uncertainty in the age dates. Results from the likelihood methods can be used with model selection criteria, such as AIC, and statistical tests relative to a Poisson null hypothesis, such as the likelihood ratio test. The Monte Carlo technique draws many samples for a specific distribution parameter(s) and repeats these draws for all possible relevant parameter values. The most likely parameters for a given distribution are the ones with the most matches to the data. The Monte Carlo technique is particularly useful when the likelihood functions or their maximum cannot be computed.

Results from mass transport deposits recorded in IODP sites U1322 and U1324 within the Ursa Basin indicate that the deposits can be described by a Poisson process, based on the center-age maximum likelihood and Monte Carlo methods, or a quasiperiodic process, based on the Type II maximum likelihood method. For the Poisson process, the return times are exponentially distributed and the hazard-rate function is a constant. The mean return time of submarine mass failures is shorter at site U1322 primarily because more mass transport events are identified over a similar period of time. Differences in estimates of the mean return time between the likelihood and Monte Carlo methods are related to how the open interval before the first event is treated. For the likelihood methods, no event is assumed to occur between a specified age S and the first event, whereas for the Monte Carlo method any event can occur before the first event. Application of likelihood and Monte Carlo techniques to the Ursa Basin data demonstrates how these methods can estimate parameters and confidence intervals for various probability models and test different hypotheses for the occurrence of submarine mass failures.

We thank Roger Urgeles for discussions regarding the age dating of mass transport deposits at the Ursa Basin IODP sites. We also gratefully acknowledge helpful and highly constructive reviews by George Hilley, Danny Brothers, and an anonymous reviewer.


The development of the likelihood function (Equation 8) is summarized here from Daley and Vere-Jones (2003) and Ogata (1978, 1999b). Consider a history of event times Ht before time t0: 0 = t0 < t1 …< tn < t in which no event can occur at the same time. For a point process such as this, the conditional probability that an event will occur in a small interval of time δ is given by:
where N(t) is a counting measure and λ(t|Ht) the conditional intensity function. In the limit as δ → 0, the conditional intensity can be thought of as the instantaneous rate of occurrence. For a finite point process on the interval (0, T), the likelihood function in terms of the conditional intensity is given by (Daley and Vere-Jones, 2003)
A renewal process is one in which the events are independent and identically distributed. In this case, the conditional intensity function depends only on the time since the last event. The conditional intensity function parameterized by vector θ is given by
where h(x) is the hazard function, f(x) is the probability density function, and F(x) is the cumulative distribution. For a finite renewal process, Daley and Vere-Jones (2003) showed that the likelihood function is given by
Finally, consider a delayed renewal process in which the history of events does not necessarily begin with t0 = 0 (Ht: 0 ≤ t1 < t2 …< tn < t). The first event is independent from, but not necessarily identically distributed with, the other events. If we designate the density and cumulative distribution functions associated with the first event by g(t) and G(t), respectively, the conditional intensity λ = g(t)/[1 – G(t)] if N(t) = 0, and Equation A3 otherwise. The likelihood function in Equation A4 is modified such that
The case where g(t1) = fθ(t1) results in the ordinary renewal process above. If the delayed renewal process is stationary, then g(t) = υθ−1[1 – Fθ(t)] (Daley and Vere-Jones, 2003), where forumla is the mean of the distribution f. The corresponding likelihood function is
For a delayed renewal process on the interval (S,T),

The log-likelihood function is then that given by Equation 8.