River effects on sea-level rise in the Río de la Plata estuary during the past century

Identifying the causes for historical sea-level changes in coastal tide-gauge records is important for constraining oceanographic, geologic, and climatic processes. The Río de la Plata estuary in South America features the longest tide-gauge records in the South Atlantic. Despite the relevance of these data for large-scale circulation and climate studies, the mechanisms underlying relative sea-level changes in this region during the past century have not been firmly established. I study annual data from tide gauges in the Río de la Plata and stream gauges along the Río Paraná and Río Uruguay to establish relationships between river streamflow and sea level over 1931–2014. Regression analysis suggests that streamflow explains 59 %±17 % of the total sea-level variance at Buenos Aires, Argentina, and 28 %±21 % at Montevideo, Uruguay (95 % confidence intervals). A long-term streamflow increase effected sea-level trends of 0.71±0.35 mm yr⁻¹ at Buenos Aires and 0.48±0.38 mm yr⁻¹ at Montevideo. More generally, sea level at Buenos Aires and Montevideo respectively rises by $(\mathrm{7.3}\pm \mathrm{1.8})\times {\mathrm{10}}^{-\mathrm{6}}$ m and $(\mathrm{4.7}\pm \mathrm{2.6})\times {\mathrm{10}}^{-\mathrm{6}}$ m per 1 m³ s⁻¹ streamflow increase. These observational results are consistent with simple theories for the coastal sea-level response to streamflow forcing, suggesting a causal relationship between streamflow and sea level mediated by ocean dynamics. Findings advance understanding of local, regional, and global sea-level changes; clarify sea-level physics; inform future projections of coastal sea level and the interpretation of satellite data and proxy reconstructions; and highlight future research directions. Specifically, local and regional river effects should be accounted for in basin-scale and global mean sea-level budgets as well as reconstructions based on sparse tide-gauge records.

1 Introduction

Tide-gauge records of relative sea level go back more than a century in some places, representing some of the longest instrumental time series of the Earth system (Hogarth, 2014; Talke et al., 2018; Woodworth et al., 2010). On long climate timescales, changes in global mean sea level are informative of global ocean warming, land ice wastage, and terrestrial water storage, whereas local and regional deviations from the global average shed light on processes including ocean dynamics and gravitation, rotation, and solid Earth deformation (Gregory et al., 2019; Horton et al., 2018; Kopp et al., 2015). Identifying the mechanisms responsible for sea-level changes observed in tide-gauge records is therefore a major goal in geophysics, oceanography, and climate science (Douglas et al., 2001; Emery and Aubrey, 1991; Lisitzin, 1974).

The nature and causes of twentieth-century sea-level changes in the South Atlantic Ocean are poorly understood compared to behavior in other ocean basins during the same time period (Dangendorf et al., 2017; Frederikse et al., 2018). This knowledge gap reflects a lack of data – the basin has few long tide-gauge records (Hamlington and Thompson, 2015; Natarov et al., 2017). Given the basin's large area (Thompson and Merrifield, 2014), the absence of long data records in the South Atlantic Ocean poses a particular challenge for estimates of global mean sea-level rise (Church and White, 2011; Dangendorf et al., 2017; Frederikse et al., 2020; Hay et al., 2015; Jevrejeva et al., 2014; Ray and Douglas, 2011), but also for our understanding of circulation and climate during the past century more generally.

