Estimation of the natural background of phosphate in a lowland river using tidal marsh sediment cores

Elevated phosphate (PO4) concentrations can harm the ecological status in water by eutrophication. In the majority of surface waters in lowland regions such as Flanders (Belgium), the local PO4 levels exceed the limits defined by environmental policy and fail to decrease, despite decreasing total phosphorus (P) emissions. In order to underpin the definition of current limits, this study was set up to identify the pre-industrial background PO4 concentration in surface water of the Scheldt River, a tidal river in Flanders. We used the sedimentary records preserved in tidal marsh sediment cores as an archive for reconstructing historical changes in surface water PO4. For sediment samples at sequential depths below the sediment surface, we dated the time of sediment deposition and analysed the extractable sediment P. The resulting time series of sediment P was linked to the time series of measured surface water-PO4 concentrations (data 1967–present). By combining those datasets, the sorption characteristics of the sediment could be described using a Langmuir-type sorption model. The calibrated sorption model allowed us to estimate a pre-industrial background surface water PO4 levels, based on deeper sediment P that stabilized at concentrations smaller than the modern. In three out of the four cores, the sediment P peaked around 1980, coinciding with the surface water PO4. The estimated pre-industrial (∼ 1800) background PO4 concentration in the Scheldt River water was 62 [57; 66 (95 % CI)] μgPO4-PL. That concentration exceeds the previously estimated natural background values in Flanders (15–35 μgTPL−1) and is about half of the prevailing limit in the Scheldt River (120 μgPO4-PL). In the 1930s, river water concentrations were estimated at 140 [128; 148] μgPO4-PL, already exceeding the current limit. The method developed here proved useful for reconstructing historical background PO4 concentrations of a lowland tidal river. A similar approach can apply to other lowland tidal rivers to provide a scientific basis for local catchment-specific PO4 backgrounds.

Abstract. Elevated phosphate (PO 4 ) concentrations can harm the ecological status in water by eutrophication. In the majority of surface waters in lowland regions such as Flanders (Belgium), the local PO 4 levels exceed the limits defined by environmental policy and fail to decrease, despite decreasing total phosphorus (P) emissions. In order to underpin the definition of current limits, this study was set up to identify the pre-industrial background PO 4 concentration in surface water of the Scheldt River, a tidal river in Flanders. We used the sedimentary records preserved in tidal marsh sediment cores as an archive for reconstructing historical changes in surface water PO 4 . For sediment samples at sequential depths below the sediment surface, we dated the time of sediment deposition and analysed the extractable sediment P. The resulting time series of sediment P was linked to the time series of measured surface water-PO 4 concentrations (data 1967-present). By combining those datasets, the sorption characteristics of the sediment could be described using a Langmuir-type sorption model. The calibrated sorption model allowed us to estimate a pre-industrial background surface water PO 4 levels, based on deeper sediment P that stabilized at concentrations smaller than the modern. In three out of the four cores, the sediment P peaked around 1980, coinciding with the surface water PO 4 . The estimated pre-industrial (∼ 1800) background PO 4 concentration in the Scheldt River water was 62 [57; 66 (95 % CI)] µg PO 4 -P L −1 .
That concentration exceeds the previously estimated natural background values in Flanders (15-35 µg TP L −1 ) and is about half of the prevailing limit in the Scheldt River (120 µg PO 4 -P L −1 ). In the 1930s, river water concentrations were estimated at 140 [128; 148] µg PO 4 -P L −1 , already exceeding the current limit. The method developed here proved useful for reconstructing historical background PO 4 concentrations of a lowland tidal river. A similar approach can apply to other lowland tidal rivers to provide a scientific basis for local catchment-specific PO 4 backgrounds.

