Vertical Stacking Statistics of Multi-facies Object-Based Models

Equations describing facies proportions and amalgamation ratios are derived for randomly placed objects belonging to two or three foreground facies embedded in a background facies, as a function of the volume fractions and object thicknesses of independent facies models combined in a stratigraphically meaningful order. The equations are validated using one-dimensional continuum models. Evaluation of the equations reveals a simple relationship between an effective facies proportion and an effective amalgamation ratio, both measured as a function only of the facies in question and the background facies. This relationship provides a firm analytical basis for applying the compression algorithm to multi-facies object-based models. A set of two-dimensional cross-sectional models illustrates the approach, which allows models to be generated with realistic object stacking characteristics defined independently for each facies in a multi-facies object-based model.


Introduction
The objectives of this work are specific and straightforward: to establish analytical expressions for certain statistics of one-dimensional vertical samples through objectbased models. The motivation for deriving precise expressions is so the stacking properties of modelled systems can be conditioned to match those of natural systems, as part of broader endeavours to build more realistic subsurface models. Classifica-B Tom Manzocchi tom.manzocchi@ucd.ie Deirdre A. Walsh deirdre.walsh@icrag-centre.org tion of sedimentary deposits into distinct facies is standard practice in both academic sedimentological studies and practical subsurface modelling, and the present work is focused on systems comprising one, two or three foreground facies embedded in a background facies. The foreground facies are composed of objects that represent discrete depositional elements. The geometries and stacking characteristics of elements within a particular facies are assumed to be similar to each other, but are different between the facies. Hence, for example, a deep-water system might contain discrete lobate and channelised elements modelled as separate facies, and the geometries and stacking characteristics of the elements will vary depending on whether the modelled volume is contained within the channel, channel-lobe transition zone or lobe of the large, fan-scale deposit (e.g. Cullis et al. 2018;Fryer and Jobe 2019).
The statistics of interest to this study are facies proportions and object amalgamation ratios. Facies proportions (defined as the fractional volumes of the model occupied by each facies) are a common measure usually specified as an input property of objectbased models. The amalgamation ratio is a measure of the connectivity of the objects and is usually an unconstrained output property of object-based models. It is typically measured at the bed scale, where it is defined as the number of bed bases in contact with an underlying bed (as opposed to shale) as a proportion of the total number of bed bases present. This definition follows Chapin et al. (1994) and was used by Stephen et al. (2001), Manzocchi et al. (2007) and Romans et al. (2009). More recently, Jobe et al. (2021) and Stanbrook and Bentley (2022) discussed practicalities of acquiring the ratio from one-dimensional data. A quite different definition of amalgamation ratio is also sometimes used in the literature (e.g. Mueller et al. 2017;Pyles et al. 2019). This alternative definition was introduced because of the difficulty of identifying individual amalgamated surfaces and is based on the thicknesses of amalgamated and non-amalgamated intervals, rather than on counts of bed bases. It is unfortunate that the same name is used for the two ratios, because although they will usually be strongly correlated, they are not the same. Considerations by Manzocchi et al. (2020, their Fig. 4) show how the former ratio might be calculated from the latter, based on an assumed bed thickness. Other measures of vertical stacking which are closely related to the amalgamation ratio have also been discussed (e.g. Funk et al. 2012;Jackson et al. 2019).
Bed-scale amalgamation ratios have been particularly well characterised in the slope and basin floor submarine channel systems of the Magallanes Basin in Chile (e.g. Romans et al. 2009;Jobe et al. 2010;Macauley and Hubbard 2013). In general, channelised systems have much larger amalgamation ratios than basin floor fan systems (Mattern 2002) for which amalgamation ratios have been measured by Manzocchi et al. (2007), Marini et al. (2015) and Pyles et al. (2019). Amalgamation ratios have also been examined at multiple, larger hierarchical levels in both lobate and channelised deep-water systems (e.g. Zhang et al. 2015;Manzocchi et al. 2020;Soni et al. 2020) with the general trend that larger scales of depositional elements tend to be more poorly amalgamated than smaller ones.
We have shown elsewhere that natural depositional systems generally have lower amalgamation ratios than are obtained in object-or pixel-based models (e.g. Walsh and Manzocchi 2021a;Manzocchi et al. in press a) and have proposed an approach (compression-based modelling) in which a geometrical transformation to a facies Cross sections through models with the same facies proportions and object shapes for two foreground facies (yellow sheets and green channels) in a black background facies. See text for discussion model results in stacking properties that cannot be obtained in conventional models but which are representative of natural systems (e.g. Manzocchi et al. 2007Manzocchi et al. , 2020Walsh and Manzocchi 2021a, b). Examples generated using cross-sectional object-based models consisting of two foreground facies with sheet-and channelshaped objects are shown in Fig. 1. Figure 1a shows a model in which the compression algorithm has not been applied, and hence the facies amalgamation ratios are unconstrained. The channels are more erosive than the sheets because they are thicker, but there are places in the model where sheets erode into the top of channels. The object stacking pattern in the model (Fig. 1a) is created by the stratigraphically appropriate facies preservation rule included in the object-based modelling algorithm and will be similar for all object-based model realisations created with these facies proportions and shapes. The other three models (Fig. 1b to d) have the same facies proportions and shapes, but in these cases the compression algorithm has been used to produce models with different staking characteristics, resulting in more poorly amalgamated sheets with unconstrained channels (Fig. 1b), more poorly amalgamated sheets and channels (Fig. 1c) and more poorly amalgamated sheets but more highly amalgamated channels (Fig. 1d). These two-dimensional models demonstrate the impact that object amalgamation ratios can have on the spatial distributions of the facies, but also highlight the most notable geometrical artefact of the compression algorithm which is the imposition of dip caused by lateral facies transitions ). This artefact is not addressed in this paper which focuses on one-dimensional model characteristics.
The transformation within the compression method applied to object-based models relies on a statement of the target values of facies proportions and amalgamation ratios, and expressions defining how these properties are related to each other in the initial object-based model prior to the transformation. This paper is concerned with estab-lishing these expressions for multi-facies models. The approach taken is to derive algebraic expressions for the stacking characteristics of idealised one-dimensional systems and to validate these expressions using numerical models. The following section defines in detail the system examined and the stacking characteristics of interest, and introduces the one-dimensional object-based modelling code used to illustrate the problem and validate the results. Sections 3, 4 and 5 develop analytical expressions for systems containing one, two and three foreground facies. The results are verified at the end of their respective sections, and Sect. 6 discusses the results and illustrates how they might be used in practice. Appendix 1 lists, and explains the logic behind, the abbreviations and symbols used throughout the paper.

