Analytical solution of folding behaviors of multi-layer viscous strata

This study extends the analytical solution of folding from a single-layer system to a multi-layer system composed of layers with varying thicknesses, viscosities and matrices, based on the stress distribution of folding with prescribed amplitude and wavelength. The parameters and factors affecting fold behavior were normalized to dimensionless quantities. In addition to single-layer systems, the influences of individual thickness, spacing and viscosity of multiple layers should be of concern. The range of matrix thickness where group folding occurs was initially identified and found to be smaller than one wavelength. Two-dimensional finite element folding simulations were also conducted to confirm the validity of the developed solutions of multi-layer folding. The 2D numerical simulation reveals when the matrix is thinner than one wavelength, the stress distribution within the matrix becomes linear and group folding occurs, and, if the matrix is thicker than one wavelength, the folding of the multi-layer system in fact is identical to independent folding of each layer alone. By means of stress superposition, this study developed theoretical solutions for multi-layer systems. The solutions developed show that the folding of viscous multi-layer and viscous matrix systems is time independent. Based on the developed solution, the effect of thickness, spacing and viscosity are then presented with the other parameters fixed. Finally, the solutions for elastic multi-layer and elastic matrix systems are also proposed. Article history: Received 5 July 2017 Revised 29 November 2017 Accepted 3 December 2017

