Modelling the impact of anthropogenic measures on saltwater intrusion in the Weser estuary

. The Weser estuary has been subject to profound changes in topography in the last hundred years through natural variations and river engineering measures, leading to strong changes in hydrodynamics. These changes are also expected to have affected the dynamics of saltwater intrusion. Using numerical modelling, we examined saltwater intrusion in the Weser estuary in four different system states (1966, 1972, 1981, 2012). Models of each system state were set up with the respective 10 topography and boundary values. We calibrated and validated each model individually to account for differences in sediments, bedforms, and the resolution of underlying bathymetric data between historical and recent system states. In simulations of one hydrological year each with realistic forcing (hindcasting study), the influence of topography is overshadowed by the effects of other factors, particularly river discharge. At times of identical discharge, results indicate a landward shift of the salinity front between 1966 and 2012. Subsequent simulations with different topographies but identical boundary conditions (scenario 15 study) confirm that topography changes in the Weser estuary affected saltwater intrusion. Solely through the topography changes, at a discharge of 300 m³ s -1 , the position of the tidally averaged and depth-averaged salinity front shifted landwards by about 2.5 km between 1972 and 1981 and by another 1 km between 1981 and 2012. These changes are significant but comparatively small, since due to seasonal variations in run-off, the tidally averaged saltwater intrusion can vary by more than 20 km. Moreover, our analysisAn analysis of the salt flux through a characteristic cross-section has shown that saltwater 20 intrusion in the Weser estuary is primarily driven by tidal pumping and only to a lesser degree due to estuarine circulation. The shift between 1972 and 1981 is likely caused by an


Introduction
Estuaries are ever-changing systems. Natural processes and anthropogenic interventions determine the topography and conditions we observe in estuaries today. Due to the high economic importance of estuaries as shipping routes, further 35 interventions to consolidate navigation channels can be expected. In order to manage adverse effects, predictability of changes associated with the engineering measures is essential. As a response to deepening measures, significant changes in hydrodynamics have been observed in estuaries. When water depth increases, the effect of bottom friction decreases. This reduces the dissipation of tidal energy, resulting in a larger tidal amplitude and often in phase shifts in the ebb and flood current durations and velocities (Winterwerp et al., 2013;. At the same time, changes in mixing processes can 40 occur, possibly affecting the length of saltwater intrusion (Grasso and Le Hir, 2019). Saltwater intrusion is directly linked to water quality and sediment transport processes (Burchard et al., 2018) and thus needs to be monitored.
The effect of topography changes on saltwater intrusion depends on the physical mechanisms that control salt transport in estuaries. Most important are the net advection, driven by freshwater discharge, tidal asymmetry effects, and the estuarine circulation, i.e. the tidally averaged estuarine exchange flow. One driver of the estuarine circulation is the combination of the 45 seaward barotropic (depth-independent) and the landward baroclinic (increasing with depth) pressure gradient, which results in a seaward flow of estuarine water near the surface and landward flow of dense water near the bed. This is known as the gravitational circulation (Geyer and MacCready, 2014). Depending on the estuary´s geometry and tidal forcing, strain-induced periodic stratification (SIPS) can occur with stratification during ebb and mixing during flood (Simpson et al., 1990). This tidal asymmetry enhances estuarine circulation (Jay and Musiak, 1994) and has been referred to as tidal straining circulation 50 (Burchard et al., 2011;Geyer and MacCready, 2014) or eddy viscosityshear covariance (ESCO, Dijkstra et al., 2017). Salt transport is also controlled by e.g. lateral momentum advection (Burchard et al., 2011) and estuarine convergence (Ianniello, 1979;Burchard et al., 2014). The strength of salt transport mechanisms depends on the water depth and it is therefore expected that channel deepening affects saltwater intrusion (Andrews et al., 2017).
Studies have been conducted in estuaries worldwide to try to quantify the impact of changes in topography, i.e., channel 55 deepening or widening, on saltwater intrusion. In recent studies, numerical models with different topographies have been used to examine mixing and transport processes in estuaries before and after implementing engineering measures. Among others, such investigations have been conducted for the San Francisco estuary, US (Andrews et al., 2017), the Hudson River estuary, US , the Danshui River estuary system, Taiwan (Liu et al., 2020), the Seine estuary, France (Grasso and Le Hir, 2019), and the Ems estuary, Germany (van Maren et al., 2015). 60 Generally, it was found that a deepening of the river channel is associated with a landward shift of the brackish water zone.
Several studies explain the increase in landward salt transport with an increase in the estuarine circulation and a decrease in vertical mixing processes, which occur due to the larger water depths (Grasso and Le Hir, 2019;Andrews et al., 2017;van Maren et al., 2015). In the case of the Danshui river estuary, some channels have been deeper in the predevelopment state, which goes along with increased gravitational circulation and a further landward limit of the saltwater intrusion in the 65 predevelopment state (Liu et al., 2020). All the other above-mentioned estuaries are deeper in the present state. Grasso and Le Hir (2019) investigated key estuarine processes in the Seine estuary in 1960 and 2010 and detected a relative increase in gravitational circulation and stratification. This caused an up-estuary shift of the salinity front and changed the response of saltwater intrusion to tides and discharge.  examined the Hudson River estuary in the present state and the predevelopment state. They detect an increase in saltwater intrusion and stratification, but only a minor change in 70 estuarine circulation and no change in the response of saltwater intrusion to river discharge. The authors conclude that dredging in the Hudson did not significantly change estuarine exchange processes.
In most of the abovementioned studies, a numerical model of the respective estuary in the present state was set up and calibrated. Different model topographies were subsequently inserted to represent earlier states of the estuary Grasso and Le Hir, 2019;Liu et al., 2020). Models representing earlier system states of the estuary were not 75 calibrated. However, differences in bed roughness are expected to occur between the system states due to sediment redistribution and changes in bed forms, which are usually not resolved in the models. Van Maren et al. (2015) calibrated models of the Ems estuary for different system states and obtained significantly larger roughness values with historical bathymetries compared to the present state. This is attributed to the observation that sediment in the Ems estuary has successively become finer in the last decades. The authors conclude that changes in bed roughness can strongly contribute to 80 shifts in hydrodynamics and transport processes. In addition to shifts in bed roughness, the resolution of data which are used to generate model topographies may differ for the system states, which can have an effect on form drag. When setting up models with different topographies, individual calibration of each model might thus be advantageous.
In the Weser estuary, no model-based examinations have yet been conducted to examine the effect of anthropogenic measures on saltwater intrusion. Instead, saltwater intrusion and influencing parameters have been tracked by measurements and it was 85 attempted to separate the impact of different factors (Krause, 1979;Grabemann et al., 1983). However, many factors affect saltwater intrusion on different time scales and it has thus not been possible to isolate the impact of topography variations (Grabemann et al., 1983). Within each tidal cycle, the position of the brackish water zone shifts by more than 15 km along the navigation channel of the Weser (Kösters et al., 2014). Changes in tidal components, meteorological conditions, long-term variations of the salinity of the North Sea, and variations in discharge also affect the position of the brackish water zone 90 (Grabemann et al., 1983). Due to variations in discharge, the tidally averaged saltwater intrusion shifts by more than 20 km, whereby the impact of discharge variations is larger for low-discharge than high-discharge conditions (Kösters et al., 2014).
This study aims to systematically quantify saltwater intrusion in a real tidal estuary in different system states and to determine to what degree anthropogenic measures, in particular deepening measures, impact saltwater intrusion in the estuary. As an example, we studied the Weser estuary in the German Bight and set up numerical model simulations in four different system 95 states (1966,1972,1981,2012). In contrast to most previous studies, we individually calibrated each model and conducted simulations with realistic and with idealized forcing. In this paper, we describe the model setup and simulation results. We discuss the methodology, the importance of calibration for the model-based examination with different model topographies, and the effect of channel deepening on processes governing saltwater intrusion.
2 Study site: the Weser estuary 100

