Critical Analysis of Mediterranean Sea Level Limit Cycles During the Messinian Salinity Crisis

The Messinian Salinity Crisis (5.97-5.33 Ma) may be one of the most significant periods of sea-level change in recent geologic history. During this period, evaporite deposition throughout the Mediterranean basin records a series of dramatic environmental changes as flow through the Strait of Gibraltar was restricted. In the first stage of evaporite deposition, cycles of gypsum appear in shallow basins on the margins of the Mediterranean. The complex environmental history giving rise to these cycles has been investigated for decades but remains controversial. Notably, whether the evaporites are connected to significant changes in Mediterranean sea-level is an open question. In one proposed model, competition between tectonic uplift and erosion at the Strait of Gibraltar gives rise to self-sustaining sea-level oscillations -- limit cycles -- which trigger evaporite deposition. Here I show that limit cycles are not a robust result of the proposed model and discuss how any oscillations produced by this model depend on an unrealistic formulation of a key model equation. First, I simplify the model equations and test whether limit cycles are produced in 64 million unique combinations of model parameters, finding oscillations in only 0.2 % of all simulations. Next, I examine the formulation of a critical model equation representing stream channel slope over the Strait of Gibraltar, concluding that a more realistic formulation would render sea-level limit cycles improbable, if not impossible, in the proposed model.


INTRODUCTION
At the end of the Miocene, a combination of environmental and tectonic processes dramatically altered the water budget and chemistry of the Mediterranean Sea, leading to the Messinian Salinity Crisis (MSC) between 5.97 and 5.33 Ma (Manzi, Gennari, F. Hilgen, Krijgsman, Lugli, Roveri, and F. J. Sierro 2013). This period is marked in the stratigraphic record by evaporite deposition throughout the basin (Hsü, Ryan, and Cita 1973). It is generally thought that the crisis unfolded in three stages (CIESM 2008;Roveri, Flecker, Krijgsman, Lofi, Lugli, Manzi, F. J. Sierro, Bertini, Camerlenghi, De Lange, Govers, F. J. Hilgen, Hübscher, P. T. Meijer, and Stoica 2014). In the first stage, cycles of gypsum and marine marl were deposited in shallow, marginal basins. In stage two, massive halite deposits formed in deeper basins. In the third and final "Lagomare" stage, large fluctuations in salinity are exhibited by gypsum-marl pairs and evaporite-free deposits.
However, fundamental questions about the causes and timing of events during the MSC remain. A precise geologic and environmental chronology is elusive because of the paucity of marine fossils in Messinian sediments and lithologies in the basin that are not well suited for accurate dating (Roveri, Flecker, Krijgsman, Lofi, Lugli, Manzi, F. J. Sierro, Bertini, Camerlenghi, De Lange, Govers, F. J. Hilgen, Hübscher, P. T. Meijer, and Stoica 2014).
In one explanation of the cycles, Garcia-Castellanos and Villaseñor (2011) proposed an elegant model where competition between tectonic uplift and erosion at the Strait of Gibraltar gives rise to oscillation of Mediterranean sea-level (MSL), repeatedly triggering gypsum deposition. Their model exhibits limit cycles (Strogatz 1994), where no oscillatory external forcing, such as Milankovitch cycles, is required to produce oscillations in Mediterranean sea level. If correct, this model would constitute an unusual and fascinating example of large-scale sea-level change that is not driven by redistribution of water between the oceans and the cryosphere.
There are geological reasons to doubt this model. While a major erosional surface is found at the top of the unit containing the gypsum beds, no subaerial erosional surfaces have been observed at the tops of individual gypsum beds and they are thought to have formed in waters shallower than 200 m (Lugli, Vinicio, Marco, and Charlotte 2010). The model, however, exhibits >400 m fluctuations in Mediterranean sea-level, considerably greater than the observations indicate.
Here I show that there are also computational and conceptual reasons to doubt the model.
In the Model Formulation section, I consolidate the original model equations into a system of two explicit, analytic, ordinary differential equations (ODEs). In the Simulations section, I describe how these equations are solved, explain my computational approach to identifying limit cycles, and show how likely limit cycles are to occur in different parameter ranges. In

