Baltic sea deep salinity: an initial and boundary value problem

Abstract

The Baltic Sea salinity is modelled using a two-box model. The simplistic approach allows for very long integrations where a large part of the phase space of the model can be probed. Particular emphasis is put on the salinity dynamics in the deeper parts of the sea and how they are affected by boundary and initial conditions. Multiple statistically steady states, corresponding to forcing from different years, are examined and the route to them through the model’s phase space is traced out. The model is forced with freshwater fluxes and sea level variations at its boundary. The respective roles of these two forcing terms is investigated using a factorization technique, and it is found the sea level variability is the dominant one for the deep salinity dynamics. The role of natural variability is also examined, and the probability for deep salinity changes for different forcing years is computed.

1 Introduction

Fjord type estuaries are water bodies, often narrow and channel like, that are separated from the ocean by a sill. The sill inhibits water exchange, especially renewal of deep water as the sill depth is often much shallower than the depth of the fjord (Howe et al., 2010). The sill is thus a natural constriction that makes it difficult for high density ocean water to penetrate into a fjord. Deep water renewal in fjords through intrusions of high density ocean water is hence often an intermittent process, and stagnation periods (periods without deep water renewal) can last many years. The limited ventilation causes environmental problems. Most notably deep water oxygen depletions (hypoxia) that makes the deeper parts of many fjords uninhabitable to higher forms of marine life (Conley et al., 2009; Meier et al., 2011; Silva and Vargas, 2014).

Fjords are remnants of prior glaciations and are thus situated in areas where precipitation and run-off typically exceeds evaporation. Many fjords can be described approximately as two layer systems, with a brackish surface layer on top of deep layer of high saline water of oceanic origin. Here we will discuss and model the Baltic Sea as such a two-layer system, but the model could easily be applied also to other fjords and many conclusions are general enough to apply to fjord type estuaries, rather than just to the Baltic Sea.

The Baltic Sea can be thought of as an oversized fjord, or a system of multiple nested fjords. It is separated from the open ocean by the narrow and shallow Danish Straits. Deep water renewal is highly intermittent and hypoxia is commonplace. Deep water renewal events in the Baltic, often called major Baltic inflows (MBIs) have gotten a lot of attention in the scientific literature owing to their importance, in particular, for alleviating hypoxia (Mohrholz et al., 2015; Mohrholz, 2018; Krapf et al., 2022). Most of these investigations have focused on the meteorological conditions needed for inflows to occur, here we refer to them as boundary conditions or simply as forcing. Our focus in the present investigation is, however, not only on the boundary conditions, but also on initial conditions. In particular on the role of the salinities in the upper and lower layer when inflow favouring boundary conditions occur.

To study the roles of initial and boundary conditions simultaneously in a statistical sense requires a very large data set and dedicated experiments. Far beyond what’s available from observations or what can be integrated with comprehensive three dimensional ocean models. For the purpose of this study we have therefore produced such a dataset using a two-box model of the Baltic Sea. Using such a crude representation of ocean dynamics comes with both perks and limitations. The main perks are the low computational cost and the high interpretability. The experiments we do here simply could not be replicated with a comprehensive model. Moreover, interpreting the results would be much harder in models with richer ocean dynamics. The main limitation is of course that it is hard, or even impossible, to know how well the response of a highly simplified system corresponds to that of the much more complex ones governing real world systems.

These caveats aside it seems evident that much of the important dynamics involved in deep water renewal can be modelled quite accurately also with simple models and also that much can be learnt from such models (Stigebrandt, 1987; Kõuts and Omstedt, 1987; Gustavsson, 2000). Our interest here is primarily in how the response of the Baltic Sea salinities to boundary conditions varies depending on the initial conditions in the basin. This is studied using integrations where boundary conditions are applied randomly and initial condition dependent distributions of their effects are estimated. Moreover, we trace out multiple statistically steady states, where a single boundary condition, corresponding to freshwater and sea level forcing from a given year, is repeated a great many times to produce a statistically steady state response.

Ultimately, our aim is to get a deeper understanding of ventilation dynamics. In particular for salinity dynamics of two layer fjord systems, which the Baltic Sea is a prominent example of. Important questions such as: what is the likelihood of natural variability giving rise to a stagnation period of length X years cannot easily be answers using either observations or comprehensive models, nor can more exotic questions about the nature of statistical steady states or the ranges of natural variability easily be examined using such data. Our approach combines the use of a classical box model framework with some diagnostics that are novel at least in the field of physical oceanography, and to our knowledge have never been applied to study deep water renewal in two layer systems.

