2.1. Scaling and dimensionless model equations

We focus on the long-wave limit and, guided by the asymptotic reductions outlined by Buckmaster et al. (Reference Buckmaster, Nachman and Ting1975) and Balmforth & Hewitt (Reference Balmforth and Hewitt2013), we set

(2.8) \begin{equation} x = \frac {{\tilde x}}{{\mathcal L}}, \quad (z,{Z},H)=\frac {({\tilde z},{\tilde Z},{\tilde H})}{{\mathcal H}}, \quad t = \frac {{\mathcal V}{\tilde t}}{{\mathcal H}}, \quad u = \frac {{\tilde u}}{{\epsilon }{\mathcal V}},\quad w = \frac {{\tilde w}}{{\mathcal V}}, \end{equation}

(2.9) \begin{equation} \tau _{xx} = \frac {{\tilde \tau _{xx}}}{{\mathcal S}}, \quad \tau _{xz} = \frac {{\tilde \tau _{xz}}}{{\epsilon }{\mathcal S}}, \quad p = \frac {{\tilde p} - \rho _{\textit{i}} g ({\tilde Z}+{\frac {1}{2}}{\tilde H}-{\tilde z})}{{\mathcal S}}. \end{equation}

Here, $\mathcal H$ denotes the initial layer thickness and $\mathcal L$ denotes a much longer horizontal lengthscale, leading to a relatively small aspect ratio,

(2.10) \begin{equation} {\epsilon } = \frac {{\mathcal H}}{{\mathcal L}} \ll 1. \end{equation}

In contrast to an extensionally dominated flow, the scalings in (2.8) significantly enhance the vertical velocity by an order of magnitude over the in-plane velocity component, to suit the bending limit, and are appropriate for the strength of gravity considered below in (2.16). The characteristic vertical velocity and initially thickness consequently set the time scale. The stress scale is set by the first line of (2.4) to be

(2.11) \begin{equation} {\mathcal S} = K \left (\frac {{\epsilon }{\mathcal V}}{{\mathcal L}}\right )^{m}\!, \end{equation}

and the shear is an order of magnitude smaller because there is zero tangential stress at the free surface.

Substituting into the dimensional equations (2.1)–(2.4), we find the dimensionless governing equations for conservation of mass, conservation of momentum and the power-law constitutive relation are

(2.12) \begin{align} {\epsilon }^2 u_x + w_z &= 0, \\[-12pt]\nonumber \end{align}

(2.13) \begin{align} \tau _{xx,x} - p_x + \tau _{xz,z} &= {\epsilon }^2{G}(1-\delta ) \left({Z} + {\frac {1}{2}} H\right)_{x}\!, \\[-12pt]\nonumber \end{align}

(2.14) \begin{align} {\epsilon }^2 \tau _{xz,x} -\tau _{xx,z} - p_z &= 0 , \\[-12pt]\nonumber \end{align}

(2.15) \begin{align} \begin{pmatrix} \tau _{xx} \\ {\epsilon }\tau _{xz} \end{pmatrix} = \big[4 u_x^2 + {\epsilon }^{-2} (u_z &+ w_x)^2\big]^{\frac {m-1}{2}} \begin{pmatrix} 2 u_x \\ {\epsilon }^{-1}(u_z + w_x) \end{pmatrix} \!, \end{align}

where we have introduced dimensionless parameters

(2.16) \begin{equation} {G} = \frac {\rho _{\textit{i}} g{\mathcal L}^2}{{\mathcal H}{\mathcal S}(1-\delta )} \quad \textrm {and} \quad \delta = 1 - \frac {\rho _{\textit{i}}}{\rho _{\textit{w}}} \gt 0 . \end{equation}

The ratio $G$ quantifies the strength of buoyancy relative to viscous stress and the choice ${G} = O(1)$ ultimately recovers the regime ${\tilde u} = O(\epsilon {\tilde w})$ assumed in our non-dimensionalisation. This is in contrast to regimes where compressive and restorative forces balance (e.g. Pfingstag et al. Reference Pfingstag, Audoly and Boudaoud2011; Ryan et al. Reference Ryan, Breward, Griffiths and Howell2025). The ratio $\delta$ quantifies the density difference between the two fluids relative to the density of the underlying fluid. We also note that the vertical buoyancy force appearing in the dimensionless force balance (2.3) has been absorbed into the pressure in the scaling (2.9).

The dimensionless boundary conditions on the stresses come from (2.5) and (2.6), and can be written as

(2.17) \begin{equation} \left . \begin{aligned} \tau _{xz} &= (\tau _{xx}-p) \left({Z}+{\frac {1}{2}} H\right)_x \\[5pt] \tau _{xx} + p &= - {\epsilon }^2 \tau _{xz} \left({Z}+{\frac {1}{2}} H\right)_x \end{aligned} \right \} \quad \textrm {on} \quad z = {Z} + {\frac {1}{2}} H \end{equation}

and

(2.18) \begin{equation} \left . \begin{aligned} \tau _{xz} &= \left [\tau _{xx}-p - {\epsilon }^2 {G} \left({Z}+{\frac {1}{2}} H - \delta H\right) \right ] \left({Z}-{\frac {1}{2}} H\right)_x \\[5pt] \tau _{xx} &+ p + {\epsilon }^2 {G} \left({Z}+{\frac {1}{2}} H - \delta H \right) = - {\epsilon }^2 \tau _{xz} \left({Z}-{\frac {1}{2}} H \right)_x \end{aligned} \right \} \quad \textrm {on} \quad z={Z}-{\frac {1}{2}} H . \end{equation}

Finally, the kinematic conditions (2.7) become

(2.19) \begin{equation} \left({Z} \pm {\frac {1}{2}} H\right)_t + {\epsilon }^2 u \left({Z} \pm {\frac {1}{2}} H \right)_x = w \quad \textrm {on} \quad z={Z} \pm {\frac {1}{2}} H. \end{equation}

2.2. Thickness-integrated equations and thin-sheet limit

We follow Howell (Reference Howell1996), Ryan et al. (Reference Ryan, Breward, Griffiths and Howell2025) and others by integrating the governing equations across the sheet thickness to derive effective bulk equations. We define a tension $\mathcal T$ and stress resultant $\mathcal Q$ by