Definition of the Problem
The one-dimensional system considered in this paper consists of up to three foreground facies (F1, F2 and F3) embedded in a background facies (F0) over a total model height (H ). For most of the derivations in this paper, H is assumed to be 1 and therefore does not generally appear in the equations. In certain places, however, it is useful to refer to H explicitly. Throughout the paper, X is used as a general facies indicator to denote any of F1, F2 or F3. Objects of facies FX have constant thickness (T X ), and the facies are ordered and named to satisfy the condition T 1 ≤ T 2 ≤ T 3 . The background facies (F0) is not modelled as discrete objects, and therefore is not associated with a particular object thickness value. The models considered in this paper are non-hierarchical, and the results could apply to any level in a depositional hierarchy. Hence, although the objects are referred to as beds throughout this paper for convenience, they could just as well represent lobes, channels or other depositional elements. The assumption of constant thickness (T X ) for objects belonging to each facies simplifies the derivations in the paper. Implications of the assumption for models with a single foreground facies are discussed by Manzocchi et al. (in press b), who show that models with broader thickness distributions have a marginally lower amalgamation ratio for the same net:gross ratio.
The code used to generate the numerical models includes many of the same assumptions as the analytical solutions derived later, and therefore a model realisation is used to illustrate the problem being addressed (Fig. 2). Initially, each facies and each bed is treated independently. Beds are assumed to be positioned randomly and are permitted to overlap with each other and with the top and bottom edges of the volume of interest (Fig. 2a). Placement of beds for each facies continues until the facies occupies a specific volume fraction (V X ) of the total model height. The random placement of beds implies that the probability that each facies is present at any location in the initial bed model is simply equal to V X and is independent of location in the column. This assumption is consistent with continuum Boolean modelling approaches (Shante and Kirkpatrick 1971;Baker et al. 2002) and is central to the analytical derivation of object stacking statistics described in this work.
The initial model (Fig. 2a) can contain multiple beds of the same or different facies at any position which must be rationalised to produce a stratigraphically meaningful stacked model. The stacked bed model (Fig. 2b) has only one facies present at any location and is created from the initial model by assuming that the preservation of overlapping beds is governed by Steno's law of superposition (Steno 1669). Hence, beds which have their tops higher in the sequence are assumed to be stratigraphically younger and to erode into stratigraphically older beds which have their upper surface further down the sequence. The beds in the numerical model are generated in a continuum, so there is never any ambiguity as to which bed is younger, and only one sequence of beds can be generated from the initial model which satisfies the law of superposition. The stacked model contains fewer beds overall than the initial model, since some are entirely eroded out of the sequence. Many of the remaining beds are partially eroded by an overlying bed, and therefore the beds in the stacked model have variable thicknesses.
The facies model (Fig. 2c) is the same as the stacked model but does not contain any information about the contacts between different beds. The most important properties of the facies model are the facies proportions (P X ), which is one of the fundamental properties defined as input in practical facies modelling exercises. It is clear from Fig. 2c that in this example, P 3 > P 2 > P 1 , despite the fact that the same volume fractions were used to generate the initial model (V 1 = V 2 = V 3 = 0.4, Fig. 2a). This is because where two facies overlap in the initial model, the thicker bed is more likely to have a higher top than the thinner bed, and therefore is more likely to be preserved given a stratigraphically sensible stacking order. Establishing a precise definition of P X as a function of V X and T X is one of the objectives of the analyses in this paper.
The version of the model shown in Fig. 2d is the same as the facies model ( Fig. 2c) except that all contacts between beds are coloured: these contacts are termed amalgamations and are labelled using the convention that FXY represents an amalgamation surface defined by an FX bed above an FY bed. Using this scheme, a non-amalgamated FX bed base that overlies F0 is termed FX0. The choice of colour used depends on the facies above the contact. Hence, blue lines represent amalgamated bases of F1 beds, red lines amalgamated bases of F2 beds, and green lines amalgamated bases of F3 beds. The black horizontal lines (Fig. 2d) represent unamalgamated bed edges where either the top of the bed is overlain by shale (F0) or the base of the bed overlies shale. Quantifying the proportions of amalgamated bed bases is an important measure of the overall three-dimensional connectivity of the facies, and if only one facies is present, the amalgamation ratio (AR) is equal to the facies proportion if all beds are the same thickness (Manzocchi et al. 2007). Establishing an expression for AR in the presence of multiple facies is the second objective of the analyses below.
Two different definitions of facies-specific AR are used in this paper: total AR and effective AR ( AR T X and AR E X , respectively). Both measures were used also by Manzocchi et al. (in press b) and are defined as follows, using the term n XY to signify the number of beds of facies X that amalgamate with a bed of facies Y, such that AR T X = n X 1 + n X 2 + n X 3 n X 0 + n X 1 + n X 2 + n X 3 ; (1) The difference between AR T X and AR E X is shown by the difference between Fig. 2d, e: the former includes all amalgamations at the bases of beds of the facies of interest, while the latter considers only amalgamations with beds of the same facies.
The amalgamation ratios used in this paper are related to, but distinct from, the embedded facies transition probabilities often used to characterise stratigraphic logs (Krumbein and Dacey 1969). Amalgamation ratios refer to the bases of individual beds and are therefore based on the information contained in the stacked bed model (Fig. 2b). Facies transition probabilities, however, refer only to facies distributions and therefore are based on the subset of information contained in the facies model (Fig. 2c). Amalgamation between beds of the same facies is not a characteristic recorded in the facies model, and n X X is undefined in the facies transition probability approach. The association between the two sets of parameters can be appreciated by stating the embedded facies transition probability ( p XY ) with the terminology used in this paper. In this way, the probability that facies FX transitions downwards to facies FY is given by This equation contains n X X twice: one explicitly, and once as one of the other terms in the denominator. These cancel out, so p XY does not depend on n X X (the equation has been written in this form so it can be generalised). Facies transition probability matrices lie at the heart of the Markov chain approach for modelling stratigraphic sections (Krumbein and Dacey 1969;Fogg 1997, 2020;Li et al. 2019), but there is no fundamental reason why amalgamation ratios (or the associated bed transition probability matrix) should not be used instead. In the Markov chain approach, a volume is populated with facies sequentially and in a specific direction, resulting in non-random sequences built as a function of input facies transition probabilities. However, this is not the approach taken in this paper which instead assumes the random bed placement method described above and calculates the resultant amalgamation ratios from the resultant sequence. The reason for taking this approach is because it is closely allied to the object-based modelling approach often used to build two-and three-dimensional facies models (Fig. 1). As discussed in the Introduction, if the amalgamation ratios of object-based models are known as a function of model input parameters, a geometrical transformation can be used to create models with specific, non-random amalgamation ratios representative of natural geological systems.