MODEL FORMULATION Original Model
The defining characteristic of the MSL limit cycle model (Garcia-Castellanos and Villaseñor 2011) is its capacity to produce oscillations without any external periodic forcing. Water is exchanged solely between the Mediterranean Sea and the ocean, without any forcing by the background climate or polar ice mass. The model includes three primary physical processes. First, the height of the sill at Gibraltar, where water flows from the Atlantic to the Mediterranean, is controlled by fluvial erosion and tectonic uplift. Second, sea level in the Mediterranean responds to a water budget including discharge over the sill, evaporation, direct precipitation, and input from continental rivers. Third, sea level in the ocean, outside the sill, changes to conserve the water lost and received by the Mediterranean.
As the original authors describe, the competition between uplift and erosion at the sill appears to give rise to an oscillatory coupling between erosion and Mediterranean sea level, possibly explaining the cyclic evaporite deposits in the first stage of the MSC. In this conception, the Strait of Gibraltar is initially open and the Mediterranean is full. As the sill is slowly uplifted, the sill depth is reduced and flow to the Mediterranean is restricted, causing MSL to drop. Next, with the equations used, the increased hydraulic head difference between the Atlantic and Mediterranean causes a nonlinearly accelerating increase in erosion at the sill. Increased erosion deepens the sill, enabling greater flow to the Mediterranean, which refills the basin and raises MSL. Then uplift continues and the cycle repeats.
More recently, Coulson, Pico, Austermann, Powell, Moucha, and Mitrovica (2019) illustrated the importance of additional sea-level physics as the Strait of Gibraltar opens and closes, extending the original model. They coupled a sophisticated sea-level model to the original erosion-uplift equations (Garcia-Castellanos and Villaseñor 2011), incorporating the effects of self-gravitation in the water bodies and crustal deformation in response to changing water load. With the additional physics, the model still exhibits limit cycles, but with a slower uplift rate at Gibraltar that is more consistent with estimates from independent geo-dynamical models (Andrews and Billen 2009;Duggen, Hoernle, and Morgan 2003;Duretz, Gerya, and May 2011;Gerya, Yuen, and Maresch 2004). The result appears to strengthen the idea of limit cycles during the first stage of the MSC, as the inclusion of well-established sea-level theory brings a key model parameter, the uplift rate, closer to expected values.