(2.20) \begin{equation} {\mathcal T} = \int _{{Z}-H/2}^{{Z}+H/2} \left (-p+\tau _{xx}\right )\ \mathrm{d} z, \quad {\mathcal Q} = \int _{{Z}-H/2}^{{Z}+H/2} \tau _{xz}\ \textrm {d} z. \end{equation}

Integrating (2.13)–(2.14) from $z = Z-H/2$ to $z=Z+H/2$ and applying the corresponding boundary conditions (2.17)–(2.19) then yields

(2.21) \begin{equation} \frac {\partial H}{\partial t} = O({\epsilon }^2),\quad \frac {\partial {\mathcal T}}{\partial x} = O({\epsilon }^2) , \quad \frac {\partial {\mathcal Q}}{\partial x} - {G} Z = 0 . \end{equation}

In view of the first relation in (2.21), we discard thickness variations and so surface gradients can be simplified to $({Z}\pm (1/2)H)_x \sim Z_x$ . The vertical force balance (2.13) also demands that $-p+\tau _{zz}$ is $O({\epsilon }^2)$ and, hence, the pressure is approximated by $p\sim \tau _{zz}=-\tau _{xx}$ .

Next, on multiplying the horizontal force balance (2.13) by $(z-Z)$ and again integrating across the sheet thickness, we find the moment equation,

(2.22) \begin{equation} {\mathcal Q} = {\mathcal M}_x + Z_x {\mathcal T} + O({\epsilon }^2), \end{equation}

where

(2.23) \begin{equation} {\mathcal M} = \int _{{Z}-H/2}^{{Z}+H/2} (\tau _{xx} - p)(z-Z) \; \textrm {d} z \approx 2\int _{{Z}-H/2}^{{Z}+H/2} \tau _{xx}(z- Z) \; \textrm {d} z \end{equation}

is the bending moment. Differentiating (2.22) with respect to $x$ and using (2.21) then yields

(2.24) \begin{equation} {\mathcal M}_{xx} + {\mathcal T} Z_{xx} + GZ = O(\epsilon ^2). \end{equation}

To achieve balance in the continuity equation (2.12) and constitutive relations (2.15), we see that the velocity field corresponds to bending without shear at leading order:

(2.25) \begin{equation} u \sim U(x,t) - (z-Z) \frac {\partial W}{\partial x} \quad \textrm {and} \quad w \sim W(x,t) \sim \frac {\partial Z}{\partial t}, \end{equation}

where $U$ and $W$ are thickness-averaged velocity components, to be determined as part of the solution. With the relatively small shear rate and stress, the constitutive law (2.15) simplifies to

(2.26) \begin{equation} \tau _{xx} \sim 2^m |u_x|^m \textrm {sgn}(u_x) , \end{equation}

which must be fed into the integrals (2.20) and (2.23) for $\mathcal T$ and $\mathcal M$ to complete the model.

2.3. Viscous model

For the special case $m=1$ , after applying velocity and constitutive equations (2.25) and (2.26), performing the integrals (2.20) and (2.23) defining $\mathcal M$ and $\mathcal T$ , and substituting into the force balances (2.21) and (2.24), we arrive at the model

(2.27) \begin{align} {\mathcal T}(t) &= 4(U_{x} + Z_{x} W_{x})= {4(U_{x} + Z_{x} Z_{{x} t})}, \\[-12pt]\nonumber \end{align}

(2.28) \begin{align} {\mathcal M}_{{x}{x}} &= - {\frac {1}{3}} W_{{x}{x}{x}{x}} = {- {\frac {1}{3}} Z_{{x} {x} {x} {x} t}} = G Z - {\mathcal T} Z_{{x}{x}} \end{align}

(cf. Buckmaster et al. Reference Buckmaster, Nachman and Ting1975; Ribe Reference Ribe2002; Pfingstag et al. Reference Pfingstag, Audoly and Boudaoud2011). We interpret (2.27) and (2.28) as a horizontal and vertical force balance, respectively. In (2.27), the horizontal and vertical motions on the right-hand side balance locally to give a spatially uniform but time-evolving horizontal tension. In (2.28), the left-hand side represents a viscous bending moment arising from the gradient of in-plane velocity in (2.25); on the right-hand side, the first term represents the stabilising effect of gravity, whereas the second term captures the effect of the thickness-averaged in-plane stress, which may be stabilising or destabilising depending on the sign of $\mathcal T$ . In general, (2.27) and (2.28) require two boundary conditions on $U$ and four boundary conditions plus an initial condition on $Z$ .

If the mean straining motion is prescribed to be

(2.29) \begin{equation} \langle U_x \rangle = \varDelta , \end{equation}

where $\langle \ldots \rangle$ denotes an average over $x$ , we then have

(2.30) \begin{equation} {\mathcal T}(t) = 4\varDelta + 2\frac {\textrm {d}}{\textrm {d} t}\big\langle Z_{x}^2 \big\rangle . \end{equation}

Equation (2.28) is now an evolution equation for the centreline $Z(x,t)$ subject to a horizontal stress (2.30). Note that nonlinearity appears through the final term on the right of (2.30), which corresponds to a global constraint on the net rate of compression. Crucially, this term takes the form of a spatial integral, so it provides a non-local constraint on the total length increase of the centreline, and transforms into a nonlinear coupling between modes under Fourier analysis. We assume that $\varDelta$ can be controlled externally, for example by imposing a relative velocity between two points in the sheet. It is possible for $\varDelta$ to be a function of time, but in this paper, we only consider the case where it is fixed and set it equal to -1 by a suitable choice of the arbitrary velocity scale in (2.8).

In §§ 3 and 4, we use the model in (2.28) and (2.30) to explore the nonlinear dynamics of viscous buckling; that is, for the special case $m=1$ . The version of the model for a power-law fluid is discussed further in § 5. We perform this analysis in spatially periodic settings for exposition, although it can be easily adapted for other scenarios. In this setting, we choose the lengthscale $\mathcal L$ so that the domain has dimensionless length  $2\pi$ .

4.1. Amplitude equations

Substituting (4.1) into (2.28) and (2.30), we find that the phases $\vartheta _{\!j}$ again remain constant, while the amplitudes $a_{\!j}$ satisfy

(4.2) \begin{equation} \frac {j^4}{3} \dot {a}_{\!j} +{\mathcal T} j^2 a_{\!j} + G a_{\!j} = 0, \end{equation}

with a corresponding net tension

