The mechanism of sapropel formation in the Mediterranean Sea: Insight from long duration box-model experiments

Periodic bottom water oxygen deficiency in the Mediterranean Sea has led to the deposition of organic rich sediments during geological history, so called sapropels. Although a mechanism linking the formation of these deposits to orbital variability has been derived from the geological record, physics-based proof is limited to snapshot and short time-slice experiments with (Oceanic) General Circulation Models. Specifically, previous modelling studies have investigated atmospheric and oceanographic equilibrium states during orbital extremes (minimum and maximum precession). In contrast, we use a concep5 tual box model that allows us to focus on the transient response of the Mediterranean Sea to orbital forcing and investigate the physical processes causing sapropel formation. The model is constrained by present day measurement data, while proxy data offers constraints on the timing of sapropels. The results demonstrate that it is possible to describe the first order aspects of sapropel formation in a conceptual box model. A systematic model analysis provides new insights on features observed in the geological record, such as timing of 10 sapropels, intra-sapropel intensity variations and interruptions. Moreover, given a scenario constrained by geological data, the model allows us to study the transient response of variables and processes that cannot be observed in the geological record. The results suggest that atmospheric temperature variability plays a key role in sapropel formation, and that the timing of the midpoint of a sapropel can shift significantly with a minor change in forcing due to nonlinearities in the system.


Background
The response of ocean circulation to changes in atmospheric forcing is an important element of the climate system. Using computer models applied to the geological past we can exploit the sedimentary record of variation in circulation for mechanistic insight. The Mediterranean Sea is of particular interest, as abundant and exceptionally well dated proxy data and present-day measurement data is available and it is a basin that displays processes such as thermohaline circulation and gateway control 20 that play a role on the global scale as well.
Presently, the Mediterranean Sea is an evaporative basin (Romanou et al., 2010) with a small annual mean heat loss to the atmosphere (Song and Yu, 2017). Water from the Atlantic flows in to the Mediterranean Sea at the Strait of Gibraltar and is then subjected to buoyancy loss due to evaporation and cooling (see Fig. 1 for a map of the Mediterranean Sea). This results Haines, 1996). During winter, in the northerly parts of the basin, situated at relatively high latitude, cold and dry winds induce a further density increase, which may lead to the formation of deep water (Schroeder et al., 2012). Specifically, deep water formation (DWF) occurs over the shallow northern Adriatic Sea (Malanotte-Rizzoli, 1991) and the Aegean Sea (Gertman et al., 2006;Roether et al., 1996) and in the form of open-ocean deep convection in the Gulf of Lion (Marshall and Schott, 1999) and the southern Adriatic Sea (Bensi et al., 2013). Dense water formed in the Adriatic and Aegean Seas, both marginal basins of 30 the Mediterranean Sea, flows out over the seafloor into the deeper parts of the main basin.
The basin's semi-enclosed nature causes the system to be very sensitive to climatic perturbations and the geological record holds an expression of this sensitivity in the form of the regular occurrence of organic rich deposits, known as sapropels (Rossignol-Strick, 1985;Rohling et al., 2015;Hilgen, 1991;Lourens et al., 1996;Cramp and O'Sullivan, 1999). Sapropels are thought to form when Nile discharge increases as a response to enhanced African summer monsoon activity during precession 35 minima (Rossignol-Strick, 1985;Rohling et al., 2015). The low density fresh water then forms a lid at the surface, stopping or reducing the strength of the overturning circulation. This commonly accepted mechanism can be nuanced by noting that the Nile does not enter the basin at a DWF site, but rather close to the location where intermediate water forms. A large part of the DWF involves this intermediate water (Schroeder et al., 2012). Reducing the density of the intermediate water implies a decrease or absence of a (positive) vertical density gradient, also diminishing or stopping the formation of deep water. In 40 contrast, runoff into the marginal basins directly affects the buoyancy at the DWF sites. For deep convection in for example the Levantine basin (which can happen with present day conditions, Gertmann et al., 1994) a decrease in surface water density directly decreases or stops DWF. With decreasing DWF, the supply of oxygen to the deep water diminishes, potentially causing anoxia and the preservation of organic matter in the Eastern Mediterranean Sea. Moreover, nutrient input increases with river outflow as well, thereby affecting primary production, export of organic carbon to the deep water and, consequently, oxygen 45 consumption (Calvert et al., 1992;De Lange and Ten Haven, 1983;Thomson et al., 1999;van Helmond et al., 2015;Weldeab et al., 2003). Sea level rise may also trigger sapropel formation (see Rohling et al., 2015, for sapropel S1), although this does not exclude monsoon intensity variability as the main cause.
In this paper we present a simple three box model of the Mediterranean Sea, which includes most elements commonly invoked to explain sapropel formation as described above. With the model we study which processes determine when and why 50 sapropels form the way they do. Our aim is to gain a new perspective on the timing of the sapropel, relative to the forcing, as a significant part of the geological time scale depends on this relation (Hilgen et al., 1995;Krijgsman et al., 1999) and views on the timing of the midpoint (the average of the top and bottom age) are contested in more recent publications (Channell et al., 2010;Westerhold et al., 2012Westerhold et al., , 2015. A low complexity model allow ::::: allows : us to perform :::: many : long runs and explore the parameter space to a much greater extent than high complexity models. ::: The :::: runs :::::::: described :: in :::: this ::::: paper ::::::: represent ::: but :: a ::::: small of the associated atmosphere box, fluxes I 1 and I 2 in Fig. 2). Note that the DWF in box 1 captures the behaviour of the marginal basins of the Eastern Mediterranean sea, but is also an approximation of the open ocean convection in the Gulf of Lion (see subsection 1.1). In the typical situation that the Mediterranean surface/intermediate water at the Strait of Gibraltar is more dense than the Atlantic water and E-P-R (the fresh water budget, where R is the total runoff) is positive, it is the outflow to the Atlantic (Qo in Fig. 2) that depends on the density difference between the adjacent water masses. The inflow 95 into the Mediterranean Sea is then the sum of the outflow to the Atlantic and the fresh water budget. The equations used in the model are further explained in subsection 2.3. In addition to the water fluxes, diffusive mixing is also included in the model. In contrast to the water fluxes, no net water transport occurs as a result of the mixing. Rather, properties are exchanged between adjacent boxes. The amount of horizontal mixing (between the upper boxes) is constant, while the vertical mixing is a function of the density difference between the boxes in question.