2 Modelling

2.1 Model description

The Baltic Sea is modelled using a simple 2-box model with three state variables (, and ), see Figure 1. The upper box corresponds to the water masses situated above 70 m and the lower box to those below. Exchange with the open ocean is modelled as a straightflow, and whether water intruding from the open ocean ends up in the upper or lower box is determined by a very simplistic plume model.

FIGURE 1

The barotrophic exchange with the Kattegatt is modelled, following Stigebrandt (1980), Omstedt (1993) and several other authors as,where is the volume transport between the Kattegatt and the Baltic Sea. , where is the haline contraction coefficient and is the sill depth of the sounds connecting the Baltic Sea and Kattegatt. is the sea level in Kategatt, that in the Baltic and is a flow resistance parameter (Stigebrandt, 1980; Mattsson, 1996). is the salinity of the Kattegatt water, while is the salinity of the upper layer water in the model. A summery of all parameter values and their units is given in Table 1.

We use as a forcing term, while is a state variable of the model that is governed bywhere is the net freshwater flux owing to evaporation, precipitation and runoff, and is the area of the basin.

The equations for and are derived from salt and mass conservation. A detailed derivation is offered in the appendix. Here we give only the resulting evolution equations. With parametrizations for the different terms the equation becomes,The equation for meanwhile is

Going more into the particulars of each term, the overbars indicates a time average over 7 days, while the time step of the model is 1 day. This time averaging is used as often times the same water runs back and forth across the entrance to the Baltic Sea, with little effect on the Baltic salinity. Using a time averaged flow in the salinity equation is thus a way of eliminating this effect from the salinity equation (i.e., to stop water slushing back and forth in the entrance area from changing its salinity instantly from to ). The effect on the salinity equations is similar to the use of a buffer zone for inflows as was implemented by Mohrholz (2018). The first term in these equations is salt diffusion between the layers, is a diffusion coefficient and the square root comes from assuming (Stigebrandt, 1987). The second term models the influence of run off, precipitation and evaporation on the salinity of the upper box.

The function, , is defined throughand is used as a continuous version of an if statement and the parameter controls the sharpness of the switch. For our chosen value, works essentially as an if statements and is applied to the sign of and the salinity differences in the model. The height over the bottom of the inflowing dense current is given by, , and controlled byThe 2/3 power is there to mimic a parabolic channel and the constant sets the height of the inflowing bottom current (a large gives a thick current). The constant is chosen so that a large inflow gets a height of about 20 m. is the height gain of the dense current that occurs owing to entrainment of water from the upper layer on transit to the deep box. It is chosen so that large inflows consist of about equal amount of Kattegat and entrained upper layer water. The salinity of the plume (bottom current), is given byThe plume model is extremely simplistic: all dense inflows entrain the same amount of upper layer water . This has the effect that weak inflows get diluted and end up in the upper layer, while strong inflows may achieve as the relative amount of Kattegatt water is much larger in those and they can consequently end up in the lower layer.

The two advective terms in the salinity equations can thus be understood as follows: the first one, term three in Equation 3, imports Kattegatt waters into the upper layer and is active when , the second advective term, term four in Equation 3 and term two in Equation 4 imports a mixture of Kattegat and upper layer waters into the deep layer and simultaneously exports an equivalent volume of deep water to the upper layer. This term is only active when , that is when the intruding waters are dense enough to replace the current deep waters in the Baltic Sea.

Taken together these equations form an admittedly very crude representation of the Baltic Sea. In essence, a considerable degree of realism is sacrificed to obtain simplicity. Nevertheless, the model captures much of the essential dynamics controlling Baltic Sea inflows as is demonstrated in Sect. 2.3, and it’s simplicity and integrations speed makes it very suitable to study the influence of boundary and initial conditions on the salinity of the Baltic Sea, which is our primary aim here.

2.2 Model forcing

The aforementioned box model has only two forcing terms: boundary conditions on freshwater and sea level denoted and respectively. The sea level forcing, , comes from sea level observations taken at the station Viken during the years 1977–2018. The mean of is set to zero each year to remove the long term trends owing to post-glacial rebound and sea level rise. The net freshwater flux is taken from the Baltic Sea Model Intercomparison Project (BMIP) (Gröger et al., 2022) for the same period. The term is derived through summation of all river runoff that ends up in the Baltic Sea and it is constrained to have a 1977-2018 average of 15,000 , a value taken from Mattsson (1996). Basin parameters and are calculated from the hypsometry used for the NEMO-Nordic model (Hordoir et al., 2019).