(4.3) \begin{equation} {\mathcal T} = 4\varDelta + 4 \frac {\textrm {d}}{\textrm {d} t}\sum _{j=1}^N j^2 a_{\!j}^2. \end{equation}

This can be rewritten more concisely as

(4.4) \begin{equation} \sum _{l=1}^N m_{\!jl} \dot {a}_l = \gamma _{\!j} a_{\!j}, \end{equation}

where the linear growth rates and mass matrix respectively are

(4.5) \begin{equation} \gamma _{\!j} = - \frac {3}{j^4} \big ( 4j^2 \varDelta + {G} \big ) \quad \textrm {and} \quad m_{\!jl} = \delta _{jl} + \frac {24 l^2}{j^2}a_{\!j}a_l . \end{equation}

For $a_{\!j} \neq 0$ , the evolution (4.2) can be rearranged into

(4.6) \begin{equation} \frac {\textrm {d}}{\textrm {d} t} \log a_{\!j} = - 3j^{-4} ({G} + j^2 {\mathcal T}\,). \end{equation}

We observe that (4.6) can be integrated for $a_{\!j}(t)$ if the integral of ${\mathcal T}(t)$ is known; this integral encodes the nonlinear coupling between modes. The system is therefore equivalent to the single ordinary differential equation (ODE)

(4.7) \begin{equation} {\mathcal T} = \dot {\varTheta } = 4 \varDelta + 4 \frac {\textrm {d}}{\textrm {d} t}\sum _{j=1}^N \big[j a_{\!j}(0) E_{\!j}(t)\big]^2 , \end{equation}

along with $\varTheta (0)=0$ and

(4.8) \begin{equation} a_{\!j}(t) = a_{\!j}(0) E_{\!j}(t), \quad E_{\!j}(t) = \exp \big [-3j^{-4}\textit{Gt} - 3 j^{-2}\varTheta (t)\big ] . \end{equation}

Integrating (4.7) now furnishes the implicit relation,

(4.9) \begin{equation} \varTheta (t) - 4 \sum _{j=1}^N \big[j a_{\!j}(0)\big]^2 \big[E_{\!j}(t)^2-1\big] = 4 \Delta t . \end{equation}

The quantity $\varTheta (t)$ can be interpreted as measuring the discrepancy between length added to the system by straining and length accommodated by wrinkles.

4.2. Viscous wrinkling cascade

Sample numerical experiments with (4.3)–(4.4) or (4.7)–(4.8) are shown in figure 2 (we use the system in (4.3)–(4.4) for relatively low numbers of modes $N\leq 20$ , integrating the ODEs in time using MATLAB’s ode23; for larger $N$ , we use the system in (4.7)–(4.8) and employ MATLAB’s ode45). In panels (a)–(c), the modes begin from the initial conditions,

(4.10) \begin{equation} a_{\!j}(0)=a_0 j^{-p}, \end{equation}

with $a_0=10^{-5}$ and $p=2$ ; in panels (d)–(f), the initial amplitudes are Gaussian random variables with a standard deviation of $a_0$ . In each case, we use $G=5$ , leading to the growth rates displayed in the inset of figure 2(c), and mode sets with $N=5$ , 10 and 20 are displayed. The solutions display interesting cascade dynamics in which the lower wavenumber modes with $j\gt 1$ (and the largest linear growth rates) appear first, but then become suppressed once higher wavenumber modes (with lower linear growth rates) amplify to a nonlinear level. When the cascade reaches the mode with highest wavenumber, that mode never becomes suppressed and continues to dictate the dynamics of the other modes via (4.6), (4.7), while itself growing algebraically as in (3.4). Before that final phase of evolution, the cascade appears to operate with a dependence of $t^{1/4}$ (the variable used for the horizontal axis).

Figure 2.

Numerical solutions of (4.4) for (a) $N=5$ , (b) $N=10$ and (c) $N=20$ , with initial conditions $a_{\!j}(0)=a_0j^{-2}$ ; $(\varDelta ,{G},a_0)=(-1,5,10^{-5})$ . Each of panels (a–f) contains three subpanels: time series of $a_{\!j}(t)$ ( $j=1,2,\ldots ,N$ , from blue to red), time series of ${\mathcal T}(t)$ (green), and a density plot of $\log _{10}|a_{\!j}|$ on the $(t,j)-$ plane, using $t^{1/4}$ for the $x\text{-}$ axes. In panels (d–f), similar solutions are presented, but using initial conditions in which the initial amplitudes are Gaussian random variables with a standard deviation of $a_0$ . The points and dashed lines show the times at which $a_{\!j}$ is maximum. The inset in panel (c) shows the wavenumbers and initial growth rates for the computations.

Figure 3.

(a) Time series of tension $\mathcal T$ and (b,c) scaled mode amplitudes $j^{-1/2}a_{\!j}(t)$ plotted against $t^{-1/4}j$ . The solutions from figures 2(b) and 2(c) with $N=10$ and 20 are shown in panel (b), and a solution with $N=200$ (from figure 4 a) is presented in panel (c). In panel (a), the tensions for all three are shown, as well as the solution with $N=5$ from figure 2(a); the dashed line in panel (a) shows the scaling $-5.9t^{-1/2}$ . In panels (b,c), modes are plotted from $j=3$ to the end of self-similar part of the cascade (a few modes short of $j=N$ ); $j$ increases as indicated by the arrow (from red to blue), and the cases $N=10$ and 20 are plotted in panel (b) as dashed and solid lines, respectively (although the modes with the same $j$ are mostly indistinguishable in the plot). The stars show the prediction from (4.17).

Figure 3. Long description

Panel A: A line graph shows the time series of tension plotted against time. The horizontal axis represents time (t) on a logarithmic scale, and the vertical axis represents tension (T(t)) on a logarithmic scale. The graph includes multiple lines representing different values of N (5, 10, 20) and a dashed line showing the scaling 5.9t^(-1/2). The lines for N = 5, N = 10, and N = 20 are shown in different colors, with N = 5 in blue, N = 10 in purple, and N = 20 in red. The tension values are negative and decay towards zero over time for all N values, with the dashed line providing a reference scaling. Panel B: A line graph shows scaled mode amplitudes plotted against t^-1/4j for N = 10 and N = 20. The horizontal represents t^-1/4j, and the vertical represents j^(-1/2)a_j(t). The graph includes multiple lines representing different values of j. The lines for N = 10 are shown as dashed lines, and the lines for N = 20 are shown as solid lines. The peaks of the lines shift to the right as j increases. Panel C: A line graph shows scaled mode amplitudes plotted against t^-1/4j for N = 200. The horizontal axis represents t^-1/4j, and the vertical axis represents j^(-1/2)a_j(t). The graph includes multiple lines representing different values of j, with j increasing from left to right. The lines for N = 200 are shown in red, with the peaks of the lines shifting to the right as j increases. Stars on the graphs indicate the prediction from equation (4.17).