100
A first order approximation of deep water oxygen concentration is included in the model to get a better understanding of when oxygen deficiency occurs. The oxygen concentration of the upper boxes is assumed to be in equilibrium with the atmosphere and is therefore constant. The oxygen concentration of the deep water (box 3) depends on the deep water fluxes, mixing and oxygen consumption. Oxygen consumption is scaled with river outflow, as a first order approximation of the nutrient input, and oxygen concentration.
In the present-day Mediterranean Sea DWF is the last step in a chain of processes (see the introduction). Our model does not include the intra-annual variability, and the basin geometry is only represented in abstract form. However, in the sense that the model does capture both the effect of salinity increase and temperature decrease on upper-water density it is expected to form a fair representation, qualitatively speaking, of the essence of the overturning circulation. To which extent this is true will 115 have to follow from more advanced models. Note that the model of Matthiesen and Haines (2003) also neglects the seasonal cycle. During winter, convection occurs (Schroeder et al., 2012)  However, deep water would not form with annual average conditions and therefore we assume perpetual winter conditions.

Surface forcing
The transient response of circulation and water properties to precession induced climate change is modelled by altering the evaporation and river outflow for each box at every time step. The analyses presented in this paper all use sine-waves to force in this paper is derived from a normalized sine-wave with a 20 kyr period, to reflect climatic precession. The amplitude and offset is then altered for evaporation and fluvial discharge in boxes 1 and 2. The phase of evaporation relative to the precession forcing is uncertain (see subsection 3.4) and is therefore varied between runs, the phase of the river discharge is kept at 0 degrees.
The fluvial discharge in box 2 (R 2 ) is interpreted as the Nile outflow and other runoff from Africa. Prior to the construction of 130 the Aswan High Dam in 1964, average Nile discharge was 2.7·10 3 m 3 /s (Rohling et al., 2015). Present day runoff from Africa is approximately 1.4 · 10 3 m 3 /s (Struglia et al., 2004). A recent modelling study (Amies et al., 2019) suggests that peak runoff from Africa may have been up to 8.8 times larger than present during sapropel S5. That :::: Note ::: that : their study does not consider changes in outflow from Europe. Fluvial discharge in to the high latitude marginal basins of the Mediterranean Sea (R 1 in the model) is presently approximately 6.7 · 10 3 m 3 /s (Struglia et al., 2004) . Increased runoff from Europe into the eastern 135 Mediterranean has been proposed as a possible source for extra fresh water during precession minima (Rossignol-Strick, 1985;Rohling et al., 2002;Scrivner et al., 2004) The current net evaporation (E-P) is approximately 0.9 m/yr (Romanou et al., 2010). During sapropel times, net evaporation is hypothesized to have decreased (Rohling, 1994), although this has not been quantified. We therefore test a broad range of net evaporation, from 0.2 to 2 m/yr to accommodate for these uncertainties.