2.3 Model hindcast

First of all, even though the runs described here are refereed to as hindcasts, the objective of the study is not to tune the model so that it accurately replicates observed salinities. A multitude of models that can do that, much more accurately than a 2-box model, are already available (Gröger et al., 2022). In fact, very little tuning has been done, and what we’re after in terms of fidelity is simply that for a given sequence of forcing, the response of the box model is in a broad sense similar to that of the real system. A demand that could be satisfied with very little tuning.

The three state variables from a 1977–2018 hindcast with the model are shown as blue lines in the top three panels of Figure 2. The box model is able to represent many features of the salinity evolution in reasonable agreement with those observed in the eastern Gotland basin of the Baltic Sea at station BY15. In particular, it is noteworthy that many of the inflow events, more precisely their related deep salinity spikes, are well captured by the model. However, we also note that there are spikes occurring in the model that do not show up in observations and vice versa. This, is of course partly owing to the simplistic representation of the physics in the model, and the very crude forcing used. However, similar biases are also seen in full 3D ocean models with much more realistic forcing and model dynamics (Hordoir et al., 2019). In section, 4 we make the case that inflow related salinity spikes have a strong sensitivity to initial conditions.

FIGURE 2

As a demonstration of this, the top three panels of Figure 2 also show a second hindcast in orange lines, identical to the first, but with its initial value reduced by 0.2 g/kg. This change clearly affects the salinity evolution, in particular, in the lower box throughout the whole integration. Note also that is almost identical between the two runs, so it is not the exchange flow with the Kattegat (i.e., the volume of the inflows) that changes. This is true also if one changes the initial condition for , which has virtually no effect on any of the three state variables beyond the first few days (not shown). To what extent this sensitivity is mirrored in more comprehensive 3D Baltic Sea models is not known, but we do know from experiments with the NEMO-Nordic model that inflow related salinity spikes can be made to disappear or to appear by changing the initial salinity conditions. Getting accurate long model reconstructions of deep Baltic Sea salinities thus appear to require not only excellent forcing and model dynamics, but also very good initial conditions. A consequence of the spin up time scale of the system being comparable to the typical length of hindcast simulations.

It is noteworthy that the salinity variations in the model hindcast are much smaller than those at the station. This is primarily a reflection of the modelled salinities being averages over the whole Baltic Sea, while the observations are at fixed depths at a single station, situated in the eastern Gotland basin of the Baltic Sea. In particular, the station BY15 is often used to study inflows since those are known to have a large effect on the deep salinities at that station. Stations further north are less affected by inflows as they are harder to reach for a dense current emanating from the Kattegatt. This is illustrated in the bottom panel of Figure 2, where the Baltic Sea mean salinity is compared to that in a recent hindcast with the 3D ocean model NEMO-Nordic (Arneborg et al., 2025). From that panel it is clear that when volume averages are compared, the interannual variability is reasonably well described and the correlation is high. The two curves in the bottom panel have a Pearson correlation coefficient of 0.8. The box model, however, has a bias in its mean salinity and a damped decadal variability. The latter can likely, at least in part, be attributed to the fact that the mean of is set to zero for all forcing years, so sea level variations on time scales longer than a year are not used to force the model. Meier H. E. M. and Kauker F. (2003) attributed the decadal salinity variability between 1902 and 1998 partly to variability in the freshwater forcing and partly to decadal variations of the large-scale sea level pressure over Scandinavia. Such sea level pressure variations would affect , but have here been filtered out.

3 Statistically steady states and the routes toward them

Alternating the forcing years naturally gives rise to variability in the state variables of the model. In this section we study the statistically steady states that the model reaches when a single forcing year is repeated a great many times. Note that a steady state is defined using the state vector of the model , as a state in which . Since our boundary conditions are never steady the model never reaches a steady state. Statistically steady states are states such that a sufficiently long time average of becomes zero. Repeating a single yearly forcing long enough puts the model in a statistically steady state. More specifically, repeating any of our forcing years long enough gives a model solution that is a periodic orbit, with a period of 1 year. Note that a periodic orbit is a statistically steady state as the average over years of , when goes through a periodic orbit with a period of 1 year. In terms of yearly means of the model finds a fixed point solution (steady state) in space after sufficiently long integration times if a single yearly forcing is repeated.