Figure 4.

Solutions to (4.8)–(4.9) for $a_{\!j}(0)=j^{-p}\times 10^{-5}$ , with $N=200$ and (a) $p=2$ , (b) $p=0$ and (c) $p=-{1}/{2}$ ; in panel (d), $a_{\!j}(0)$ are chosen as Gaussian random variables with standard deviation $10^{-5}$ ( $G=5$ , $\varDelta =-1$ ). The top and middle rows show mode amplitudes $a_{\!j}(t)$ as a time series (lighter red lines) or density plot, respectively; the bottom row shows the time series of $\varPsi (t)=-\varTheta (t)/\sqrt {\textit{Gt}}$ . The thicker (blue) line in the top row displays the prediction in (4.24), and superposed on the density plots along the middle row is the line (4.23) (dashed red); both are calculated using the computed $\varPsi (t)$ . Following figure 2, the top two rows are plotted against $t^{1/4}$ ; the last row shows a logarithmic scale for time. The stars in panels (a,b) and the dotted line in panel (c) indicate the predictions from (4.29); the lighter lines and stars show results with $p=\{-{1}/{2},0,1,2,3\}$ . The lighter grey lines in panel (d) present twenty-five more realisations of the same randomly initialised computation.

Note that the net compression, $-{\mathcal T}(t)$ , falls each time a mode reaches its nonlinear stage and suppresses its predecessor, as shown in more detail in figure 3(a). For large times, the descending steps in $-{\mathcal T}(t)$ converge towards a scaling ${\mathcal T}\sim T t^{-1/2}$ for some constant $T$ . Similarly, for the more regular cascades in figures 2(b) and 2(c), we obtain a suggestive collapse of the solutions if we rescale each mode amplitude $a_{\!j}(t)$ by $j^{-1/2}$ and then plot the results against $t^{-1/4}j$ ; see figure 3(b). In this panel, the modes at the beginning and the end of the cascade are not shown (we plot modes from $j=3$ up to $j=N-2$ ). The rescalings collapse the mode amplitudes towards a common curve, suggesting that the cascade dynamics is self-similar, at least until the highest mode begins its final algebraic runaway.

The finer features of the cascade are illustrated further in figure 4, which shows computations with a rather larger number of modes ( $N=200$ ) and different sets of initial conditions. In particular, figure 4(a–c) displays computations using the initial conditions in (4.10) with $a_0=10^{-5}$ and $p=2$ , 0 and $-1/2$ ; figure 4(d) displays another case with a set of randomly initialised modes. The tension and rescaled mode amplitudes for the case with $p=2$ are also included in figure 3. Altogether, the various sets of computations with different initialisations demonstrate how the cascade is relatively insensitive to initial conditions. The random cases lead to less regular-looking cascades.

Similarly, the dynamics is not sensitive to varying the parameter $G$ (provided it is not set to zero): choosing other values leads to no significant differences in the qualitative features of the cascade. The only difference with lower values of $G$ concerns the first Fourier mode, which is damped for the value of $G$ used in figures 2, 3 and 4 (see the inset of figure 2 c). Regardless of nonlinear mode interactions, this mode always decays for $G=5$ . For sufficiently small $G$ , however, the first mode becomes unstable and can then grow in the initial stages of the cascade. Conversely, with higher values for $G$ , more of the lower modes become stable, demanding that the cascade begins at higher wavenumber.

Figure 5 shows a reconstruction of the deflection implied by the evolving mode amplitudes of figure 2. The figure presents the deflection as a density over the $(t,x)$ -plane and illustrates the spatio-temporal pattern of the buckling dynamics. The cascade in mode amplitude corresponds to a continued wrinkling of the buckling pattern. In particular, the transitions between modes with adjacent wavenumbers correspond to the spontaneous emergences of new wrinkles.

Figure 5.

Reconstruction of the deflection $Z(x,t)$ for the modal solutions shown in figure 2, with phases chosen randomly. The red lines show the loci of the nodes.

Figure 5. Long description

A heat map showing deflection patterns with nodes marked by red lines. The heat map represents the reconstruction of deflection for modal solutions with randomly chosen phases. The overall structure is a grid layout with six panels arranged in two rows and three columns. Each panel displays a heat map with a color scale indicating deflection values. The color scale ranges from -20 to 20, with darker colors representing lower values and lighter colors representing higher values. The horizontal axis is labeled with t to the power of 1/4, and the vertical axis is labeled with x. The red lines in each panel indicate the loci of the nodes. Panel (a) shows a heat map with fewer nodes. Panel (b) displays a more organized pattern with nodes forming distinct lines. Panel (c) shows a highly organized pattern with nodes forming parallel lines. Panel (d) has a similar pattern to panel (a). Panel (e) shows a pattern similar to panel (b) with nodes forming distinct lines. Panel (f) displays a highly organized pattern similar to panel (c) with nodes forming parallel lines. The overall trend across the panels shows an increase in the organization and regularity of the deflection patterns as the value of t to the power of 1/4 increases.

4.3. Two-mode interaction

To understand the cascade dynamics in more detail, we consider two interacting modes, labelled $j$ and $l$ , assuming all others have too low an amplitude to be relevant. Rewriting (4.4), we find the mode amplitudes satisfy

(4.11) \begin{equation} \begin{pmatrix} \dot a_{\!j} \\[3pt] \dot a_l \end{pmatrix} = \big(1+ 24 a_{\!j}^2 + 24 a_l^2 \big)^{-1} \begin{pmatrix} \big[\gamma _{\!j} - 72 {G} j^{-2} \big(j^{-2}-l^{-2}\big) a_l^2 \big]a_{\!j} \\[9pt] \big[\gamma _l + 72 {G} l^{-2} \big(j^{-2}-l^{-2}\big) a_{\!j}^2 \big]a_l\end{pmatrix}\!. \end{equation}