Model equations and parameters 140
Here we first discuss the flux equations resulting from the model set-up and assumptions described above, followed by the equations used to integrate all flux-equations into a fully functioning model. All parameters are given in Table 1. We use a matrix vector representation to calculate the temperatures, salinities, densities, oxygen concentration of the next time step. The (water and heat) flux magnitudes and mixing intensities define the elements of the matrices used for these calculations. This same matrix-vector representation could be used for an arbitrary configuration (and number) of boxes, to represent different 145 oceanographic settings. Observational and modelling studies (Herrmann et al., 2008;Schroeder et al., 2012) have shown that during colder winters, more deep water is formed. Hence, it makes sense that the magnitude of the vertical, downward fluxes (D 1 and D 2 , see Eqs. 2 and 3) depends on oceanographic (and thereby indirectly also atmospheric) conditions. The most simple way of implementing this behaviour on a yearly resolution, is to assume a linear relationship between the density difference and flux magnitude (similar to Matthiesen and Haines, 2003). When the density of the overlying water mass is smaller than 150 that of the deep water, the water column is stratified and no vertical flux exists. To clip negative components of a flux to 0, we use the form F j,i = max(0, a), where a is the flux in question. We therefore define the following mathematical operator: The proportionality of DWF to surface to deep water density difference is determined by an efficiency constant, c 13 and c 23 for D 1 and D 2 respectively. The magnitude of these constants is chosen in such a way that a realistic deep water flux occurs at a 155 present-day density difference. In the current circulation deep convection in the Levantine basin (represented by D2) does not 5 occur every year (Gertman et al, 1994;Pinardi et al., 2015), making it difficult to determine c 23 empirically. By assuming that the DWF process in box 2 is the same as in box 1 (i.e. linearly dependent on the vertical density difference), c 23 can be taken as 4 times larger than c 13 , proportional to the difference in surface area of boxes 1 and 2. We therefore define the DWF fluxes in Eqs. 2 and 3, where ρ 1 , ρ 2 and ρ 3 are the densities of boxes 1 to 3 respectively.
At the Strait of Gibraltar, the exchange has two components from which the in-and outflow is calculated (see equations below): a density driven flux Q o (Eq. 4) and a compensating flux Q i (Eq. 5). The magnitude of Q o has a square-root relation 165 to the horizontal density difference at the strait (where ρ 0 is the density of the Atlantic ocean), in accordance with Bryden and Kinder (1991) . Theoretically, this flux should be able to change direction, when the density difference changes sign. We therefore multiply the square-root of the absolute value of the density difference with the sign of the density difference. Note that the direction of the fluxes (i.e. whether it goes in or out of the Mediterranean Sea) is determined in Eqs. 11 and 10. The strait efficiency c 20 (the coefficient of proportionality between volume transport and density difference) is again calibrated on 170 present-day conditions (Schroeder et al., 2012;Jordà et al., 2017;Hayes et al., 2019). The compensating flux Q i can then be calculated as the difference of Q o and the total freshwater budget of the Mediterranean Sea, to allow for conservation of volume.
The equations above describe all fluxes driven by gradients. By combining these fluxes with the surface forcing, we can derive the other fluxes by assuming constant box volume. The next set of equations (Eq. 6-16) define elements of a matrix F, representing all water fluxes.
180 F 2,A2 = e 2 (7) Mixing has a major impact on oceanic circulation, and must therefore be included in the model. Unlike the water fluxes described above, mixing does not cause a net water transport between boxes, but rather an exchange of properties (salt, heat and oxygen). In the model, we distinguish between horizontal and vertical mixing. Horizontal mixing, between boxes 1 and 2, depends on a fixed length scale over which mixing occurs and diffusivity (see Eq. 17). Vertical mixing (see Eq. 18 and 19) depends on the density difference between the boxes in question, where a larger density gradient causes more mixing. d 1 , d 2

205
and d 3 in Eq 18 and 19 are the depths of boxes 1, 2 and 3 respectively and A 1 and A 2 the surface areas of boxes 1 and 2. Thereby the diffusivity of vertical mixing effectively depends on the density difference. When the water column is stratified, mixing does not stop completely, but rather decreases to a background level, representing the internal waves and other disturbances. In the model this is included by clipping the vertical mixing to a fixed level (k bg ) when the density gradient becomes very small or negative. Equations 17-19 define elements of a matrix M.
In our model, heat exchange with the atmophere is represented by a relaxation to a prescribed air temperature (e.g. Ashkenazy et al., 2012). When one uses a prescribed heat flux instead :: In : a :::::::: previous :::::: version :: of :::: the ::::: model ::: we :::: used :: a ::::::: constant :::: flux ::: (of :: 5 :::::: W/m 2 ), in general similar results are obtained. However, when the fresh water budget of the margins approaches zero and the circulation (almost) stops, the results are not realistic. In this situation the margins become almost completely isolated from the rest of the basin, causing a massive temperature drop that does not stop until the circulation starts again. In reality this 220 temperature drop would be limited by the atmospheric temperature, something the relaxation representation does capture.
We thus multiply the temperature difference between the atmosphere and the water by a relaxation parameter c A in W/(m 2 · K). The value of this parameter is chosen such that at present day temperatures, a heat loss to the atmosphere of approximately 5 W/m 2 occurs, in accordance with Song and Yu (2017) and Schroeder et al. (2012). In the matrix-vector representation the two relaxation boundary conditions correspond, upon the necessary conversion, to two elements of a matrix H.

225
In a previous version of the model we used a constant flux (of 5 W/m 2 ). With a present-day circulation this gives similar results, however, when the fresh water budget of the margins approaches zero, and the circulation (almost) stops, the results are not realistic. In this situation the margins become almost completely isolated from the rest of the basin, causing a massive temperature drop that doesn't stop until the circulation starts again. In reality this temperature drop would be limited by the atmospheric temperature, hence we use a relaxation.
The model parameters (excluding surface forcing) used in all runs are given in Table 1. Initial temperatures are set to 16 • C and salinities to 37 for all dynamic boxes. They have no effect on the outcome of the model runs after spin-up. With the strait efficiency used in this paper, the model has a typical equilibrium time of less then 1000 years, while a spin-up of 20 kyr is removed from the output. The temperature in boxes 1 and 2 relaxes to both the Atlantic temperature and the air temperatures, , T A2 therefore effectively set the temperature range of the dynamic boxes. The winter air 255 temperatures are within the range given by the Naval Oceanography Command (1987): 10 • C in the northwest to 15 − 16 • C in the southeast. The river inflow also affects temperature, but has a much smaller impact due to the relatively small amount of water (2 orders of magnitude smaller than the Atlantic exchange). The air temperatures are chosen as winter values, since average air temperatures do not result in a realistic atmospheric heat loss and DWF. The temperature of the river water does not have a large influence on the model outcome.