This is illustrated in Figure 3, where and are shown at the end of a 2000 years long simulation with repeated yearly forcing from either 1987 or 1993. These particular years are chosen because the year 1987 has a very low probability of having an inflow, while 1993 has a very high probability, see Sect. 4. Nevertheless, the periodic orbits these two simulations reach are qualitatively similar. Both feature deep inflows, two of them occur in each year with the 1987 forcing and a single one with 1993 forcing. Thus, in a sense, the forcing does not control the variability as much as it does the mean state. The mean state is much saltier with 1993 forcing than with 1987 forcing. The periodic orbits are reached when the deep salinity is just high enough that the yearly deep salinity increase owing to inflows exactly matches the yearly diffusive salt loss from the lower layer to the upper one. For this to work, years with a high inflow probability must reach a state with high mean salinity, while the opposite is true about years with a low inflow probability. Note also that it is not the exchange flow, , that takes a long time to equilibrate, but the salinities in the different boxes. The time rate of change of owing to inflows depends essentially on , in which equilibrates very fast and much slower. The coupling between having a large volume inflow and a large salinity change is thus quite weak.

FIGURE 3

The paths toward the periodic orbits from different initial conditions are shown in Figure 4. Here the forcing from year 1990 is used, but the figure is similar for all forcing years (not shown). The figure is made by starting simulations from initial conditions at a given radii in space away from the periodic orbit. That is, from circles in space around the origin taken to be the periodic orbit corresponding to the 1990 boundary conditions. Points on this circle with negative salinities and with are excluded.

FIGURE 4

It is clear from the figure that the model prefers to approach its periodic orbit from two main directions. When is too large it tends to first equilibrate and then enter a state of no inflows, with slowly decreasing . When both and are too low the tendency is instead to increase them in fixed proportions. A third interesting branch occurs when is too high, but is not too far from its preferred value for the given forcing. Here we see overshoot, the reason is that this is an inflow spike favouring situation: high upper layer salinity gives a small dilution owing to entrainment which leads to high values that increases the likelihood that the dense current ends up in the lower box rather than in the upper one.

The time it takes the model to reach its periodic orbit is illustrated in colours in the top two panels. For the most part the model reaches its periodic orbit within a few decades, in agreement with spin-up time scales reported for a range of Baltic Sea models (Meier M. and Kauker F., 2003; Feistel et al., 2006; Omstedt and Hansson, 2006). An exception to this, best seen in the middle panel, is when the model’s initial condition has a much too high . Under such circumstances the equilibration occurs through diffusive salt loss from the lower layer, a process that can take several centuries for the range of initial values used in this experiment.

The lower panel of Figure 4 shows that there is a near linear relationship between and in the periodic orbits for the different forcing years. Since the difference between and determines the diffusive flux between the two layers. This also shows that the salt increase in the lower layer owing to inflows during a year is similar for all forcing years. Again suggesting that the forcing controls the mean state in a sense more than it does the variability.

4 Long integrations with random forcing years

In this section we consider a very long run with randomly drawn boundary conditions. An integration of length 164,000 years is performed where the forcing year is drawn randomly from the boundary conditions for the 1977-2018 period. That is, at the start of every new year a random forcing year is chosen so that a random sequence of forcing years is used to force the model. Essentially, this is done to mimic the randomness of the atmospheric forcing driving the Baltic Sea. In particular, we wish to study the range of the models response to a single forcing year that depends on what forcing years preceded it. That is, the sensitivity to initial conditions. In total each forcing year is run around 3,900 times, a consequence of them being given equal probability of occurrence. Figure 5 shows an attractor in space for this simulation. The colouring shows frequency with which the model visits different regions of the phase space. The attractor is truncated in two senses: the third state variable is omitted and the finite length of the run suggests that some low probability areas of the phase space that can be visited have not yet been so.

FIGURE 5

It is evident from the figure that there is a positive correlation between and , as can be expected given that straitflow dynamics control the salinity import into both layers. It is also clear that with the given forcing years, the model has a preferred state with and slightly in excess 8.4 and 12.2 g/kg respectively, but with a considerable range. The standard deviations of the daily and values are 0.12 and 0.07 g/kg respectively. The and values of the periodic orbit reached with the repeated 1987 forcing in Figure 3 is well within the range of those manifested in this experiment, while those for repeated 1993 conditions are far outside. That is, conditions such as those seen in the repeated 1993 integration belong to a very low probability part of the phase space that is not visited in this 164,000 years long integration.