A phase portrait for this system, with $G=5$ and $(j,l)=(2,3)$ , is shown in figure 6. This phase portrait illustrates how mode $j$ grows first, but is then replaced by mode $l$ (when $j\lt l$ ), regardless of initial condition, provided both modes begin at relatively low amplitude. We rationalise this dynamics in more details as follows.

Figure 6.

A phase portrait for ${G}=5$ and $(j,l)=(2,3)$ , with time increasing along the arrows. Diagonal lines tracking up and right illustrate modes growing at comparable rates; horizontal lines tracking left to right illustrate $a_3$ growing while $a_2$ stays approximately fixed (because $|a_3| \ll |a_2|$ ), and vertical lines tracking from high to low illustrate $a_2$ decaying while $a_3$ maintains its size. The vertical dashed line indicates where $\dot a_{\!j}$ reverses sign.

If $j\lt l$ and $\gamma _{\!j}\gt \gamma _l$ , mode $j$ will first grow exponentially until it reaches the nonlinear state with

(4.12) \begin{equation} a_{\!j} \dot a_{\!j} \sim \frac {\gamma _{\!j}}{24}. \end{equation}

At that stage, the growth rate of the other mode, still of low amplitude, becomes modified so that

(4.13) \begin{equation} \dot a_l \sim 3{G} l^{-2} (j^{-2} - l^{-2}) a_l \gt 0. \end{equation}

That is, mode $l$ remains unstable, but the effective growth rate is reduced. There are therefore no barriers to the continued growth of mode $l$ until it, too, reaches the nonlinear regime.

Once $a_l$ becomes sufficiently large, that mode exerts its influence on mode $j$ . Indeed, it is clear from (4.11) that $\dot a_l$ is always positive (when $a_l\gt 0$ ), whereas $\dot a_{\!j}$ reverses sign when

(4.14) \begin{equation} \gamma _{\!j} \lt 72 {G} j^{-2} (j^{-2}-l^{-2}) a_l^2 . \end{equation}

Therefore, the growth of mode $l$ cannot be stopped and it must eventually take control over mode $j$ . Once mode $l$ achieves control and sufficiently suppresses the amplitude of mode $j$ , a new nonlinear state emerges with

(4.15) \begin{equation} a_l \dot a_l \sim \frac {\gamma _l}{24} \quad \textrm {and} \quad \dot a_{\!j} \sim - 3{G} j^{-2} (j^{-2} - l^{-2}) a_{\!j} \lt 0. \end{equation}

In other words, mode $j$ enters an exponential decline once controlled by mode $l$ .

Figure 7.

Mode amplitudes for $j=2,3,4$ (red to blue) for the numerical solution with $N=5$ also shown in figure 2(a). For $j=3$ and 4, the dotted lines show the initial, linear, exponential growth and the dots show the nonlinear phases with $a_{\!j}=\sqrt {\gamma _{\!j} t/12}$ ; the dots for $j=2$ plot (3.4). The dashed triangles labelled $\{j,l\}$ show the expected exponential trends when the growth of mode $j$ is nonlinearly controlled by mode $l$ .

In figure 7, the predictions in (4.13) and (4.15) are compared with the numerical solution for $N=5$ from figure 2. Note that the interactions exposed in (4.13) and (4.15) rely on the gravity term $G$ . Indeed, without this term, there is no cascade and the longest mode dominates, as found previously by Buckmaster et al. (Reference Buckmaster, Nachman and Ting1975). Also, when another mode $n$ enters the fray, with $\gamma _n\lt \gamma _l$ , this mode continues to grow (with a growth rate reduced by mode $l$ ) regardless of the interaction between $j$ and $l$ . The subsequent interaction between $n$ and $l$ then eventually switches control from $l$ to $n$ . This switch, and the ensuing decline of $l$ , does not, however, permit the re-emergence of mode $j$ , because mode $n$ , now in the nonlinear regime with $ a_n \dot a_n \sim \gamma _n/24$ , continues the suppression of $j$ in place of $l$ (with an even stronger damping rate).

Despite the success of the comparison in figure 7, the situation becomes murkier further down the cascade: it is clear from figures 2(c) and 2(f) that interactions for higher mode number involve more than two modes, and the nonlinear phase in (4.12) becomes less prominent. In fact, given the exponential growth predicted by (4.13), the threshold in (4.14) suggests that mode $a_{\!j}$ reaches its maximum for $t\sim O(j^5 \log j)$ when $j\gg 1$ . This estimate differs from the observation in figure 2 that the cascade operates with a time-dependence in which $j\sim t^{1/4}$ .

4.4. A self-similar wrinkling cascade?

The numerical computations in figures 2, 3 and 4 suggest that the cascade may be self-similar. To examine this possibility in more detail, we make a continuum approximation of the problem in which we replace the sum over the discrete wavenumber $j$ by an integral and set $a_{\!j}(t)=\alpha (j,t)$ . From (4.3) and (4.8), we then find

(4.16) \begin{equation} {\mathcal T} = 4\varDelta + 4 \frac {\textrm {d}}{\textrm {d} t} \int _1^{N} j^2 \alpha ^2 \,\textrm {d} j, \quad \alpha (j,t) = \alpha (j,0) e^{-3(\textit{Gt}j^{-4}+\varTheta j^{-2})} . \end{equation}

The structure of these equations, and in particular the factors $tj^{-4}$ and $\varTheta j^{-2}$ in the argument of the exponential, suggests a self-similar solution of the form

(4.17) \begin{equation} \begin{aligned} \alpha (j,t) &\sim j^{1/2} A_0 e^{-3(Gq^{-4}-2\varGamma q^{-2})} , \quad {\mathcal T} = - \frac {\varGamma }{\sqrt {t}}, \quad q = \frac {j}{t^{1/4}}, \\ -\varDelta & \sim A_0^2 \int _{q_{\textit{min}}}^{q_{\textit{max}}} q^3 e^{-6(Gq^{-4}-2\varGamma q^{-2})}\, \textrm {d} q +O(t^{-1/2}) , \end{aligned} \end{equation}

where $\varGamma$ and $A_0$ are constants, $(q_{\textit{min}},q_{\textit{max}})=t^{-1/4}(1,N)$ , and the factor $j^{1/2}$ in the amplitude of $\alpha$ is required for the integral equation in (4.16) to transform into the version in (4.17), in which the integrand is independent of $t\gg 1$ at leading order.