260
The volumes of the boxes are calculated from the depth and surface area for all dynamic boxes, where V , A and d are all vectors with three elements: Except for the lack of a flux from box 3 to box 1, water can flow between all boxes in both directions. In the model, three (boxes n and n − 1). The matrix G describes the water fluxes for the calculation of the new temperature, where I is the identity matrix and l the unit vector.
The matrix P describes the water fluxes for the calculation of the new salinity, where J is a matrix of ones. The only difference 275 with G being that evaporation is excluded (since evaporated water does not contain salt).
The matrix N describes the mixing fluxes for the calculation of both the new temperature and salinity.
W is a vector so that W i = 1 Vi . Similar to F, the heat fluxes (equations 20 and 21) are placed in H, which is of the same size 280 as F and where all undefined elements are zero. We use a time step, dt, of 1 year, unless noted otherwise. Then the change in temperature for each time step equals: and for salinity: The density for the next time step is calculated from the temperature and salinity using the EOS80 formula (on Oceanographic Tables, 1986). Note that vectors V , S, T and ρ include both static and dynamic boxes. The deep water oxygen concentration of the next time step is similarly: Note that the oxygen concentration is only calculated for the deep water (making O and O consumption effectively scalars) and

Statistical analysis
One of the results of the model is that slight variations of forcing parameters can cause significantly different sapropel duration and timing. We therefore introduce a statistical test to determine the magnitude thereof, given the uncertainty of each of the 295 forcing variables. With eleven forcing parameters (the phase of evaporation and minima and maxima of R 1 , R 2 , T A1 , T A2 and evaporation) it is not feasible to calculate all permutations at a meaningful resolution. We therefore randomly pick and run 200 permutations (fewer permutations would produce unreliable results), given the uncertainty of each parameter , and calculate the 1σ and minimum and maximum values of the resulting oxygen concentrations per time step (see Fig. A1 for an example).
During testing, we can thereby visualize much more of the parameter space than when doing individual runs.

300
3 Analysis and results

Reference experiment
In the reference experiment the sine functions for the forcing are calibrated such that the precession maximum corresponds to present-day values-given that the orbital configuration is close to a precession maximum today. The curves are shown in falling at T=0 and T=20 kyr. The precession minimum sits at T=10 kyr. Nile outflow (R 2 ) increases from 5 · 10 3 to 3 · 10 4 m 3 /s, while river outflow from Europe only increases from 5 · 10 3 to 1.2 · 10 4 m 3 /s and evaporation decreases from 0.9 to 0.75 m/yr. While quantitative reconstructions of fluvial discharge and evaporation during sapropel formation are not available, these minimum and maximum values are in agreement with Marzocchi et al. (2015). All other parameters, found in Table 1, 310 do not vary with time. Table 2 shows the forcing parameters that are the same for all presented runs and Table 3 gives the forcing parameters that are changed between runs. ::: As :::::::: explained :: in ::: the ::::::::::: Introduction ::: we :::::: choose :::: these ::::::::: particular :::: runs ::::::: because ::: they :::::::: illustrate :::: well ::: the :::::::: sensitivity ::: of ::: the ::::: model :: to ::: the :::::: various :::::::::: controlling :::::: factors. : From Time=0 towards the precession minimum, the river outflow increases and, as a result, salinities decrease, as shown in for a decrease in temperature at the margins in the interval surrounding the precession minimum (we will come back to this below). As temperature does not change much, density variability ( Fig. 3D) is largely determined by changes in salinity. The dip in marginal temperature has an opposite effect on density compared to the salinity fluctuation, consequently the decrease in surface to deep density gradient is relatively small, and the marginal density does not drop below the deep water density.
Nevertheless, the decrease in the vertical density difference causes a decrease in DWF (D 1 Fig. 3E, also see equation 2).

Addition of atmospheric temperature variability
As described in the introduction, temperature variation due to precession likely also affected buoyancy loss. In order to examine this aspect, we run the model with a 3 • C temperature increase at the precession minimum relative to the precession maximum.