A recent study brings together available tide-gauge records along with other data, proxies, and models to quantify rates and mechanisms of twentieth-century South Atlantic sea-level change (Frederikse et al., 2021). Those authors determine that sea level in the South Atlantic rose about 0.3 mm yr⁻¹ faster than the rate of global mean sea-level rise owing to a combination of ocean dynamics and gravitational, rotational, and deformational effects from contemporary mass redistribution. Importantly, their estimate of twentieth-century sea-level rise over the South Atlantic rests heavily on a handful of long tide-gauge records in and around the Río de la Plata, which feature large sea-level trends that have been reported on previously (Aubrey et al., 1988; Brandani et al., 1985; D'Onofrio et al., 2008; Dennis et al., 1995; Douglas, 1997, 2001, 2008; Emery and Aubrey, 1991; Fiore et al., 2009; Isla, 2008; Lanfredi et al., 1988, 1998; Melini et al., 2004; Pousa et al., 2007; Verocai et al., 2016).

The Río de la Plata is a long, broad, shallow salt-wedge estuary that widens from ∼50 to ∼250 km and deepens from ∼5 to ∼20 m between Buenos Aires, Argentina, and Punta del Este, Uruguay, before emptying out onto the shelf (Guerrero et al., 1997; Verocai et al., 2016; Figs. 1, 2). The estuary is typified by a strong salinity and turbidity front at Barra del Indio Shoal between Punta Piedras, Argentina, and Montevideo, Uruguay, with fresher, more turbid waters upstream to the northwest, and saltier, less turbid waters downstream to the southeast (Acha et al., 2018; Guerrero et al., 1997; Moreira and Simionato, 2019). These features, and the region's hydrography and ecology generally, are strongly shaped by the situation of the estuary at the confluence of the Río Paraná and Río Uruguay, which are two of the world's largest rivers by streamflow and drainage.

Figure 1 Study region. Color shading is bathymetry (m) from the GEBCO 2021 grid (GEBCO Compilation Group, 2021). Note the nonlinear color scale. Red symbols locate tide gauges. The green star is the river mouth, selected as the confluence of the Río Paraná and Río Uruguay near Isla Oyarvide. Black dots identify other locations referenced in the text. The inset shows the study area in a global context.

Figure 2 Black curves illustrate the (a) average depth and (b) width of the Río de la Plata as a function of distance along the estuary from the river mouth based on the GEBCO 2021 grid (GEBCO Compilation Group, 2021). Values are determined by identifying all marine grid cells (depths <0) in successive 5 km increments from the river mouth. The average depth is computed as the arithmetic mean of all grid cell depths, and the width is defined as the maximum distance between the marine grid cells within the given 5 km increment. Dark blue curves and light blue shading represent best estimates and 95 % confidence intervals, respectively, of exponentials fit to the black curves using ordinary least squares. To account for residual autocorrelation, the uncertainties are based on the effective degrees of freedom assuming residuals are described by an order-1 autoregressive model.

Streamflow into the Río de la Plata increased over the past century (Dai, 2016; Dai and Trenberth, 2002; Dai et al., 2009; see Fig. 3). However, the possible influence of the increased streamflow on multidecadal and centennial sea-level trends has not been considered. Discussions of the connection between streamflow and regional sea level are mostly qualitative and center on interannual variability at Buenos Aires in relation to El Niño; for example, precipitation over the Plata Basin, streamflow of the Río Paraná and Río Uruguay, and sea level at Buenos Aires tend to increase in succession during El Niño events (Douglas, 2001; Frederikse et al., 2021; Isla, 2008; Meccia et al., 2009; Papadopoulous and Tsimplis, 2006; Raicich, 2008; Santamaria-Aguilar et al., 2017; Thompson et al., 2016; Verocai et al., 2016). Douglas (2001) and Thompson et al. (2016) argue that sea-level trends calculated from the Buenos Aires tide gauge are effected by sea-level variability during the 1982–1983 El Niño. While both studies relate this variability to river effects, Douglas (2001) favors an interpretation in terms of ocean dynamics, whereas Thompson et al. (2016) appeal to gravitational, rotational, and deformational effects. Aubrey et al. (1988) and Melini et al. (2004) give alternative interpretations of regional tide-gauge trends generally in terms of continental crustal rifting as well as subsidence and the sea-level response to the 1960 Valdivia earthquake, respectively. Therefore, it remains unclear what processes mediate the relationship between streamflow and sea level, how these two variables are related more broadly as a function of time, and whether such considerations are relevant for interpreting long-term sea-level trends. A dedicated comparison of long stream- and tide-gauge records that provides a physical interpretation and establishes causality is needed.

Figure 3 Yearly river-gauge streamflow records (Table 1). The thick blue Río de la Plata time series is the sum of the thick red Río Paraná time series from Timbúes and the thick yellow composite Río Uruguay time series. Thin time series show data from individual gauges.

Did streamflow effect long-term sea-level trends at tide gauges in the Río de la Plata? If so, what processes were involved? To answer these questions, I apply statistical analyses to annual data from stream gauges and tide gauges over the past century, and I formulate simple theories based on ocean dynamics to interpret the results. I conclude that local estuarine and coastal ocean dynamics forced by changes in streamflow had an important impact on twentieth-century sea-level rise in the Río de la Plata. Once adjusted for these effects and background late Holocene rates, both of which contribute negligibly to changes in global ocean water volume, the tide gauges show trends more in line with contemporary estimates of twentieth-century global mean sea-level rise (Dangendorf et al., 2017; Frederikse et al., 2020; Hay et al., 2015). The remainder of this paper is structured as follows: in Sect. 2, I describe the datasets; I report on results of the observational analysis, which involves correlation and regression methods applied to the data, in Sect. 3; in Sect. 4, I develop simple analytical models of the sea-level response to streamflow forcing to interpret observational results from Sect. 3; finally, I conclude with a summary and discussion in Sect. 5.

2 Data

2.1 Streamflow

I use yearly streamflow records from the Global Streamflow Indices and Metadata Archive (GSIM; Do et al., 2018; Gudmundsson et al., 2018a). The GSIM database gives data from 3 stream gauges along the Río Paraná and 12 from the Río Uruguay (Table 1; Fig. 3). To estimate Río de la Plata streamflow, I combine data from the two rivers. Records from the Río Paraná are long and complete. Therefore, I use the time series from Timbúes, which spans 1905–2014 and has the largest gauged area. Data from the Río Uruguay are shorter and more gappy; for example, the station with the largest drainage, Aporte Salto Grande, only gives data for 2012–2016. Since drainage area and mean streamflow are strongly correlated across stream gauges along this river (Pearson correlation coefficient >0.99; Table 1), I create a composite streamflow time series for the Río Uruguay covering 1931–2016 by averaging the available records after scaling each station's time series by the ratio of the total drainage area to the drainage monitored by that particular gauge. Summing the Río Paraná data at Timbúes and the composite Río Uruguay record gives a complete time series of Río de la Plata streamflow for 1931–2014 (Fig. 3).

Table 1 GSIM river-gauge records (Do et al., 2018; Gudmundsson et al., 2018a; Fig. 3). Long and Lat are degrees west longitude and south latitude, respectively. Completeness is the percentage of years during the span featuring data. Area is the gauged drainage area. Mean flow is the time mean streamflow over the record length.

2.2 Relative sea level

I use annual relative sea-level records from the Permanent Service for Mean Sea Level (PSMSL; Holgate et al., 2013; PSMSL, 2022). The PSMSL database extracted on 21 March 2022 provides long (>50-year) time series reduced to a common datum for six tide gauges from three regions in and around the Río de la Plata: Buenos Aires and Palermo towards the head of the estuary in Argentina, Montevideo and La Paloma near the mouth of the estuary along the coast of Uruguay to the north, and Mar del Plata and Quequén outside the estuary along coastal Argentina to the south (Table 2; Figs. 1, 4). To extend record length and to reduce dimensionality and errors, I average adjacent pairs of tide-gauge records relative to their common period, creating longer virtual station records (Dangendorf et al., 2017; Frederikse et al., 2021; Jevrejeva et al., 2014) at Buenos Aires (1905–2019), Montevideo (1938–2018), and Mar del Plata (1918–2019). For each station, I interrogate the period of overlap between virtual station and stream-gauge data.

Table 2 PSMSL tide-gauge records (Holgate et al., 2013; Figs. 1, 4). Long and Lat are degrees west longitude and south latitude, respectively. Completeness is the percentage of years during the span that feature data.

2.3 Late Holocene trends

To distinguish late Holocene trends related to background geological processes from modern rates of change due to ocean circulation and climate in the tide-gauge records, I use proxy reconstructions of relative sea level from Santa Catarina, Brazil, compiled by Milne et al. (2005) and originally reported by Angulo et al. (1999) based on Vermetid snails (Table 3). These mollusks are sea-level indicators because they grow formations between the infra- and mid-littoral zones, so formations fossilized in growth position are informative of low water (Laborel, 1986). Applying Bayesian linear regression to the data and accounting for the relative sea level and age errors, I determine a relative sea-level trend during the past 2000 years of $-\mathrm{0.54}\pm \mathrm{0.32}$ mm yr⁻¹ (95 % posterior credible interval); the Bayesian model is detailed in Appendix A. This negative rate of change arises from ocean siphoning and continental levering (Mitrovica and Milne, 2002), and past modeling studies of the glacial isostatic adjustment process report similar rates over the past few millennia (Caron et al., 2018; Peltier, 2004).

Table 3 Proxy sea-level reconstructions for the past 2 millennia from Santa Catarina (Milne et al., 2005). Milne et al. (2005) give calendar ages as min–max ranges, which I take to be 95 % confidence intervals. I take the center point as the best estimate and one-quarter of the range as 1 standard error. I also assume that sea-level errors given by Milne et al. (2005) correspond to 2 standard errors.

3 Results

Mean Río de la Plata streamflow is $(\mathrm{2.2}\pm \mathrm{0.1})\times {\mathrm{10}}^{\mathrm{4}}$ m³ s⁻¹ (Fig. 3), which is one of the largest river flows in the world, and consistent with values in past studies (Guerrero et al., 1997). Unless otherwise indicated, ± values identify 95 % bootstrap confidence intervals. The record standard deviation of $(\mathrm{4.8}\pm \mathrm{1.0})\times {\mathrm{10}}^{\mathrm{3}}$ m³ s⁻¹ quantifies variability across interannual to multidecadal timescales, including a long-term trend of 96±37 m³ s⁻¹ yr⁻¹, which has been reported on previously (Dai, 2016; Dai et al., 2009). Interannual variations in streamflow partly correspond to the El Niño–Southern Oscillation (ENSO); the correlation coefficient between streamflow and the Niño 3.4 Index (Rayner et al., 2003) is 0.33±0.16, and peak streamflow occurred during the 1982–1983 and 1997–1998 El Niños. Such relationships between streamflow and ENSO have been extensively documented (Berri et al., 2002; Cardoso and Silva Dias, 2006; Depetris et al., 1996; Grimm et al., 1998; Robertson and Mechoso, 1998; Ropelewski and Halpert, 1987). Also apparent is a regime shift from the late 1960s to early 1980s when streamflow increased substantially. This transition has been ascribed to increased precipitation and decreased evaporation over the drainage basin due to changes in land use, deforestation, and large-scale climate modes (Lawrence and Vandecar, 2015; Medvigy et al., 2011).

The virtual station data similarly show that relative sea level varies over all periods (Fig. 4). These records also exhibit spatial structure. Detrended series at Buenos Aires and Montevideo are significantly correlated with one another (correlation coefficient 0.44±0.20), but neither is correlated with the detrended record at Mar del Plata (coefficients $-\mathrm{0.01}\pm \mathrm{0.20}$ and 0.18±0.28, respectively). While the time series at Mar del Plata is uncorrelated with ENSO (correlation coefficient 0.12±0.17 with Niño 3.4), the records from Buenos Aires and Montevideo both show correlation with ENSO (coefficients 0.26±0.19 and 0.25±0.20 with Niño 3.4, respectively). These results are consistent with past studies (Douglas, 2001; Papadopoulous and Tsimplis, 2006; Raicich, 2008; Verocai et al., 2016) and suggest that there are processes that drive common sea-level changes at Buenos Aires and Montevideo but which do not affect sea level along Mar del Plata. Considering the longest timescales, I compute a long-term rate of change at Buenos Aires of 1.46±0.36 mm yr⁻¹ based on ordinary least squares linear regression, which is larger than the trends of 1.03±0.53 and 1.00±0.35 mm yr⁻¹ obtained for Montevideo and Mar del Plata, respectively (Fig. 7). These values agree with previous studies of regional sea-level rise cited in the Introduction. After adjusting for a late Holocene rate (Sect. 2.3; Fig. 5), I find an average sea-level trend across virtual stations of 1.70±0.40 mm yr⁻¹, which is faster than modern estimates of twentieth-century global mean sea-level rise referenced earlier and similar to conclusions from Frederikse et al. (2021).

Figure 4 Yearly tide-gauge relative sea-level records (Fig. 1, Table 2). Virtual station time series are shown as thick lines, and individual tide-gauge records are shown as thin lines. The time series are shifted vertically by an arbitrary amount for ease of visualization.

Figure 5 Proxy sea-level reconstructions (orange) and Bayesian linear regression (blue). Orange shading identifies best estimates plus and minus twice the standard errors. Blue shading corresponds to 95 % posterior credible intervals. The Bayesian model is detailed in the Appendix.

Streamflow explains a substantial portion of the sea-level variation at Buenos Aires and to a lesser extent Montevideo, and it largely accounts for the apparent faster-than-global rate of regional sea-level rise (Figs. 6–8). To quantify the influence of streamflow on sea level, I evaluate a multiple linear regression model at each virtual station with sea level as the dependent variable and streamflow, time, and unity as the independent variables. The streamflow regressor explains 59±17 %, 28±21 %, and $-\mathrm{6}\pm \mathrm{9}\mathit{\%}$ of the sea-level variance at Buenos Aires, Montevideo, and Mar del Plata, respectively (Fig. 6). This suggests that streamflow generally has more of an influence on sea level closer to the mouths of the Río Paraná and Río Uruguay. Regression coefficients between streamflow and sea level for Buenos Aires, Montevideo, and Mar del Plata are $(\mathrm{7.3}\pm \mathrm{1.8})\times {\mathrm{10}}^{-\mathrm{6}}$, $(\mathrm{4.7}\pm \mathrm{2.6})\times {\mathrm{10}}^{-\mathrm{6}}$, and $(-\mathrm{1.1}\pm \mathrm{1.6})\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s, respectively (Fig. 7). This structure shows that sea level is more sensitive to streamflow closer to the mouths of the rivers. Finally, linear trends computed from the virtual station time series from this regression model are 0.75±0.34, 0.56±0.58, and 1.11±0.37 mm yr⁻¹ at Buenos Aires, Montevideo, and Mar del Plata, respectively (Fig. 7). Compared to trends reported in the last paragraph, this implies that streamflow effected sea-level rates of 0.71±0.35, 0.48±0.38, and $-\mathrm{0.11}\pm \mathrm{0.17}$ mm yr⁻¹ at the respective virtual stations (Fig. 7). Averaging the streamflow-corrected sea-level trends and adjusting for the background geologic rate, I obtain a mean rate of 1.34±0.40 mm yr⁻¹, which is more in line with recent global mean sea-level trends for the past century from Hay et al. (2015), Dangendorf et al. (2017), and Frederikse et al. (2020).

Figure 6 Scatter plots comparing yearly average Río de la Plata streamflow (horizontal axes) and relative sea level (vertical axes) at Buenos Aires (blue), Montevideo (orange), and Mar del Plata (yellow). Sea-level values from the different sites are shifted vertically by an arbitrary amount for ease of visualization.

Figure 7 (a) Regression coefficients between sea level and streamflow found empirically from linear regression (blue) and predicted theoretically from ocean dynamics (orange). (b) Trend computed from tide gauges without (blue) and with (orange) adjusting for river effects. (c) Sea-level trend due to streamflow found empirically from linear regression (blue) and predicted theoretically from ocean dynamics given the streamflow trend (orange). To evaluate predicted values at Buenos Aires, I use a value of x=65 km from the source in Eq. (13).

To establish the robustness of these results, I also consider alternative regression models and analysis approaches. First, I evaluate the same regression model but using ridge regression. This is meant to account for collinearity between predictors (e.g., the linear trend in streamflow). Results obtained for a wide range of ridge parameter values are essentially identical to the results found from ordinary least squares discussed earlier (not shown). From this, I conclude that the model is well posed and that collinearity between streamflow and time does not present a serious issue. Second, I evaluate the same regression model using ordinary least squares but considering sea-level and streamflow data with ENSO effects removed prior to analysis. I remove ENSO effects by regressing the quantity of interest against the Niño 3.4 Index and its Hilbert transform to capture arbitrary phase relationships between quantities. If river effects on sea level are restricted to ENSO events, then results from this analysis should give no meaningful relationship between sea level and streamflow. However, in this analysis, I find very similar regression coefficients between sea level and streamflow ($(\mathrm{6.9}\pm \mathrm{1.7})\times {\mathrm{10}}^{-\mathrm{6}}$ at Buenos Aires; $(\mathrm{3.9}\pm \mathrm{2.6})\times {\mathrm{10}}^{-\mathrm{6}}$ at Montevideo; $(-\mathrm{1.3}\pm \mathrm{1.8})\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s at Mar del Plata) as well as sea-level variance explained by streamflow (55±18 % at Buenos Aires; 21±19 % at Montevideo; $-\mathrm{7}\pm \mathrm{10}\mathit{\%}$ at Mar del Plata) as previously when I did not remove ENSO effects prior to analysis. From this, I conclude that river effects on sea level in the Río de la Plata are not restricted to ENSO events, which have been the focus of past studies cited above, but are rather more general. This conclusion is consistent with results from wavelet coherence analysis (Grinsted et al., 2004). Streamflow is significantly coherent with sea level at Buenos Aires across most time periods and frequency bands resolved by the data, as well as with sea level at Montevideo for particular periods and frequencies (e.g., decadal scales generally, interannual scales during the 1990s and 2000s), but it is mostly incoherent with Mar del Plata sea level (Fig. 8).

Figure 8 Magnitude-squared wavelet coherence between streamflow and sea level at (a) Buenos Aires, (b) Montevideo, and (c) Mar del Plata. Solid black lines identify where values are significant at the 68 % confidence level, and white dashed lines mark the cone of influence. Statistical significance is determined from 1000 simulations with phase-scrambled versions of the observational data (Theiler et al., 1992). Values are based on the analytic Morlet wavelet.

4 Interpretation

Findings in the preceding section are based on correlation and regression analysis. They do not necessarily demonstrate that streamflow and coastal sea level are causally connected. To provide physical interpretation and establish causality, I develop simple theories for the relationship between streamflow and coastal sea level based on ocean dynamics in Sect. 4.1 and 4.2, and I compare model predictions to observational results in Sect. 4.3. To facilitate model development, I provide a schematic illustration of the region identifying key model quantities in Fig. 9.

Figure 9 Schematic of the study region with key model quantities identified. The top shows a plan view of the region. The bottom shows cross sections at various locations in and around the Río de la Plata. Locations a and b are upstream near Buenos Aires, where the barotropic theory developed in Sect. 4.1 applies. Location c is downstream of Montevideo near La Paloma, where the baroclinic theory developed in Sect. 4.2 applies. At the bottom, yellow ⊗ identifies flow into the page, blue indicates fresher less dense water, and green denotes saltier more dense water. Other symbols and quantities are as defined in the text of Sect. 4.1 and 4.2. Illustration by Natalie Renier, WHOI.

4.1 Theory for Buenos Aires

Around Buenos Aires and Palermo, the Río de la Plata is relatively shallow, narrow, and fresh (Guerrero et al., 1997; Fig. 9). To model sea level in this region, I use the following conservation laws:

$$\begin{array}{}\text{(1)}& {\displaystyle }{u}_x+v_y+w_z=\mathrm{0},\text{(2)}& {\displaystyle }{p}_z=-{\mathit{\rho }}_{\mathrm{f}}g,\text{(3)}& {\displaystyle }\mathrm{0}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{f}}}}{p}_x+\mathit{\nu }{u}_{zz}.\end{array}$$