Introduction
Excess phosphorus (P) concentrations in surface waters is a global problem (Azevedo et al., 2015;Dodds and Smith, 2016;Elser et al., 2007). Eutrophication by excess nutrients, including phosphate (PO 4 ) and nitrogen (N), can lead to hypoxia, acidification, and harmful algal blooms (Azevedo et al., 2015;Correll, 1998;Watson et al., 2018). Therefore, limiting P concentrations in the surface water is crucial to ensure a good ecological status. Lowland river systems are at higher risk for eutrophication than upland streams (Watson et al., 2018). As a result, eutrophication of lowland rivers is on the international agenda (Jarvie et al., 2006;Mainstone and Parr, 764 F. Lauryssen et al.: Estimation of the natural background of phosphate in a lowland river thophosphate (PO 4 ), almost identical to the reactive P determined by a colour reaction. Other P forms present in surface water include organic P fractions and P adsorbed to mineral colloids. Total P (TP) refers to all P forms together. The environmental limits for P are expressed as reactive P (equated to PO 4 -P limits), TP limits, or both.
The lowland rivers of densely populated regions do not achieve good water quality mainly due to the excess of nutrients (Bitschofsky and Nausch, 2019;Huet, 1990;Van Der Molen et al., 1998;Van Puijenbroek et al., 2014;Rönspieß et al., 2020;Schulz and Herzog, 2004). For example, since 2004, the average PO 4 concentration in Flemish waterways has stabilized at 290 µg PO 4 -P L −1 . That concentration is well above the limits varying between 70-140 µg PO 4 -P L −1 for different river types (Smolders et al., 2017;VMM, 2018). Despite the current net-zero P balance in agricultural soils in that region, the situation has not improved since 2004 (VMM, 2017). The question arises when or even whether these limits can be achieved.
Since 2000, the European Union has regulated surface water quality with the Water Framework Directive 2000/60/EC (WFD), which does not prescribe limits but provides a framework for local regulations. The WFD has identified a high ecological status of a river if nutrient conditions remain within the range associated with undisturbed conditions, i.e natural background levels (European Commission, 2000). However, the definition of the natural background has been subject to debate for many river basins (Matschullat et al., 2000;van Raaphorst et al., 2000). The natural background can be defined as the situation (concentration range) found in the environment without any human activity, reflecting only natural geochemical processes (Laane, 1992;Reimann and Garrett, 2005). This definition implies that concentrations must be estimated before human activity, which is not always feasible. Therefore, a pre-industrial background can be defined instead, inferred from samples dating before the industrial revolution (Reimann and Garrett, 2005). The preindustrial background can logically be affected by anthropogenic processes. Alternatively, natural background concentrations can be estimated by sampling regions with an expected minimal anthropogenic influence, i.e. reference lakes and rivers (Cardoso et al., 2007). Natural background concentrations have been established for different chemical elements in rivers in Europe; however that did not include TP or PO 4 -P (Salminen et al., 2005).
For P in densely populated regions such as Flanders, the natural background concentrations can only be inferred indirectly. Natural background TP concentrations for Flanders were set at 15-35 µg P L −1 , inferred from reference lakes sampled in central and Baltic states in Europe (Cardoso et al., 2007). However, the natural background PO 4 -P concentrations of lowland rivers may be higher than those of upland rivers. For example, diatom assemblages revealed natural eutrophic conditions in the Spree River in Germany with TP concentrations of 80 µg L −1 , compared to recent data of 120 µg TP L −1 (Schönfelder and Steinberg, 2004;Zak et al., 2006).
Lowland rivers can be significantly loaded with P derived from the riverbed sediment, defined as internal loading (Froelich, 1988). The internal loading of P better explains temporal and spatial trends in river P concentrations than the P emissions of the lowland river system of Flanders (Smolders et al., 2017). Biogeochemical processes in the sediment explain that result. Ferric iron (Fe(III)) and aluminium oxyhydroxides have a high affinity for PO 4 anions and limit the PO 4 in solution (Borggaard et al., 1990;Holtan et al., 1988). However, anoxic conditions lead to the reductive dissolution of those Fe(III) minerals, releasing the associated P to the overlying water when the sediment is strongly reduced (Baken et al., 2015;van Dael et al., 2020). The small rivers and ditches in lowland regions have a low water velocity and are nutrient-rich, leading to anoxic conditions during summer and autumn. Those anoxic events explain the typical summer peaks in PO 4 in small rivers. Moreover, regional differences in sediment P / Fe concentration ratios explain regional differences in surface water PO 4 concentrations (Smolders et al., 2017).
Sediment analysis is linked to surface water P and can be a valuable tool for assessing historical river water quality. In surface waters, sediments can serve as a sink or a source of PO 4 , depending on the sediment surface chemistry and water concentrations (Froelich, 1988;House and Denison, 1998;Simpson et al., 2021;van der Zee et al., 2007). For example, P storage on fine bed sediments can amount to 60 % of a catchment nutrient budget (Ballantine et al., 2009;Svendsen and Kronvang, 1993). The essential processes for PO 4 are adsorption and desorption from Fe oxyhydroxides, present in the suspended matter or bed sediments (Froelich, 1988;van Raaphorst and Kloosterhuis, 1994;van der Zee et al., 2007). Sediment P concentrations can likely predict surface water P concentrations and have been relevant for the long-term reconstruction of P in the environment (Wang et al., 2009;Zhou et al., 2005). For example, Boyle et al. (2015) used P profiles from lake sediments in the UK to infer catchment P inputs over the last 10 000 years and linked that to the historical evolution in population density. Similarly, banded iron formations in deep oceanic waters allowed the inference of oceanic P concentrations of over 2 billion years ago (Bjerrum and Canfield, 2002). Likewise, the sediments deposited by rivers or lakes react with surface water PO 4 and are deposited in regularly flooded areas. Thus, those sediments can serve as an archive for reconstructing historical P emissions trends and provide useful information on historical P concentrations in adjacent water bodies (Birch et al., 2008).
In lowland rivers with tidal influence, like the Scheldt, vegetated tidal marshes develop along the river banks. Tidal marshes directly adjacent to tidal rivers are regularly flooded during high tides: these flooding events deposit sediments and associated elements like P on the densely vegetated marsh surface (Friedrichs and Perry, 2001;De Swart and Zimmerman, 2009;Temmerman et al., 2004a). The net accumulation of sediments increases the elevation of tidal marshes over time (Temmerman et al., 2003a). Therefore, researchers have used tidal marshes as sediment archives of deposited substances other than P, such as organic carbon (Van de Broek et al., 2019) and silicon (Struyf et al., 2007). However, it remains to be investigated to what extent P concentrations measured in tidal marsh sediment archives can be used to reconstruct historical changes in PO 4 concentrations in the adjacent estuary.
This study tested and evaluated a methodology to estimate the pre-industrial background water PO 4 concentrations based on the analysis of tidal marsh sediment. Those sediments had been deposited over multiple centuries on the banks of Flanders' largest tidal river, the Scheldt. Using sediment analysis and a sorption model provided the first estimate of pre-industrial PO 4 levels in a large lowland river. First, we described the tidal marsh sediment sorption characteristics by linking the P concentration of tidal marsh sediments to historical measurements of PO 4 in the Scheldt River water. Those sorption characteristics allowed an estimation of historical river water PO 4 concentrations. This estimate was based on an analysis of sediments deposited in the 1800s or before industrialization. The underlying assumption is that sediment P remains immobile and that the sediment's depth profile reflects the historical trend of PO 4 in the Scheldt River. Accordingly, we argue that the sediment P composition in deeper sediment layers of tidal marshes provides an archive of the historic PO 4 concentration of the adjacent river. A database containing measurements of the PO 4 concentration in the Scheldt River's surface water (1967current) verified this assumption. This study hypothesizes that the previously estimated natural background P of this major lowland river is larger than that estimated earlier for lakes (15-35 µg P L −1 ).

Study area
Freshwater tidal marshes were sampled at four locations along the Scheldt River ( Fig. 1 and Table S1 in the Supplement). The Scheldt estuary is located in northern Belgium and the southwestern Netherlands and flows into the North Sea. The river basin of the Scheldt covers a large part of Flanders (71 %) and the adjacent region of northern France; the total catchment area is approximately 22 000 km 2 . The population living in the river basin is about 10 million (Meire et al., 2005). The tidal wave extends from the mouth (Vlissingen) to 160 km upstream near Ghent, where sluices stop the tidal wave. The estuary's freshwater tidal zone reaches from Ghent to Rupelmonde (Fig. 1). This research focused on freshwater tidal marshes, i.e. located in this freshwater tidal zone of the estuary. Brackish waters experience the mixing of seawater, making it difficult to distinguish the anthropogenic sources from seawater influence. Furthermore, saltwater in the North Sea has PO 4 concentrations about a factor of 10 lower than fresh water in the Scheldt River (Burson et al., 2016). This research was focused on freshwater lowland river systems and the human influence on the P concentrations; saltwater environments were beyond the scope of this study.
Sediment accreting in tidal marshes originates from the deposition of riverine suspended matter, including inorganic mineral sediment and organic matter (Callaway et al., 1996). We discriminate between old and young tidal marshes, hereafter referred to as marshes. Old marshes have a higher elevation compared to young marshes. As a general mechanism, young marsh surfaces accumulate sediments quickly and increase their elevation asymptotically up to an equilibrium level, which is around the mean high water level (MHWL) (Pethick, 1981;Temmerman et al., 2003a). Temmerman et al. (2003a) defined an old marsh as visible on topographic maps of Ferraris (1774-1777), so it was formed before the 19th century (Temmerman et al., 2003a). Young marshes in the Scheldt estuary were formed more recently, by the natural establishment of pioneer marsh vegetation on formerly bare tidal mudflats, generally after 1944. During the last decades, the young marshes had a surface elevation below MHWL. As a result, young marshes experienced more frequent inundations and had larger sediment accretion rates than old marshes. For example, between 1931 and 1951, young marshes accumulated at rates of 1.6 to 3.2 cm yr −1 . In contrast, the elevation of old marshes was very close to the yearly MHWL increase rate of 0.3 to 0.6 cm yr −1 in the western Scheldt (Temmerman et al., 2003a).
This study analysed depth profiles of sediment cores originating from tidal marshes along the freshwater Scheldt River. The analysis contained two old and two young marshes, Old1, Old2, Young1, and Young2 (locations indicated in 766 F. Lauryssen et al.: Estimation of the natural background of phosphate in a lowland river Fig. 1). The coordinates of sampling locations can be found in Van de Broek et al. (2018 and Sect. SI.I in the Supplement. Marshes Old1 and Young1 originated from the tidal marsh named the Notelaer, Old 2 from Grembergen, and Young2 from Mariekerke. Eight cores were analysed: three replicate cores for both sites Old1 and Young 1 and one core for Old2 and Young2.

PO 4 concentration in surface waters
The Flanders Marine Institute (IMIS) provided surface water phosphate (PO 4 ) data measured colourimetrically on a filtered water sample and total phosphorus (TP) by acid digestion and a segmented flow analyser. Concentration data of PO 4 in Scheldt River were available from 1967 to 2019, originating from different sources and compiled by the research programme Onderzoek Milieu Effecten Sigmaplan (OMES). The OMES programme did additional quality controls on the data (ECOBE -UA; The Flemish Waterway, 2019). The different sources are described in Sect. SI.V in the Supplement (De Pauw, 2007;ECOBE -UAntwerpen, 2007; Institute for Hygiene en Epidemiology (IHE), 2007; ECOBE- UA and De Vlaamse Waterweg, 2016;Van Meel, 1958).
The open-source software R (R Core Team, 2020) was used to compile all available datasets for PO 4 closest to the study sites (Temse) and to calculate annual means by averaging all observations within a year. The annual means of PO 4 were used to visualize the evolution of PO 4 in the Scheldt River (Fig. 2). The emissions of P mainly originate from point sources due to domestic loading (Billen et al., 2005). As a result, the increasing surface water P concentration between 1950 and 1975 can be related to the rise in the number of households connected to sewer systems. At first, no wastewater treatment was in place, resulting in a sharp increase in nutrient loads to the river. However, since 1985, wastewater treatment has significantly improved the situation (Billen et al., 2005).

Sediment sampling
The sediment samples used here had been collected during a previous study about carbon sequestration in tidal marsh sediments in the Scheldt estuary ( Van de Broek et al., 2018. Collection of undisturbed sediment profiles on the tidal marshes took place between July and September 2016 (Old1, Young1, Old2, Young2; Fig. 1). Undisturbed sediment cores were taken at each sampling location using a gauge auger (0.06 m diameter). The cores were divided into subsamples with a 0.03 m interval. The sediment samples were dried at a maximum temperature of 50 • C for 48 h, crushed, and sieved to a < 2 mm grain size. Macroscopic vegetation residues were removed manually using tweezers (Van de Broek et al., 2018). Bulk density, grain size distribution, and organic carbon (OC) content were analysed by Van de Broek  2018). We refer to Van de Broek et al. (2016 for further sample collection and processing information.

Sediment analysis
The dried sediment samples were analysed for oxalateextractable P, Fe, Al, and Mn (P ox , Fe ox , Al ox , Mn ox ; Schwertmann, 1964). The preparation of extraction solution and dilutions were made with ultrapure water (Milli-Q ® ), and all glassware was acid-soaked overnight in a 1 % HCl acid bath to prevent P contamination. That acid oxalate extractant, a mixture of ammonium oxalate (0.2 M) and oxalic acid (0.2 M) at pH = 3, targets poorly crystalline oxyhydroxides of Fe, Al, and Mn and the associated P (Schwertmann, 1964). Those poorly crystalline oxyhydroxides are the most reactive due to their large specific surface area (Hiemstra et al., 2010). The extraction was done with 1 g of dry sediment in 50 mL extraction solution over 2 h in an end-overend shaker at 20 • C (26 rpm). The suspension was filtered through a 0.45 µm membrane filter (CHROMAFIL ® Xtra PET -45/25). Analytical blanks, internal reference samples, and duplicate samples were included in every batch to ensure the analysis's quality, purity, and reproducibility. The extracts were diluted 20 times and measured by inductively coupled plasma optical emission spectrometry (ICP-OES). The degree of P saturation (DPS; %) was calculated as in Eq. (1). The DPS represents the extractable (P ox ) ratio to the P sorption capacity of the sediment. This P sorption capacity is estimated as half of the sum of oxalate Fe ox and Al ox , because not all the Fe and Al in the soil are available for P sorption with Fe ox , Al ox , and P ox in molar units.
The DPS is expressed as a percentage and can be interpreted as the ratio of sorption sites on the sediment occupied by P. Previous research used the DPS to identify agricultural areas sensitive to phosphate leaching and showed a good correlation with pore water P concentrations (Breeuwsma et al., 1995;Lookman et al., 1995;Schoumans and Chardon, 2015;Schoumans and Groenendijk, 2000;van der Zee, 1988). The DPS relation for porewater-soil systems was developed and verified by van der Zee et al. (1990). Lexmond et al. (1982) illustrated that the maximal sorbed P was about half the pool available after a long-term precipitation experiment. Therefore, the factor 0.5 is an empirical value representing the soil's sorption capacity. The parameter α primarily affects the maximum sorption capacity. So they set α at 0.5 ± 0.1. However, even among soils, this parameter varied between 0.3 and 0.6 (Lexmond et al., 1982). For this research, low background concentrations are most important, so maximal sorption occurring at high PO 4 concentrations is less relevant.

Age-depth model
The sediment analysis and the surface water PO 4 data had to be linked by a corresponding date and location to fit a sorption model. Therefore, an age-depth model was used to calculate the time since deposition of each sediment sample. Temmerman et al. (2004b, a) developed a time-stepping marsh sedimentation model (MARSED). That model estimates sediment deposition rates and the resulting evolution of the tidal marsh elevation in the Scheldt estuary. Hence, we could use MARSED to determine the time since deposition of sediments throughout the sampled sediment profiles. The MARSED model simulates the tidal supply of suspended sediments and the settling to the marsh surface during tidal inundation cycles integrated over the years. The model was calibrated and validated against measured sediment deposition rates on the Scheldt estuary tidal marshes from 1945 until 2002 (Temmerman et al., 2003a(Temmerman et al., , 2004b. The empirical data on sediment deposition rates were derived from radiometric and palaeoenvironmental dating of sediment cores at the exact locations sampled for the present study (Temmerman et al., 2004a, b). For our current study, we extrapolated the MARSED model simulations of sediment accretion from 2002 until 2016, the sampling date of the sediment cores ( Van de Broek et al., 2018). However, simulations overestimated the observed marsh surface elevation in 2016 by 25 cm for sampling location Old1, 29 cm for Young1, 19 cm for Old2, and 8 cm for Young2 (observed by RTK GPS surveying; Van de Broek et al., 2018;Poppelmonde, 2017). The MARSED model was initially designed to simulate the overall sediment accretion and surface elevation changes in tidal marshes in response to sea-level rise scenarios, for which those errors were acceptable. In contrast, the most important was the time of sediment deposition throughout the sediment profile for the present study. Therefore, the original age-depth relation calculated by MARSED was recalibrated using observed age-depth points. The observed age-depth points originated from GPS measurements of marsh elevation in 2016 (M. Van de Broek, unpublished data) and previously published radiometric and palaeoenvironmental dating (Temmerman et al., 2003(Temmerman et al., , 2004. This rescaling procedure is explained in the Supplement (Figs. S1-S4 in the Supplement).
An approximate extrapolation procedure was used to estimate the sediment deposition time from depths below the oldest measured age-depth points (mentioned in the previous sentence). The observed age-depth points were available from 1958 for sampling site Old1, 1947for Young1, 1963 for Old2, and 1968 for Young2 (Temmerman et al., 2004). This extrapolation procedure could only be applied for old marshes, which were defined as marshes that existed at least since the end of the 18th century (Temmerman et al., 2003a(Temmerman et al., , 2004. Two sediment cores originated from old tidal marshes (Old1 and Old2). Based on observed age-depth points, it has been proven that the old marshes reached equilibrium with the MHWL before 1944. After 1944, the old marshes have built up their elevation at a rate comparable to local MHWL rise (Temmerman et al., 2003a). Here, we assumed that also between 1800 and 1944 these old marshes accreted at a rate comparable to the MHWL rise.
Historical tide gauge data of MHWL rise were available from 1901 for site Old1 and 1930 for site Old2 (Schelde-Monitor Team and VNSC, 2020; Temmerman et al., 2003a), and linear regression of the MHWL against time was used to estimate the marsh surface elevation before 1944 (Figs. S6 and S7 in the Supplement). However, the dating accuracy will be lower going further back in time. Furthermore, such extrapolation to earlier dates is not appropriate for young marshes, as they were only formed after 1950 by pioneer vegetation establishment on formerly bare mudflats (Temmerman et al., 2003a(Temmerman et al., , 2004. Those mudflat sediment profiles do not have continuous sedimentary records as tidal marshes and are likely to be disturbed by erosion and sedimentation alternations (Belliard et al., 2019). Therefore, the sediment deposition time could not be extrapolated for the young marsh sampling locations.
2.6 Relating surface water PO 4 with sediment P: the sediment-water model The age-depth model and linear regression of MHWL provided a deposition year for each sediment sample. Thereby, the dataset of water PO 4 between 1967-2016 was linked to the sediment DPS for each core. The resulting dataset contained all available surface water PO 4 readings between 1967 and 2016, closest to the tidal marshes in Temse (n = 1932) and a corresponding DPS value. The DPS value of a sediment sample originates from a specific layer of one sediment core or a mean DPS of the replicate sediment samples. This dataset allowed us to fit a sorption model further termed the 768 F. Lauryssen et al.: Estimation of the natural background of phosphate in a lowland river sediment-water model. Schoumans and Groenendijk (2000) presented a Langmuir-type sorption model to predict PO 4 concentration leaching from a soil layer based on the DPS Eq. (2).
with [PO 4 ] phosphate concentration in kg L −1 , K the sorption constant (L kg −1 ), and DPS (degree of P saturation; %). This model adequately described P sorption in soil across a wide range of pH values, including the Scheldt River pH (Schoumans and Groenendijk, 2000;. The model relies on surface complexation between PO 4 and Fe and Al oxyhydroxides in the sediment. That complexation is determined by a chemical equilibrium between solid (adsorbed) and dissolved PO 4 phase (Warrinnier et al., 2019). The parameter K of existing soil models has been calibrated for the soil-pore water system, and the sediment-water parameter (K) is unlikely equal. Therefore, the model was calibrated by fitting parameter K (Eq. 2) on sediment DPS measurements and recent Scheldt water PO 4 measurements. As a result, the fitted K value is adapted to the local geochemistry of tidal marsh sediments and the surface water. We explored 16 different scenarios to fit the sedimentwater model Eq. (2). These scenarios illustrate the statistical uncertainty surrounding the estimated PO 4 concentrations. The model was fitted separately for each site sediment core or on the combined replicate cores for Old1 and Young1 (Sect. SI.VI in the Supplement). Every sediment sample had between one and three replicates, depending on the depth and the site. The average value of these replicates was used or the individual replicate's DPS values. One sediment sample covered several deposition years, so multiple PO 4 observations corresponded with each sediment sample. Again, the average of all corresponding PO 4 readings was taken, or all available values were used separately. The combination of mean or individual DPS and PO 4 resulted in 16 models (Table S2 in the Supplement). For each of these, the parameter K was fitted by non-linear least-squares regression with R using the Rstudio interface (R Core Team, 2020; RStudio Team, 2015).

Evaluation of model performance
The predictions of the sediment-water model were evaluated based on several parameters: the residual standard error (RSE), the Nash-Sutcliffe model efficiency (E), and the measured surface water PO 4 plotted against predicted PO 4 between 2007 and 2016 (Table S2 and Fig. S10 in the Supplement). Additionally, the percentage bias (PBIAS) was calculated for data points between 2007 and 2016. The PBIAS measures the average tendency of the simulated data to be larger or smaller than their observed counterparts. That difference is expressed as a percentage of deviation from the observations (Moriasi et al., 1983;Eq. 3). The predictions of recent years are interesting to evaluate the model's per-formance for two reasons. First, the most recent surface water PO 4 concentrations are relatively low and more representative of background concentrations. Second, the monitoring data have a high temporal resolution, and the age-depth model is more accurate at shallow depths.

History of surface water PO 4 concentrations
The Scheldt PO 4 concentrations varied greatly over the past decades, with the peak in surface water PO 4 concentrations between 1975 and 1985 (Fig. 2). In Temse, the annual mean concentrations rose from 410 µg PO 4 -P L −1 in 1967 and peaked in 1980 with 1570 µg PO 4 -P L −1 . Between 1990 and 2003 concentrations decreased and stabilized between 160 and 200 µg PO 4 -P L −1 in Temse. The current PO 4 levels are a factor of 2 lower than in 1967 and almost a factor of 10 lower than the peak in 1980 ( Fig. 2 and Table 1).

Sediment cores
The P ox in the sediments ranged between 370 and 13 000 mg P kg −1 , while the DPS ranged between 13 % and 94 % (Table 1). In all soil cores starting at the surface, the DPS and P ox increased with depth and peaked at about 0.5 m below the surface (Figs. S7 and S8 in the Supplement). In deeper (> 1.0 m) sediment layers, P ox and DPS decreased and stabilized for Old1, Young1, and Young2 (Table 1). Overall, the P ox increased by an average factor of 3.5 between the surface and the maximum concentrations ( Fig. S8 and Table 1). The sediments with these peak DPS values were deposited between 1960 and 1985 in three of the four sediment cores (Fig. 2). Only the core Old2 peaked earlier (ca. 1940-1950). Most importantly for this work, DPS for Old1 showed an apparent stabilization in deeper or older layers, which indicated undisturbed sediment layers (Figs. 3 and S7).
Within the first metre, Fe ox was stable in the three soil cores (Old1, Young1, Old2) with concentrations around 20 000 mg kg −1 , except for Young2 for which Fe ox was a factor of 2 larger (Fig. S9 in the Supplement). For Young1 and Young2, Fe ox concentration decreased at depths > 1 m. For Old1, Fe ox showed a steady decline from 20 000 mg kg −1 at the surface to 10 000 mg kg −1 at the bottom of the profile (Fig. S9). The Al ox concentrations showed a similar trend as the P ox concentrations, with an initial increase followed by a decrease with depth. The strong correlations of Al ox and Fe ox with P ox (r Al = 0.73 and r Fe = 0.65) illustrate the positive effect of Fe and Al oxyhydroxides on P sorption.  Figure 3. The degree of phosphate saturation (DPS) timeline based on four tidal marsh sediment sites. Each dot represents a sediment analysis. The year assigned to each sediment analysis was calculated with the age-depth model. Before 1930, no model dates were available. Therefore a linear regression of the MHWL was used to extrapolate the dates for the old marshes. Dates before 1930 are increasingly uncertain going further back in time. For young marshes, such extrapolation was not possible. The points before the formation of the marshes are indicated with a question mark.

Sediment core selection
Under the assumption that PO 4 does not migrate, the tidal marsh sediment cores can provide an archive for river water PO 4 . Considering P migration, it is crucial to evaluate the distance from a creek within the tidal marsh. That distance is essential for two reasons. First, within 10 to 20 m of the creeks, the groundwater table fluctuates largely with the tides, which can induce vertical P migration (Van Putte et al., 2020). Secondly, sediment accretion is more difficult to predict at closer distances to the creeks and can affect the age-depth relation (Temmerman et al., 2003b). The distance from the sediment cores to the nearest creek was 21 m for Old1, 56 m for Young1, 35 m for Young2, and 5 m for Old2.
The assumption that PO 4 does not migrate may be most violated at Old2 and Young2. The profile of Old2 indicated P migration because it had a peak of P ox at an earlier date (1950) than was expected from surface water data (1980) (Fig. 2). Consequently, Old2 was not taken up to interpret the relation between DPS and PO 4 . For core Young2, deeper sediment layers had a larger DPS than the surface layers (Table 2). Additionally, the age estimation of sediments older than 1968 was impossible due to this tidal marsh's young age. Furthermore, Fe ox concentrations were a factor of 2 larger than the other cores (Fig. S9) and a factor of 2 larger than the average sediment Fe concentration of the upper Scheldt basin (VMM, 2019). The local enrichment in iron lowers the DPS values and makes the core less representative of the average situation in the Scheldt. These observations made Young2 inappropriate to fit the relation between DPS and PO 4 .
The two remaining soil cores, Old1 and Young1, originated from the same tidal marsh area named "the Notelaer", located near the city of Temse (Fig. 1). That marsh has been the topic of multiple studies on sediment accretion (Temmer -Table 2. The predicted concentrations of phosphate (PO 4 -P µg L −1 ) in the Scheldt River based on the degree of phosphate saturation (DPS) in the sediment layers of marsh Old1, dating back to 1800 (pre-industrial), where DPS values stabilized with depth at 20 %. The predicted concentration dated to 1930, where DPS stabilized at 36 %. The PBIAS is the mean difference of simulated and observed data between 2007 and 2016, expressed as a percentage of the observed data. The conversion of DPS to river phosphate concentration is based on the association of DPS with PO 4 -P. That association was calibrated to data from 1967-2016, thereby using different sediment-water models; the details of models are in Table S2. Model 3b (in italics) is proposed as the most accurate (see text).  (Fig. 2). However, the core Old1 was taken up for the model fitting because it dates back to 1800 at the deepest levels and is essential to predict the background. Furthermore, the DPS concentrations stabilized before 1920, indicating that P has not migrated to these depths, making it suitable for background prediction. These observations suggested a wellpreserved P ox and DPS profile, essential for the DPS-PO 4 relation. Therefore, Old1 and Young1 are considered the best profiles for applying the sediment-water model and interpretation of background concentrations.

Sediment-water model fit
The sediment-water model Eq.
(2) was fitted on DPS-PO 4 data from the different sediment cores (Table S2). Two observations were omitted because the DPS values were too large (0.99-1.02) and produced artefacts in the results. The Nash-Sutcliffe model efficiency (E) ranged between 0.04 and 0.85 depending on the input data (Table S2; Nash and Sutcliffe, 1970). The sediment-water model was fitted on each core's data separately and combined data from Old1 and Young1, as they came from the same tidal marsh location. The models fitted on data from sites Old2 and Young2 were not considered as migration likely affected those cores (see Sect. 3.3). The models fitted on an average DPS (across replicates) associated with individual PO 4 readings were considered most suitable (models 1b, 2b, 3b; Table S2). A single sediment sample analysis represents an average P signal over the sediment's deposition period. However, the age-depth relation can vary slightly due to the marsh surface elevation variation. By taking an average DPS from replicate cores, the variation in the independent variable was reduced. Furthermore, the prediction error increased in most models by relating individual rather than mean DPS values with individual PO 4 measurements (Table S2). Models using unique DPS associated with single PO 4 data duplicated or even triplicated the PO 4 data, artificially creating more degrees of freedom (model 1c, 2c, 3c; Table S2). Using mean PO 4 values artificially reduced the degrees of freedom, compromising the model predictions by increasing RSE and widening confidence intervals (models 1a, 2a, 3a; Table S2). The fitted parameter K (L kg −1 ) ranged between 1.0 × 10 6 and 5.4 × 10 6 for the different input datasets, with the 95 % confidence intervals ranging between 0.8×10 6 and 7.2×10 6 . The variation in parameter K for the various input datasets was larger than the individual confidence limit variation (Table S2). Thus, the uncertainty was more pronounced due to the variability in sediment samples than due to the model fit.

Model performance
The sediment-water model performance was evaluated with several parameters (RSE, PBIAS, and E) and the actual PO 4 by predicted plots of the PO 4 concentrations over the last decade. Those recent PO 4 concentrations are more comparable to the background (Fig. S10 in the Supplement and Table 2; mean Temse [2007][2008][2009][2010][2011][2012][2013][2014][2015][2016] = 170 µg PO 4 -P L −1 ). Model 3b was considered the most suitable for predicting background concentrations. The PBIAS was the lowest for recent observations for model 3b. The average tendency of simulated data compared to the observations was only −14.9 %, which is within the acceptable range of ±25 % (Moriasi et al., 1983). Model 2b had an underestimation of observed data of more than 28 %, and model 1b overestimated recent observations by more than 60 %. Such an overestimation is unwanted for calculating the background, and therefore, both were considered unsuitable ( Table 2). The actual PO 4 by predicted plots illustrates a similar message (Fig. S10). Based on these observations, model 3b was con- sidered the best model, although the residual standard error (RSE) was lower for model 2b (Table S2). The selected model 3b successfully reconstructs the rise and fall in surface water PO 4 concentrations based on the sediment characteristic DPS (Fig. 4).
Maxima of monitored and predicted PO 4 concentrations coincide in time and have a similar size (Fig. 4). For example, in 1973, the average PO 4 concentration predicted by the model was 1200 µg PO 4 -P µg L −1 , and measured concentrations were on average 1300 µg PO 4 -P µg L −1 . The maximal predicted PO 4 concentration was 2200 µg PO 4 -P L −1 , while the maximal observed concentration was 3000 µg PO 4 -P L −1 . Predictions for recent years are within 15 % of the observed data (e.g. 2015: model: 133 µg PO 4 -P L −1 , measured 155 µg PO 4 -P L −1 ). Between 1940 and 1990, the modelled PO 4 concentrations had more variation. Likewise, monitored PO 4 data are spread more between 1967 and 1990 (Fig. 2). Before 1930, modelled PO 4 concentrations stabilized at levels below current observations (Fig. 4).

Estimating background PO 4 concentrations in the Scheldt River
The deepest sediment layers are most suitable for predicting background PO 4 concentrations of the Scheldt River water. These layers are the oldest and expected to have experienced the lowest impact of P additions from anthropogenic sources. The Old1 marsh site was appropriate for this purpose as it developed before 1774, before the industrial revolution in Belgium. The average DPS for the bottom sediments, dated between 1800 and 1840, was 20 % for core Old1 (Table 1

concentrations
This work presents a novel approach to reconstruct background surface water PO 4 concentration in a tidal river using the DPS of adjacent tidal marsh sediments. The background concentration is essential in the context of developing local nutrient limits. The predicted pre-industrial background concentration (62 µg PO 4 -P L −1 ; Table 2) is about half of the current surface limit of the Scheldt (120 µg PO 4 -P L −1 ; Flemish Government, 1995). Remarkably, the predicted background concentrations are about a factor of 2 larger than the background estimates of lake waters for Flanders today (15-35 µg PO 4 -P L −1 ; Cardoso et al., 2007). That pre-industrial PO 4 concentration is about 3 times lower than the current concentration in the Scheldt. For example, between 2007 and 2016, the mean PO 4 concentration of the Scheldt in Temse was 170 µg PO 4 -P L −1 . However, in the 1930s, the concentration was estimated at 140 µg PO 4 -P L −1 and larger than current limits, at a time before widespread connection to sewer systems, P-loaded detergents, and application of mineral fertilizers.
Those results suggest that the sediment internal loading triggered by summer anoxia in lowland rivers contributes to larger PO 4 concentrations than estimated before (see introduction). The summer PO 4 peak lasts about 5 months per year in Flanders and largely affects the rivers' mean P concentrations (Smolders et al., 2017). Summer anoxia can occur in eutrophic lakes or sometimes in oligotrophic brown water lakes (Nürnberg, 1995). Additionally, lowland rivers in Flanders are primarily groundwater-fed, and 73 % of streamflow can be attributed to base flow. The groundwater in Belgium has a median P concentration between 150-320 µg P L −1 (Edmunds and Shand, 2009). Therefore, groundwater feeding the river waters logically affects the river P concentrations. In contrast, primarily rain-fed lakes will have lower P concentrations, with rain P ranging between 1.5 and 120 µg P L −1 (Migon and Sandroni, 1999).

Limitations of the model
Care needs to be taken with background extrapolations to ensure that post-depositional processes have not modified the biogeochemical patterns and that the area represents the area of interest (Reimann and Garrett, 2005). Several factors can obscure the reconstructed background concentrations. First, vertical migration of P can enrich deeper sediment layers, causing an overestimation of the background. Second, the sediment profiles at the tidal marshes are almost permanently saturated, so the intrusion of P-rich groundwater could affect the P concentrations in the tidal marsh sediment. Moreover, depending on the tidal marsh elevation, periodic flooding occurs at an approximate range of 300-350 inundations per year (Temmerman et al., 2003b). These conditions could favour P migration due to the reductive dissolution of Fe (oxy)hydroxides (Baken et al., 2015;van Dael et al., 2020).
Two cores with indications of PO 4 migration were removed from the analysis to address the issue (Old2 and Young2). These cores were identified by the DPS age profile and considering the distance from the nearby creeks (Figs. 2 and 1). Additionally, the DPS levels of the deepest sediment layers were compared with layers at the surface. The surface layers had lower DPS levels than the deepest layers for one core (Young 2). The two remaining cores (Old1, Young1) had lower DPS levels in deeper sediment layers (Fig. S7). More importantly, the modelled peak in PO 4 concentrations based on the cores Old1 and Young1 were found within 2 years of the monitored peak and had a similar magnitude (Fig. 4). The coinciding peaks illustrate little migration of PO 4 in Old1 and Young1, thereby justifying these cores as an archive for water PO 4 .
The limited migration is also logical: at the average DPS of 90 % in sediment showing at the peak, the sorption models predict that the solid-liquid P concentration ratio is 2900 L kg −1 with the average K value of models of Table 2. That solid-liquid ratio can be converted to a dimensionless retardation factor representing the ratio of the distance migrated by the PO 4 compared to the distance travelled by percolating water. For example, the retardation was calculated to be 7500 with a bulk density (ρb) of 1.3, porosity (θ ) of 0.5, and a net vertical annual water percolation of about 2 m. That retardation corresponds to a net vertical P migration rate of 2.5 cm over 100 years, i.e. vanishingly small (calculation details not shown).
Secondly, there is uncertainty on the age-depth estimation of the sampled sediment profiles. The age-depth model is expected to be most reliable for the Young1 sediment core, as it is based on a fitting of a modelled age-depth relation to four observed age-depth points, while we only had two observed age-depth points available for the other cores (Temmerman et al., 2004a). Additionally, observed age-depth points were not older than 1944. Hence, the extrapolation of the agedepth model to periods before the older available age-depth points is increasingly uncertain.

Pre-industrial and natural background values
The population increase between 1800 and 1930 can provide a first very crude estimate of the population-DPS relation in the Scheldt basin. In 1800 the population in Belgium was around 3 million. Later, in 1930, this number had more than doubled to 7 million (Vanhaute, 2003). A linear relation between both suggests that the DPS is 8 % for the pristine pre-anthropogenic environment, corresponding with a PO 4 concentration of 19-41 µg PO 4 -P L −1 , i.e. close to what researchers have indicated for pristine lakes. Such predictions need to be corroborated with older sediment observations and other archaeological information. The Scheldt River is logically more aerated than smaller lowland rivers where summer anoxia is naturally more present; i.e. the pristine PO 4 -P values will be higher.

Conclusions
Our study illustrated that tidal marsh sediments could evaluate pre-industrial background PO 4 concentrations of the freshwater rivers like the Scheldt River. A sediment assessment can record time-integrated environmental events, providing useful spatial and temporal information. Our data estimated the pre-industrial background concentration at 62 µg PO 4 -P L −1 [95 % CI (57; 66)], about half of the environmental limits set for surface waters in Flanders and neighbouring countries. Around 1930, the PO 4 levels were only about 20 % lower than today, which is a remarkably large concentration at a time before the massive application of mineral fertilizers, with lower population density and limited connection to sewer systems. The current PO 4 concentrations decreased 10 times from the peak found 40 years ago, reflecting wastewater treatment efforts and reducing diffuse P emission. It is also clear from this study that the pristine, preanthropogenic PO 4 -P concentrations in the Scheldt River are well below the current ambient ones.
Data availability. The Supplement provides the sediment data analysis and age-depth model results in csv format. In addition, results of surface water data are available upon request at the IMIS (Flanders Marine Institute).
Author contributions. FL, ES, and PC designed the research. FL conducted the investigation process and developed the methodology under supervision of ES. MVDB carried out the fieldwork and conceptualized the use of the samples. ST provided the methodology for the age-depth model and software. TM validated the use of the surface water data. EVM and FL placed the results in perspective with historical data. All the authors contributed to discussion and data interpretations, review, and editing of the work.
Competing interests. The contact author has declared that neither they nor their co-authors have any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. The results of this research greatly depended on the data collected by the OMES-monitoring and The Flemish Waterway. Many years of intensive data collection and quality assessment of the Scheldt River resulted in a unique and valuable phosphate time series. We have the utmost respect for their work and are thankful we could apply the dataset for this research. We acknowledge Dries Grauwels and Kristin Coorevits for technical assistance. We recognize the efforts from the anonymous reviewers for their constructive comments on the work, which improved the quality of the result. Finally, thanks to the Scheldt for providing this beautiful sediment archive to travel back in time and explore environmental history.
Financial support. This research has been supported by the Research Foundation -Flanders (FWO), project identifier: G089319N.
Review statement. This paper was edited by Aninda Mazumdar and reviewed by two anonymous referees.