It is said that the deformation of viscous material can be a time-dependent behavior.Biot (1961), Currie et al. (1962), and Jeng and Huang (2008) revealed that the deformation behavior of a viscous single-layer buckling fold embedded in a viscous matrix (V-V combination) and an elastic single-layer buckling fold embedded in an elastic matrix (E-E combination) are both time independent based on a theoretical solution.In addition to that, the theoretical folding solutions were verified by means of two-dimensional finite element method (2D FEM) numerical simulations.The geometry of the analysis model and boundary constraints are similar to those adopted by previous research (Zhang et al. 1996(Zhang et al. , 2000;;Mancktelow 1999;Jeng et al. 2002;Huang et al. 2010).The competent layer enclosed by matrix on both sides, and a "roller" constraint.If a perfectly straight layer without any perturbation is compressed, the layer is shortened without buckling.A perturbation is therefore required to induce folding in a numerical analysis and some initial geometric configuration different from a perfect-straight layer was generally adopted for this purpose (Zhang et al. 1996(Zhang et al. , 2000;;Mancktelow 1999;Jeng et al. 2002;Huang et al. 2010).The adopted perturbation, the end-rotational method, involves the imposition of a boundary rotation with an angle, small in magnitude, at one end of the competent layer.By employing this end-rotation, the shortened layer is applied a very small moment, which accounts for the subsequent buckle folding.The perturbation will be applied to the boundary at the beginning of the analysis, the system is then shortened to the pre-specified strain in one increment within the first step, and finally the system buckles, owing to this shortening and boundary perturbation, in later increments and steps.
Though single-layer folds is admittedly common; multi-layer folds by no means complex phenomena both in field and in reality.Therefore, the primary goal of this study is to explore theoretical folding solutions extending from a single-layer to a multi-layer system along with the associated folding behaviors.To better characterize the folding behavior of a multi-layer system compared with a singlelayer system, the following fundamental factors should be taken into consideration: (1) the spacing between the layers of a group of folds; (2) the stiffness contrast of each competent layer with regard to the matrix; and (3) the more complex compositions of different competent layers such as the thickness, spacing and viscosity.
To decipher the folding behavior of a viscous multilayer system, the matrix repulsion caused by the folding is initially studied to obtain the governing equation for the multi-layer strata.Once the theoretical solution of the governing equation is determined, the relationship and the effects of aforementioned factors are then discussed, e.g., layer thickness, spacing, viscosity, and competence contrast.Subsequently, 2D FEM-based numerical simulations are conducted in order to verify the validity of the theoretical solutions.

MATRIX REPULSION OF MULTI-LAYER FOLD SYSTEM
Biot (1961), Currie et al. (1962), and Jeng and Huang (2008) showed that the matrix repulsion caused by competent layer folding is an essential component to constructing the governing equation of the fold and also plays a major role in rendering the theoretical solution.In this study, the repulsion interaction between the matrix and the competent layers in the buckled fold is initially studied (Fig. 1) and is accordingly used to construct the governing equation.The obtained governing equation is then solved to determine the folding behavior of multi-layer folds.

Stress Distribution of Sine Wave Deformation in Infinite Field
The viscous matrix repulsion caused by competent layer sine-wave folding in a half-space infinite field can be written as (Biot 1961) where is the frequency of the fold and Y = Y(x, t) is the waveform function that describes the outlook of the deformed stratum in the y-direction.The subscript s represents a single-layer system.All the symbols are listed and defined in Table 1.
In the case of single-layer fold, the competent layer is embraced by two half-infinite matrices.The stress caused by matrix repulsion in a distance y regarding the competent layer can be written as (Biot 1937) Figure 2 illustrates the relationship between the normalized stress ( q ) versus the normalized distance (y/L), where L is the wavelength of the fold based on Eq. (2).However, by using a two-dimensional FEM-based numerical simulation, the stress field of the matrix is also obtained and found to be concordant with the prediction of Eq. (2) as depicted in Fig. 2. It should be noted that when y/L approaches one, i.e., the distance to the fold equals one wavelength, the strain/ stress caused by folding approaches zero (Fig. 2).

Stress Field Within Matrix of a Multi-Layer Folding System
When the competent layers are close to each other, i.e., within the range of one wavelength, they are mutually affected, resulting in group folding such that all of the layers in the system have the same wavelength and form a harmonic fold.In this situation, the frequency L 2 m r = can be imposed onto Eq. ( 2), where L then becomes the wavelength of the harmonic fold.Using Eq. ( 2) and employing superposition of the stress fields of two adjacent layers, the normalized repulsion ( q ) acting on the surface of competent layer is then related to S/L as Adequateness of this superposing has been verified using numerical analyses, and is to be presented in a later section.By concurrently using a 2D FEM simulation, we establish models comprised of six competent layers with viscosities of 2 × 10 20 Pa•s and 2 × 10 17 Pa•s along with a contrast ratio (R v ) of 1000 for the competent layer and the matrix, respectively, where the R v expresses the viscous compe- versus S/L obtained from both Eq.(3) and the numerical simulation and indicates good agreement.Equation (3) offers a good prediction of q , , y m o s v versus the S/L relationship notably, even when the number of layers is finite.
The stress distribution within the layers can actually be nonlinear, and a direct superposition of the stress intuitively yielded by different layers may not necessarily be applicable.We must examine why and how the two approaches render the almost identical results shown in Fig. 3. Figure 4 illustrates the superposed stress field in a matrix where the distribution of S/L varies from 2.0 -0.4, as obtained from Eq. (3).The stress at the middle of the matrix vanishes and shows its maximum magnitude on the two surfaces of the competent layer.As indicated in Fig. 3, when the thickness of the matrix S approaches 1.0 L or greater, the repulsion stress q o, m approaches that of a single-layer system q o, s .Specifically, when S ≥ L, the folding behavior of the multi-layer system converges to an individual single-layer system such that group folding (or harmonic folding) can only occur when S < L. Furthermore, when S < L, the stress inside the matrix exhibits a linear-like distribution along the y-direction and vanishes at the center of the matrix (Fig. 4b).The predictions of stress distribution shown in Fig. 4 have been compared to the results of the 2D FEM simulation and were found to have good agreement.Although the stress distribution along the y-direction can be either nonlinear (Fig. 4a) or linear (Fig. 4b), the good agreement suggests that direct superposition is admissible and attributable to the linear nature of the material itself.This comparison thus assures the validity of Eq. ( 3), which will then be adopted to construct the governing equation for group folding in a later section.
Consequently, when the matrix thickness exceeds one wavelength, a multi-layer system can be simplified as a single-layer competent folding; this concurs with the onewavelength zone of contact concept proposed by Ramsay and Huber (1987) from field observations and regressions.Furthermore, when the spacing is less than one wavelength, the repulsion within the matrix rapidly decreases with its thickness and eventually vanishes to zero with a matrix thickness of zero [Fig.3 and Eq. ( 3)].

Theoretical Solutions
In considering a multi-layer system having equally spaced competent layers with the same viscosity and thickness, an induced repulsion can be correlated with the thickness of the matrix (S/L) as shown by Eq. ( 3).
If such a system possesses infinite identical competent and matrix units, then the problem can be simplified by looking at one unit of the system.The basic unit is composed of Matrix repulsion on surface of competent layer, load applied by the matrix in y-direction q o , s Matrix repulsion (load) applied by the matrix in a single-layer system q o , m Matrix repulsion (load) applied by the matrix in a multi-layer system q i Shared load applied by the matrix to the i-th layer in y-direction  one competent layer and two half-thicknesses of the matrix (total matrix thickness being S) situated on both sides of the competent layer, as illustrated in Fig. 5.The entire pile of the multi-layer fold may be considered as a permutation and combination of the units.The deformation of each unit during the folding process is identical due to kinematic similarity.The surface repulsion for an individual competent layer (q o, m ) can be obtained by combining Eqs. ( 3) and (1) as follows: Similar to the process used by Biot (1961) and Jeng and Huang (2008), the relationship of q o, m with the lateral thrust (f, illustrated in Fig. 5) at the impending buckling can be established as where f is Knowing the basic waveform function of the competent layer fold to be ( , ) ()sin Y x t A t x x B f m = (Jeng and Huang 2008) and considering a constant lateral compression strain rate, then imposing Eqs. ( 4) and (6) onto Eq. ( 5) provides sin sin s in by reforming the equation, we obtain where x B f is the lateral accumulative compressive strain of the layer-matrix unit.With the application of the frequency L 2 m r = to Eq. ( 9), the theoretical solution of the viscous multi-layer fold is thus obtained as follows: We can then further normalize the wavelength with the thickness of the competent layer h.Normalized wavelength l = L/h and normalized matrix thickness s = S/h, are obtained as such, therefore s/l = S/L.Equation (10) then becomes The normalized dominant wavelength (l d, m ) can be obtained by taking l 0 where l 2 6 , o ds 3 r h h = is the dominant wavelength in a single-layer fold system (Biot 1961).
Finally, Eq. ( 12) can be reformed as the purely viscous multi-layer folding is thus time independent just as the purely viscous single-layer folding.
(2) As revealed by Eq. ( 12), the dominant wavelength of a multi-layer system (l d, m ) fold is related to either S/L or s/l in addition to the contrast ratio ( R v o h h = ).
(3) Using the viscosity contrast ratio R v = 1000 as an example, Fig. 6 illustrates the curves of the folding behavior according to different S/L values.It can be seen that when S/L = 1, the multi-layer folding behavior performs almost the same as a single-layer fold.By decreasing the matrix thickness (S), the required folding strain decreas-es and the dominant wavelength enlarges.Finally, when S/L = 0.1 the multi-layer folding behavior of approaches that of a matrix-free folding material.(4) In Fig. 7, the results from the numerical simulation agree well with the theoretical solution, as shown in Fig. 7.However, some minor discordance could be detected under the condition of small matrix thicknesses (S/L).Such discrepancy is possibly result from the fact that a numerical model can only construct a finite length, whereas a simulated layer can have infinite length.(5) Figure 7 further demonstrates the relationship of Fig. 6.The result of the example R v = 1000 by Eq. ( 10) which shows the influence of matrix thickness (S/L).The uppermost curve is the waveform curve of a single-layer system with infinite matrix and the lowest is the waveform curve of a single-layer system with no matrix around itself.When S/L = 1, the waveform curve of a multi-layer system is close to a single-layer system.When S/L is rather small, say 0.1, the multi-layer system folding is close to a "no matrix" folding situation.Fig. 7.The influence of matrix thickness on the resulting dominant wavelength.The prediction of the analytical solution [Eq. ( 13)] compares well with the results of 2D FEM simulation.A simplified curve is also proposed by Eq. ( 14).The subscripts d, m, and s denote "dominant wavelength", "multi-layer system", and "single-layer system", respectively.It shows that closer competent layers lead to greater dominant wavelengths.l d, m /l d, s vs. s/l d, m by means of Eq. ( 13), in which the influence of R v vanishes.In the condition of s/l d, m > 0.4, the values of l d, m /l d, s are all close to one.In other words, when the matrix thickness against the wavelength exceeds 0.4, the wavelength of the multi-layer fold is actually similar to the single-layer dominant wavelength; within the condition of gradually diminishing of matrix thicknesses (s/l d, m < 0.4), the wavelength of the multilayer fold increases rapidly.(6) If the competent layers have the same thickness and spacing, then the dominant wavelength of the multilayer fold (l d, m ) is directly correlated with the contrast of viscosity between the competent layer and the matrix [Eqs.( 12) and ( 13)].Equation ( 13) may be simply regressed from Fig. 7 This regression function is suitable for use under the condition of s/l d, m < 1.0, as indicated by Fig. 7.If the wavelength measured in this field is close to L d, m and together with the measured S and h (for the uniform thickness multilayer system), then l d, s can be more easily accessed using Eq. ( 14) and ( ) R v o h h = can then be inferred.

Multi-Layer System with Different Viscosities
Consider the multi-layer system shown in Fig. 8, where the viscosity of the i-th competent layer is i h and all of the layers have same thickness of h, adjacent with a thin matrix between the layers.Recall from Eq. (3) and Fig. 3 that the repulsion of the matrix rapidly decreases with its thickness and eventually vanishes to zero with a matrix thickness of zero.We therefore set the matrix thickness to be very thin (close to zero) so that there is no extra repulsion inside the grouped layers and repulsion only exists on the two outer-most surfaces of the grouped layer.
All of the strata equally deform with a constant lateral strain rate x f o so that the compressive stress in the i-th layer fold should be f 4 o and the matrix repulsion distributed for the i-th layer is q i .The governing equation in the i layer is obtained as The total repulsion can be obtained by summing the q i values: The viscous matrix repulsion, located on the two outermost surfaces of the group layer, is q t Y 4 and supposing that the deformation of the strata is compressed under a constant lateral strain rate, we then have Reformatting and simplifying the equation gives us In applying L 2 m r = to Eq. ( 19), the theoretical solution of an assemblage of this type of geometry and physical property multi-layer strata should be Applying the normalized wavelength l = L/h into Eq.( 20), the theoretical solution is obtained as As such, the folding behavior of multi-layer system is related to the matrix viscosity o h and the summation of the layer viscosity  21).
In the condition of L 2 m r = , the normalized dominant wavelength of the multi-layer strata l d, m is Considering the entire set of multi-layer folding strata as a single-layer (group layer) fold, the higher the value of equivalent viscosity ( ) will be for the system.Based on Eqs. ( 21) and ( 22), it is thus shown that for a pile of strata with different i h values, the folding behavior is not affected by the strain rate and again remains time independent.
As derived from Eq. ( 21), with the increment of total viscosity of the strata pile ( to initiate the buckle fold decreases.Similarly, as hinted in Eq. ( 22), the normalized dominant wavelength (l d, m ) then increases.

Folding Behavior Analysis
To explore the effect of total viscosity on a multi-layer fold, the model is evaluated by three different sets of viscosities in the competent layers.Each set is composed of five layers with a constant thickness and an identical matrix The equivalent total viscosity of these three sets therefore are 4 × 10 19 , 1 × 10 20 , and 4 × 10 20 (Pa•s), respectively.Figure 9 illustrates three different folding curves obtained from these three different sets of viscosity according to Eq. ( 22).
If the viscosity of each competent layer is identical, then Eq. ( 22) converges into l n 2 6 , where n is the number of layers.This result is identical to the analytical theory of multi-layer folds as proposed by Biot (1965).

Multi-Layer System with Different Thicknesses
Due to environmental perturbation and variation, it is unnecessary for all competent layers to form with identical thicknesses.Thus, a pile of strata with different thicknesses will appear.In a single-layer fold, the wavelength proportionally increases with the increment of the thickness of competent layer.After normalizing the wavelength with the layer thickness, it has been shown that the normalized wavelength only relates to the contrast ratio of strata; the effect of the layer thickness vanishes (Jeng and Huang 2008).However, the thickness of each individual layer may differ in a multi-layer fold.The effect of varying thicknesses associated with multi-layer folding must be taken into consideration.This section targets only at the evaluation of the thickness factors of the competent layers.The competent and matrix layers with identical viscosities and thin matrices are presumed on this premise.
The analyzed model is shown in Fig. 10, with the thickness of the i-th layer assigned as h i therein, and its associated inertia moment therefore be I h The strata are laterally deformed with an identical lateral strain rate such that all strata receive the same accumulated lateral compressive strain.
The resulting governing equation then takes the form of: The viscous matrix repulsion is q t Y 4 . Assuming deformation of the strata is compressed in a constant lateral strain rate with the lateral compression force of 4 o and substituting the basic waveform ( , ) ()sin This can be reorganized and simplified to to Eq. ( 26), the analytical solution of thickness-varying multi-layer strata should be Under the condition of l 0 Because the thickness of all layers is no longer the same, we do not have l d, m for this case.
Based on Eq. ( 27), it has been shown that the folding behavior of the thickness-varying viscous multi-layer strata is not only affected by the resistance contrast [compared to Eqs. ( 12) and ( 22)] but is also correlated to the term of h i thickness.The thickness of each bed is set at 0.5, 1, and 1.5 m, and the total thickness is set at 2.5, 5, and 7.5 m.Under this situation, the total thicknesses ( h i i n 1 = / ) and the sums of the cubic power of thickness ( h i i n 3 1 = / ) of the multi-layer folds are different.Figure 11 illustrates the curves of the three corresponding fold behaviors, where the wavelength increases with total thickness.Based on Eq. ( 28), a thicker stratum creates a greater inertia moment in the buckle multilayer strata (related with h i i n 3 1 = / ), and the wavelength of the fold is thus longer.The result of the theoretical solution fits well with the numerical simulation, as shown in Fig. 11.
In the case of identical bed thicknesses, the wavelength can be reformed as L h n 2 6 where n is the number of competent layers), which again is identical with the solution proposed by Biot (1965).

Same Bulk Thickness by Varying Individual Bed Thickness
Equation ( 28) indicates that the wavelength of a harmonic fold related with h i i n 3 1 = / .In order to evaluate this ef- fect, three models are thus set in a condition of maintaining the same bulk thickness yet with five varied bed thicknesses.The three models are as follows: Fold 1: 1, 1, 1, 1, 1; Fold  2: 0 .75,0.75,2,0.75,0.75;Fold 3: 0.5,0.5,3,0.5,0.5 (m), with an identical bulk thickness of 5 m.The maximum thicknesses of the three sets are 1, 2, and 3 m, respectively.
As illustrated in Fig. 12, the result demonstrates that the greater the value of h i i n 3 1 = / , the longer the wavelength and the greater the folding-required strain prior to fold initiation.Even though the total thickness is identical, the total inertia moment increases with the maximum thickness of the individual bed.Furthermore, the strata require more lateral compression x B f in order to generate sufficient strain to buckle the strata, as shown in Fig. 12.
In summary, given the same bulk thickness, a thicker maximum bed thickness results in a greater inertia moment of the strata, and the folding-required lateral compression and the wavelength are thus larger.These results have been verified with numerical simulation and compare well with each other, as indicated in Fig. 12.

Competent Layers with Different Viscosities and Thicknesses
Environmental nature would vary with time resulting in the changes of the rock unit in viscosity and thickness.This section seeks to develop the theoretical solution in general condition with a system, whereby the layers have different viscosities and thicknesses with a thin matrix in between.Other conditions, including fixed matrix and very thin matrix viscosities, are chosen as before.
The respective viscosity and thickness of the i-th layer are i h and h i , and the resultant inertia moment then is I h The buckle strata are assumed to have same lateral compressive strain rate, so the lateral compressive stress in each layer is 4 , where the i-th layer shares the matrix repulsion by q i .The governing equation thus becomes The matrix repulsion is q with the lateral compressive stress of the competent layer 4 by imposing all of the parameters, we then have After reformation and simplification, then for Eq. ( 32); then the analytical solution of the varying physical property and different thickness of the multi-layer fold system is Based on this theoretical solution, the folding behavior of multi-layer strata is now correlated with the sum of the product of the viscosity with the cubic power of the layer thickness.Equation ( 34) serves as the general solution for competent layers with different thicknesses and viscosities.

APPLICATIONS -SOLUTIONS FOR E-E MULTI-LAYER SYSTEM
The governing equation for the folding of an elasticmatrix elastic-competent multi-layer system can be established and solved using processes similar to those previously mentioned.A model with the following conditions is solved as a demonstration: (1) Under plane strain conditions, the equivalent elastic modulus is (2) The matrix is very thin (S close to zero) and is identical for every layer.
(3) The thickness of each competent layer is identical, h.(4) The pile of the strata is composed of a competent bed and matrix with an infinite group of layers.(5) Using stress superposition, we can determine the repulsion in the form of Combining Eq. ( 35) with the lateral forces, the governing equation of the elastic fold can be constructed as This solution has a similar format compared with the case of purely viscous multi-layer fold.Likewise, the folding behavior is time independent.By changing the parameter set ( , ) o h h to be ( , ) E Eo in Eqs. ( 1) -( 34), it can be shown that elastic and viscous material folds demonstrate similar folding behaviors.

CONCLUSION
(1) A two-dimensional numerical fold simulation of a multilayer system was initially conducted to explore the actual stress distribution within a matrix.When the matrix is thinner than one wavelength, the stress distribution within the matrix becomes linear (Fig. 4) and group folding occurs.However, if the matrix is thicker than one wavelength, the folding of the multi-layer system in fact is identical to independent folding of each layer alone.In other words, group or harmonic fold is most likely to occur the condition that the matrix is "thin enough", which can be defined as "smaller than one wavelength".(2) By means of stress superposition, this study developed theoretical solutions for multi-layer systems in which the viscosity and the thickness of the layers were set to be identical or different.
(3) The theoretical solution allows direct study of the effect of the matrix thickness S in case of S/L value greater than one, the matrix repulsion is similar to a single-layer fold system, matching with the concept of one wavelength zone of contact proposed by Ramsay and Huber (1987).In the case of S/L being less than one, the matrix repulsion decreases with decreasing spacing (Fig. 3).
The wavelength of the multi-layer fold thus increases as a consequence.(4) In the case of competent layers with varying viscosity, the folding behavior correlates with the sum of viscosity   22)], the longer the wavelength (Fig. 9).This result is identical with the analytical solution proposed by Biot (1965).
(5) In the case of competent layers with varying thickness, the folding behavior correlates with the sum of the cubic thickness [ h i i n 3 1 = / , Eq. ( 28)].The higher the sum, the longer the wavelength.If the bulk thickness is identical, the folding behavior will be dominated by the thickest layer (Fig. 12).