Here u, v, and w are velocities in along-estuary (x), across-estuary (y), and vertical (z) directions, respectively, p is hydrostatic pressure, ρ_f is a reference freshwater density, g is acceleration due to gravity, ν is kinematic viscosity, and x, y, and z subscripts are spatial derivatives (Fig. 9). Equations (1) and (2) are familiar forms of the continuity equation and hydrostatic balance (Gill, 1982). Equation (3) specifies along-estuary momentum conservation in terms of a balance between pressure gradient and viscous forces; it omits the time tendency given the long periods under consideration, and it also neglects nonlinear advection and Coriolis acceleration under the assumptions of a small Reynolds number and large Ekman number, which are reasonable given the spatial scales of the problem.

Integrating Eq. (1) over the depth H(x) and width W(x) of the estuary (Fig. 9), applying kinematic boundary conditions at the bottom and along the sides, and ignoring the time tendency gives

$$\begin{array}{}\text{(4)}& {\left(\langle \stackrel{\mathrm{‾}}{u}\rangle WH\right)}_x=\mathrm{0},\end{array}$$

where the overbar and bracket are depth and across-estuary average, respectively. Integrating Eq. (2) vertically, substituting into Eq. (3), and averaging over depth and width yields

$$\begin{array}{}\text{(5)}& \mathrm{0}=-g\langle \mathit{\zeta }{\rangle }_x-{\displaystyle \frac{C_{\mathrm{d}}U}{H}}\langle \stackrel{\mathrm{‾}}{u}\rangle ,\end{array}$$