Stacking Properties of Single-Facies Models
The properties of models containing only one foreground facies (F1) are simple. It is obvious that P 1 = V 1 since no other objects are present. Less obvious is the fact that AR E1 = V 1 . The reason for this was explained by Walsh and Manzocchi (2021a) as follows: if all beds are of equal thickness, then it is impossible for a higher bed to erode the base of a lower bed. Therefore, all bed bases are preserved in the stacked bed model. Since the beds are located randomly, so too are the bed bases. Hence, there is a probability of exactly V 1 that each bed is positioned on top of (and therefore is amalgamated with) a lower bed, and a probability of exactly (1 − V 1 ) that it is not amalgamated. Therefore, if there are n 1 beds altogether, the expected numbers of F11 and F10 bed bases are and These can be combined using either Eq. 1 or 2, since the different versions of AR are the same when only one facies is present. This yields Bed stacking can be considered in a bit more detail to derive the distribution of bed thicknesses in the stacked model. This is necessary for the solutions involving two or three foreground facies derived in the following sections. Shante and Kirkpatrick (1971) showed that the expected total number of beds in the system (n 1 ) is given by This is because, if a bed is added to a system that already has a facies proportion P 1 , then the facies proportion will increase, on average, to P 1 + T 1 H (1 − P 1 ). Expressing this as a differential equation gives which, when integrated, gives Eq. 7. The bed bases are located at random, so it must be the case that the distances between bed bases follow a negative exponential distribution. When these distances are less than the original thickness of the beds (T 1 ), the distance must represent the thicknesses of beds which have been partially eroded and are capped by an amalgamation surface (i.e. t 1 < T 1 ). It is impossible for the bed thickness to exceed T 1 , and therefore, distances greater than T 1 represent beds that are preserved in their entirety (t 1 = T 1 ), plus an overlying region of F0. Hence, the distribution of preserved F1 bed thicknesses is where n t1 is the number of beds of thickness in the range t 1 to t 1 + dt 1 . The ideal distribution (Eq. 9) is compared to that obtained in three realisations of two systems (Fig. 3). Differences between realisations reflect the small size of the systems considered in this illustration, but they clearly follow the expected distribution.

Introduction
When two foreground facies (F1 and F2) potentially with different thicknesses are present (T 2 ≥ T 1 ), there is the potential to entirely erode out F1 beds and, as discussed in the Introduction, we expect that P 2 ≥ P 1 even if V 2 = V 1 , because thicker beds are more likely to be erosive. Calculation of facies proportions and amalgamation ratios are therefore much more complicated than when only one facies is present.
The initial model ( Fig. 4a) can be converted to independent single-facies stacked models in which P 1 = V 1 and P 2 = V 2 (Fig. 4b, c). Since the location of beds is entirely random, there is a probability equal to V 1 V 2 that any position is occupied by beds belonging to both F1 and F2. The particular bed retained at these locations in the combined stacked model (Fig. 4d) is a function of the locations of the two bed tops Stacked bed models of systems with T 1 = 1m, H = 100m and a V 1 = 0.3, b V 1 = 0.7. Blue lines show amalgamation surfaces. c, d Analytically derived cumulative distribution (black) and cumulative distribution of three realisations (coloured lines) of resultant bed thickness (t) for the models shown in a and b, respectively and will be calculated later. For now, it is simply expressed as a facies preservation probability ( p 1|12 ) that F1 rather than F2 is preserved at locations where both are present. It follows, therefore, that the facies proportions in the combined stacked model (Fig. 4d) are and  Table 1 and text for discussion The amalgamation ratios are calculated as a function of different types of bed bases (Eqs. 1, 2). An F1 bed base in a single-facies stacked model can be one of two types, F11 or F10, with the quantities of each given by Eqs. 4 and 5. In positions where a bed base in a single-facies model overlies a bed from the other facies, the surfaces can either be preserved, changed to a different type of surface, or eroded out of the system altogether. The possible outcomes for these surfaces are summarised in Table 1. The outcome probabilities in Table 1 are expressed in the form p X X XY |XY , which means that an FXX surface from a single-facies model is preserved as an FXY surface in the two-facies model when the surface is directly above a position occupied by beds of both facies X and Y in the independent single-facies stacked models. In these expressions, a hyphen (−) is used to denote erosion, so p X 0−|Y refers to the probability that an FX0 surface is eroded out of the system where an FY bed overlaps it. Further details of the precise logic applied to defining the symbols are given in Appendix 1. These surface outcome probabilities ( p O B ; Table 1) apply only where the bed base coincides with a position occupied by the other facies. This coincidence is defined by a surface overlap probability ( p O A ) which is equal to the volume fraction of the other facies. If the other facies is not present at this locality, the base must be preserved. So, for example, an F11 base cannot be transformed into anything else where there is nothing else to interfere with it. Hence, p 1111|1 = 1 and p 2222|2 = 1. Table 1 omits terms for p 22−|1 and p 20−|1 , which would both be associated with an F1 bed eroding out an F2 bed base. This is impossible (i.e. has a probability of zero), since the facies are defined to satisfy the condition that T 1 ≤ T 2 , and a bed cannot erode down to the base of a thicker (or equal thickness) bed.
The surface outcome probabilities listed in Table 1 (which will be calculated in Sect. 4.3) can be combined with the number of bed bases present in the single-facies models (Eqs. 4 and 5) and the surface overlap probabilities ( p O A ) to calculate the number of bed bases of different types (n XY ) which are required to calculate the amalgamation ratios (Eqs. 1 to 3). For example, where the first term (n 1 V 1 ) is the number of F11 surfaces in the single-facies model (Fig. 4b, Eq. 5), the second term (1 − V 2 ) is the proportion of these surfaces that do not coincide with F2 in the other single-facies model (Fig. 4c) and therefore cannot be altered, and the third term (V 2 p 1111|12 ) gives the probable proportion of these contacts that could be altered but which are not. Performing a similar calculation for the other bed contact types, we obtain and The discussion above and Eqs. 1, 2 and 10 to 17 show how the facies proportions and amalgamation ratios can be calculated in a two-facies model. The equations contain various probabilities that define the likelihood that different facies, or types of bed bases, are preserved when the stack model (Fig. 4d) is created from the single-facies models (Fig. 4b, c). The following two sub-sections address these probabilities within the context of the facies proportions and amalgamation ratio calculations.