370
For atmospheric box A 1 the temperature increases from 12 to 15 • C, and for box A 1 the temperature increases from 10 to 13 • C. Both air temperature curves are described by sine waves, as shown in Fig. 4C. We decide to maintain a constant temperature difference between the two atmospheric boxes as there is insufficient evidence for other options. All other parameters are set as described in the reference run.
The overall behaviour of the model is similar to that in the reference run, except that the temperatures of all boxes are now 375 higher during the interval surrounding the precession minimum (Fig. 4C). We still observe a minor decrease in marginal water temperatures at the precession minimum (cf. Fig. 3C), albeit much smaller than in the reference run, since it is now imposed on top of the trend caused by the changing atmospheric temperature. The net effect of a homogeneous basin wide temperature increase during the precession minimum is a further decrease in DWF during this time interval. We find a sapropel from t=8084 years to 10970 years, which therefore lasts 2013 years and the midpoint leads the precession minimum by 473 years (see significantly affects the duration of sapropel conditions in the model. Since both observational and modelling studies find this temperature variability (Marzocchi et al., 2015;Herbert et al., 2015), it will be included in all following model runs.

Nonlinear behaviour
Next, we explore the effect of a transition to and from a time interval with a positive freshwater budget. Whether or not the freshwater budget of the Mediterranean Sea becomes positive during sapropel formation has been widely debated (Rohling,395 1994, and references therein). Although our model cannot directly prove whether or not this has happened, it does allow us to study what the implications for the water properties and circulation would be, which should help in recognising the expression 13 of a budget switch in the geological record. First we consider a scenario where only the freshwater budget of the margins becomes positive, in a subsequent run we force the model in such a way that the freshwater budget of the entire basin changes sign.

400
To have the freshwater budget of the margins become positive, the maximum of R 1 is increased from 1.2 · 10 3 m 3 /s to 1.4 · 10 4 m 3 /s (Fig. 5A), all other parameters are kept the same as in the temperature variability run (Fig. 4) : , :: as ::::: shown :: in ::::: table   : 3 :: in ::: the :::: row :: of :::::: "fwb1 :::: run". At a similar timing as the dip in temperature observed in the reference run, we now see a very large decrease in salinity at the margins from 9 to 13 kyr (see Fig. 5B). During this interval we observe that temperatures at the margins approach the temperature of the overlying atmospheric box, while deep water temperatures approach those found in 405 the open Mediterranean (see Fig. 5C). All are an expression of the disappearance of DWF at the margin (see Fig. 5E
430 Salinities (Fig. 6B) decrease in response to the decrease in net evaporation. When the freshwater budget reverses, the exchange with the Atlantic decreases, causing less relatively saline water to flow into the upper boxes. Consequently, the salinity of the upper boxes further decreases. The deep water salinity only begins to decrease more when DWF at the margin starts again. When the freshwater budget becomes negative again, the salinities abruptly increase and then follow the freshwater budget more or less linearly. The main features of the temperature curves (Fig. 6C) are caused by the same events that are 435 described above for the salinity variability, although temperature is also affected by heat loss to the atmosphere. Consequently, the same main features can be identified, with the difference that 1) the temperature of the upper boxes follows the air temperature curves, and 2) the amplitude is smaller, because the heat exchange with the atmosphere acts as negative feedback. The changes in densities are predominantly determined by salinity, as the changes in temperature are relatively small in this run.
Reversing the freshwater budget also causes the density difference between the Atlantic and open Mediterranean sur-440 face/intermediate box to change sign. Consequently, the density driven flow goes from the Atlantic to the Mediterranean, instead of the other way around. In Fig. 6E this is represented by the flux becoming negative. Note that this shift occurs almost instantaneously. In this run we find a very sharp termination of the sapropel, followed by a brief period with lower oxygen concentration (as shown in Fig. 6). This is caused by a peak

Phase of evaporation 450
Recent modelling studies (Marzocchi, 2016) have shown that while evaporation and river outflow are both forced by precession, they may have a different phase relation to their forcing. Runoff and evaporation are the only transient forcings in the presented model runs, therefore shifting runoff for example 2 kyrs forward in time gives the exact same wave shape as shifting evaporation 2 kyrs backwards in time. The only difference would be that the waveform would be shifted by 4 kyrs. Since we are primarily interested in the transient response rather than the absolute timing, we only investigate the effect of the phase of evaporation.

455
To assess the effect of the phase of evaporation on sapropel formation, we calculate the sapropel midpoint and duration for a set of runs, with varying evaporation phase (all other parameters remaining unaltered between runs). Apart from the phase of the evaporation forcing, the model is forced exactly the same as in the atmospheric temperature variability experiment (as described in subsection 3.2).
As shown in Fig. 7, we find a maximum in sapropel duration when evaporation is almost in phase with precession. This is to 460 be expected, as minimum evaporation then coincides with maximum river outflow. Similarly, a minimum in sapropel duration is found when evaporation is almost exactly in anti-phase with precession. In between these peaks the sapropel duration as a function of the phase of evaporation is described by a cosine. The sapropel midpoint is found to vary significantly (up to thousands of years with a large amplitude of evaporation variability) when varying the phase of the evaporation forcing. The shift in sapropel midpoint relative to the precession minimum is at a maximum when evaporation lags or leads approximately anti-phase (i.e. a shift of 0 or 10 kyr), the 5 kyr lead/lag falls right in between these points. The timing of the midpoint as a function of the phase of evaporation in between these extremes is described by a nearly perfect sine wave. Note that the peaks are not exactly at -5 and and 5 kyr, but slightly shifted, this is likely a result of the equilibrium time of the system. This experiment, combined with the systematic testing of the parameter space, highlights that although the exact timing and 470 duration depend on the exact forcing, the minimum in deep water oxygen concentration always occurs close to the precession minimum and the model response is always quasi-linear, as long as fresh water budgets are not reversed. We also find that the magnitude of the effect of the phase of evaporation on sapropel timing and duration depends on the amplitude of the evaporation variability (not shown). This makes sense, as the changes in circulation largely respond to freshwater budgets (the only difference between river inflow and evaporation being their respective temperatures) and the amplitude of evaporation 475 variability scales linearly with its impact on the variability of the freshwater budgets.
When systematically varying the components of the water budget within the limits mentioned in model setup, we find that in the regimes where the freshwater budget of (part of) the basin changes sign, sapropels are cut short considerably. When performing the same analyses described above, but now using the forcing of the first run in subsection 3.3, we find that this causes the midpoint of the sapropel to occur prior to the precession minimum (Fig. 8). Note that runs with multiple sapropel 480 intervals cannot be described as having a single midpoint or duration.

Model convergence
The time step of one year naturally results from the concept that deep water forms during winter storms, making it the highest resolution possible as long as seasonal variability is not included. From a purely mathematical perspective, however, the time 485 resolution should not affect the outcome significantly, as long as a sufficiently small time step is used to prevent aliasing. We tested this by varying the temporal resolution given a certain forcing. We find no significantly different results with a time step below 10 years. With a time step larger than 10 years aliasing occurs. We conclude that a time step of 1 year is sufficient for the analyses in this study.