Simplified Model Equations
The original model has four dependent variables: ocean sea-level, Gibraltar sill height, western Mediterranean sea level, and eastern Mediterranean sea level (Garcia-Castellanos and Villaseñor 2011). First, to consolidate and simplify the system, I assume that the Mediterranean behaves as a single basin, setting aside the separate treatment of the western and eastern regions. The Sicily sill, which separates the eastern and western basins, may have played a role when MSL was low (Just, Hübscher, Betzler, Lüdmann, and Reicherter 2011). However, in the original model, oscillations are shown to occur almost entirely above the Sicily sill level (-430 m), so the exclusion of the sill should not alter the dynamics of the entire system. This assumption will be further discussed in the final section.
where A m is the Mediterranean surface area, z m is MSL, and parameters c 1 , c 2 , α 1 , and α 2 and are found by fitting equation (1) to present-day hypsometry in Figure 2 of (P. Meijer and Krijgsman 2005). The values of these parameters are shown in Table 1. Figure 1 shows the result of the fit and the residual, which is less than 1 % to a depth of at least 1200 m.
Equation (1) is simpler than linear interpolation and a slightly better representation of a smoothly varying surface area. It also has the advantage that it can be used to formulate an expression for the ocean sea level, as I explain next.
The model assumes all water lost from the Mediterranean is instantaneously received by Equation (2)  Integrating equation (2), using equation (1) for A m , produces an expression for the ocean level as a function of the Mediterranean level, whereż s is the time derivative of the sill height (dz s /dt), U is the uplift rate at the sill, k b is an erodibility coefficient, τ is the shear stress exerted by the flowing water, τ c is the critical shear stress, and a is an erosion exponent. The "max" operator prevents erosion from occurring when τ < τ c . The shear stress τ is a function of the channel's depth and slope, where ρ is water's density, g is gravitational acceleration, z o − z s is the approximate depth of the channel, and S is the slope of the water surface. Garcia-Castellanos and Villaseñor (2011) compute the slope using where L is a constant length of 100 km, representing the approximate half-width of the Betic-Rifean orogen (Garcia-Castellanos and Villaseñor 2011). Equation (6) is meant to approximate the mean channel slope. Plugging equations (6) and (5) into equation (4) yieldsż a single expression for the rate of change of the sill height that depends only on z s and z m because z o is a function of z m defined by equation (3).
The second component of the model is the level of the Mediterranean. Changes in MSL are governed by input and removal of water, where P is direct precipitation into the Mediterranean, E is evaporation, R accounts for input from continental rivers, and Q is input from the ocean. By design, discharge over the sill only occurs from the Atlantic into the Mediterranean, without any return flow. The Mediterranean area A m is computed with equation (1). Garcia-Castellanos and Villaseñor (2011) compute Q with a simple geometric relationship, where W is the width of the channel flowing over the sill, z o − z s is the channel depth, and V is the flow velocity. The velocity is represented by Manning's formula, where n is a roughness coefficient, R h is the hydraulic radius, and S is again the slope.
Because the channel depth is expected to be considerably smaller than the channel width, the hydraulic radius is approximated by the channel depth, R h = z o −z s . Garcia-Castellanos and Villaseñor (2011) choose a form for the channel width W that accounts for the effect of uplift (Turowski, Lague, and Hovius 2007), where C w is an empirical constant. For clarity, I let By plugging equations (11) and (10) into equation (9), then plugging the result into equation (8), the expression forż m becomeṡ Here the "max" operator handles cases where the sill becomes higher than the ocean. When this occurs, the operator yields zero and flow from the ocean is shut off, disconnecting the Mediterranean and allowing the sill to rise indefinitely. Equation (13)  • Equation (1)   Next, I integrate the model for wide ranges of six key parameters. The "tested range" column of Table 2  If the first check is satisfied, the second check is performed. If the second check is also satisfied, the model is assumed to be stationary and integration stops. Fixed points are located using multivariate Newton's method and a finite difference approximation of the system's Jacobian. Performing this check in two stages prevents the initiation of Newton's method when the solution is not near a fixed point. This stringent, two-stage test is designed to prevent any oscillatory solutions from being erroneously identified as stable ones.    (4) and (5). It also controls discharge into the Mediterranean through equations (9) and (10). In this section, I analyze the choice of equation (6) for the channel slope and its consequences for model oscillation.

The Slope Equation
Garcia-Castellanos and Villaseñor (2011) use equation (6)  average slope of the newly exposed terrain is 1 %, the horizontal length of the channel has increased 10 km. Figure 5 shows a schematic of this relationship between MSL, the channel length L, and the average channel slope S. More realistically, L and z m would covary and the slope would be nearly constant for small changes in z m .
To calculate the true average slope along the channel, detailed bathymetry would be required. We do not know the exact bathymetry of the Mediterranean during the MSC, but we can consider modern data for intuition. Figure 6 shows modern bathymetry at the Strait of Gibraltar and into the western Mediterranean basin. If modern bathymetry is any guide, there is little reason to expect a constant channel length of 100 km as MSL varies over hundreds of meters. In fact, the average slope would probably decrease when MSL drops, not increase, because deeper parts of the strait and nearby basin are flatter.
Further, the nonuniform slope through the modern strait raises doubts about whether the average slope is most applicable when MSL varies hundreds of meters. Although the average slope is simple and convenient, the local slope near the sill may be different than the average slope when MSL drops significantly. We do not know the shape of the sill during the Messinian, but we can consider the modern configuration again as an example. If MSL dropped below the Camarinal Sill today, erosion on the sill would likely be controlled by the higher local slope, at least until erosion significantly modified the sill's profile.
In this scenario, where erosion occurs primarily over steeper terrain near the sill, the slope could theoretically be computed over a fixed length L and depend on the height of the sill as it rises and falls. This conception is illustrated in Figure 7. However, this representation would require an equation for the slope that depends on sill height, z s , and equation (6) does not depend on the sill height in any way. It depends only on z o and z m , with constant L. Therefore, equation (6) does not represent an average local slope on the sill being slowly modified by uplift.
To summarize, if equation (6) is meant to represent the average slope between the Atlantic and Mediterranean sea levels, it should account for changing channel length L instead of assuming L is constant. There is no reason to expect Mediterranean sea level to rise and fall hundreds of meters without any lateral movement of the shoreline. If instead, equation (6) is meant to represent a more local slope near the sill as it is uplifted and eroded, it fails to capture this process because it does not depend on the sill height.