Facies Proportions
The key to the facies proportion calculations (Eqs. 10 and 11) is the probability term p 1|12 , which expresses the probability that F1 is preserved where both F1 and F2 are present. To solve p 1|12 , it is necessary to consider location as a function of distance below the top of a bed irrespective of whether or not this bed top has been eroded. This requires consideration of three length terms. The first two are T X and t X , which are the thickness of the initial input beds and of the resultant, potentially eroded, beds in the single-facies stacked models, respectively. The third length term (z X ) is defined as a distance below the original bed top (Fig. 5). The approach taken is to calculate the total thickness of preserved facies X material present at a particular distance (i.e. between z X and z X + dz; Fig. 5) below bed tops in the single-facies stacked models (e.g. Fig. 4b, c). This thickness (which is termed Lz X and which increases with increasing z X since higher portions of beds are more likely to be eroded) is then used to derive the probabilities required. For every bed of preserved thickness t X ≥ T X − z X , there is a thickness dz that contributes to Lz X (Fig. 5). Hence, where n t X is the number of beds of thickness t X which is given by Eq. 7. Replacing this in Eq. 18 and solving, we obtain Figure 5 illustrates a situation in which both an F1 and an F2 bed are present within the distance z X to z X + dz. If z 1 > z 2 , F1 will be preserved and F2 eroded, and if z 1 < z 2 (the case illustrated in Fig. 5), F2 will be preserved and F2 eroded. The case z 1 = z 2 does not need to be considered, since the beds are considered within a continuum in which it is vanishingly probable that z 1 = z 2 . For a particular value of z 1 , the probability that F1 is preserved is given by Replacing Eq. 19 into Eq. 20, we have To get the full probability p 1|12 for all values of z 1 requires integration over the full possible range of z 1 . Hence, which ultimately gives Applying this in Eqs. 10 and 11 gives the facies proportions P 1 and P 2 as a function of expected number of beds (n 1 and n 2 , given by Eq. 7) and the system input parameters T 1 , T 2 , V 1 and V 2 . Note that n 1 and n 2 are themselves functions of T 1 , T 2 , V 1 and V 2 , and so could have been eliminated as terms in this expression. However, we have Fig. 6 a An F1 bed of original thickness T 1 is partially eroded to leave a thickness t 1 . The dashed and solid blue lines represent the original top and base of the overlying eroding bed; the dashed and solid black lines are the same for the lower, eroded bed. b An F2 bed also occupies the same locality as the amalgamation surface (within the shaded region dz). Whether the F11 surface (the solid blue line) is eroded out, remains an F11 surface or is transformed into a F12 surface depends on the locations of the two F1 beds and the F2 bed found it more convenient here and elsewhere to keep n 1 and n 2 as separate terms in the equations.

Amalgamation Ratios
Calculation of the amalgamation ratios requires consideration of different probabilities for four different types of surface in the single-facies models (Table 1). These are considered in turn below, in each case by first considering a constant value of z 1 or z 2 (Figs. 6, 7) and then integrating over the full range of possible values of z 1 or z 2 to provide the final probabilities required. Figure 6 considers the possible outcome of an F11 surface (Fig. 6a) when its position coincides with an F2 bed (Fig. 6b). There are three possible outcomes. If z 2 > T 1 , it will be eroded out of the system because all of the upper F1 bed, including its amalgamated base, will be replaced by the F2 bed. If T 1 > z 2 > z 1 , it will be transformed into an F12 surface because the top of the lower F1 bed, but not the base of the upper F1 bed, will be replaced by F2. If z 1 > z 2 , it remains an F11 surface since the lower F1 bed erodes out this portion of the F2 bed. These three cases are defined by the probabilities p 11−|12 , p 1112|12 and p 1111|12 , respectively (Table 1). The first of these does not need to be calculated, as it does not feature in any of the equations used to determine the number of different bed base types (Equations 12-17), but the other two do. Fig. 7 a An F2 bed of original thickness T 2 is partially eroded to leave a thickness t 2 . The dashed and solid red lines represents the original top and base of the overlying eroding bed, and the dashed and solid black lines the same for the lower, eroded bed. b An F1 bed also occupies the same locality as the amalgamation surface (within the shaded region dz). Whether the F22 surface (the solid red line) remains an F22 surface or is transformed into a F21 surface depends on whether the original top of the F2 bed is above or below the original top of the F1 bed. The surface cannot be eroded out, since the top of the upper F2 bed must be above the top of the F1 bed A similar procedure to the one applied to derive p 1|12|z 1 (Sect. 4.2) yields