where ζ is ocean dynamic sea level, C_d is a drag coefficient, and U is a reference velocity scale. To obtain Eq. (5), I assumed that the ζ slope across the estuary is linear and that

$$\begin{array}{}\text{(6)}& v{u}_z=C_{\mathrm{d}}U\stackrel{\mathrm{‾}}{u}\end{array}$$

along the bottom. To solve Eqs. (4) and (5) for 〈ζ〉, I specify that along-estuary transport equals the streamflow q at the origin,

$$\begin{array}{}\text{(7)}& \langle \stackrel{\mathrm{‾}}{u}\rangle WH=q\mathrm{at}x=\mathrm{0},\end{array}$$

and that 〈ζ〉 vanishes far from the source:

$$\begin{array}{}\text{(8)}& \underset{x\to \mathrm{\infty }}{lim}\langle \mathit{\zeta }\rangle =\mathrm{0}.\end{array}$$

Combining Eqs. (4) and (7), substituting for $\langle \stackrel{\mathrm{‾}}{u}\rangle$ in Eq. (5), integrating along the estuary from x to ∞, and applying the boundary condition from Eq. (8) gives

$$\begin{array}{}\text{(9)}& \langle \mathit{\zeta }\rangle ={\displaystyle \frac{C_{\mathrm{d}}Uq}{g}}\underset{x}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\mathrm{1}}{H^{\mathrm{2}}W}}\mathrm{d}{x}^{\prime }\end{array}$$

for arbitrary depth and width profiles. For an estuary with exponential width and depth (Fig. 2),

$$\begin{array}{}\text{(10)}& {\displaystyle }W=W_{\mathrm{0}}\exp \left(x/L_W\right),\text{(11)}& {\displaystyle }H=H_{\mathrm{0}}\exp \left(x/L_H\right),\end{array}$$

where W₀ and H₀ are initial values and L_W and L_H are length scales, the solution to Eq. (9) is

$$\begin{array}{}\text{(12)}& \langle \mathit{\zeta }\rangle ={\left({\displaystyle \frac{\mathrm{2}}{L_H}}+{\displaystyle \frac{\mathrm{1}}{L_W}}\right)}^{-\mathrm{1}}{\displaystyle \frac{C_{\mathrm{d}}Uq}{g{H}^{\mathrm{2}}W}}.\end{array}$$

The 〈ζ〉 response is linear in q and controlled by friction and the geometry of the estuary; it is larger for stronger friction C_dU, narrower initial width W₀, shallower initial depth H₀, longer width and depth scales L_W and L_H, and decays rapidly with distance from the origin (Fig. 10).

Regression coefficients computed between sea-level and streamflow data (Fig. 7) can be understood as approximate observational estimates of the derivative of the former with respect to the latter. From Eq. (12), it follows that

$$\begin{array}{}\text{(13)}& \langle \mathit{\zeta }{\rangle }_q={\left({\displaystyle \frac{\mathrm{2}}{L_H}}+{\displaystyle \frac{\mathrm{1}}{L_W}}\right)}^{-\mathrm{1}}{\displaystyle \frac{C_{\mathrm{d}}U}{g{H}^{\mathrm{2}}W}}.\end{array}$$

Below, I numerically evaluate Eq. (13) and compare the values to the empirically determined regression coefficients to test whether the theory is consistent with the observations.

Figure 10 Sea-level response 〈ζ〉 to streamflow forcing q described by Eq. (12) as a function of distance along the estuary away from the mouth of the rivers for different values of (a) streamflow q, (b) friction C_dU, (c) depth length scale L_H, (d) initial depth H₀, (e) width length scale L_W, and (f) initial width W₀. Default values are $q=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ m³ s⁻¹, C_dU=0.001 m s⁻¹, L_H=150 km, H₀=2 m, L_W=150 km, and W₀=30 km.

4.2 Theory for Montevideo

The solution for Buenos Aires (Eq. 12) is not applicable to Montevideo. The estuary becomes wider, deeper, and more saline by this point (Guerrero et al., 1997; Figs. 1, 2, 9), so stratification and rotation effects cannot be neglected as they were previously. I develop a theory for the ζ response at Montevideo based on past studies of bottom-advected (slope-controlled) plumes (Chapman and Lentz, 1994; Lentz and Helfrich, 2002; Yankovsky and Chapman, 1997). I take x, y, and z to be the alongshore, across-shore, and vertical coordinates, respectively. As a mental model, I envision a narrow alongshore jet over a sloping bottom H(y) in thermal wind balance with a sharp density front at a location y_p offshore (Fig. 9). I imagine the jet transport includes both the fresh river water and salty ocean water brought into the plume by turbulent mixing. These features are represented by the following governing equations:

