Travertine crystal growth ripples record the hydraulic history of ancient Rome’s Anio Novus aqueduct

Travertine crystal growth ripples are used to reconstruct the early hydraulic history of the Anio Novus aqueduct of ancient Rome. These crystalline morphologies deposited within the aqueduct channel record the hydraulic history of gravity-driven turbulent flow at the time of Roman operation. The wavelength, amplitude, and steepness of these travertine crystal growth ripples indicate that large-scale sustained aqueduct flows scaled directly with the thickness of the aqueous viscous sublayer. Resulting critical shear Reynolds numbers are comparable with those reconstructed from heat/mass transfer crystalline ripples formed in other natural and engineered environments. This includes sediment transport in rivers, lakes, and oceans, chemical precipitation and dissolution in caves, and melting and freezing in ice. Where flow depth and perimeter could be reconstructed from the distribution and stratigraphy of the travertine within the Anio Novus aqueduct, flow velocity and rate have been quantified by deriving roughness-flow relationships that are independent of water temperature. More generally, under conditions of near-constant water temperature and kinematic viscosity within the Anio Novus aqueduct channel, the travertine crystal growth ripple wavelengths increased with decreasing flow velocity, indicating that systematic changes took place in flow rate during travertine deposition. This study establishes that travertine crystal growth ripples such as those preserved in the Anio Novus provide a sensitive record of past hydraulic conditions, which can be similarly reconstructed from travertine deposited in other ancient water conveyance and storage systems around the world.

Initial irregularities have a strong and persistent effect on ripple development. Ripples of larger amplitude and wavelength have initially grown up from visible protrusions (defects) from the normally smooth channel floor. An example is the ripple in Supplementary Figure 1B between the upstream portion of the RF 9 m sample and the left-most green trough line extending from the bottom to the top of the deposit. This ripple has the greatest amplitude and asymmetry, as well as more downstream-transported SiO2 lee sands, than any others observed at Roma Vecchia, which is manifested throughout the entire vertical stratigraphic sequence of the rippled travertine sample (Supplementary Fig. 1). A ripple of seemingly longer wavelength (almost as long as the sample itself), but smaller amplitude, is visible immediately downstream.
Another example is the ripple behaviour created by the slope of the unconformity at RF 140 m, which is distributed throughout the entire stratigraphic section ( Supplementary Fig. 1A), with a particularly large and steep lee face of the ripples.
Downstream instability (i.e., migration of troughs downstream over time with upsection accumulation) at the largest ripple scale generally occurs when wavelength is increasing, and upstream when the large-scale wavelength is decreasing (Fig. 6). RF 140 m is a particularly clear example of different wavelength scales migrating in concert with each other and with the instability. In this sample, the two largest wavelength scales and the amplitude increase when the instability at the largest scale is in the downstream direction and decrease when the instability is upstream.

S3.2 Quantitative Characterization of Travertine Ripples
At times during the quantitative characterization of the up-section variation of the travertine ripples at Roma Vecchia (Fig. 6), the chosen confidence interval excluded some spectral peaks that corresponded to wavelengths clearly visible during ground truthing ("false negatives" shown in Fig. 6 as grey triangles). Amplitude tends to move in the same direction as wavelength (i.e., as amplitude increases, wavelength does too, and vice versa). The changes in wavelength are of greater magnitude than those in amplitude, however, steepness (2a/l) develops in the opposite direction to wavelength. In general, asymmetric ripples, with long, gentle convex stoss slopes and short, steep lee slopes, are more common where wavelength is systematically increasing (Fig. 6). Where it could be measured on the wall and vault at Empiglione Bridge, the mean wavelength at each point decreases with increasing height (y) above the floor ( Supplementary   Fig. S3). It is difficult to tell whether the relationship is linear, exponential etc. from 3 points.

Supplementary
From the non-boundary averaged versions of Equations 1 and 3, we would expect shear velocity and stress to be inversely proportional to wavelength. Shear stress decreases linearly with increasing y under uniform flow, however . This is the same as wavelength, rather than being inverse. This finding needs checking on more, and better preserved, rippled channels. If it is found to be generally true, then it deserves further investigation.