490
All models require assumptions and simplifications to be made, as they are by design a simplified version of (part of) a system.
Simple box models, such as the model presented in this paper, aim to identify the smallest subset of processes that can describe a certain phenomenon. As such, this model represents a generic semi-enclosed basin, given that no specific geometry is included.
This implies that by altering the parameter values and in some cases the strait exchange equations, the model can easily be adapted to other semi-enclosed basins, such as the Black Sea.

495
In our model we parametrize intra-annual variability. Including seasonality would require separate intermediate water boxes (increasing complexity), while for the oxygenation of deep water only the amount of DWF and mixing with the overlying water mass is truly relevant. Furthermore, we would have to make assumptions regarding the annual variability of the forcing parameters (river outflow and evaporation), which are not well constrained for geological history. We therefore decided to parameterize the seasonal variability, by calculating a yearly averaged DWF flux based on winter temperatures. This allows us 500 to study the fundamental mechanisms of sapropel formation. OGCMs would be more appropriate to study the role of seasonal variation.
We do not include separate boxes for the Eastern and Western Mediterranean in our model. The aim of a conceptual model is to capture the first order aspects of a process with a minimal setup. The current setup does this. Incorporating separate subbasins would imply doubling the number of boxes and would also double the number of forcing parameters and equations, all 505 of which add uncertainty to the model (quantitative reconstructions do not exist for most of these parameters). Moreover, the complexity quickly increases, making it much harder to test and describe the parameter space and identify key mechanisms.
While this could be tested in a future study, we find it important to first understand the behaviour of the a semi-enclosed basin with a gateway before studying what is essentially a second order system. We expect that the effect in the Eastern basin is similar to the difference between a first order and second order filter: a larger shift of the midpoint (larger group delay), and 510 likely a higher sensitivity in the Eastern basin. Any resonances in the system are also expected to become more prominent (since the resonant frequencies are now amplified twice).
The model forcing used in this study is chosen to reflect either the variability envisaged in the commonly accepted mechanism (as sketched in subsection 1.1), or oceanographic and climatic variability deduced from modelling studies and the geological record as accurately as possible. Other processes such as the North Atlantic oscillation and solar activity are not taken into 515 account, because they are not thought to be of first order importance for sapropel formation, as described in Rossignol-Strick (1985) and Rohling et al. (2015), for example. While these processes likely influence sapropel formation, they are unlikely to be essential. Besides precession, obliquity also affects sapropel formation, but it is not our aim to reconstruct the exact conditions during specific time intervals. For an individual sapropel, adding an obliquity component would effectively slightly modulate the frequency and amplitude of the forcing. Since the model is not very sensitive to the exact frequency of the forcing, and we 520 already extensively tested the parameter space in terms of amplitude, a simple (20 kyr) sine wave suffices as forcing. Also note that since obliquity does not have a harmonic relation with precession, the modulation would not have the same effect on every sapropel. It likely affects the thickness of a sapropel for example, but the effect may work both ways when comparing different sapropels.
The model output comprises an average value for deep water oxygen, as the deep water is a single box. In reality, however, have to be included to specifically study bottom water oxygenation. We expect that the main difference with a biogeochemical 530 model is that in our model river input directly affects oxygen consumption, while the surface/intermediate boxes would act as a reservoir for nutrients (with their own feedbacks) in the biogeochemical model, thereby delaying the response to changes in river input. Note finally that in reality DWF occurs following two different mechanisms (open ocean convection and mixing of the water at margins during winter storms, which then cascades to the deep basin, see subsection 1.1), as well as in multiple sub-basins that each have their own conditions. The regime in Fig. 5 relies on the freshwater budget of the margins changing 535 sign. In reality there are many different marginal water masses in the sub-basins, rather than one single "margin". This makes it likely that the freshwater budget becoming positive in any one of these sub-basins will have similar consequences for the circulation. Since the freshwater budgets of these basins are independent, it would be possible to drastically alter the circulation multiple times during a single precession cycle. Presently, the Adriatic sea has a positive freshwater budget (Raicich, 1996), and the Aegean sea is known to have had a positive fresh water budget in the past (Zervakis et al., 2004).

540
The simplicity of the model makes it especially suitable for describing transient, nonlinear behaviour, allowing for the identification of crucial mechanisms. More complex models, while providing other benefits, are generally too difficult to interpret on this level, or do not allow for runs of sufficient length to study the transient response over a full precession cycle.
The presented model runs give an overview of the behaviour of the model. When systematically testing the parameter space, we find that this behaviour largely depends on general trends and reversal of fresh water budgets, rather than specific forcing 545 or parameter choices. This makes the results of the study much more robust and meaningful.
Exchange with the Black Sea also affects circulation in the Mediterranean Sea (Soulet et al., 2013, show increased runoff during HS1). For sapropels during which there is exchange through the Bosphorus strait, the exchange is not constant through time, and also depends on the inflow of Mediterranean water into the Black Sea. Opening or closing of the strait prior to, or during, a sapropel may impact the circulation. When the sill becomes deep enough to allow for a two layer exchange, a large 550 amount of saline water would flow into the Black Sea (following the same principle as at the Gibraltar Strait), thereby causing extra, relatively fresh, water to flow out into the Mediterranean Sea. During some sapropels, the strait may have been closed.
Furthermore, there is very little data regarding the exchange opening and closing of the Bosphorus Strait prior to the most recent opening (approximately 11 ka), perhaps with the exception of the Pontian (which is beyond the scope of this paper).
Consequently, we do not include exchange with the Black Sea. For the cases where there was a steady outflow of fresh water 555 (or an exchange that can be parametrized as such) this could indeed be seen as an extra fresh water source for the margins. We have already tested this effect by varying the river outflow into the margins.
Melt water pulses likely affect sapropel formation, but we do not consider them to be of first order importance. During many sapropels, melt water pulses did not occur. In future applications of this model where a specific sapropel/interval is studied, drivers such as melt water pulses should be included.