Figure 6 shows the probability density for yearly salinity changes in the lower layer for different years. For each forcing year a distribution of of the last day of the year minus of the first day of the year is thus computed. It is clear from the figure that a vast majority of years are neither salinity spike or non-salinity spike years, but can be either depending on the years that preceded it. Another noteworthy observations is that these distributions are all based on around 3,900 years (the amount of time each year is drawn), however, we know that all these forcing years if repeated long enough give rise to a periodic orbit. That is, each year would have a distribution that contains zero on the x-axis in Figure 6 if the experiment was ran for an infinite time. However, in our 164,000 years long experiment there are distributions, like those for years: 1991,1993 and 2017 that have all their probability on the positive side of the x-axis. Meaning that exceptionally high salinities are needed to not get a salinity spike during those years. These are thus years with a very high probability of spikes, and one would have to start from very unlikely initial conditions to not get an increasing salinity during such years. The strongest spike year in this simulation is 1993, which is also one of the strongest in our observed history. Another noteworthy feature is that 2014, which is the year in our observed history with the largest deep salinity change at BY15, is not forcing wise one of years with the most spike favouring conditions. Again suggesting that forcing alone is a poor determinant of deep salinity changes, as initial conditions also play a considerable role. Another interesting year is 2003, which in reality was a year with a large salinity spike, but in our model the mode of its distribution is actually situated on negative salinity changes. Our forcing is almost certainly too crude to accurately represent all conditions that makes a year spike favouring, but it also seems likely that the observed salinity spike of 2003 could have been considerably smaller had the year not been preceded by a stagnation period.

FIGURE 6

The long integration can also be used to estimate probabilities for having stagnations periods (here defined as periods with little or no salinity increase in the deep water) of different lengths. Figure 7 shows the yearly probability for having stagnations period of length up to 30 years. Two different thresholds on salinity changes in the deep water are used, , and . With a stagnation period is defined as a consecutive period where never increases, while for we allow the salinity difference between the beginning and the end of any year within a stagnation period by a maximum of 0.01 g/kg. The latter threshold is more in line with how stagnation periods are defined in the real Baltic Sea, where small increases in salinity are usually not taken to mean a stagnation period has ended. Indeed, small inflow related salinity increases are rather common as is evidenced in Figure 2.

FIGURE 7

The threshold has considerable influence on the stagnation period probabilities, especially for longer stagnation periods. Having a decade long stagnation period, for example, has a yearly probability of only 0.01 with the threshold and one of 0.17 with the threshold. It is not possible to compare these probabilities to similar ones for the real Baltic Sea since observational series are much to short for such statistics to be computed. Moreover, picking an appropriate thresholds for model to observational comparisons is also challenging. One notable comparison one can make, however, is to the work of Schimanke and Meier (2016). Schimanke and Meier (2016) examined an 850 years long integration with a 3D ocean model, forced with atmospheric variables from a dynamically downscaled global climate model experiment intended to mimic natural variability, for periods of prolonged salinity decrease at BY15. They found that a 10 year period with a salinity decrease comparable to that which occurred between 1983 and 1992 has a return period of about 100 years. This is the same return period we find for a 10 year stagnation period with the condition. However, one should not overemphasize this similarity given the considerable differences both in the experiments and the criteria for stagnation periods used in the current study and that by Schimanke and Meier (2016). The general shape of these curves with its rapidly decreasing yearly probability with increasing period length, however, is almost certainly a property not just of models, but also of the real Baltic Sea. In the present model experiment, the yearly probability of having a 30 years stagnation period starting is only 0.004 with the threshold, while with the threshold no such period occurs in the whole integration.

5 Forcing factorization

In this section we probe the respective role of the two forcing terms using a factorization technique developed by Stein and Alpert (1993). The method separates the role of the forcing terms into their direct effects and their interactions. In our case where we have only two forcing terms (freshwater and sea level) we get a total of three factors: one directly related to , one directly related to and a interaction term between and .

The function we factorize (here called ) is the probability density for yearly salinity changes in the lower layer for different years (i.e., the same as shown in Figure 6). The effects of the forcing terms are separated through experiments where the time dependent forcing is exchanged with climatological daily averages over the 1977-2018 period. The full factorization system is given belowThe hat operator (ˆ) signifies factorized variables. That is, is the probability density for yearly salinity changes in the lower layer for a run forced with both and , while shows the contribution to that is only due to the interaction between and . , since the run where both forcings are climatological reduce to a periodic orbit (i.e., all probability is situated in a very small interval around zero), but is left in the equations for completeness. To sum it up, is the part of that is only due to sea level variability in the Kattegatt, is the part of only owing to variability in the freshwater forcing and is the interaction between the two forcing terms.