S4.1 Heat/Mass Transfer Crystalline Bedforms
It has been demonstrated by previous studies that convective heat/mass transfer crystalline bedforms begin at points of surface roughness (defects) that locally affect flow and hence mass transfer 6,[17][18][19] , which is consistent with our observations of the rippled-marked aqueduct travertine at Roma Vecchia as outlined above. Crystalline bedforms precipitated from chemically supersaturated water could result from several different mechanisms that are capable of causing increased precipitation at or near these defect irregularities 20,21 . The first is driven by the reduction of the depth of the flow as it passes over a raised surface defect 9 . This shallowing causes increased CO2 degassing and increased velocity due to the Bernoulli effect, which leads to thinning of the boundary layer and increased diffusion. Another mechanism is the deformation of concentration gradients in the viscous sublayer of the precipitating waters due to flow disruption caused by the surface defect. In Roman aqueducts such as the Anio Novus, great care was generally taken during Roman construction to smooth the mortar lining of the channel. As a result, the depth of flow within aqueduct channels was significantly larger than the height of any mortar surface defects and shallowing would have had a minimal effect.
Therefore, the presence of well-preserved travertine ripples in the absence of significant shallowing suggests that the viscous sublayer was the most important hydraulic factor within

S4.2 Instability in Anio Novus Travertine
One bedform characteristic which could be used to reconstruct the flow conditions that formed the bedforms is instability. For crystallization from a laminar falling film, such as stalactites (Camporeale and Ridolfi, 2012a) and icicles [30][31][32] , upstream instability is dominant, although earlier theoretical calculations predicted downstream instability in icicles 33 , and the rate (or angle) of upstream instability is inversely proportional to Reynolds Number (10 -3 < Re < 10 -1 ).
Maximum mass transfer is generally on the stoss side of bedforms 17 , probably because this is generally the location of flow reattachment. Instability in precipitation bedforms is thought to occur by adding to and extending crests upstream towards this location of maximum mass where J ( is the equivalent roughness height, which is the sum of J 6 the roughness component  Supplementary Table S2.

S4.4 Temperature Variation in the Anio Novus
Heating as the water descended more than 300 m 13  Regarding flowrate, if ' % decreases, then, by Equation 2, in a taller-than-wide rectangular channel, the cross-sectional flow area % must also decrease. When both ) 4 and % decrease, then by Equation 7, the flow rate 5 must also decrease.
Thus, observed wavelength variation at Roma Vecchia likely has a significantly greater effect on ) 4 calculated by Equation 6 than variations in 2 * (or temperature variation, as discussed above).

S4.6 Likely Time and Impact of Flow Rate Changes of Units 1 and 2 at Roma Vecchia.
In all three locations at Roma Vecchia, the travertine is deposited on an apparently pristine mortar floor surface (the t0 surface; Figs. 4, 5) and there is no evidence of rebuilding of the channel or its mortar lining after the aqueduct's initial construction 12 . Given the hardness and strongly cemented nature of the travertine, it is highly unlikely that floor travertine could have been removed during maintenance at all three Roma Vecchia without causing visible mortar damage 40 . This implies that Units 1 and 2 were deposited directly on the mortar soon after the 52 CE construction and opening of the Anio Novus.
At the best-preserved location (RF 140 m), Units 1 and 2 represent the lowest 8 cm of a 27-cm-thick floor deposit 16 . Assuming the Anio Novus operated for the longest possible time period (748 years: 52 -800 CE: see section 4: Regional Setting), this 27 cm-thick accumulation implies a minimum mean depositional rate of 0.36 mm/year. If this mean rate holds for the first 8 cm, then that 8 cm was deposited over a maximum period of 222 years. However, the actual period is likely to significantly shorter, since the minimum rate assumes neither maintenance removal of travertine, which clearly occurred as they are visible at unconformities 40 , nor periods where the aqueduct was out of use. Therefore, the mean depositional rate was likely significantly higher, especially when compared to the extremely high depositional rates observed in other natural travertine deposits 9,41 .
The effects of the reduced flow rate inferred from the qualitative relationships are clear.
Since the Anio Novus was the highest aqueduct in the Aniene valley, there would have been greater demand for its water upstream of Roma Vecchia, i.e. outside the city. More users stood to be affected by these reductions in its flowrate than by those of the other aqueducts. The wellknown increased demand for perishable and luxury goods fuelled by Rome's increasing wealth and population during the first centuries BCE and CE would have increased pressure on the Anio Novus. Much land close to the city was given over to the production of high value-density products such as flowers, luxury meats, dairy and vegetables, which could not be imported from further afield 42 . Many of these products required irrigation to realise the greatest profit and irrigation technology was installed in Rome's hinterland for this purpose. Since private water concessions were given out by the emperor at this time 10 , the wealthy and elite villa owners involved likely used their influence to gain water concessions for their fields (and elaborate gardens). These concessions would have been drawn above all from the Anio Novus, whose lower-quality water Frontinus considered more suitable for irrigation that that of the spring-fed Aqua Claudia and Aqua Marcia. The concessions would have come, however, at the expense of Rome's large urban population. Tiber River lead contamination shows that the city's water distribution system, although probably contracting, was still near its peak at the time of the two flow reductions noted above 43 .