Fig. 1 .
Fig.1.Schematic diagram and the mechanics of the multi-layered fold.(L: wavelength, S: spacing of the competent layer, equal to the thickness of the matrix, h: thickness of the competent layer, f: lateral compressive stress, xy: 2D coordinate system, q o : matrix repulsion applied to the layer.)

Fig. 2 .
Fig. 2. Stress distribution ( , y s v ) within matrix, caused by folding of a single competent layer embedded in infinite matrix.The subscript s denotes "single-layer system."Note the coincidence between the analytical solution and the numerical simulation.h = 2 × 10 20 Pa•s, o h = 2 × 10 17 Pa•s, contrast ratio R v = 1000.
Fig. 5. Schematic diagram shows the regular pile of multi-layered folds; the system may simplified by one rock unit.

19 )Fig. 8 .
Fig. 8. Schematic diagram illustrates the analyzed model: competent layers with the same thickness, but with different viscosities closely spaced.

Fig. 9 .
Fig. 9. Waveform curves yielded from models with different sets of viscosities.The analytical predictions compare well with 2D FEM simulation by checking the dominant wavelengths.
Fig. 10.Schematic diagram illustrates the analyzed model: competent layers with different thicknesses, but with same viscosity closely spaced.

Fig. 11 .
Fig. 11.The waveform curves of the models with different total thickness of the competent layers.The analytical predictions compare well with 2D FEM simulation by checking the dominant wavelengths.