Figure 3 offers numerical evidence in support of the solution (4.17), with $\varGamma \approx 5.9$ chosen to fit the trend of the computed tension and $A_0$ fixed using (4.17 b) with $(q_{\textit{min}},q_{\textit{max}})=(0.1,3)$ . Nevertheless, the solution (4.17) corresponds to the specific initial condition, $a_{\!j}(0)=j^{1/2}A_0$ , which is not equivalent to those used in figures 2 and 4. Therefore, if these solutions do become truly self-similar, it can only be in a large-time asymptotic sense.

Two further observations obscure this interpretation, however. First, the cascade with $N=200$ and $p=2$ , shown in figures 3 and 4, displays a convergence towards a solution like that in (4.17), but also then slowly drifts beyond. This drift is particularly notable in the compensated plot of $\varPsi (t)=-\varTheta (t)/\sqrt {\textit{Gt}}$ shown in figure 4(a), which should converge to a constant if ${\mathcal T}\propto t^{-1/2}$ . The cascade beginning with (4.10) and $p=0$ in figure 4(b) suffers a similar fate, as do further solutions with other values of $p$ (shown by lighter lines in the figure); only for that with $p=-1/2$ in panel (c) does the net tension appear to approach the correct limiting behaviour.

A second observation regards the awkward normalisation of $\alpha$ by $j^{1/2}A_0$ in (4.17): the limits $(q_{\textit{min}},q_{\textit{max}})$ in the integral constraint are now time-independent, implying that $A_0$ cannot be constant, and taking the limit $q_{\textit{max}}\to \infty$ leaves the integral divergent. Only if one selects $q_{\textit{min}}\ll 1$ and $q_{\textit{max}}$ to lie over a range of order one values (which is approximately $[1,20]$ for the parameter settings of figure 3) does this integral converge to a finite result that is approximately independent of the value chosen. In fact, because $q_{max } = t^{-1/4}N$ , the normalisation implied by (4.17) must inevitably break down for later times. This breakdown is nothing more than the termination of the viscous cascade, arising when $q_{\textit{max}}=t^{-{1}/{4}}N$ reaches the range of $q$ where the integrand of (4.17 b) is $O(1)$ . That is, from figure 3, for $q\approx 6/5$ . This implies that the cascade terminates for $t^{{1}/{4}}\gtrsim 5N/6$ , in comfortable agreement with figure 2.

Overall, the lack of clear convergence to self-similar form and the issues with the normalisation integral in (4.17 b) both indicate that the cascade dynamics cannot be truly self similar. We expose this feature more clearly with the deeper analysis outlined in § 4.5.

4.5. General long-time dynamics

To extract the long-time dynamics explicitly for more general initial conditions, we consider the alternative formulation of the problem in (4.8)–(4.9), and again approximate the sum as an integral. We adopt the initial condition

(4.18) \begin{equation} a_{\!j}(0) = a_0 \varUpsilon (j) , \end{equation}

where $a_0$ is again the overall magnitude and $\varUpsilon (j)$ describes some order-one spectral form (for our previous power law in (4.10), $\varUpsilon =j^{-p}$ ). After defining the new variables,

(4.19) \begin{equation} v = \frac {\sqrt {\textit{Gt}}}{j^2\varPsi }, \quad \varPsi (t) = - \frac {\varTheta }{\sqrt {\textit{Gt}}} , \end{equation}

the mode amplitudes become

(4.20) \begin{equation} a_{\!j} = a_0 \varUpsilon (j) e^{-3(j^{-4}\textit{Gt} + j^{-2}\varTheta )} = a_0 \varUpsilon (j) e^{3\varPsi ^2 v(1-v)}, \end{equation}

and we arrive at the implicit equation,

(4.21) \begin{equation} \sqrt {\frac {G}{t}} \varPsi + \frac {2 a_0^2 G}{(\textit{Gt})^{1/4} \varPsi ^{3/2}} \int _{v_-}^{v_+} \big[e^{ 6 \varPsi ^2 v (1-v)}-1\big] \; \frac {\varUpsilon ^2\,\textrm {d} v}{v^{5/2}} = 4 |\varDelta | , \end{equation}

in which

(4.22) \begin{equation} (v_-,v_+)=\frac {\sqrt {\textit{Gt}}}{\varPsi }(N^{-2},1) \quad \textrm {and} \quad \varUpsilon =\varUpsilon \left (\mbox{$\dfrac {(\textit{Gt})^{1/4}}{v^2\sqrt \varPsi }$}\right )\!. \end{equation}

The factors $E_{\!j}=e^{3\varPsi ^2 v(1-v)}$ in (4.20) and $E_{\!j}^2-1=e^{ 6 \varPsi ^2 v (1-v)}-1$ in (4.21) both become exponentially localised to $v={1}/{2}$ when $\varPsi \gg 1$ . This location in $v$ corresponds to a wavenumber of

(4.23) \begin{equation} j = j_*(t) = (\textit{Gt})^{1/4}\sqrt {\frac {2}{\varPsi }} = \sqrt {-\frac {2\textit{Gt}}{\varTheta }} . \end{equation}

As long as $1\lt j_*\lt N$ , mode $j$ therefore reaches a peak at the moment that $j=j_*(t)$ , and the maximum amplitude is

(4.24) \begin{equation} \textrm {max}(a_{\!j}) \equiv a_{\textit{max}}(t) = a_0 \varUpsilon (j_*) \exp \left ( {\frac {3}{4}} \varPsi ^2 \right ) = a_0 \varUpsilon (j_*) \exp \left ( \frac {3\varTheta ^2}{4 \textit{Gt}} \right )\!. \end{equation}

The predictions for the wavenumber $j_*$ and maximum amplitude from (4.23) and (4.24) respectively are included in the top two rows of figure 4 for the power-law spectrum in (4.10). Evidently, the values of $\varPsi (t)$ achieved in the computations (which are shown in the lowest panels) are sufficiently large to ensure that $E_{\!j}$ becomes exponentially localised, thereby building up power around the locus $j=j_*(t)$ in the $(j,t)\text{-}$ diagrams and placing an envelope on the combined mode amplitudes. The generation of a cascade therefore hinges on the dominance of this exponential factor; for the cascade to become nearly self-similar, with $j_*\sim t^{1/4}$ , $\varPsi (t)$ must depend only weakly on time, as we argue next.