Figure 8 shows the factorization for four different years, chosen because they depict different conditions. The year 2017 was chosen because it depicts somewhat “standard” conditions. Here we find that the role of the interaction term is essentially to negate the effect of the freshwater forcing, while having a modest effect in the range of where the sea level forcing is effective. The interaction term in fact negates the effect of the freshwater forcing for all years (not shown), demonstrating that the sea level forcing is the dominant term. However, there are years where the interaction term gives a significant perturbation also in the range of where the sea level forcing is effective. Two such examples are shown in years 1998 and 2002. In both cases, one of the two forcing terms favours a deep salinity spike while the other does not, and in both cases we see that while is still much closer to than to , is significantly perturbed away from in the direction of . Showing that also the freshwater forcing has a non-negligible effect on the probability of having a deep salinity spike. The year 1993 was chosen just to show what an extreme outlier is was in terms of sea level forcing. There are other years that also feature very long tailed distributions, but none of them also have their mode very far to the right.

FIGURE 8

Figure 9 shows how the factorization for the forcing year 2017 changes under parameter variations. Six different variations are considered. The diffusion coefficient and the flow resistance parameter are either doubled or halved and in one experiment both and are doubled. In addition, an experiment is performed where . That is, in addition to drawing a random forcing for each year, each drawn year’s is also perturbed by a number drawn from a uniform distribution between . This perturbation to introduces additional variability in , which has the effect of greatly spreading out the distributions of . In this experiment we also see a new peak in and at values close to . This peak is a consequence of the distribution being bounded in its left tail by the diffusive salt loss to the upper layer. That is, drawing a large could give rise to a large if other conditions also favour that. However, drawing a small does not lead to a small , since years with small are years when no inflows end up in the lower box. is then a consequence of diffusion between the boxes, which is only weakly dependent on . The effect of doubling (halving) is similar to that from halving (doubling) , which can be seen both in their respective panels and in the similarity between the factorization for the case and the original 2017 factorisation shown in Figure 8. The effects of these two perturbations is straightforward to understand. Increasing increases the diffusive salt loss from the lower box to the upper, enabling larger salinity spikes. The effect of halving is to increase the variance in , which similarly enables larger salinity spikes in the lower box.

FIGURE 9

6 Discussion & conclusion

A very simple two box model was applied to study Baltic Sea salinity dynamics. In particular, focus is put on how inflow related salinity spikes depend on initial and boundary conditions. The amplitude of these salinity spikes depend not only on the magnitude of the inflow (i.e., the volume transport into the Baltic Sea), but also on the salinities in the two box system. While that is quite elementary, it is worth mentioning and the fact that a large inflow transport does not always lead to a large deep salinity spikes have also been discussed in relation to observationally based time series (Mohrholz, 2018). In reality, there is a further complication that is not represented in our model, namely, that inflowing volume and inflowing salt mass are related by a time dependent salinity, not a constant as used in this model. Loptien et al. (2025) makes this point in a recent article, where they come to the conclusions that MBIs come in two different flavours depending on the atmospheric conditions that forced them, and that the salinity of the inflow is significantly higher in one flavour than in the other. Regardless of that simplification and many others, one may, rather, generally conclude that a deep salinity spike implies an inflow, but not that an inflow implies a deep salinity spike. Moreover, the amplitude of the spike is strongly dependent on the initial conditions, and is thus not, in itself, a good proxy for the salt transport of the inflow.

What’s most novel in this study is that the simplicity of the model allows us to study multiple statistically steady states and a very large range of natural variability, something that cannot be done with either observations or comprehensive 3D ocean models. That is, large parts of the state space of the model can be probed and very unlikely states can be examined. An interesting finding relating to this is that the periodic orbits reached when a single forcing year is repeated long enough are dynamically similar, but have very different mean salinities. Essentially, they all find a deep salinity just low enough that the yearly deep inflow related salinity increase matches the yearly diffusive salt loss to the upper layer. It is possible that parameter or parametrisation variations could give rise to longer than yearly orbits, however, we found no evidence of such in this study. Another interesting behaviour is the models preference to approach its periodic orbit from two different directions. If the initial is too high the model tends to first fix to the level of the orbit and then gradually approach the final value from a fixed value. This route is consequence of inflowing plumes not reaching the deep box when is very high. When the initial is instead too low the tendency is instead to approach the periodic orbit along a strait line in space. A consequence of inflows entering the lower box and lower box water simultaneously being heaved into the upper box.