Outcomes for F11 Surfaces
and Integrating over the full possible range of z 1 values gives the probabilities and Note that p 1111|12 (Eq. 26) is identical to p 1|12 (Eq. 23), but p 1112|12 is not the same as p 2|12 (which was not calculated explicitly in Sect. 4.2 since p 2|12 = 1 − p 1|12 ). This is because a surface can be eroded entirely out of the system, so there are three possible outcomes for each surface. But one or other facies must be preserved, so there are only two possible outcomes for facies proportions.

Outcomes for F10 Surfaces
There are only two possible outcomes for unamalgamated F1 bed bases: they either become F12 surfaces, or they are eroded out of the system. Only the first of these outcomes feature in Eqs. 13 to 17, so only p 1012|2 needs to be calculated. Transformation of the bed base to an F12 surface requires z 2 < T 1 . Therefore, which gives

Outcomes for F22 Surfaces
For the F22 surfaces, two probabilities are important ( p 2222|12 and p 2221|12 ), but only one needs to be calculated since p 2222|12 = (1− p 2221|12 ). The criterion differentiating between these possibilities is whether or not z 1 > z 2 (Fig. 7). Considering the case where z 1 > z 2 , for a particular value of z 2 we have This integrates across all values of z 2 to give Additionally, as discussed,

Outcomes for F20 Surfaces
F20 bed bases cannot be eroded out of the system. Therefore, the only thing that can happen to them if their position in the F2 facies model (Fig. 4c) coincides with F1 beds in the F1 facies model (Fig. 4b) is that they are transformed into a F21 contact. Therefore,

Validation of the Two-Facies Results
Results are validated by running a simple experiment consisting of 27 (i.e. 3 3 ) onedimensional continuum models, comprising all possible combinations of three settings of three variables. Variables are V 1 and V 2 , for which settings of [0.15, 0.5, 0.85] are used, and T 2 using values of [1.0m, 1.5m, 5.0m]. T 1 is kept constant at 1.0m, and model resolution is defined by setting the model height (H ) to a value consistent with an expectation of at least 100 beds for each facies. Hence, H is different in different models and is calculated by taking the model height to be the minimum of H 1 and H 2 , where (following Eq. 7) these facies-specific heights are given by H X = − 100T X ln(1−V X ) . Example models are shown at a scale of H /10 in Fig. 8, and results for 100 realisations of all systems in Fig. 9. These graphs show unbiased 1:1 relationships between the analytically derived expected values and modelled values of facies proportions (Fig. 9a, d) and amalgamation ratios (Fig. 9b, c, e, f), and therefore confirm the equations. Fig. 9 Comparison between the analytically derived facies proportions and amalgamation ratios, and those obtained in the numerical models, for the two facies cases. Error bars represent +/-one standard deviation around the mean value obtained in 100 realisations. See text for discussion

Stacking Properties of Models Containing Three Facies
The solution for three facies follows the same basic considerations as for two facies. The facies are ordered such that T 1 ≤ T 2 ≤ T 3 , and the objective is to derive the facies proportions (P X ) and amalgamation ratios ( AR T X , AR E X ) as a function of these thicknesses, and of the three input volume fractions V 1 , V 2 and V 3 .

Facies Proportions
The probability that all three facies overlie each other at any location is V 1 V 2 V 3 , and that two facies (FX and FY) but not the third (FW) are present at any location is V X V Y (1 − V W ). Hence, the facies proportions in a combined three-facies model (e.g. Fig. 2c) are given by and The probability p 1|12 has already been calculated (Eq. 23), and p 1|13 and p 2|23 are derived in exactly the same way. They are and The probabilities p 1|123 and p 2|123 are a little more complicated since the locations of beds of three different facies need to be considered. However, the same basic approach (Sect. 4.2) is applied to give and V 1 e (n 1 +n 2 +n 3 )T 1 − 1 n 1 + n 2 + n 3 − e (n 1 +n 2 )T 1 − 1 n 1 + n 2 − e (n 2 +n 3 )T 1 − 1 n 2 + n 3 + e n 2 T 1 − 1 n 2 + e (n 2 +n 3 )T 2 − e (n 2 +n 3 )T 1 n 2 + n 3 − e n 2 T 2 − e n 2 T 1 n 2 . (40)

Bed Base Types
When three rather than two facies are considered, the type and origin of possible bed bases increases significantly. Table 2 lists all possible outcomes of the two types of bed bases present in the stacked single-facies F1 models when they are merged to form the combined three-facies model. Similar tables have been constructed for F2 and F3 bed bases but are not shown since only one table is required to explain the general procedure followed. The probabilities of the different outcomes for the F11 and F10 bed bases are rationalised in Table 2 using two different terms. The first term is the probability that the location of a particular type of bed base in the single-facies Eroded p 10−|23 model is coincident with the presence of other specific facies in the other singlefacies models, and is termed p O A . The second term is the probability that a particular surface type results from the interactions of facies present, and is termed p O B . There are often several different combinations of situations that can lead to the formation of a particular surface type, and the total number of surfaces is calculated by summing these outcomes. For example, there are four ways of creating F13 surfaces (Table 2), and the total expected number of F13 surfaces is given by The shorter terms in this expression ( p 11|13 , p 11|123 , p 10|3 and p 10|23 ) are the overlap probabilities generalised as p O A , and the longer terms ( p 1113|13 , p 1113|123 , p 1013|3 and p 1013|23 ) are the outcome probabilities generalised as p O B (Table 2). These two types of probability ( p O A , and p O B ) are dealt with in the next two sub-sections.