$$\begin{array}{}\text{(14)}& {\displaystyle }fu=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}}}{p}_y,\text{(15)}& {\displaystyle }{p}_z=-\mathit{\rho }g,\text{(16)}& {\displaystyle }Q=q+E,\text{(17)}& {\displaystyle \frac{Q}{H}}\underset{-H}{\overset{\mathrm{0}}{\int }}\mathit{\rho }(y,z)\mathrm{d}z=q{\mathit{\rho }}_{\mathrm{f}}+E{\mathit{\rho }}_{\mathrm{0}},\end{array}$$

where f=2Ωsin ϕ is the Coriolis frequency for Earth rotation rate Ω and latitude ϕ, ρ₀ is an ambient ocean density, Q is volume transport of the vertically sheared geostrophic jet, and E is entrainment flux. Equations (14) and (15) are geostrophic and hydrostatic balances, respectively. Equation (16) is a form of the continuity equation, which states that volume is conserved within the jet. Density conservation in Eq. (17) is equivalent to steady-state heat and salt conservation for a linear equation of state¹. Boundary conditions are that alongshore velocity vanishes everywhere along the bottom and that velocity shear is zero at the foot of the front (Chapman and Lentz, 1994; Lentz and Helfrich, 2002; Yankovsky and Chapman, 1997)

$$\begin{array}{}\text{(18)}& {\displaystyle }u=\mathrm{0}\mathrm{at}z=-H\left(y\right),\forall y\text{(19)}& {\displaystyle }{u}_z=\mathrm{0}\mathrm{at}y=y_p,z=-H\left(y_p\right)\doteq -H_p\end{array}$$

A solution to Eqs. (14)–(19) is obtained by giving a functional form to the density field. I picture an infinitely narrow front, with ambient ocean density everywhere offshore, and a mixture of fresh river water and salty ocean water onshore of the front, which I model as (Fig. 11a, b)

$$\begin{array}{}\text{(20)}& \mathit{\rho }(y,z)={\mathit{\rho }}_{\mathrm{0}}+{\displaystyle \frac{{\mathit{\rho }}^{\prime }}{H_p}}\left(z+H_p\right)\left[\mathcal{H}\left(y_p-y\right)-\mathrm{1}\right],\end{array}$$

where ρ^′ is a density increment and ℋ is the Heaviside step function. The alongshore velocity field in thermal wind balance with this density structure, obtained by cross-differentiating Eqs. (14) and (15) and then integrating vertically subject to the boundary conditions, is

$$\begin{array}{}\text{(21)}& u(y,z)=-{\displaystyle \frac{g{\mathit{\rho }}^{\prime }}{\mathrm{2}{\mathit{\rho }}_{\mathrm{0}}f{H}_p}}{\left(z+H_p\right)}^{\mathrm{2}}\mathit{\delta }\left(y_p-y\right),\end{array}$$

where δ is the Dirac delta (Fig. 11c).

Figure 11 Idealized (a) density structure onshore of the front (Eq. 20), (b) density structure offshore of the front (Eq. 21), and (c) velocity structure within the front (Eq. 21) as a function of depth. Sea-level response ζ described by Eq. (28) as a function of (d) streamflow q, (e) latitude ϕ, and (f) ambient ocean density ρ₀. Default values: $q=\mathrm{2}\times {\mathrm{10}}^{\mathrm{4}}$ m³ s⁻¹, $\mathit{\varphi }=\mathrm{35}{}^{\circ }$, ρ₀=1030 kg m⁻³.

To obtain the sea-level solution corresponding to Eq. (21), I integrate geostrophic balance at the surface,

$$\begin{array}{}\text{(22)}& fu=-g{\mathit{\zeta }}_y,\end{array}$$

over all offshore locations, which gives

$$\begin{array}{}\text{(23)}& \mathit{\zeta }={\displaystyle \frac{{\mathit{\rho }}^{\prime }{H}_p}{\mathrm{2}{\mathit{\rho }}_{\mathrm{0}}}}\left[\mathrm{1}-\mathcal{H}\left(y_p-y\right)\right].\end{array}$$

That is, ζ takes on a constant value of ${\mathit{\rho }}^{\prime }{H}_p/\mathrm{2}{\mathit{\rho }}_{\mathrm{0}}$ onshore of the front, experiences a step change at the front, and vanishes offshore of the front. The ζ solution can be written more explicitly in terms of streamflow q and river and ocean densities ρ_f and ρ₀ as follows. First, I express Q in terms of q and density. Given Eq. (20), the vertically averaged density within the front is

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\mathrm{1}}{H}}\underset{-H}{\overset{\mathrm{0}}{\int }}\mathit{\rho }(y_p,z)\mathrm{d}z={\mathit{\rho }}_{\mathrm{0}}-{\displaystyle \frac{{\mathit{\rho }}^{\prime }}{\mathrm{4}}},\end{array}$$

which, substituting into Eq. (17) and combining with Eq. (16) to eliminate E, implies

$$\begin{array}{}\text{(25)}& Q={\displaystyle \frac{\mathrm{4}q\left({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{f}}\right)}{{\mathit{\rho }}^{\prime }}},\end{array}$$

which is analogous to a form of Knudsen's hydrographical theorem (Dyer, 1997). Second, I solve for H_p in terms of Q and density. Integrating both sides of Eq. (21) over all depths and offshore locations and rearranging gives

$$\begin{array}{}\text{(26)}& Q=-{\displaystyle \frac{g{\mathit{\rho }}^{\prime }{H}_p^{\mathrm{2}}}{\mathrm{6}{\mathit{\rho }}_{\mathrm{0}}f}},\end{array}$$

or, after rearranging and solving for H_p (and recalling that f<0 in the Southern Hemisphere),