The condition $1\lt j_*\lt N$ , or $v_-\lt {{1}/{2}}\lt v_+$ , also ensures that the maximum arises within the range of the integral in (4.21). In this circumstance and provided that the variation of $\varUpsilon (j)$ does not modify the exponential localisation of $E_{\!j}^2-1$ to $v= 1/2$ , we may apply Laplace’s method (Bender & Orszag Reference Bender and Orszag2013) to find

(4.25) \begin{equation} \varPsi \sqrt {\frac {G}{t}} + \frac {8a_0^2 G}{\varPsi ^{5/2}(\textit{Gt})^{1/4}} e^{3\varPsi ^2/2} [\varUpsilon (j_*)]^2 \sqrt {{\frac {1}{3}}\pi } \sim 4 |\varDelta | . \end{equation}

Assuming that the first term on the left is exponentially smaller than the second, we conclude

(4.26) \begin{equation} \varPsi ^{-5/2} e^{3\varPsi ^2/2} [\varUpsilon (j_*)]^2 \sim \frac {|\varDelta |}{2a_0^2G} \sqrt {\frac 3\pi } (\textit{Gt})^{1/4}. \end{equation}

Given the exponential $e^{3\varPsi ^2/2}$ on the left-hand side, we anticipate a logarithmic dependence of $\varPsi ^2$ on time. However, without making a choice for $\varUpsilon (j)$ , extracting further detail is difficult because the unknown function $\varPsi (t)$ permeates the formulae. In particular, it is challenging to confirm that $\varPsi ^2$ does indeed become large and verify the validity of all the various approximations that have been made in arriving at (4.26). That said, the presence of the factor $a_0^{-2}t^{1/4}$ on the right-hand side leads one to expect that this limit may be achieved for either low-amplitude initial conditions, long times or a conspiracy between both.

When $\varPsi (t)$ does turn out to depend logarithmically on time, $j=j_*$ is close to a multiple of $t^{1/4}$ and a viscous wrinkling cascade with a roughly self-similar form is expected as long as

(4.27) \begin{equation} 1 \lt (\textit{Gt})^{1/4} \sqrt {\frac {2}{\varPsi }} \lt N . \end{equation}

The failure of either of these conditions implies that Laplace’s method no longer provides the dominant contribution to the integral in (4.21). The first condition refers to the initiation of the cascade, which must arise for sufficiently late times; the failure of the second, for $N\sim (\textit{Gt})^{1/4}\sqrt {2/\varPsi }$ , improves our earlier estimate for the time at which the cascade becomes arrested.

For the special case of the power-law initial spectrum (4.10), the reduction can be more consistently traced out (see Appendix A) and we find

(4.28) \begin{equation} \varPsi ^{p-5/2} e^{3\varPsi ^2/2} \sim \frac {2^{p-1} |\varDelta |}{a_0^2 G} \sqrt {\frac {3}{\pi }} (\textit{Gt})^{(1+2p)/4} = B . \end{equation}

Hence, after taking a logarithm and iterating once,

(4.29) \begin{equation} \varPsi \sim \sqrt {{\frac {2}{3}} \log \left [B \left (\frac 23 \log B \right )^{(5-2p)/4} \right ]} ,\end{equation}

which predicts the limiting long-time dependence of $\varPsi$ , as illustrated in figures 4(a) and 4(b) for $p=2$ and 0. For $p=- 1/2$ , $\varPsi$ is actually constant and approximately equal to $4.0$ , in agreement with the long-time numerical solution in figure 4(c). Note that the cases with $p\lt 3/2$ require more careful attention to justify the application of Laplace’s method; see Appendix A.

Finally, with (4.28) in hand, we can establish the specific conditions that lead to a viscous cascade for the power-law initial condition (4.10). When $p\gt {{1}/{2}}$ , $\varPsi ^2\gg 1$ is guaranteed for either long enough time $t\gg 1$ (provided (4.27) is not violated) or sufficiently low-amplitude initial condition $a_0\ll 1$ . When $p\lt {{1}/{2}}$ , $\varPsi ^2\gg 1$ can be achieved when $a_0\ll 1$ (cf. Appendix A). For the computations shown in figures 2 and 4, the initial condition has sufficiently low amplitude for $\varPsi ^2$ to attain values of approximately 16, and so a viscous cascade is observed.

5.1. Linear dynamics and preliminary single-mode buckling states

For $|Z|\ll 1$ , we enter the linear buckling limit, in which the bending moment remains small (of order $W$ ), whereas the tension remains order one and is given by the mean straining motion. Hence, $\varXi \gg 1$ and

(5.10) \begin{align} {\mathcal T} &\to 2^{m+1} |\varDelta |^m \textrm {sgn} (\varDelta ) \equiv 2 \varSigma , \\[-12pt]\nonumber \end{align}

(5.11) \begin{align} {\mathcal M} &\to -{\frac {1}{3}} m |2\varDelta |^{m-1} W_{{x}{x}} \equiv -{\frac {1}{3}} m \mu W_{{x}{x}}. \end{align}

Plugging these results into (5.6)–(5.8), and again restricting the solution (3.1) to the first Fourier mode with $Z=A_1e^{ix}+c.c.$ and $A_1=a_1=a(t)$ gives

(5.12) \begin{equation} \dot a = - \frac {3}{m\mu }\left ( {G} + 2 \varSigma \right ) a \equiv \gamma a . \end{equation}

Figure 9.

Single-mode buckling solution of (5.6)–(5.8) for $(m,\varDelta,{G},N)=(1/3,-1,3/2,256)$ . The first mode, $a_1(t)=|A_{\!j}(t)|$ , is initialised with an amplitude of $10^{-5}$ ; the higher modes have zero initial amplitude ( $a_{j\gt 1}=0$ ). Shown are time series of $a_1(t)$ and ${\mathcal T}(t)$ . The dotted line is the next strongest mode amplitude, with $j=3$ , and the dashed lines show ${\mathcal T}=-G$ and (5.14). The dash-dotted lines indicate the expected linear, exponential growth.