Geomorphology and hydrology
The Weser estuary, located in Northern Germany, is of high ecological and economic importance. It is divided into two sections: the Lower Weser and the Outer Weser (see Fig. 1). The Lower Weser stretches from the tidal weir at Bremen-Hemelingen (km 0, the tidal limit) to Bremerhaven (km 66.7). The funnel-shaped Outer Weser starts from Bremerhaven and opens towards the North Sea (Weser-km 126). The estuary is an important shipping channel, providing access to the container 105 terminal of Bremerhaven, the port of Bremen, and several smaller ports. In 1941-2015, the mean annual discharge, measured at station Intschede, was 321 m³ s -1 with high seasonal variability (NLWKN, 2018). In the same period, the mean low discharge  The semidiurnal tidal wave from the North Sea propagates through the Weser estuary in about three hours. At the northern 110 part of the Outer Weser, the tidal range is 2.8 m. On its way through the estuary, it increases to 3.8 m near Bremerhaven due to estuarine convergence, slightly decreases between Weser-km 50 and 30, and increases to 4.1 m near Bremen. In most stretches of the navigation channel of the Lower and Outer Weser, ebb currents are slightly stronger than flood currents. Flood currents dominate in some stretches, such as Weser-km 80-95 (Lange et al., 2008).
Before mixing with seawater, the river Weser has an initial salinity of about 0.5 psu as a caustic potash solution is discharged 115 further upstream. The mean position of the 2 psu isohaline is located at Weser-km 45 at the reversal from flood to ebb (high water slack) and Weser-km 60 at the reversal from ebb to flood (low water slack) (Lange et al., 2008).