Describing nonlinear relationships and transient response
The occurrence of sapropels is often considered from a binary perspective: a sediment is either a sapropel, or it is not. The dominant forcing mechanism (astronomic variability), however, can be easily described by a combination of a limited number of sine-waves: the resonant frequencies of the planetary bodies in our solar system (for example Laskar, 1988). For a single sapropel, only climatic precession-a nearly a perfect sine-is considered to be of first order importance in controlling bottom 565 water oxygenation (Rossignol-Strick, 1985), with the rest of the orbital configuration mostly modulating the effect of preces-sion. If a model strives to capture the current hypothesis of sapropel formation, starting from astronomic variability, it must therefore transform a sine wave into something that is not a sine wave, requiring the model to be nonlinear. Our model allows for such behaviour. Even when considering intra-sapropel variability, thereby surpassing a binary approach, the sapropel record is clearly not sinusoidal (see for example Grant et al., 2016;Dirksen et al., 2019).

570
One of the main research questions of this study is when sapropel formation occurs. In a linear system, one would simply calculate the phase of the output with respect to the input. However, as the output is no longer linearly related to the input, this is not possible. A simple threshold analysis is not ideal either, as the cut-off level can have major impact on both timing and duration, while a clear definition is not readily available. Furthermore, even when the threshold is defined, this method would not be usable for sapropelic marls, which are thought to be the result of the same process, but do not share the same chemical 575 composition. We partly avoid this problem by not only considering the midpoint of the sapropel (when assuming a certain oxygen threshold, see subsection 4.5), but also the full wave form (e.g. which intervals could be sapropelic with a slightly different forcing). In the sedimentary record this is generally not possible, since the non-sapropelic intervals do not record all parameters and are often bioturbated. So while our approach cannot be applied to the sedimentary record, it does give insight into the factors that influence sapropel timing. Even this definition of sapropel timing becomes problematic when one or more 580 interruptions occur, since in that case there is more than one midpoint.
We find that, when using realistic model forcing, stable sapropel conditions do not occur. Even when using constant forcing, a permanent complete stop of DWF either does not occur, or only under very specific conditions. Note that in the Black Sea permanent stratification does appear to occur, but here, the positive freshwater budget allows some of the water flowing into the Black Sea at the Bosphorus Strait to sink to the deep water (Bogdanova, 1963), keeping it relatively saline and dense. However, 585 our results indicate that sapropel conditions can occur transiently without a positive freshwater budget, with realistic forcing.
We therefore conclude that studying the oceanographic state during sapropel conditions by modelling steady-state conditions with a stratified water column results in a very limited understanding of sapropel formation.