Implications for Oscillation
Whichever vision of the average slope ( Figure 5 or 7) is preferred, equation (6) is not sufficiently representative. This is an important concern, as the choice of equation (6) is critical to the model's capacity to oscillate. For oscillation to occur, there must be feedback that increases erosion when the Mediterranean level drops. In this model, a decrease in z m causes an increase in the slope S, strengthening erosion and opening the channel to flood the Mediterranean and raise z m again. The dependence of S on z m creates the necessary feedback. Without this dependence, the feedback from z m toż s is broken.
Crucially, it is the proposed form of equation (6), with average slope calculated using fixed L, that introduces the dependence of S on z m . Because changes in z o are much smaller than changes in z m (< 1 %), the slope in equation (6) is approximately a linear function of z m . When z m drops, S increases proportionally. However, as discussed above, the mean slope would probably be constant or decrease, not increase, if the change in L and basin bathymetry are accounted for. Alternatively, if the local slope on the sill is more important than the mean slope, S would not respond to the value of z m at all. In both cases, a more realistic representation of the slope is likely to render limit cycles impossible in this simple model because it would remove the proportional relationship between S and z m . When MSL drops, there would not be increased erosion on the sill to deepen the channel, refill the Mediterranean, and generate limit cycles.
In the Revisiting Channel Slope section, I explained how equation (6) is not physically representative of an average slope between the Atlantic and Mediterranean. To calculate the average channel slope, it assumes the horizontal position of the Mediterranean shoreline is fixed, even as the Mediterranean sea level varies tens or hundreds of meters. It is also not representative of the scenario where erosion is controlled by the local slope on the sill as the sill is uplifted and eroded. This is simply because equation (6) is completely independent of the sill height z s . In either case, the channel slope should not significantly increase when z m decreases. The rate of erosion on the sill need not depend on Mediterranean sea level because the shore is downstream of the erosion process. This is a critical problem because the unrepresentative dependence of S on z m introduced by equation (6) (2019) show, future work must also account for the gravitational and isostatic effects of the changing water loads.
In conclusion, the proposed model only rarely exhibits limit cycles, and limit cycles require parameter values considerably higher than those originally reported. However, the proposed model includes an unphysical representation of the channel slope. This is im-  (Hunter 2007). All code used for the analysis and plots in this paper is permanently archived in Zenodo (Baum 2021a) and available with documentation at github.com/markmbaum/messinian-salinity-crisis.
Parameter Fit Value c 1 2.068 × 10 12 m 2 α 1 2754 m c 2 4.035 × 10 11 m 2 α 2 127.5 m Table 1: Parameter values for equation (1), the surface area of the Mediterranean, and equation (3), the ocean height.  (1), and the residual. The fit is accurate to better than 1 % down to at least 1200 m.  Simulation Outcome Percent disconnection 20.9 % oscillation 0.2 % fixed point 78.9 % Table 3: Outcomes of the 64 million simulations. "Disconnection" refers to solutions where zs > zo, the Strait of Gibraltar closes, and the Mediterranean becomes completely disconnected from the Atlantic. "Fixed point" refers to solutions that arrive at stable fixed points with an open channel between the Mediterranean and Atlantic, as in Figure 2.  Table 2. As the legend indicates, the black line shows the sill level, the light blue line shows the Mediterranean level, and the dark blue line shows the ocean level.  Table  2). Notably, there is a complete lack of oscillatory solutions near the reference values of τc and k b .