The very long integration with random forcing is intended to mimic natural variability. In the real Baltic Sea there has been discussions about whether recent stagnation periods are the result of natural variability or anthropogenic changes (e.g., Stigebrandt, 1992; Almroth-Rosell et al., 2021; Lehmann et al., 2022). Given the simplicity of the model used here one cannot give a definite answer to that based on our results. However, what one can say is that the system we model is very capable of putting itself into stagnation periods just through natural variability and that decade long stagnation periods are not terribly uncommon, which has also been shown by Schimanke and Meier (2016). Another key finding is that we can characterise the different forcing years according to their propensity to give inflow related salinity spikes. A surprising result of that investigation is that most forcing years can give rise to a very wide range of lower layer salinity changes, depending on the initial conditions in the basin. Only a few are almost always giving salinity spikes, and even those years show a wide range for the amplitude of the spike. The forcing factorisation experiments further showed that sea level forcing is the primary driver of the significant variability in the deep salinity spikes, while the freshwater forcing has a more modest effect. This holds true also under parameter variations. The largest spike in recorded history occurred in 2014 (Mohrholz, 2018), which in our analysis is a year with a clear spike preference, but not one of the most extreme. The factorisation for 2014 (not shown) also suggests that this was a year whose propensity to give a large salinity spike was lessened by the freshwater forcing. Another example (also not shown) is 2003, where a large salinity spike occurred at BY15, but not in our hindcast. Our factorization showed that this year had a highly spike favouring freshwater forcing, but a sea level forcing that did not favour spikes.

To what degree this characterisation of forcing years is applicable to the real Baltic is, of course, questionable owing to the simplicity of the dynamics represented by our model and the simplified forcing. We have, for example, not employed explicit wind forcing. Nevertheless, we have seen in experiments tailored to tune the comprehensive NEMO-Nordic model (Hordoir et al., 2019), where the 2014 inflow was ran starting from different initial conditions, that the magnitude of the salinity spikes vary considerably there as well and for some initial conditions does not appear at all. Suggesting that a similar sensitivity to initial conditions is not just a model artefact, but a property also of the real system. A property that makes long term salinity projections in the Baltic Sea a difficult undertaking, and warrants a discussion about spin-up strategies. In particular, in the context of model intercomparision projects, such as BMIP for the Baltic Sea (Gröger et al., 2022).

References

• 1

  Almroth-RosellE.HanssonM.Vã¤liG.EilolaK.AnderssonP.et al (2021). A regime shift toward a more anoxic environment in a eutrophic sea in northern Europe. Front. Mar. Sci.8. 10.3389/fmars.2021.799936

  • CrossRef
  • Google Scholar

• 2

  ArneborgL.HieronymusM.PembertonP.LiuY.FredrikssonS. T. (2025). Response of a semi-enclosed sea to perturbed freshwater and open ocean salinity forcing. EGUsphere2025, 1–25. 10.5194/egusphere-2025-4735

  • CrossRef
  • Google Scholar

• 3

  ConleyD. J.BjörckS.BonsdorffE.CarstensenJ.DestouniG.GustafssonB. G.et al (2009). Hypoxia-related processes in the baltic sea. Environ. Sci. & Technol.43, 3412–3420. 10.1021/es802762a

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  FeistelR.NauschG.HagenE. (2006). Unusual Baltic inflow activity in 2002-2003 and varying deep-water properties. Oceanologia48, 21–35.

  • Google Scholar

• 5

  GrögerM.PlackeM.MeierH. E. M.BörgelF.BrunnabendS.-E.DutheilC.et al (2022). The baltic sea model intercomparison project (bmip) – a platform for model development, evaluation, and uncertainty assessment. Geosci. Model Dev.15, 8613–8638. 10.5194/gmd-15-8613-2022

  • CrossRef
  • Google Scholar

• 6

  GustavssonB. G. (2000). Time-dependent modeling of the baltic entrance area. 2. Water and salt exchange of the baltic sea. Estuaries23, 253–266. 10.2307/1352831

  • CrossRef
  • Google Scholar

• 7

  HordoirR.AxellL.HöglundA.DieterichC.FransnerF.GrögerM.et al (2019). Nemo-nordic 1.0: a NEMO-based ocean model for the Baltic and north seas – research and operational applications. Geosci. Model Dev.12, 363–386. 10.5194/gmd-12-363-2019

  • CrossRef
  • Google Scholar

• 8

  HoweJ. A.AustinW. E. N.ForwickM.PaetzelM.HarlandR.CageA. G. (2010). “Fjord systems and archives: a review,” in Fjord systems and archives (geological society of London).

  • Google Scholar

• 9

  KõutsT.OmstedtA. (1987). Water cooling in the entrance of the baltic sea. Tellus A Dyn. Meteorology Oceanogr.39, 254. 10.3402/tellusa.v39i3.11758

  • CrossRef
  • Google Scholar