Comparison with geological data and other models
Comparing model results to geological data is most effective when an accurate age model is available for the geological 590 data.We will therefore only consider the 5 youngest sapropels in this paper. The most recent sapropel (S1) is thought to have been triggered by sea level rise, which in turn resulted in a connection between the Black Sea and the Mediterranean Sea (Rohling et al., 2015). As both sea level variability and exchange through the Bosphorus Strait are beyond the scope of this paper, sapropel S1 is not suitable for comparison. We therefore focus on sapropels S3, S4, and S5 in the rest of the section.
In the run with variable air temperature (Fig. 4), the modelled sapropel duration and timing is within dating uncertainty with 595 what has been found for sapropel S3 and S4 in core LC21 (Grant et al., 2016). Note that the same study finds that sapropels S1 and S5 lag precession by 2.1-3.3 kyr. This suggests that our model is capable of capturing the most relevant mechanisms for S3 and S4, but that other features not included in the model affected the timing of S1 and S5. For S5 there is evidence suggesting that the Black Sea reconnected to the Mediterranean Sea within the dating uncertainty of the onset of sapropel S5 (Grant et al., 2016(Grant et al., , 2012Wegwerth et al., 2014) It should be noted that while the model often shows a midpoint lag (relative to the insolation minimum) of a few hundred years, uncertainties related to radiometric dating methods are often larger. However, we find that midpoint lag becomes larger with decreasing strait efficiency, implying that during times with low sea level or otherwise restricted exchange, the lag might become very relevant. A prime example of such a case would be the Messinian Salinity crisis and the surrounding intervals (Topper and Meijer, 2015). Moreover, as shown in Figures 5 and 6, alternative regimes (where the freshwater budget of either 605 part of the basin, or the entire basin changes sign) can shift the midpoint of the sapropel considerably, as shown in Fig. 8.
Unlike the relatively minor shifts in midpoint resulting from only changing the phase of evaporation, the resulting difference in sapropel timing is sufficiently large to be detectable in the geological record.
De Lange et al. (2008) find that the freshening of the surface waters starts earlier and lasts longer than the suboxic bottom water conditions during S1. Our model also shows this behaviour in all of the presented runs (most notably in figures 3, 4 and 610   6). This makes sense as the DWF does not stop completely when the surface water starts to become less saline; it only reduces slowly. Furthermore, the oxygen has to be depleted for suboxic condition to occur, this is limited by oxygen consumption, and further slowed down by vertical mixing and DWF. How much longer the period of reduced sea surface salinity lasts compared to the period of suboxic bottom water likely depends on the exact location and water depth. The model is therefore mostly useful to gain insight into the mechanisms, rather than in the exact timing.

615
The regime described in Fig. 6 can show a sudden termination of the sapropel, which is similar to that seen in records of sapropel S5 (Dirksen et al., 2019 Sapropel interruptions commonly occur in the stratigraphic record (for example in S5 in core LC21, Rohling et al., 2006Rohling et al., , 2015. With slightly different settings, the sharp peak in deep water oxygen in Fig. 6E can be made to occur earlier and less intense, resulting in an interrupted sapropel. The model therefore suggests that such interruptions can occur without further external forcing. This hypothesis could be tested in the stratigraphic record by looking for evidence of a reversed freshwater budget of (part of) the basin during such interrupted sapropels, and constructing high resolution intra-sapropel age models to 630 assess the relative timing of the relevant features compared to insolation.
Each sapropel in the geological record is different, this already becomes apparent when considering the most recent 6: S1 has a different timing and is likely related to sea level variability (Grant et al., 2016;Hennekam, 2015;Grimm et al., 2015). S2 is not found at all. S4 has and interruptions in core LC21 (Grant et al., 2012), and the high resolution records of for example trace metals show very different characteristics. In the same core, the timing of the midpoint of S3 and S4 compared to insolation is the same, while that of S1 and S5 are different. S5 is much longer than all of the previous sapropels, does not have any burn down at at least one location (Dirksen et al., 2019), has an interruption in other locations (Rohling et al., 2006) and again shows generally different characteristics in different cores (Dirksen et al., 2019;Grant et al., 2016;Rohling et al., 2006). S6 again looks very different.
We conclude that both in the geological record and in the model, a typical sapropel does not exist. The timing and mechanisms 640 involved may differ considerably between sapropels and locations, as highlighted by our model results. Moreover, we find that an increase in freshwater budget alone is not sufficient to describe all key aspects of sapropel formation. An increase in atmospheric temperatures during the precession minimum (as observed in data and modelling studies, Marzocchi et al., 2015;Herbert et al., 2015) directly affects buoyancy loss during the interval in which sapropels form. This makes atmospheric temperature variability an integral feature of the system. Without it unrealistic evaporation or river outflow is needed to result 645 in sufficiently long sapropels. Our model results support this hypothesis.

Conclusions
The analysis presented in this paper illustrates that relatively simple models can give new, fundamental insights into the physical processes driving sapropel formation. The timing of sapropels relative to insolation has been widely studied in the sedimentary record. On the basis of our novel long-duration experiments we find that the timing of sapropels is very sensitive to the 650 exact climatologic and oceanographic conditions. The nonlinear response to insolation forcing implies that the sapropel record does not have a linear phase relation with insolation. The strongly nonlinear regimes in our model highlight that the midpoint of a sapropel (the average of the ages of the top and bottom of the sapropel) can be shifted significantly with a minor change in forcing, by cutting it short with a sudden termination, while during the rest of the precession cycle the response can be very similar to the nearly linear regime presented in the reference experiment. Our model results suggest that an 655 increase in freshwater input alone, as in the general hypothesis for sapropel formation, does not provide a sufficient reduction in buoyancy loss to form sapropels as they are found in the geological record. We propose that precession controlled atmospheric temperature variability also plays a key role in the process of sapropel formation.

21
Appendix A: Appendix A Figure   Competing interests. The authors have no competing interests.

Response to reviewer
We thank the reviewer for his/her thoughtful comments.
Main Comment: Comment: Line 268: Please provide justification for the assigned values of OCR and OCO .

Response:
Primarily, these oxygen consumption parameters are chosen to achieve present day deep water oxygen concentrations (close to the values found in Powley et al., 2016) under present day conditions of forcing. Furthermore, we assume that anoxic conditions can be reached in a few hundred years after stopping the circulation, in accordance with Bianchi et al. (2006). These two conditions constrain the two oxygen consumption parameters to their chosen values. We added this to line 246-249 in the revised manuscript.

Comment:
Line 336: I still feel that the mention of runs are still coming put of the blue. I appreciate the authors views on not putting the parameters into the methods but maybe a sentence or even half a sentence in the introduction can be said saying something along the lines that multiple simulations are run to test sensitivity etc.

Response:
We changed: "A low complexity model allow us to perform long runs and explore the parameter space to a much greater extent than high complexity models." (Lines 54-55 of the Introduction).
To: "A low complexity model allows us to perform many long runs and explore the parameter space to a much greater extent than high complexity models. The runs described in this paper represent but a small fraction of the experiments that have been conducted and were chosen to give an overview of the behaviour of the model." And added to the start of the Analysis and Results ("line 336"): Table 2 shows the forcing parameters that are the same for all presented runs and Table 3 gives the forcing parameters that are changed between runs. As explained in the Introduction we choose these particular runs because they illustrate well the sensitivity of the model to the various controlling factors. Figure 1: is it possible to add a scale for bathymetry? Response: Yes, we now added a scale for bathymetry. We now mention these acronyms in the text.

Other changes:
Two minor textual errors have been corrected (lines 83 and 133) . The description of an earlier heat exchange formula was repeated, we removed the repetition (lines 226-230). The units of c 13 and c 23 in Table 1 have been corrected (this does not affect the results).