Number of Bed Bases and Overlap Probabilities
There are two types of bed bases in the single-facies stacked models: bases that overlie another bed of the same facies (FXX) and bases that overlie shale (FX0). The expected numbers of these are given by and where n X (the number of beds) is given by Eq. 7. The overlap probabilities ( p O A ) define the probability that one of these bed bases in the single-facies model lies in the same position as one or more different facies in the other single-facies models. The are six probabilities for each facies (e.g. Table 2) which are straightforward to calculate and are given by the general equations and where X, Y and W refer to the three foreground facies, and 0 to the background one.

Surface Outcome Probabilities
Table 2 lists all possible surface outcome probabilities ( p O B ) associated with F1 beds, but not all need to be calculated, since some are associated with outcomes that are not included in the amalgamation ratio calculations (e.g. p 10−|23 , associated with eroding out a bed based), and others can be combined to simplify the equations. For example, the term p 1113|123 (present in Eq. 41 to derive n 13 ) can be combined with p 1113|123 into a different form which is easier to solve. Hence, the number of bed bases of different types in the combined facies models is calculated using the following nine equations.
These equations contain significantly fewer surface outcome probabilities than are listed in Table 2 for F1, and in the equivalent lists for F2 and F3. Many of those that remain define particular cases of a common general condition. For example, P 2222|12 and P 3333|23 both represent the probability that a thicker bed (an F2 or F3 bed, respectively) is preserved in the merged models when the amalgamation surface overlaps a bed from a thinner facies (F1 or F2, respectively). Hence, both probabilities can be expressed using the general term P B B B B|AB where A and B can be any of facies 1, 2 or 3, provided T B ≥ T A . Note that this term, and the other similar ones discussed below, use A and B as general facies identifiers rather than X and Y which are used elsewhere. This is to make clear that facies X and Y can be any of facies 1, 2 or 3, but facies A and B are constrained by the requirement that T B ≥ T A . The various surface outcome probabilities are calculated in the following sub-sections.

Outcomes for Amalgamated Bed Bases in the Presence of One Other Facies
The first set of conditions considered is when an amalgamated bed base in the singlefacies model overlaps with a bed of a different facies in one of the other single-facies models. There are four basic groups of situation for this, which have all been discussed already for the two-facies models in Sect. 4. Preservation of an amalgamation to a thinner bed in the presence of a thicker bed ( p A A A A|AB ) has been calculated for the specific case of p 1111|12 (Eq. 26), but also applies to p 1111|13 and p 2222|23 . The general form is Probabilities p 1112|12 , p 1113|23 and p 2223|23 refer to the case where a bed from a thinner facies beneath an amalgamation surface is replaced by a bed from a thicker facies. They are of the form given by Eq. 27 as Probabilities p 2221|12 , p 33331|13 and p 3332|23 refer to cases where a bed from a thicker facies beneath an amalgamation surface is replaced by a bed from a thinner facies. p 2221|12 is given in Eq. 31, and the general form is Finally, p B B B B|AB refers to the case where the thicker bed beneath an amalgamation surface is preserved in the presence of a bed from a thinner facies. Since it is impossible for a thinner bed to erode out the base of a thicker bed (i.e. p B B−|AB = 0), we have

Outcomes for Non-Amalgamated Bed Bases in the Presence of One Other Facies
Equations 47 to 55 contain two forms of outcome for the case where a nonamalgamated bed base in the single-facies model is converted to an amalgamated base in the presence of one other facies. Probability p 1012|2 was derived in the two-facies case (Eq. 29), and p 1013|3 and p 2023|3 are of the same form. The general expression for these is In general, p A0 AB|B < 1, since it is possible for the thicker FB bed to entirely erode out the thinner FA bed above the surface (i.e. p A0−|B > 0). The other form refers to the case resulting in a bed from a thicker facies overlying a thinner one ( p B0B A|A ). In this case, p B0−|A = 0, and therefore p B0B A|A = 1. (61)

Outcomes for Bed Bases in the Presence of Two Other Facies
The final eight surface outcome probabilities in Eqs. 47 to 55 are not of a form that has already been considered for the two-facies case, since they involve interactions between beds of all three facies. They can, however, be calculated using similar considerations to those applied for the two-facies cases (Figs. 6,7,Sect. 4.3), and therefore their derivations are not given, and only the final equations are reported here. Unlike those described above, these probabilities cannot be expressed in general forms since the relative thicknesses of all three beds are important. However, they can be discussed in three related groups. The first group considers the probabilities that unamalgamated bed bases in the single-facies model are eroded when they overlie locations occupied by beds in both other single-facies models. These probabilities are and Although p 10−|23 and p 20−|13 express similar conditions, the latter is a simpler equation since only one facies (F3) can erode the F20 surface, but the F10 surface can be eroded by both F2 and F3. Since neither F1 nor F2 can erode out an F30 surface, p 30−|12 = 0 and therefore does not appear in Eq. 54. The second group of probabilities define conditions where amalgamated bed bases in the single-facies models are preserved where they overlie locations occupied by beds in both other single-facies models. It turns out that for F1 and F2, these surface outcome probabilities are the same as the facies preservation probabilities given in Eqs. 39 and 40. Hence, and It follows that The final group of probabilities assesses whether amalgamated bed bases in the single-facies models are eroded out of the system. These turn out to be the same as the probabilities of erosion of the equivalent unamalgamated surfaces (Eqs. 62 and 63). Hence, p 33−|123 = 0, and