S4.7 Climate-induced Variation in Aniene River Discharge
The Aniene River could probably have always supplied more water to the Anio Novus aqueduct than it could carry (2 m 3 /s). For instance, in 1991, baseline flows in the Aniene River near the Anio Novus' entrance were always at least 5 -8 m 3 /s, with storms elevating the flow rate above 40 m 3 /s at times 44 . The modern Aniene River has hydroelectric dams along its course, which is similar to the situation in antiquity when there were three dams at or above the intake to the Anio Novus, one of which fed the Anio Novus aqueduct 10,45 . The rainfall proxy evidence does not suggest that the Tiber basin was extremely dry during the 1 st and 2 nd century CE 46 . Hence it appears that flow in that period was comparable to the modern day and would rarely, if ever, be low enough to affect flows in the Anio Novus. This is especially likely when considering that Roman aqueduct intakes from rivers often incorporated a weir that diverted flow into the aqueduct 47 , especially in times of low flow.

S5.1 Aqueduct travertine Sampling Ssites
Galleria Egidio is located on the Tivoli loop, where little of the channel has been found, surveyed or accessed. The available evidence is largely consistent with a uniform flow assumption, however. The elevation of Galleria Egidio was not recorded during previous topographic surveys 48 , but the mean gradient over the reach containing Galleria Egidio is between 1 m/km and 1.5 m/km. This is one of the lowest gradients found along the entirety of the Anio Novus aqueduct. Given that a continual downhill gradient is required to prevent stagnation, there is little leeway for localized gradient variation along this reach or at Galleria To minimize the effects of minor damage to the channel, resulting in possible unevenness of the cross section being displayed, average eastings and northings of the interior surface of each wall (both shown in Supplementary Fig. S5) were taken and the distance between them calculated as a representation of the standard channel width.

S5.3 Aqueduct Travertine Ripple Characterization
MIPAV's livewire function produced specific coordinates for the horizons at irregular x intervals (where x is the coordinate in the downstream direction), permitting the Lomb-Scargle algorithm to be used on unevenly sampled datums. The Lomb-Scargle algorithm with WOSA was implemented using the REDFIT program 51 via the RED2CON interface (www.geo.unibremen.de/geomod/staff/mschulz) between MATLAB and REDFIT. The highest wavenumber analysed was always the Nyquist wavenumber (which is equal to half the average sampling rate and is implemented by setting the "hifac" parameter to 1). The oversampling factor was set to 4.0. There was no prescribed value for r (the average autocorrelation coefficient), except where problems with the significance testing in REDFIT required it, as suggested by the creators of REDFIT 52 . The Fourier transforms sometimes returned significant wavelengths longer than those recognizable during ground truthing.
Four amplitude proxies were tested against ground-truthed amplitude (which was the mean of 2 measurements at each horizon). These were calculated in the following way. When exported from MIPAV, the digitised ripple horizons were located entirely above the x-axis (i.e. all y-values were positive). The difference between the largest and smallest y-values (i.e. the vertical range) was divided by two to give the Vertical Range (no baseline correction) proxy (shown as blue triangles in Supplementary Fig. S6). The horizons were then baseline corrected 53 . First the horizon was centred around the x-axis by subtracting the value of a 4thorder polynomial (fitted to the horizon) where its second differential was equal to zero. A fitted sinusoidal baseline was then subtracted, to remove a ripple scale larger than the sample, if one existed. The amplitude of this sinusoidal baseline was another (shown as red triangles in Supplementary Fig. S6). The vertical range was then calculated again to give the Vertical Range (baseline corrected) proxy (shown as black triangles in Supplementary Fig. S6). Then mean of the absolute value of the horizon was calculated to determine the mean displacement proxy (shown as grey triangles in Supplementary Fig. S6).