Historical development
The topography of the Weser estuary has been strongly altered through river deepening and correction measures since the end of the 19 th century. The main objective of these measures was to establish and maintain a continuous shipping route with 120 sufficient depth and width for ships of increasing size to pass through (Wetzel, 1988).
Since 1960, three major deepening measures have been conducted. The Outer Weser was modified in 1969-1971 to create a continuous depth of nautical chart datum 'Seekartennull' (SKN) -12 m and in 1998-1999 to increase the depth to SKN -14 m (Lange et al., 2008;Wetzel, 1988). Engineering measures at the Lower Weser were conducted in 1973-1978. Thereby, the stretch between Brake (Weser-km 39) and Bremen (Weser-km 0) was deepened to nautical chart datum SKN -9 m, and the 125 stretch between Bremerhaven and Nordenham to SKN -11 m (Lange et al., 2008;Wetzel, 1988   conducted in the navigation channel (Lange et al., 2008). The bathymetry height along the navigation channel in the four system states is depicted in Fig 2a.  130

Methods
We built numerical models of the Weser estuary in four system states: 1966, 1972, 1981, and 2012. A simulation of one hydrological year was conducted for each model with realistic forcing to hindcast hydrodynamics and saltwater intrusion ('hindcast study'). Results from previous calibration runs with respective topographies were used as initial conditions and a spin-up time of four weeks was applied. Additionally, simulations of each system state, but with identical forcing, were 135 conducted and analysed for one spring-neap cycle (after six weeks spin-up time) so that the impact of topography changes and roughness changes could be individually evaluated ('scenario study').

Description of the numerical model
Numerical simulations are based on UnTRIM² (Casulli, 2008). UnTRIM² solves the Reynolds-averaged Navier Stokes equations, the continuity equation, and transport equations based on a semi-implicit finite volume / finite difference approach 140 to calculate current velocities, surface elevations, and tracer concentrations (Casulli and Walters, 2000). The equations are solved on a horizontally unstructured grid with vertical z-layers of 1 m thickness. The model considers 3D hydrodynamics, daily freshwater discharge from the Weser, and the transport of salinity and heat. Vertical turbulent mixing was estimated with a two-equation k-ε model; horizontal mixing was modelled with a constant viscosity (0.1 m² s -1 ) and diffusivity (0.1 m² s -1 ) in all simulations. UnTRIM² was coupled with the sediment transport model SediMorph (BAW, 2002), but only to calculate 145 bottom roughness. Sediment transport was neglected in this study. The model is well-suited to simulate flows in estuaries, as the algorithm accurately preserves mass where wetting and drying occurs (Casulli, 2008). For example, it has been used to simulate flows in the San Francisco Bay (Andrews et al., 2017), the German Bight (Hagen et al., 2021b;Rasquin et al., 2020), and the Weser Estuary (Kösters et al., 2014).

Model topographies 150
The model area includes the Lower Weser, the Outer Weser, and the adjacent Jade (see Fig. 1). Tributaries of the Weser and their discharge are not represented, as bathymetric and hydrographic data of the tributaries are not available for all system states. The same unstructured orthogonal calculation grid was used in all models. It contains cells of different sizes, arranged to increase the resolution in areas of interest. Cell spacing in the navigation channel is 50-250 m.
Model topographies were created by interpolating respective topographical data on the calculation grid, smoothening 155 transitions, and correcting essential landscape features. The accuracy of the data sets differs depending on the quality and density of the original data. The model topography of 2012 is based on airborne laser scanning (ALS) for flat areas, multibeam echo-sounding for deep channels, and single-beam echo-sounding for shallow areas. Older bathymetries are mostly based on single-beam echo-sounding measurements. For most of the Outer Weser, topographical maps at a scale of 1:20.000 containing tracks with 70-200 m distance were used. For the Lower Weser, measurements along cross-sections with 125 m longitudinal 160 spacing were used. Topographical measurement data were not available for all regions for each time period. For example, historical data for the northern Outer Weser were only available from nautical charts of 1970 with lower data density. With temporal and spatial interpolation, considering deepening measures, all data were compiled to historical topographies of the Lower and Outer Weser (see Fig. 1, hatched area) in different historical system states (BAW, 2020(BAW, , 2021. We interpolated the bathymetric data from 1966, 1972, and 1981 on the computational grid preserving the wet volume, which ensures that the 165 water volume in the model agrees exactly with the average depths of the bathymetric data. We then filled the remaining area of the North Sea and Jade with bathymetric data from 2012. Subsequently, we corrected important landmarks and structures such as summer dikes, side channels, constructions, and the transition between Outer Weser (historical bathymetric data) and Jade (bathymetric data of 2012) in the case of the model topographies 1966, 1972, and 1981. Due to the differences in temporal and spatial data coverage, the quality and level of detail of the model topographies differs. 170 For example, bedforms such as dunes landwards of Weser-km 55 (Lange et al., 2008) are represented in the model topography of 2012 as variations in depths (see Fig. 2c). Historical recordings were less detailed, i.e. model topographies of 1966, 1972, and 1981 do generally not contain information on bedforms. Hence, the older topographies are smoother than 2012 (see Fig. 2b). We account for the different resolutions of the underlying bathymetric data and their effect on form drag by individually calibrating the models of each system state with the respective topographies (see Sect. 3.6). 175

Initial sediment distribution
We prescribed an identical initial sediment distribution in all models based on diverse sediment samples as described in Milbradt et al. (2015). The sediment distribution represents the system state 2012; data from previous years were additionally used at locations with insufficient data coverage. The distribution comprises eight fractions of sediments: very coarse sand, coarse sand, medium sand, fine sand, very fine sand, coarse silt, medium silt, and fine silt. We use the same initial sediment 180 distribution in all models as historical data to suitably represent changes in sediment distribution over time were not available.
In the models, the initial sediment distribution is used for roughness calculation based on sediment type and bedform prediction (see Sect. 3.6).

Model forcing
As boundary conditions, we prescribed time series for salinity, temperature, and water level at the seaward boundary to the 185 North Sea. We also provided salinity, temperature, and discharge of the Weser at the landward boundary. The hindcast study aimed to represent the system states as realistically as possible. Measurement values were retrieved from the Waterways and Shipping Authorities Bremerhaven, Bremen and Hannoversch Münden and time series were constructed by linear interpolation between observations. Measured values were available for most parameters, but not for water level and salinity at the seaward model boundary in the North Sea in 1966Sea in , 1972Sea in , and 1981 and for salinity and temperature at the landward boundary in 1966. 190 For the water level however, historical records for tidal high and low water levels were available for the station 'Leuchtturm Alte Weser' (ALW, see Fig. 1) close to the model boundary. Thus, we generated a synthetic time series of 1965-2012 by reconstructing the astronomical tide at station ALW, fitting the signal to measured high water and low water values and inducing a phase shift and amplitude amplification to account for the distance of the station to the model boundary. Salinities at four positions along the open boundary were approximated utilising neural networks of two layers (55 and 11 neurons) with 195 a Levenberg-Marquardt algorithm. For each position, 100 networks were trained with salinity values of 1996-2016, which were extracted from a validated model of the German bight, the Easy GSH model (Hagen et al., 2021a). As reference data, the tidal range and tidal mean water level at station ALW, salinity records at Helgoland station (Wiltshire et al., 2008), and the discharge of the Weser and Elbe were used. Subsequently, the correlation of network results and the extracted target salinity values were calculated and the best network for each position was selected. With the selected networks, salinity at each position 200 was predicted for all system states . For evaluation of the performance of neural networks, a skill estimator was applied according to Murphy (1988), implemented in the form suggested by Ralston et al. (2010).
where Sriver is salinity and Q is discharge. For temperature, measured values of 1968 were used instead of 1966, as the variations over the years were assumed to be similar and temperature is only of secondary importance in this study.
In contrast to the hindcast study, we used identical boundary values in all simulations for the scenario study. The synthetic 220 time series of the hindcast model 2012 was used for water level at the open boundary. For all other parameters, cross-scenario boundary values were generated; representing an average of the four system states (see Table 1).

Analysis methods
Model results were analysed by calculating tidal characteristic values, which describe the tidal curve and facilitate 230 characterization of the system's behaviour and comparison between simulations (Lang, 2003). Thus, for example, the tidal range and the minimum, mean, and maximum salinity per tidal cycle were calculated in each scenario and compared.
In addition, the position of the brackish water zone was calculated for all model results by determining the position of the 5 psu and 20 psu isohaline along the navigation channel (see Fig. 1) based on the tidal mean of the vertically averaged salinity. The saltwater intrusion length was defined as the distance from the estuary mouth (Weser-km 126.2, see Fig. 1) to the 5 psu 235 isohaline along the navigation channel. In order to identify processes driving the saltwater intrusion the following analysis methods have been applied.

Regression analysis
The effect of discharge on the saltwater intrusion length was quantified by means of a regression analysis based on data from the hindcast simulations. Saltwater intrusion was calculated for each tidal cycle within the respective hydrological years and 240 discharge conditions were assigned. Data points were then sorted into 50 m³ s -1 bins of discharge between 150-400 m³ s -1 and 100 m³ s -1 bins between 400-1,100 m³ s -1 and averaged to reduce the effect of additional factors and outliers. If fewer than 10 entries were available for a category, these entries were excluded. Based on the resulting average intrusion length for the specified discharge categories, multiple nonlinear regression with interaction effects was performed with the Levenberg-Marquardt algorithm to obtain trends of the form 245 * = , with saltwater intrusion length L * , factor k, discharge Q, and exponent m. The analysis was repeated without data aggregation to examine the method's validity. 250 We decomposed the total salt flux Ftot into three components: barotropic flux Fbf, estuarine exchange Fexf, and tidal pumping flux Ftp. The method is based on the decomposition analysis of sediment fluxes of Becherer et al. (2016). We focused our analysis on the brackish water zone and have chosenapply their method to analyse salt flux through a cross-section in the Lower Weser (Weser-km 60, see Fig. 1) to carry out1), located in the brackish water zone of the estuary. We conducted the 255 salt flux decomposition. We for equidistant points along the cross-section (distance of 25m) and integrated the results over the whole width of the cross-section. In a first step, we interpolated our data on a σ-layer grid with M=50 equidistant layers.
Then we calculated the temporal average of salinity s and horizontalcurrent velocity in the main flow direction u at σ-layers Where X is the respective quantity, h is the height of the σ-layer, and T is the analysis period. Temporal averages are represented by angle brackets. Becherer et al. use a moving average window to track effects of variable forcing on the flux components.
As the simulations are run with constant river discharge as forcing, and only contain variations in water level from the North 265 Sea (see Table 1, scenario study), we can simply average over the entire analysis period of the scenarios (one spring-neap cycle). We calculated deviations from the mean, indicated by curly brackets, with 270 Vertical averages were calculated according to From the temporal and vertical averages and deviations, we determined the salt flux components and integrated over the cross-275 section's width W.
The barotropic flux Fbf contains residual barotropic flows due to river discharge, barotropic ebb-flood asymmetries, and wind stress. The vertical exchange flow Fexf indicates the estuarine circulation and the tidal pumping flux Ftp indicates upstream transport by the correlation of salinity and velocity fluctuations (Becherer et al., 2016). 285

Potential energy anomaly
According to Simpson , the potential energy anomaly Φ is the energy required to homogenise the water column with given density stratification instantaneously. Burchard and Hofmeister define it as where H is the mean water depth, η is the water level elevation, D = η+H is the actual water depth, g is the gravitational acceleration, and ρ is the density. Thus, an increase in potential energy anomaly indicates an increase in stratification assuming a constant water level..

Calibration and validation 295
Models of each system state were calibrated in two steps by adjustment of the roughness settings: first, tuning of a form roughness predictor, and second, definition of an additional roughness in the Lower Weser. The simulation time for each system state was four weeks. Bottom roughness in our models is described by Nikuradse´s effective roughness coefficient ks, which relates to the bed roughness length z0 with z0 ≃ ks/30 (Malcherek, 2010). In nature, roughness can be seen as the sum of grain and form related roughness. The sediment transport model SediMorph, which is coupled with the hydrodynamic model 300 UnTRIM², calculates the grain roughness ks,grain at each element from the prescribed sediment distribution with ks,grain = 3dm, where dm is the mean grain size. Form roughness is estimated at each time step of the simulation based on sediment properties, water depth, and velocity, according to van Rijn (2007). Calibration factors, which determine the prevalence and importance of bedforms (ripples, mega ripples, dunes), were varied in the calibration to optimise the agreement between simulated and observed water levels, if available (system state 2012), high water and low water values (all scenarios), and the tidal range (all 305 scenarios). The effect of dunes could not be reliably modelled with the predictor but had to be prescribed as an additional roughness landwards of Weser-km 55. The additional roughness also compensates model effects, such as due to the omission of the Weser's tributaries. Tributaries typically dampen the tidal wave through low water depths and high energy dissipation rates. The additional roughness increases from Weser-km 55 towards Weser-km 26 and decreases afterwards again. In this way, damping of the tidal wave could be well represented. Roughness as calculated by the numerical model (effective bed 310 roughness) along the navigation channel is described in Table 2. For each model, different roughness settings were determined.
While the estimated form roughness is in the range of 0.01 to 0.04 m, the maximum additional roughness ranges from 0.22 m (system state 2012) to 0.36 m (system state 1981). The mean RMSE of tidal range in scenarios 1966, 1972, and 1981 improved by almost 50 % after adjustment of the roughness settings by individual calibration (see Table 2).
In other estuaries, changes in sediments and bedforms were observed over time. For example, dredging induced fine sediment 315 accumulation in the Ems (Winterwerp et al., 2013;van Maren et al., 2015) and reduced drag associated with bedforms in Columbia river (Jay et al., 2011), both contributing to lower friction and an amplification of tides. There are no data to assess whether there have been comparable changes in sediment inventory in the Weser estuary. But the height of bedforms will almost certainly have changed with changing water depths. Such changes are not represented in our models, which could explain to some extent differences between scenarios. 320 In addition, effects induced by the different resolutions of the original bathymetric data (see Sect. 3.2) are balanced out through roughness calibration. To some extent, dunes are depicted in the numerical gridsif present in the bathymetric data. In the examined years prior to 2012 bathymetric data almost did not resolve dunes. This has an effect on depth variation in the model topographies and the associated form drag. The lower depth variation in the older model topographies is compensated by defining a higher additional roughness. The largest additional roughness was required in the model of system state 1981 (see 325 Table 1 Fig. 1). An overall bias of 5 cm indicated a slight, consistent overestimation of water levels. A 345 visual inspection of observed and modelled water levels suggested a correct phasing of tides, which is corroborated by the low RMSE. Saltwater intrusion was slightly overestimated, with higher computed values than measured values in Lower Weser (see Fig. 3b). The intratidal variation in Weser-km 60-100 was slightly underestimated, with higher minimum salinity and lower maximum salinity per tide. The mean RMSE in hindcast model 2012 was 1.4 psu for mean salinity, averaged over all stations with available measured values. Overall, the magnitude and dynamics of saltwater intrusion were reproduced well in 350 the hindcast model 2012. It is expected that the model performance will be similar for the other system states.

Natural variability of saltwater intrusion in the Weser estuary
In order to give an idea of temporal variation of salinity in the Weser estuary measurement data of the hydrological year 2012 are shown in Fig. 4. In each semi-diurnal tidal cycle, the 5 psu isohaline moves up-and downstream by almost 20 km (see 355 Fig. 4a). It is located about 2 to 3 km further landward during spring tide compared to neap tide. The mean position of the 5 psu isohaline per semi-diurnal tidal cycle shifts by more than 30 km within the year. A major seaward shift occurs in December and January, where the isohaline moves from Weser-km 40 to Weser-km 73. This is induced mainly by an increase in discharge from 100 m³ s -1 to 1000 m³ s -1 .

Figure 4. Intratidal and seasonal variation of saltwater intrusion in the Weser estuary. The water level at station Leuchtturm Alte Weser (ALW), discharge at station Intschede (INT), and salinity along the navigation channel of the Weser estuary are displayed on different time scales. Salinity values were interpolated linearly between measurement stations with available measurements along the navigation channel. The 5 psu isohaline is displayed in white. On top, salinity variation within one day, September 1 st , 2012, is depicted (a). Below, variation within the hydrological year 2012 of the mean salinity per tidal cycle is depicted (b). We applied a moving median filter to the water level data in (b) to show the seasonal variation (filter size 12h 25min).
Hydrodynamics and saltwater intrusion in the Weser estuary were simulated in hindcast models of 1966, 1972, 1981, and 2012. These models contain respective model topographies and realistic forcing (see Table 1) and they were individually calibrated (see Sect. 3.6). The mean tidal range in the estuary is up to 0.6 m larger in 1981 and 2012 than in 1966 and 1972 (see Fig. 3a). This reflects the impact of the deepening measures in Lower Weser in 1973-1978, which had a strong effect on tidal dynamics. Saltwater intrudes further in 1972 and 2012 compared to 1966 and 1981 (see Fig. 5). In 1972 and 2012, the 365 5 psu isohaline is positioned about 10 km more landward. The reason for these large differences is particularly high discharge values in 1966 and 1981, which were on average about twice as high compared to the other system states (see Table 3). Not all variation can be explained with the discharge conditions. Compared to 1972, there was higher discharge in 2012; however, the position of the salinity front is about 3 km more landward. In this case, other influencing factors outweigh the impact of discharge. The mean tidal level at station ALW in system state 2012 was considerably higher compared to 1972, indicating 370 different meteorological conditions and reflecting an increase in mean sea level (Wahl et al., 2013). Moreover, the tidal range at station ALW was larger in 1981 and 2012, possibly due to topography changes (Hubert et al., 2021) and in response to the nodal tide. Along with these factors, changes in topography might affect the differences in the saltwater intrusion.

Identifier
Mean salinity

Effect of discharge on saltwater intrusion in different system states
The sensitivity of saltwater intrusion to discharge can be described with analytical expressions, which depend on dominant 380 salt flux mechanisms and estuary shape. With a regression analysis (see Sect. 3.5.1), we examined the influence of discharge on the saltwater intrusion length in all hindcast simulations. Amongst the simulations, discharge and other parameters differ (realistic forcing, see Table 1). Two versions of the analysis were conducted: with and without pre-processing by data aggregation. Data from the hindcast simulation of 1972 were excluded because the range of discharge values was insufficient for the analysis. 385 In the scenarios representing 1966, 1981, and 2012, we identified a clear relationship between discharge and saltwater intrusion length according to Eq. (43), with L~Q -0.15 . In system state 2012, the saltwater intrusion length increases by approximately 9 km if the discharge decreases from 1000 m³ s -1 to 380 m³ s -1 or from 380 m³ s -1 to 150 m³ s -1 . The trend lines for 1966, 1981, and 2012 have comparable gradients in the logarithmic plot (comparable exponents m). A significant difference occurs between the exponent in HY1966 compared to the exponent in HY1981 (p<0.001) and HY2012 (p=0.005) in the regression analysis 390 with non-aggregated datapoints (Fig. 6b).
The shift in trend lines indicates that saltwater intrusion increased between 1966 and 1981 and between 1981 and 2012 for similar discharge conditions (see Table 4). For example, with 300 m³ s -1 discharge, the saltwater intrusion increased by approximately 2.4 km between 1966 and 1981 and another 1.8 km between 1981 and 2012. This trend could be linked to differences in topography, i.e. the deepening measures conducted in-between the system states. 395  4. The average landward shift of the 5 psu isohaline in the hindcast simulations, based on regression analysis (Fig. 6b). 1966-19811981-2012Overall shift 1966

Impact of topography changes on saltwater intrusion 400
The effect of topography on saltwater intrusion was further examined in simulations with identical boundary values but different model topographies (scenario study). One set of simulations (1966,1972,1981,2012) Fig. 7b). According to that, if the navigation channel is deepened by 1 m in Weser-km 40-110, the position of the 5 psu isohaline shifts landwards by about 2 km. Deepening measures in the Outer Weser seem to induce an overall extension of the brackish water zone. Between scenario 1966Between scenario and 1972Between scenario and between scenario 1981 and 2012, when deepening measures in the Outer Weser were conducted, the length of the brackish water zone increased by 1.1 km and 1.8 km. In contrast, an extension of only 0.5 km occurred between 1972 and 1981 (see Fig. 7a). An extension of the brackish water zone means that salinity increases more gradually towards the sea.

Relevance of roughness calibration for the examination of saltwater intrusion in different scenarios
To evaluate the impact of the roughness calibration on the results, we compared simulations with different model topographies  Fig. 7a). However, simulations with identical roughness settings do not adequately reproduce tidal energy, e.g. the tidal range in 1981 is larger than in 2012 425 contrary to observations. The erroneously increased tidal range exceeds the effect of changes in other processes, so that the limit of saltwater intrusion shifts seawards, and not landwards, between 1981 and 2012.

Impact of topography changes on tidal effects and estuarine circulationsalt transport processes
Topography changes can influence tides in the estuary, e. g. the tidal range may increase due to increasing depths. For the Weser estuary, this effect is confirmed in the simulations with identical forcing. For example, the tidal range at Weser-km 60 430 increases from 3.65 m in 1966, 3.67 m in 1972, 3.86 m in 1981, to 3.88 m in 2012. This affects saltwater intrusion. Along with tidal effects, other salt transport processes such as the gravitationalestuarine circulation can change through topography changes. In order to disentangle individual processes, the salt flux through a characteristic cross-section was analysed in scenario simulations of 1966, 1972, 1981, and 2012 were analysed.. The scenarios contain identical boundary values but different topographies and roughness settings. (see Table 1, scenario study). 435 For a better understanding of the salt transport processes, weWe decomposed the total salt flux through a cross-section in the Lower Weser (Weser-km 60, see Fig. 1) over the period of one spring-neap cycle for each scenario into advectivebarotropic flux, tidal pumping flux, and the estuarine exchange flow (see Sect. 3.5.2). Results of all scenarios (see Fig. 8 Note that the total salt flux contains the effect of the river discharge. The river discharge is in all scenarios constant with a salinity of 1.2 psu (see Table 1, scenario study), contributing to the total salt flux with a seaward salt flux of 300 m³ s -1 * 450 1.2 psu = 360 psu m³ s -1 . The remaining total flux reflects variability in salt storage, which is expected as the simulation is not in a steady state, i.e. the water level at the North Sea boundary contains both tidal components and surgeNote that the scenarios contain constant river discharge forcing but spring-neap tidal variations. Thus, the position of the brackish water zone varies on the spring-neap cycle but does not show a river discharge driven shift in position. One would expect the total salt flux to be Results indicate that the seaward volume transport is in all scenarios higher than the landward volume transport, reflecting the river discharge input (see Fig. 9). The seaward salt transport is 2 to 4 % higher than the landward salt transport. This corresponds, when adding Qinsin and Qoutsout, to a residual salt transport (total salt flux) of 388 (1966), 359 (1972), 418 (1981) and 369 psu m³s -1 (2012). This approximately matches the total salt fluxes determined in the salt flux decomposition. 465 Landward and seaward volume transports increase between 1972 and 1981, probably due to deepening measures in Lower Weser and the associated increase of the tidal volume. However, the ratio of landward to seaward volume transport slightly increases between 1972 and 1981 from 0.91 (1966 and 1972) to 0.92 (1981 and 2012). As salt transport in the Lower Weser is primarily governed by tidal processes (see Sect. 4.6.1), this induces an increase in salt transport and salinity between 1972 and 1981. An additional increase in salt transport and salinity occurs between 1981 and 2012; this cannot be explained by the 470 volume transports.

Discussion
Saltwater intrusion is subject to high natural variability and salinity measurements of the past decades are not always available or exhaustive. This makes it difficult to quantify the effect of anthropogenic measures based on measurements only. With numerical simulations, we analysed the impact of anthropogenic measures on saltwater intrusion in two ways: First, we evaluated saltwater intrusion in hindcast-models of the hydrological year 1966, 1972, 1981 and 2012 and analysed it as a 485 function of discharge. Second, we evaluated saltwater intrusion in a scenario study with identical boundary values.
In the model setup, a key challenge was the lower frequency and coverage of historical measurements compared to the present state, i.e. for water level, salinity, and temperature (see Sect. 3.4). This can lead to differences in model accuracy and compromise the comparability of model results. In our case, a detailed validation of the historical models was not possible due to the lack of historical salinity measurements. With a systematic approach and equal treatment of the respective system states, 490 wherever possible, we created consistency and minimized bias. For example, we generated boundary values of all models with identical methods and calibrated each model individually.
In our study, an individual calibration of roughness parameters proved essential (see Sect. 3.6). Without calibration, changes in sediments and bedforms and effects of different resolutions of bathymetric data strongly influenced model results and lead to false conclusions (i.e. seaward shift instead of landward shift due to topography changes between 1981 and 2012). It has to 495 be noted that calibration is only possible if observational data are available to support a robust calibration. Recalibration with imprecise data of low resolution might even negatively affect model results. In addition, it might not always be necessary, depending on the representation of roughness and sediments, and the resolutions of bathymetric data and grid. However, we recommend that a thorough evaluation of bathymetric data and, if needed and possible, calibration of each model should be conducted when results of models with different topographies are to be compared. 500 The analysis of hindcast results showed that, at similar discharge conditions, there is a landward shift in the position of the brackish water zone between 1966-1981and 1981-2012. The exact distance of the shift depends on the discharge condition (see Table 4). When comparing saltwater intrusion at times with a discharge of 300 m³ s -1 , there is a 2.4 km shift between 1966 and 1981 and a 1.8 km shift between 1981 and 2012. In the scenario study with identical boundary conditions, i.e. Q=300 m³ s -1 , the 5 psu isohaline shifts by 2.5 km between 1972 and 1981 and by 1 km between 1981 and 2012 505 (see Sect. 4.4). Both analyses show that anthropogenic measures in the Weser estuary led to an increase in the saltwater intrusion. This is in agreement with findings for other estuaries, where the navigation channel was also deepened (Andrews et al., 2017;Grasso and Le Hir, 2019;. In the hindcast study, additional factors such as tides and salinity in the North Sea influence saltwater intrusion. In addition, the evaluation of saltwater intrusion depending on discharge does not account for adaptation time to changing discharge conditions. Nevertheless, results of both the hindcast study and the 510 scenario study indicate an overall shift in the same range, 3.5-4 km. Measures in the Lower Weser, which were conducted in 1973-1978, had the highest impact. However, the exact extent of saltwater intrusion will also be influenced by factors not included in the model such as changes in tributaries, side channels, and the construction of waterways infrastructure. The influence of discharge on saltwater intrusion follows the relationship L~Q -0.15 in our simulations (see Sect. 4.3). Thus, changes in discharge have a higher impact for low discharge than high discharge conditions, in agreement with previous studies 515 (Kösters et al., 2014).  examined the relationship between discharge and saltwater intrusion in the Hudson river in idealized scenarios with different discharge conditions. For both, bathymetries representing the present state and the predevelopment state, the sensitivity of the saltwater intrusion to discharge follows L~Q -0.28 . The large difference from the sensitivity found in our study could be related to differences in bathymetry and transport processes between the estuaries.
While saltwater intrusion into the Hudson is dominated by estuarine exchange flow, we show that the Weser is dominated by 520 tidal pumping (see Sect. 4.6.1). Following , a change in the sensitivity of saltwater intrusion to discharge could indicate a change in salt flux mechanisms. We detect a significant, but very small difference in the sensitivity in 1966 compared to 1981 and 2012 (see Sect. 4.3). However, the exact exponent depends on the definition of the saltwater intrusion length. Moreover, not only discharge, but also other influencing factors vary in the examined simulations and lag times in the adjustment of saltwater intrusion to changed conditions were not considered in our analysis. Therefore, the small 525 difference in sensitivity to discharge between 1966 and 1981 cannot be clearly attributed to a change in salt flux mechanisms.
The impact of anthropogenic measures on processes governing saltwater intrusion was further examined based on results of the scenario study (idealized forcing) (see Sect. 4.6). For all scenarios, a decomposition of the salt flux through a cross-section in the Lower Weser was conducted according to Becherer et al. (2016). The salt flux was divided into contributions from advectivebarotropic flux, tidal pumping, and the estuarine exchange flow. Results indicate that salt flux in the Weser estuary 530 is dominated by tidal pumping contributions and advectivebarotropic flux. This also explains the result that saltwater intrusion is further landward during spring tides compared to neap tides (see Sect. 4.1). In contrast, saltwater intrusion in estuaries dominated by estuarine exchange flow decreases during spring tides. Grasso and Le Hir (2019) found that the effect of the spring-neap cycle on saltwater intrusion can change if the governing processes in the estuary change. This aspect could be further examined in subsequent studies. A more thoroughextensive study disentangling the role of different processes on 535 seasonal timescales in Weser estuary is presently carried out by Becker (personal communication, July 4, 2022). The estuarine exchange flow through the cross-section strongly increases between 1972 and 1981, while tidal pumping increases between 1981 and 2012. In both cases, this is outbalanced by a simultaneous increase of the barotropic flux. A possible interpretation could be that deepening measures in the Weser estuary locally induce an increase in estuarine exchange flow and globally induce an increase in tidal pumping. In addition, the shape of the estuary and the depth distribution are complex and the impact 540 of a deepening measure depends on its location, meaning a deepening of the outer estuary as between 1981 and 2012 has a different impact than a deepening of the inner estuary as between 1972 and 1981. This aspect could be further explored in future studies, which should also include the analysis at different cross-sections and an analysis of lateral processes. In the outer Weser estuary, which is characterized by wide cross-sections and a two-channel system (Gundlach et al., 2021), lateral variations could significantly influence salt transport processes (Wei et al., 2022). In addition, other decomposition methods 545 could be tested (e.g. Dyer, 1974;Dronkers and van de Kreeke, 1986;Park and James, 1990;Hughes and Rattray, 1980). Additionally, the exchange flow through the same cross-section was analysed according to Knudsen , as described in Burchard et al. . Results confirm a small increase in landward and seaward volume transports (and thus in exchange flow). Representative inflow and outflow salinities are increasing substantially over the years, such that also the salt transports exhibit a strong increase, especially between 1972 and 1981. The residual salt transports determined in the analysis approximately match the 550 total salt fluxes determined in the salt flux decomposition.
An analysis of the potential energy anomaly (PEA) on a transect along the navigation channel of the Weser further indicates that stratification increases between 1972 and 1981 in the Lower Weser due to the deepening of the Lower Weser in [1973][1974][1975][1976][1977][1978]. It also shows the high variability of tidal mixing along the longitudinal profile, which points the interaction of global trends and local effects. In case of the PEA, local effects could comprise e.g. an increased stratification at especially deep 555 stretches. It has to be noted that the location of the transect is identical in all scenarios and there is deviation from the thalweg on some stretches in case of the historical system states, which could affect the results. These limitations might be overcome by aggregating data over the entire navigation channel rather than analysing data along one transect..

Conclusions
Saltwater intrusion in the Weser estuary is highly variable and dependent on natural forcing factors, i.e. river discharge, as 560 well as anthropogenic impacts, i.e. channel deepening. A systematic study of the effect of topography changes on saltwater intrusion was carried out to disentangle these overlapping influencing factors. This study used numerical models of four system states (1966, 1972, 1981, and 2012) with respective topographies to hindcast the historical development of saltwater intrusion (hindcast study) and examine the effect of anthropogenic measures (scenario study). The models were individually calibrated and validated to account for differences in sediments, bedforms, and the resolutions of underlying bathymetric data. 565 In the hindcast simulations, the influence of topography is overshadowed by the effect of other factors, particularly discharge.
The salinity front is about 10 km more landward in 1972 and 2012 compared to the other system states, as the average discharge is almost double in these years. However, at similar discharge, a landward shift of the salinity front is indicated between 1966 and 1981 and between 1981 and 2012. This trend is confirmed in the scenario study. With identical boundary values (i.e. Q=300 m³ s -1 ), the salinity front shifts landwards by about 2.5 km between scenario 1972 and 1981 and by about 1 km between 570 scenario 1981 and 2012. It has to be noted that the exact distance by which the salinity front shifts depends on the discharge due to the nonlinear relation between discharge and saltwater intrusion. Analyses of the salt flux components, the exchange flow, and stratification indicate that the processes governing saltwater intrusion changed only by through a small amount in response to anthropogenic measures. However,cross-section in the Lower Weser, an increase in estuarine circulation is detected in response to deepening measures in the Lower indicate that the Weser in 1973-1978 and an increase in estuary is 575 dominated by barotropic flux and tidal pumping in response to deepening in the Outer Weser in 1998-1999. , while the estuarine exchange flow is of secondary importance in all scenarios.