• 10

  KrapfK.NaumannM.DutheilC.MeierH. E. M. (2022). Investigating hypoxic and euxinic area changes based on various datasets from the baltic sea. Front. Mar. Sci.9, 823476. 10.3389/fmars.2022.823476

  • CrossRef
  • Google Scholar

• 11

  LehmannA.MyrbergK.PostP.ChubarenkoI.DailidieneI.HinrichsenH.-H.et al (2022). Salinity dynamics of the baltic sea. Earth Syst. Dyn.13, 373–392. 10.5194/esd-13-373-2022

  • CrossRef
  • Google Scholar

• 12

  LoptienU.RenzM.DietzeH. (2025). Major baltic inflows come in different flavours. Commun. Earth Environ.6, 232. 10.1038/s43247-025-02209-0

  • CrossRef
  • Google Scholar

• 13

  MattssonJ. (1996). Some comments on the barotropic flow through the danish straits and the division of the flow between the belt sea and the Öresund. Tellus A Dyn. Meteorology Oceanogr.48, 456–464. 10.3402/tellusa.v48i3.12071

  • CrossRef
  • Google Scholar

• 14

  MeierH. E. M.KaukerF. (2003). Modeling decadal variability of the baltic sea: 2. role of freshwater inflow and large-scale atmospheric circulation for salinity. J. Geophys. Res. Oceans108. 10.1029/2003JC001799

  • CrossRef
  • Google Scholar

• 15

  MeierM.KaukerF. (2003). Sensitivity of the Baltic sea salinity to the freshwater supply. Clim. Res.73, 231–242. 10.3354/cr024231

  • CrossRef
  • Google Scholar

• 16

  MeierH. E. M.AnderssonH. C.EilolaK.GustafssonB. G.KuznetsovI.Mã¼ller-KarulisB.et al (2011). Hypoxia in future climates: a model ensemble study for the baltic sea. Geophys. Res. Lett.38. 10.1029/2011GL049929

  • CrossRef
  • Google Scholar

• 17

  MohrholzV. (2018). Major Baltic inflow statistics †revised. Front. Mar. Sci.5, 384. 10.3389/fmars.2018.00384

  • CrossRef
  • Google Scholar

• 18

  MohrholzV.NaumannM.NauschG.KrügerS.GräweU. (2015). Fresh oxygen for the baltic sea †an exceptional saline inflow after a decade of stagnation. J. Mar. Syst.148, 152–166. 10.1016/j.jmarsys.2015.03.005

  • CrossRef
  • Google Scholar

• 19

  OmstedtA. (1993). Deep water exchange in the baltic proper. Tellus A Dyn. Meteorology Oceanogr.45, 311–324. 10.3402/tellusa.v45i4.14895

  • CrossRef
  • Google Scholar

• 20

  OmstedtA.HanssonD. (2006). The baltic sea ocean climate system memory and response to changes in the water and heat balance components. Cont. Shelf Res.26, 236–251. 10.1016/j.csr.2005.11.003

  • CrossRef
  • Google Scholar

• 21

  SchimankeS.MeierH. E. M. (2016). Decadal-to-centennial variability of salinity in the Baltic sea. J. Clim.29, 7173–7188. 10.1175/JCLI-D-15-0443.1

  • CrossRef
  • Google Scholar

• 22

  SilvaN.VargasC. A. (2014). Hypoxia in chilean patagonian fjords. Prog. Oceanogr.129, 62–74. 10.1016/j.pocean.2014.05.016

  • CrossRef
  • Google Scholar

• 23

  SteinU.AlpertP. (1993). Factor separation in numerical simulations. J. Atmos. Sci.50, 2107–2115. 10.1175/1520-0469(1993)050<2107:fsins>2.0.co;2

  • CrossRef
  • Google Scholar

• 24

  StigebrandtA. (1980). Barotropic and baroclinic response of a semi-enclosed basin to barotropic forcing from the sea. Boston, MA: Springer US, 141–164. 10.1007/978-1-4613-3105-6_5

  • CrossRef
  • Google Scholar

• 25

  StigebrandtA. (1987). A model for the vertial circulation of the baltic deep water. J. Phys. Oceanogr.17, 1772–1785. 10.1175/1520-0485(1987)017<1772:amftvc>2.0.co;2

  • CrossRef
  • Google Scholar

• 26

  StigebrandtA. (1992). Bridge-induced flow reduction in sea straits with reference to effects of a planned bridge across Öresund. Ambio21, 130–134.

  • Google Scholar