Figure 9 displays a solution of the initial-value problem for (5.6)–(5.8) with $m={1}/{3}$ . In this computation, the first mode is initialised with an amplitude of $10^{-5}$ ; the higher modes (up to $j=N=256$ ) all have zero amplitude. Although the nonlinearity in the constitutive model immediately drives growth in the higher modes, the solution is dominated by the first mode over the time interval displayed in the figure. That mode grows exponentially at early times, following the prediction in (5.12). Once the $j=1$ mode amplitude reaches order-one values, however, nonlinear effects slow down the growth of the mode. At this stage, $\varXi$ had declined towards order-one values, and the nonlinearity of the constitutive model in (5.7) and (5.8) becomes relevant.

For the later times plotted in figure 9, the buckles again enter into a tension-dominated regime in which ${\mathcal T}\sim - G$ and the mode amplitude grows algebraically with time. This growth implies that $|W_{{x}{x}}|$ remains relatively small even though $Z(x,t)$ is now order one, so that $\varXi$ mostly remains larger than unity. Hence,

(5.13) \begin{equation} {\mathcal T} \sim 2^{m+1} |(U_{x}+Z_{x} W_{x})|^{m-1} (U_{x} + Z_{x} W_{x}) \sim - {G} , \end{equation}

which, in combination with a spatial average, gives

(5.14) \begin{equation} a \sim t^{\frac 12} \left [ -\varDelta - {\frac {1}{2}} \left ( {\frac {1}{2}}{G}\right )^{\frac 1m}\right ]^{\frac 12}. \end{equation}

The predictions in (5.13) and (5.14) are included in figure 9.

At even later times, the excitation of higher modes ensures that the tension-dominated, single-mode state (5.14) cannot survive indefinitely. Indeed, the amplitude of the next strongest mode, with $j=3$ and also displayed in figure 9, grows faster than $a_1(t)$ , suggesting an impending nonlinear cascade as in the viscous buckling problem.

5.2. Modified cascades

Three sample solutions to (5.6)–(5.8) with $(m,G,N)=(1/3,3/2,256)$ are displayed in figure 10(a–c). The first case corresponds to the continuation of the computation in figure 9 with the single-mode initial condition. The other two cases show examples with different initialisations. For the second, the amplitudes of the first five modes are all set to $a_0=10^{-5}$ (and the others zeroed), complementing the viscous solution in figure 2(a). In the last case, we use the initial condition in (4.10) with $p=2$ ; cf. figure 4(a).

Figure 10.

Multi-mode buckling solutions of (5.6)–(5.8) for $({m},\varDelta ,{G},N)=(1/3,-1,3/2,256)$ . The initialisations are (a,d) $a_1(0)=10^{-5}$ and $a_{j\gt 1}(0)=0$ , (b,e) $a_{j\leq 5}(0)=10^{-5}$ and $a_{j\gt 5}(0)=0$ , and (c,f) $a_{\!j}(0)=j^{-2}\times 10^{-5}$ . Shown are time series of $a_{\!j}(t)$ ( $j=1,2,\ldots ,N$ ) and density plots of $\log _{10}a_{\!j}(t)$ on the $(t,j)\text{-}$ plane, split into (a–c) earlier and (e–f) later sections. The (red) points and dashed lines show the times at which $a_{\!j}$ reaches a maximum exceeding one. The corresponding tensions ${\mathcal T}(t)$ are plotted in panel (g) (with the solutions in panels (a,d), (b,e) and (c,f) shown by solid, dotted and dashed lines, respectively; the last two cases are hard to distinguish in the plot); the triangle illustrates the scaling $t^{-1/2}$ .

Figure 10. Long description

The image contains multiple heatmaps and line graphs depicting the buckling behavior of a thin sheet of power-law fluid under compression. Panel A: A heatmap shows the time series of mode j versus t to the power of 1/4, with a color scale indicating the magnitude of a(t). Red lines and stars mark times when a(t) reaches a maximum exceeding one. Panel B: Similar to Panel A but with different initial conditions. Panel C: Another heatmap with different initial conditions. Panel D: A heatmap showing a broader range of modes j versus t to the power of 1/4, with red lines and stars marking significant times. Panel E: Similar to Panel D but with different initial conditions. Panel F: Another heatmap with different initial conditions. Panel G: A line graph plots the tension T versus t on a logarithmic scale, with solid, dotted, and dashed lines representing different cases. The triangle illustrates the scaling of -1/2.

As in the viscous problem, for all three examples in figure 10(a–c), the uncontrolled growth of the higher modes eventually triggers a nonlinear buckling cascade. However, the dynamics is less regular, with power building up more quickly in the higher modes as a result of the nonlinear bending moment and tension; see the continuations of the computations to longer times in figure 10(d–f). After reaching their peaks and declining, the lower modes also do not become completely suppressed. Consequently, power is no longer concentrated onto a small subset of the modes, but spreads more evenly across the entire spectrum at later times. All the while, the net compression is slowly relieved, with $\mathcal T$ again following a trend close to $t^{-1/2}$ (figure 10 g). The decay of the tension, coupled with the emergence of order-one bending moments, implies that $\varXi$ now decreases well below unity except near the nodes of $M(x,t)$ . In other words, during the cascade, the rate of change of the curvature dominates in-plane straining.

Note that the relatively sudden emergence of power in the higher modes for $t\approx 0.9$ in figures 10(a) and 10(b) corresponds to the moment that $\varXi$ first falls below unity in these computations. In the first example, the breach through this threshold does not immediately lead to the excitation of all the higher modes, however, because the nonlinearity in (5.7) is cubic when $m=1/3$ , and transfers power only through odd $j$ by itself. The second two examples in figure 10 also display similar dynamics for the cascade, despite their different initial conditions. Evidently, the cascade of power from the lower modes proceeds in largely the same fashion, washing out the differences in initialisation.

Figure 11.

Reconstruction of the deflection $Z(x,t)$ for the solutions shown in (a) figures 10(a) and 10(d) and (b) figures 10(b) and 10(e). The red lines show the loci of the nodes.

A reconstruction of the deflections for the first two solutions is presented in figure 11 (the third solution in figure 10 generates a very similar spatio-temporal buckling pattern to that in panel b). As for the viscous buckling problem, the nonlinear cascade corresponds to protracted wrinkling of the layer. However, because the lower modes now fail to decay, the buckling pattern retains longer-scale spatial modulation. For the symmetrical solutions in figure 11, the result is that the wrinkles become localised to half of the domain.