Verification of the Equations
Equations 1, 2 and 34-68 provide an analytical solution for the stacking statistics of models containing three foreground facies. The results are verified in the same way as for the two-facies results (Sect. 4.4), by conducting a high-resolution numerical experiment and comparing the stacking properties obtained in numerical models with those calculated using the equations. In this case, 100 realisations of 243 (i.e. 3 5 ) models are examined, comprising all possible combinations of three settings of five variables (Table 3). The models are higher resolution than those used for the two-facies cases, since when three facies are present, each bed can create a wider variety of surface types which sometimes have only a very low probability of occurring. Hence, model height (H ) is defined on the basis that at least 300 beds are expected for each facies. The median model (V 1 = V 2 = V 3 = 0.5, T 2 = 1.5m, T 3 = 2.25m, H = 974 m) is shown as an example (Fig. 10). Different positions are shown from this single realisation, highlighting the wide variability in stacking patterns and frequencies of bed base types which can be present within logs that sample only a few tens of each bed type. Even with samples of at least 300 beds of each facies, there is still quite a wide range of values obtained in the experiment, particularly for AR E (Fig. 11b, e, h), as reflected by the length of the error bars. However, the scatter is unbiased, and all parameters are estimated accurately by the equations. Therefore, we conclude from this comprehensive experiment that the analytical expressions for facies proportions and amalgamation ratios are robust.

Application of the Results
The principal results of this paper are the equations defining the stacking statistics of systems containing two or three foreground facies, presented in Sects. 4 and 5.
The reason for deriving these equations is to provide a firm analytical foundation for applying the compression method to facies models containing several foreground facies. The compression method allows object-based facies models to be built in which amalgamation ratio and facies proportions are defined independently and honoured in the output model (e.g. Walsh and Manzocchi 2021a, b). The compression method in general, and the specific results derived in this paper, are applicable to either high or low net:gross systems. However, the approach is particularly relevant for geological systems in which readily available facies modelling approaches are unable to create models with geologically realistic amalgamation ratios (Manzocchi et al. in press b). An example of such a system is deep-water turbidite deposits in which coupled deposition of shale drapes over sand lobes can result in systems with low AR at high NTG. The compression algorithm is a two-step process. In step 1 a conventional facies model is created with initial facies proportions and object thicknesses (P I X and T I X , respectively) set to values that will ensure the model has the target amalgamation ratios. In step 2, the thicknesses of model cells are expanded by facies-specific expansion factors (E X for the foreground facies and E 0 for the background facies) which transform the initial thicknesses and facies proportions to their target values without changing the amalgamation ratios. The result is a model with independently defined AR X , P X and T X values.
If a model contains only one foreground facies (F1), the facies proportion (P X ) is equal to the net:gross ratio (N T G) by definition. There is only one definition of AR This relationship, which is also derived in the present paper but expressed using different symbols (Eq. 6), provides the basis of the geometrical transformation underlying the compression algorithm when applied to a model containing a single foreground facies (Manzocchi et al. 2007;Walsh and Manzocchi 2021a).
The single-facies algorithm is parameterised using a compression factor (c F ) which represents the ratio between the amount by which the thickness of shale and sandstone cells must be modified during the geometrical transformation associated with the compression algorithm (Manzocchi et al. 2007). Hence, Fig. 11 Comparison between the analytically derived facies proportions and amalgamation ratios, and those obtained in a numerical model, for the three facies cases. Error bars represent +/-one standard deviation around the mean value obtained in 100 realisations. See text for discussion AR is linked to NTG through If c F = 1, then AR = N T G, which, as discussed, is the case in a constant-thickness object-based model. However, natural depositional systems often have c F << 1, with values of 0.05 to 0.3 typical (e.g. Manzocchi et al. 2007Manzocchi et al. , 2020Romans et al. 2009;Zhang et al. 2015;Pyles et al. 2019;Mueller et al. 2017;Soni et al. 2020). In these cases, application of the compression method using P I 1 = AR in step 1 would result in a geometrically transformed object-based model with the correct stacking characteristics.
The simple relationship (Eq. 69) arises because bed bases cannot be eroded in models of constant-thickness beds. Therefore, the equivalent relationship for multi-facies models applies only to the thickest facies present (F3 in the three-facies models), since the thinner F1 or F2 bed bases can be eroded, which will influence their amalgamation statistics. Since all F3 bed bases are preserved in the model, the probability that one of them overlies any one of F1, F2 or F3 is given by the overall net:gross ratio. This means that which has been verified by examining the results of the experiment described in Sect. 5.3. However, no similar, simple relationships exist for AR T 1 and AR T 2 . Therefore, although it is possible to create multi-facies compression-based models which honour user-defined input values of AR T X (see Manzocchi et al. in press b for an example with six facies distributed in pairs over two hierarchical levels), it is much simpler to use the alternative definition of amalgamation ratio ( AR E X ) when parameterising systems for use in compression-based modelling. This approach is developed below. Evaluation of stacking property values calculated from the equations for multifacies models shows that the effective amalgamation ratio ( AR E X ) is equal to an effective facies proportion (P E X ) which is defined on the basis of only the foreground facies in question and the background facies. Hence, where P E X = P X P X + P 0 .
This relationship is confirmed by the results of the two-and three-facies experiments (Fig. 12) and is compatible with the known single-facies relationship (Eq. 69), since in that case, P E X = P X = N T G.
Equation 73 provides the anchor needed to apply the compression method to multiple-facies models. In multiple-facies models, facies-specific compression factors (c F X ) are related to P E X and AR E X through the expression  (Table 4). In each case, AR E X has been calculated as a function Fig. 12 Cross-plots of AR E X versus P E X for a the two-facies models and b the three-facies models of P E X and c F X using Eq. 75: the combinations of P E X and AR E X values defining each facies are shown overlying the corresponding c F X curves in Fig. 13c. Although c F X is a very useful parameter for understanding the geometry of compression-based facies models and relating them to natural systems, it is not used directly in the compression algorithm which instead uses facies-specific expansion factors (E X , see Walsh and Manzocchi 2021a for details). The four E X values can be derived as a function of the three c F X and P X values (Table 4) by solving the equations and Once the expansion factors are determined, the thicknesses of each facies, and the facies proportions in the initial step 1 model (T I X and P I X ), can be determined from the expressions and P I X = P X /E X . Although the target thickness value for each facies is identical in all six cases, they all have different values of T I X (Table 4), illustrating the dependency of E X on the properties of all the facies in the model (Eqs. 76 to 79). Facies proportions are the user-defined input required by many object-based modelling software packages, in which case Eqs. 80 and 81 provide all the information Table 4 Parameters used in the six compression-based cases examined in These dimensions are arbitrary, and the systems are not based on any specific example needed to construct the step 1 model. However, the object-based models used in this paper are built from facies-independent volume fractions (V X ) rather than from merged facies proportions, and therefore, an extra step is required in which Eqs. 34 to 36 are used to determine the initial volume fractions (V I X ) as a function of the initial facies proportions and thicknesses (P I X and T I X ) given by Eqs. 80 and 81. This inversion is done numerically. Figure 13d compares P I X and AR E X values obtained in continuum onedimensional models (examples of which are shown in Fig. 13a) which have a similar resolution to those used previously (e.g. Fig. 10). The mean P I X and AR E X values match accurately the target values (the large circles) with similar variability in AR E X between realisations as observed previously (e.g. Fig. 11b, c, h). These onedimensional models therefore validate the approach described above (Eqs. 73 to 81).
The two-dimensional models (Fig. 13b) were constructed using the implementation of compression-based modelling discussed in detail by Manzocchi et al. (in press b). These models are more representative of the model resolution that might be used in practice, and differ from the one-dimensional models in at least three important ways. First, at 50 m thickness, the models are only a few times thicker than the thickest object present (T I 3 = 10 m in case 1; Table 4), and therefore any single realisation is unlikely to be representative. Second, unlike the one-dimensional models where objects are placed in a continuum, these two-dimensional models are gridded. The grid cells in these examples were scaled to be six times thinner than the thinnest objects in the initial model, which means that the models have a consistent resolution, but the presence of a grid will introduce a systematic bias in facies connectivity and hence amalgamation ratios (e.g. Koza et al. 2014). Third, the objects have tapered shapes and terminate laterally within the model area. This leads to the systematic bias in P X of compression models discussed by Walsh and Manzocchi (2021a), who proposed and applied a third step in the algorithm to correct for it (this step has not been applied here). These biases in P I X and AR E X account for the small differences between the target values and the mean values obtained in 100 realisations of the two-dimensional models, which nonetheless show a reasonably close correspondence between target and model stacking characteristics (Fig. 13e). Hence, the new expressions for amalgamation ratios and facies proportions derived in this paper provide the analytical foundations for applying the compression method to multi-facies object-based models.