$$\begin{array}{}\text{(27)}& H_p={\left(-{\displaystyle \frac{\mathrm{6}Qf{\mathit{\rho }}_{\mathrm{0}}}{g{\mathit{\rho }}^{\prime }}}\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Finally, I substitute Eq. (25) for Q in Eq. (27), insert the resulting expression for H_p in Eq. (23), and cancel common terms to give

$$\begin{array}{}\text{(28)}& \mathit{\zeta }={\left[-{\displaystyle \frac{\mathrm{6}fq\left({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{f}}\right)}{{\mathit{\rho }}_{\mathrm{0}}g}}\right]}^{\mathrm{1}/\mathrm{2}}\left[\mathrm{1}-\mathcal{H}\left(y_p-y\right)\right].\end{array}$$

The ζ response is nonlinear in q and controlled by stratification and rotation; it is larger for higher latitude, stronger streamflow, and a sharper density contrast (Fig. 11d–f). While there is no alongshore dependence in Eq. (28), it assumes that the location of interest is downstream in the far field of the river mouth. Given Eq. (28), the derivative of ζ with respect to q, which can be evaluated numerically and compared to regression coefficients from observations, is

$$\begin{array}{}\text{(29)}& {\mathit{\zeta }}_q={\left[-{\displaystyle \frac{\mathrm{3}f\left({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{f}}\right)}{\mathrm{2}{\mathit{\rho }}_{\mathrm{0}}gq}}\right]}^{\mathrm{1}/\mathrm{2}}\left[\mathrm{1}-\mathcal{H}\left(y_p-y\right)\right].\end{array}$$

4.3 Model–data comparison

To test whether empirical results from Sect. 3 are consistent with theories developed in Sect. 4.1 and 4.2, I evaluate Eq. (13) for Buenos Aires and Eq. (29) for Montevideo using parameter values in Table 4 and then compare the predictions to the observed values (Fig. 7). See Appendix B for a discussion of the numerical values of the model parameters in Table 4. Equation (13) gives a theoretical regression coefficient between streamflow and sea level for Buenos Aires of $(\mathrm{7.0}\pm \mathrm{4.0})\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s, where the error bar reflects uncertainties in the parameter values (Table 4). Multiplying this coefficient by the long-term trend in streamflow estimated earlier (96±37 m³ s⁻¹ yr⁻¹), I obtain an expected sea-level trend at Buenos Aires due to streamflow of 0.68±0.47 mm yr⁻¹. These theoretical estimates agree with the coefficient of $(\mathrm{7.3}\pm \mathrm{1.8})\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s and the streamflow-driven sea-level trend of 0.71±0.35 mm yr⁻¹ found earlier from regression analysis of observed streamflow and sea level at Buenos Aires (Fig. 7). Following the same approach and evaluating Eq. (29), I find a theoretical regression coefficient of $(\mathrm{4.0}\pm \mathrm{0.1})\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s and an anticipated sea-level trend forced by streamflow of 0.41±0.19 mm yr⁻¹ for Montevideo. Again, these values from first principles are consistent with the regression coefficient of $(\mathrm{4.8}\pm \mathrm{2.7})\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s and the streamflow-induced sea-level trend of 0.48±0.38 mm yr⁻¹ found from the observational data (Fig. 7). The consistency between theory and observation suggests that the statistical connections found earlier between measured streamflow and sea level at Buenos Aires and Montevideo identify cause-and-effect relationships, which are consistent with the physics prescribed above.

Table 4 Parameter values used to evaluate Eqs. (13) and (29). See Appendix B for a detailed discussion on how the numerical values of the model parameters were selected.

The lack of a significant relation between streamflow and sea level in Mar del Plata in the data (Figs. 6–8) is also consistent with the theories developed in Sect. 4.1 and 4.2. The response described by Eq. (12) imagines a rapid decay away from the rivers. Indeed, given its strong exponential dependence, the sea-level response predicted by this theory is vanishingly small at Mar del Plata (Figs. 1, 10). The response described by Eq. (28) envisions coastal sea level coupled to a buoyant longshore current in the sense of coastal waves: counterclockwise along the Uruguay coast and then equatorward along the Brazil coast (Piola et al., 2005). In other words, given this mechanism, Mar del Plata is not downstream of the Río de la Plata, and hence no signals are communicated between the two locations according to these physics.

Figure 12 (a) Correlation coefficient and (b) regression coefficient (m m⁻³ s) between annual streamflow in the Río de la Plata (Fig. 3) and sea surface height anomaly from along-track satellite altimetry data (Birol et al., 2017) during 1993–2014 over the study regions. Values are only shown where correlation coefficients are positive at the 95 % confidence level determined through bootstrapping. Contours identify the 20, 50, 100, 200, and 500 m isobaths.

5 Conclusions

The Río de la Plata estuary in South America features the longest tide-gauge records in the South Atlantic Ocean (Figs. 1, 2). However, the causes of long-term relative sea-level changes in this region have not been firmly established. I interrogated data (Figs. 3–5) and developed theories (Figs. 8–11) to argue for cause-and-effect relationships between low-frequency streamflow and sea-level changes in the Río de la Plata over 1931–2014 (Figs. 6, 7). Streamflow forcing generally explained one-half of the sea-level variance on interannual and longer timescales observed at Buenos Aires and one-quarter of the sea-level variance at Montevideo over the study period. Specifically, a trend in streamflow of ∼100 m³ s⁻¹ yr⁻¹ during the past century caused sea level to rise at rates of ∼0.7 mm yr⁻¹ at Buenos Aires and ∼0.5 mm yr⁻¹ at Montevideo. These findings advance understanding of local, regional, and global sea-level changes; clarify basic sea-level physics; inform future projections of coastal sea-level change as well as the interpretation of satellite data and proxy reconstructions; and highlight future research directions.

This paper complements past tide-gauge studies on mean sea-level changes in the Río de la Plata on interannual to centennial timescales (e.g., Aubrey et al., 1988; Brandani et al., 1985; D'Onofrio et al., 2008; Dennis et al., 1995; Douglas, 1997, 2001, 2008; Emery and Aubrey, 1991; Fiore et al., 2009; Frederikse et al., 2021; Isla, 2008; Lanfredi et al., 1998; Meccia et al., 2009; Melini et al., 2004; Papadopoulous and Tsimplis, 2006; Pousa et al., 2007; Raicich, 2008; Santamaria-Aguilar et al., 2017; Thompson et al., 2016; Verocai et al., 2016). Previous authors have established that streamflow and sea level in the Río de la Plata covary on interannual timescales during ENSO events, but they do not identify the causal mechanisms responsible for the observed statistical correlations, nor do they consider how these two variables correspond more generally on longer timescales. My paper builds on their foundation by showing that river effects on sea level are not restricted to ENSO events in particular, but are also apparent more generally at multidecadal and centennial periods, and by identifying ocean dynamic mechanisms that mediate the relationship between streamflow and sea level. These results corroborate the hypothesis due to Douglas (2001) that interannual sea-level variation at Buenos Aires over the 1982–1983 El Niño can be understood in terms of ocean dynamic processes, but they do not necessarily falsify suggestions that contemporary gravitational, rotational, and deformational effects also played a role (Isla, 2008; Thompson et al., 2016). Likewise, while they suggest that streamflow changes contributed importantly to long-term sea-level rise observed at Buenos Aires and Montevideo, these results do not rule out the possibility that other geophysical processes also effected regional sea-level trends (Melini et al., 2004; Aubrey et al., 1988).

My results have implications for twentieth-century global sea-level reconstructions and budgets (e.g., Church and White, 2011; Dangendorf et al., 2017; Frederikse et al., 2018, 2020, 2021; Hamlington and Thompson, 2015; Hay et al., 2015; Jevrejeva et al., 2014; Natarov et al., 2017; Ray and Douglas, 2011; Thompson and Merrifield, 2014; Thompson et al., 2016). The streamflow-driven sea-level effects highlighted here are local to regional in scale; they do not contribute meaningfully to sea-level changes on basin or global scales. Hence, such river effects on tide gauges in the Río de la Plata should be removed prior to analysis if the data are used in large-scale circulation and climate studies, lest this local or regional “noise” alias onto the basin or global “signal” of interest (e.g., Papadopoulous and Tsimplis, 2006; Thompson et al., 2016). Given the heavy weight placed on tide gauges from the Río de la Plata, streamflow-driven ocean dynamics could contribute to the lack of sea-level-budget closure and faster-than-global trends across the South Atlantic during the twentieth century found by Frederikse et al. (2018, 2021). Since tide-gauge records in and around the Río de la Plata are the main (if not sole) data constraint in the South Atlantic prior to 1950 in twentieth-century global mean sea-level reconstructions (Fig. 1b in Hamlington and Thompson, 2015; Fig. S1a in Dangendorf et al., 2017), it would be informative to estimate twentieth-century global mean sea-level rise from tide-gauge records adjusted for river effects, which are typically not considered in global budgets and reconstructions.

Theories developed here (Eqs. 13 and 29) clarify relationships between streamflow and coastal sea level, the physics of which have not been well understood (Durand et al., 2019). Piecuch et al. (2018a) formulate a theory for the far-field coastal sea-level response to buoyant river discharge in the limit of a pure surface-advected plume (their Eqs. 5 and 6). This study improves upon their work in two ways. First, I developed a barotropic theory for the sea-level response within an estuary (Eq. 13), where frictional effects and the shape of coastlines and bathymetry are important. Second, I formulated a far-field theory for the coastal sea-level adjustment in the alternative limit of a purely bottom-advected (or slope-controlled) plume (Eq. 29), which is more germane to the problem at hand. These new theories allow the relationship between sea level and river discharge to be studied in a wider range of settings. In a future study, I plan to develop a more general far-field theory for the buoyancy-driven sea-level response to an intermediate buoyant plume that falls between the extremes of a surface-advected plume and a bottom-advected plume (Yankovsky and Chapman, 1997; Lentz and Helfrich, 2002).

I demonstrated that the sea-level response to buoyant coastal discharge can depend sensitively on density gradients over short scales and the geometry of coastlines and bathymetry. With some exceptions (Haarsma et al., 2016), the current generation of coupled models used for climate projections is too coarsely resolved to represent such features (Holt et al., 2017). Theories developed here may be helpful in this regard. Equations (13) and (29) may be instructive for obtaining basic scales and magnitudes of future coastal sea-level changes due to streamflow, assuming that the details of coastlines and bathymetry are known and given projected changes in continental freshwater runoff into the coastal ocean. In the future, as global climate models with improved representation of the coastal ocean and shelf seas become more widely available, numerical experiments could be performed to broadly test analytical predictions made here regarding relationships between sea level and streamflow.

Due to my focus on long-term trends, I interrogated sea-level records from tide gauges. However, streamflow-driven sea-level changes are also apparent in data from other observing systems, including satellite altimetry. Comparing annual streamflow and sea surface height anomaly from along-track altimetry over 1993–2014 (Birol et al., 2017), I observe a region of significant correlation between the two variables extending broadly over the Uruguay coast from Montevideo past La Paloma towards Brazil and onshore of the ∼100 m isobath (Fig. 12a; see Fig. 1). The shape of the region mirrors the structure of low-salinity water near the mouth of the estuary (e.g., Piola et al., 2005). Regression coefficients obtained between Río de la Plata streamflow and sea surface height anomaly are consistent with theoretical expectations: more upstream in the estuary, values are $\lesssim \mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m m⁻³ s, similar to predictions from barotropic theory developed in Sect. 4.1 (Eq. 13), whereas values downstream in the far field are $\sim \mathrm{4}\times {\mathrm{10}}^{-\mathrm{6}}$ m m⁻³ s, consistent with values anticipated from the baroclinic theory from Sect. 4.2 (Eq. 29) (Fig. 12b; see Fig. 7). The offshore extent of the region of significant correlation between streamflow and sea surface height also corroborates basic theoretical expectations: for strong slope control and large river discharge, the offshore and vertical scales of a buoyant coastal plume are expected to be ∼100 km and ∼100 m, respectively (e.g., Yankovsky and Chapman, 1997).

Findings here may have implications for proxy reconstructions of late Holocene sea level from natural archives, which have temporal resolution of decades to centuries (e.g., Kemp et al., 2009; Khan et al., 2019). Whereas past studies reason that river effects contribute to sea-level variability on interannual and shorter timescales (e.g., Durand et al., 2019; Woodworth et al., 2019), I showed that streamflow changes can be an important driver of sea-level changes over multidecadal and longer periods. This result has (at least) two important implications for proxy reconstructions. First, it implies that river effects may be important to consider when interpreting proxy sea-level reconstructions from large rivers or estuaries (e.g., Gerlach et al., 2017; Kemp et al., 2018). Second, it suggests that proxy sea-level reconstructions produced from strategic locations may inform past changes in streamflow and thus complement estimates from more traditional archives like tree rings (e.g., Margolis et al., 2011; Devineni et al., 2013).

Other major rivers including the Mississippi, Yenisey, and Lena have also undergone significant streamflow trends in the past century (e.g., Dai, 2016; Dai and Trenberth, 2002; Dai et al., 2009). However, the effect of these historical changes in streamflow on long-term sea-level change has not been considered. Future studies should take advantage of the growing number of available runoff and streamflow datasets (e.g., Do et al., 2018; Gudmundsson et al., 2018a; Tsujino et al., 2018) to test the analytical models developed here and observationally constrain river effects on historical sea-level rise more globally, which could inform studies of ocean circulation and climate change.

Appendix A: Bayesian hierarchical model

I apply Bayesian linear regression to proxy reconstructions from Milne et al. (2005) to quantify late Holocene rates of sea-level change. Bayesian linear regression is chosen over more traditional approaches like least squares or maximum likelihood because Bayesian methods provide a more transparent means for incorporating data errors into the formal uncertainty quantification. I design the Bayesian hierarchical model following similar algorithms developed in past studies (Ashe et al., 2019; Cahill et al., 2015, 2016; Walker et al., 2020). The model used here is essentially the time component of the space–time model from Piecuch et al. (2018b). While I give a brief description for the sake of completeness, readers are referred to Piecuch et al. (2018b) for a more detailed presentation.

Temporal Bayesian hierarchical models comprise three levels: a process level that prescribes the temporal evolution of the sea-level process, a data level that codifies the relationship between the uncertain proxy reconstructions and the sea-level process, and a parameter level at which prior constraints are specified.

For the process level, I model sea level $\mathit{y}={\left[y_{\mathrm{1}},y_{\mathrm{2}},\mathrm{\dots },y_n\right]}^{\mathsf{T}}$ as a linear function of time $\mathit{x}={\left[x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots },x_n\right]}^{\mathsf{T}}$ according to

$$\begin{array}{}\text{(A1)}& y_k\sim \mathcal{N}\left(\mathit{\alpha }{x}_k+\mathit{\beta },{\mathit{\gamma }}^{\mathrm{2}}\right),k\in \left[\mathrm{1},n\right],\end{array}$$

where ∼ means “is distributed as”, 𝒩(a,b²) is the normal distribution with mean a and variance b², and α, β, and γ² are uncertain slope, intercept, and residual variance parameters, respectively. For the data level, I represent the proxy reconstructions of relative sea level $\mathit{z}={\left[z_{\mathrm{1}},z_{\mathrm{2}},\mathrm{\dots },z_n\right]}^T$ and age $\mathit{w}={\left[w_{\mathrm{1}},w_{\mathrm{2}},\mathrm{\dots },w_n\right]}^T$ as noisy versions of the respective processes, viz.,

$$\begin{array}{}\text{(A2)}& {\displaystyle }{z}_k\sim \mathcal{N}\left(y_k,{\mathit{\delta }}_k^{\mathrm{2}}\right),\text{(A3)}& {\displaystyle }{w}_k\sim \mathcal{N}\left(x_k,{\mathit{ϵ}}_k^{\mathrm{2}}\right),\end{array}$$

where ${\mathit{\delta }}_k^{\mathrm{2}}$ and ${\mathit{ϵ}}_k^{\mathrm{2}}$ are the data error variances, which are provided (Table 3). To close the model, I assume normal priors for α and β and an inverse-gamma prior for γ²,

$$\begin{array}{}\text{(A4)}& {\displaystyle }\mathit{\alpha }\sim \mathcal{N}\left(\stackrel{\mathrm{̃}}{\mathit{\mu }},{\stackrel{\mathrm{̃}}{\mathit{\kappa }}}^{\mathrm{2}}\right),\text{(A5)}& {\displaystyle }\mathit{\beta }\sim \mathcal{N}\left(\stackrel{\mathrm{̃}}{\mathit{\eta }},{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}\right),\text{(A6)}& {\displaystyle }{\mathit{\gamma }}^{\mathrm{2}}\sim {\mathcal{G}}^{-\mathrm{1}}\left(\stackrel{\mathrm{̃}}{\mathit{\xi }},\stackrel{\mathrm{̃}}{\mathit{\chi }}\right),\end{array}$$

where tildes identify fixed hyperparameters (see below for numerical values).

Given Bayes' rule and the model equations, I assume the posterior distribution is

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& p\left(\mathit{y},\mathit{x},\mathit{\alpha },\mathit{\beta },{\mathit{\gamma }}^{\mathrm{2}}|\mathit{z},\mathit{w}\right)\propto p\left(\mathit{\alpha }\right)p\left(\mathit{\beta }\right)p\left({\mathit{\gamma }}^{\mathrm{2}}\right)\\ & \prod _{k=\mathrm{1}}^n\left[p\left(z_k|y_k\right)p\left(w_k|x_k\right)p\left(y_k|x_k,\mathit{\alpha },\mathit{\beta },{\mathit{\gamma }}^{\mathrm{2}}\right)\right],\end{array}\end{array}$$

where p is probability, | is conditionality, and ∝ indicates “proportional to”. To evaluate the posterior, I use a Gibbs sampler (Gelman et al., 2013), evaluating the full posteriors (Wikle and Berliner, 2007):

$$\begin{array}{}\text{(A8)}& {\displaystyle }\begin{array}{rl}& \mathit{\alpha }|\cdot \sim \mathcal{N}({\left[{\stackrel{\mathrm{̃}}{\mathit{\kappa }}}^{-\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\sum _{k=\mathrm{1}}^n{x}_k^{\mathrm{2}}\right]}^{-\mathrm{1}}\\ & \left[{\stackrel{\mathrm{̃}}{\mathit{\kappa }}}^{-\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\mu }}+{\mathit{\gamma }}^{-\mathrm{2}}\sum _{k=\mathrm{1}}^n{x}_k\left\{y_k-\mathit{\beta }\right\}\right],\\ & {\left[{\stackrel{\mathrm{̃}}{\mathit{\kappa }}}^{-\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\sum _{k=\mathrm{1}}^n{x}_k^{\mathrm{2}}\right]}^{-\mathrm{1}}),\end{array}\text{(A9)}& {\displaystyle }\begin{array}{rl}& \mathit{\beta }|\cdot \sim \mathcal{N}({\left[{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{-\mathrm{2}}+n{\mathit{\gamma }}^{-\mathrm{2}}\right]}^{-\mathrm{1}}\\ & \left[{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{-\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\eta }}+{\mathit{\gamma }}^{-\mathrm{2}}\sum _{k=\mathrm{1}}^n\left\{y_k-\mathit{\alpha }{x}_k\right\}\right],\\ & [{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{-\mathrm{2}}+n{\mathit{\gamma }}^{-\mathrm{2}}{]}^{-\mathrm{1}}),\end{array}\text{(A10)}& {\displaystyle }{\mathit{\gamma }}^{\mathrm{2}}|\cdot \sim {\mathcal{G}}^{-\mathrm{1}}\left(\stackrel{\mathrm{̃}}{\mathit{\xi }}+{\displaystyle \frac{n}{\mathrm{2}}},\stackrel{\mathrm{̃}}{\mathit{\chi }}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{k=\mathrm{1}}^n{\left[y_k-\mathit{\alpha }{x}_k-\mathit{\beta }\right]}^{\mathrm{2}}\right),\text{(A11)}& {\displaystyle }\begin{array}{rl}& y_k|\cdot \sim \mathcal{N}({\left[{\mathit{\delta }}_k^{-\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\right]}^{-\mathrm{1}}\\ & \left[{\mathit{\delta }}_k^{-\mathrm{2}}{z}_k+{\mathit{\gamma }}^{-\mathrm{2}}\left\{\mathit{\alpha }{x}_k+\mathit{\beta }\right\}\right],\\ & {\left[{\mathit{\delta }}_k^{-\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\right]}^{-\mathrm{1}}),\end{array}\text{(A12)}& {\displaystyle }\begin{array}{rl}& x_k|\cdot \sim \mathcal{N}({\left[{\mathit{ϵ}}_k^{-\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\gamma }}^{-\mathrm{2}}\right]}^{-\mathrm{1}}\\ & \left[{\mathit{ϵ}}_k^{-\mathrm{2}}{w}_k+{\mathit{\gamma }}^{-\mathrm{2}}\mathit{\alpha }\left\{y_k-\mathit{\beta }\right\}\right],\\ & {\left[{\mathit{ϵ}}_k^{-\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}{\mathit{\gamma }}^{-\mathrm{2}}\right]}^{-\mathrm{1}}),\end{array}\end{array}$$

where $|\cdot$ is conditionality on all other processes, parameters, and data. I set weak, uninformative priors ($\stackrel{\mathrm{̃}}{\mathit{\mu }}=\mathrm{0}$ mm yr⁻¹, ${\stackrel{\mathrm{̃}}{\mathit{\kappa }}}^{\mathrm{2}}=\mathrm{0.001}$ mm² yr⁻², $\stackrel{\mathrm{̃}}{\mathit{\eta }}=\mathrm{0}$ m, ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}=\mathrm{100}$ m², $\stackrel{\mathrm{̃}}{\mathit{\xi }}=\mathrm{0.5}$, $\stackrel{\mathrm{̃}}{\mathit{\chi }}=\mathrm{0.02}$ m²). I discard 1000 burn-in draws to eliminate startup transients. I reduce autocorrelation of the samples by keeping only every 10th draw of the subsequent 10 000 iterations of the Gibbs sampler. This gives a 1000-member ensemble of posterior estimates for y, x, α, β, and γ². Figure 5 shows summary statistics for the posterior solution of αx+β for x from 500 BCE to present.

Appendix B: Model parameters

To evaluate Eqs. (13) and (29), numerical values need to be assigned to the various model parameters. Here I detail the rationale behind the values tabulated in Table 4.

I chose standard reference values for freshwater density ρ_f, seawater density ρ₀, gravitational acceleration g, and the Coriolis parameter at the latitude of the Río de la Plata.

Ranges on the initial estuary width W₀ and depth H₀ as well as the width and depth scales L_W and L_H from Table 4 correspond to best estimates plus and minus 2 standard errors on these parameter values as determined by fitting exponentials in the form of Eqs. (10) and (11) to the bathymetry data in Fig. 2 using nonlinear least squares as described in the Fig. 2 caption.

The q value is the time mean of the Río de la Plata streamflow time series in Fig. 3 plus and minus twice its standard error.

Typical values for C_d range from 0.001 to 0.003 (Adcroft et al., 2018, 2019). I selected a middle-of-the-road value of 0.002.

The velocity scale U parameterizes the influence of unresolved processes and is typically selected to represent tidal motions (Adcroft et al., 2019). Regional tidal current amplitudes are 0.5–0.8 m s⁻¹ (O'Connor, 1991; Piedra-Cueva and Fossati, 2007). Multiplying by a factor of $\mathrm{2}/\mathit{\pi }$, the average amplitude of a sine wave, gives the range of 0.3–0.5 m s⁻¹ used here. This is larger than the background mean flow from river discharge, $\langle \stackrel{\mathrm{‾}}{u}\rangle =q/H\left(x\right)W\left(x\right)$, which is 0.12 m s⁻¹ at Buenos Aires (x=65 km), using the q, H₀, W₀, L_H, and L_W values in Table 4 discussed earlier.