Summary and Conclusions
This paper has developed equations for the expected facies proportions and amalgamation ratios of object-based models consisting of one, two or three foreground facies of constant-thickness beds embedded in an additional background facies. Expressions have been derived for the probabilities of facies preservation given a stratigraphically meaningful stacking of beds, and for the total numbers of bed bases which are unamalgamated or have different types of amalgamation surface. The results have been validated by comparing their predictions against the same statistics measured in multiple realisations of large, one-dimensional models built using the same assumptions that were made in the analytical solutions.
The objective of the work has been to provide a firm theoretical basis for applying the compression-based modelling approach to systems containing multiple foreground facies. The most important result for this is the demonstration that P E X = AR E X . This relationship underlies a relatively simple procedure for determining the parameters required to construct compression-based models including multiple facies parameterised interpedently. The procedure has been illustrated using a set of wellparameterised two-dimensional compression-based models containing a background and three foreground facies. Although the paper has only considered cases with one, two and three foreground facies, it is likely that P E X = AR E X irrespective of the number of facies present, and therefore that the procedure developed can be applied more generally.
• Total amalgamation ratio: AR T X . This applies to a foreground facies.
• Effective amalgamation ratio: AR E X . This applies to a foreground facies.
• Bed thickness (m): T X . This is specific for individual foreground facies.
• Bed thickness in the single-facies stacked model (m): t X . This is ≤ T X , since the bed may or may not be partially eroded. • Facies volume fraction: V X .
• Facies proportion in the combined stacked model: P X .
• Effective facies proportion in the combined stacked model: P E X .
• Facies transition probability: p XY . The probability that facies X transitions downwards into facies Y. • Facies preservation probability: p X |XY , p X |XY W . This defines the probability that the facies in the single-facies model is preserved in the combined stack model given an overlap with one or more other facies. The single character to the left of | indicates which facies is preserved. The two or three characters to the right of | are a list of the facies present at this location across the full set of single-facies stacked models. • Surface overlap probability. This defines the probability that a bed base in the single-facies model is in the same location as beds of other facies. The general term used for this is p O A , with specific probabilities having forms like p X X|XY W or p X 0|Y . The two characters to the left of | indicate a surface type in the singlefacies stacked model (FXX and FX0 in these examples), and the one, two or three characters to the right of | list the facies that are present in the other single-facies stacked models directly beneath this surface. • Surface outcome probability. This defines the probability of a specific outcome for a surface in the single-facies model, if a specific overlap is present. Outcomes include preservation of the surface, transformation of it into a different type of surface, or erosion of the surface out of the system. The general term used for this is p O B , and specific probabilities have forms like p X X XY |XY W or p X 0−|Y . There are either four characters to the left of |, or two characters followed by "-". The first two characters indicate a surface type in the single-facies stacked model. The third and fourth characters indicate the outcome surface type in the combined model. Replacement of these by "-" indicates that the surface is eroded out and does not exist in the combined model. The one, two or three characters to the right of | list the facies that are present in the single-facies stacked models directly beneath this surface. • Initial bed thickness in step 1 of the compression algorithm (m): T I X .
• Initial facies volume fraction in step 1 of the compression algorithm: V I X .
• Initial facies proportion in step 1 of the compression algorithm: P I X .
• Compression factor: c F X . This is facies-specific for all foreground facies. • Expansion factor: E X . This is facies-specific for all facies (foreground and background).