Bridge Pier Scour in Complex Environments: The Case of Chacao Channel in Chile

: Chacao channel bridge is located in a tidal channel with highly-energetic hydrodynamics conditions and signiﬁcant erodibility potential. Once ﬁnished, this 2.5 km long cable-stayed bridge will be the largest in South-America. Here we report an integrated procedure to estimate scour around two of its three towers, both located on a relatively complex but different soil matrices. A high-resolution hydrodynamic model based on the Reynolds-averaged Navier–Stokes equations (RANS), physical tests of in situ soil samples in a Rotating Erosion Testing Apparatus (RETA) and empirical formulas for scour estimation are combined to provide a reliable estimation of scour depth under a periodic tidal ebb-ﬂow regime. The relatively homogeneous soil material at the North Tower shows a high susceptibility to hydrodynamic erosion, which is estimated with SRICOS methodology. The Central Tower, in contrast, needs a combined approach based on the current state of the rock, information collected from underwater explorations and theoretical progress made about rock scour in order to reduce the uncertainty of the soils’ substrate. This study reveals that scour estimation for engineering design purposes in complex soils can be achieved with a joined vision of different disciplines and modelling tools for minimizing the uncertainty.


Introduction
Chacao Channel bridge is one of the most important civil works built in recent decades in Chile and, once in operation, will be the longest cable-stayed bridge in South-America [1][2][3].Once the concession contract expires, the bridge will be managed by the Ministry of Public Works [4].The 2.5 km long bridge will span between continental Chile and Chiloé Island (Figure 1a), sites which are currently connected by a system of ferries along a highly energetic tidal channel (Figure 1b).Travel times will be reduced in nearly 27 min once the bridge is in operation and the connection with the island will be guaranteed during storms, which currently cause operational downtime of ferries.The bridge (Figure 1c) consists of four traffic lanes and will cross a channel of up to 100 m deep.At its narrowest 2.2 km wide section, Chacao Channel exhibits an isolated pinnacle, Remolinos Rock (Figure 1d), which separates the stream in two smaller waterways [5] and where maximum speeds may reach about 6 m/s during spring tides.Remolinos Rock provides a shallow basis for its Central Tower (CT), currently under construction (Figure 1e) [6].The North Tower (NT), conversely, will be located in a shallow base near the northern coast (Figure 1f).The design life of the structure is T =100 years [1].The hydrodynamics of the nearly 26 km long Chacao Channel, which connects the Pacific Ocean with the Chilean Inland Sea (Figure 1a), is strongly dominated by differences in tide levels and phases between its two ends.At its easternmost side, the channel flows into a complex system of channels extending from 41°S to 47°S, which in turn is divided in two main water bodies with different yet correlated tidal patterns located north and south of the Guafo Mouth [7].The northern extreme of the Chilean Inland Sea is characterized by a strong amplification of the semidiurnal tidal constituents, explained by resonant effects of tides entering through Guafo Mouth [5,7].As a consequence of this amplification, tidal ranges of up to 7 m are observed in Puerto Montt during spring tides, largely exceeding the 2 m range found in the open ocean.
Simulating currents in Chacao Channel is a challenging task, mainly because the hydrodynamic system is fairly complex.Tides are the most important driving force within the channel, when contrasted to meteorological effects and river discharges [8].Several measurements performed in Chacao Channel, e.g., ICUATRO-COWI [6], and tide levels from the national tide gauge network have been used to calibrate hydrodynamic models within the Chilean Inland Sea [9][10][11][12], in which tides were prescribed as boundary conditions in relatively small domains.Winckler et al. [13] improved the modelling approach by using a much larger domain where the tidal potential of the main constituents was used in both the momentum equations and boundary conditions far from the continental shelf.The boundary conditions of the high-resolution hydrodynamic models used in the present study for the CT and NT are based on such model.
Environmental conditions in the site are particularly difficult including, aside from the high currents mentioned above, a complex soil structure and high seismicity.These conditions have prompted the introduction of innovative engineering methodologies or adaptation of existing ones in order to reduce the uncertainty in the project [2,3,6,14].Still, significant gaps remain in design conditions, particularly those related to the long-term fluid structure interaction and the estimation of scour [15,16].To this concern, design guidelines to estimate scour are readily available in non-cohesive granular soils [17][18][19][20], but they are scarce and subject to large uncertainty in cohesive, cementated or rocky materials [21][22][23][24][25][26][27].
This study deals with the scour estimation of the soil substrates supporting the NT and CT of the bridge, which geometry is shown in Figure 2.Both structures are located within a hydrodynamically complex zone, dominated by intense bidirectional tidal currents which will persistently act during the lifetime of the project.Soil conditions are particularly complex as well: While NT's substrate is formed by highly cohesive materials, CT's substrate in Remolinos Rock has a very heterogeneous lithological composition, whose origin is not well known [1].This rock shows some faces offering little to null erodibility, whereas others can be eroded by simply touching them [28].These characteristics introduce high uncertainties in the estimation of scour depth, thus forcing designers to adapt already existing methods.This study, supported in a close public-private collaboration, discusses the strengths, limitations and uncertainties of the methods used in the design, together with the complex environmental conditions judged atypical for this kind of studies.

Modelling Flow Fields
The numerical simulation of the flow, particularly in the vicinity the bridge, was calculated with OpenFoam [30], which solves the Reynolds-Averaged Navier-Stokes equations [31]: where i, j = 1, 2, 3 and u i are the components of the velocity field in rectangular coordinates denoted by x j , P the hydrodynamic pressure, g i the components of gravity (g 1 = g 2 = 0; g 3 = 9.8 ms −2 ), ν the fluid viscosity (=10 −6 m 2 s −1 ) and ν t the turbulent viscosity.This last parameter is a function of the turbulent kinetic energy of the flow k and the and turbulence specific dissipation rate ω [31] both of which are estimated from the k − ω SST turbulence closure model of an incompressible flow [30][31][32][33]: where ūj denote the time-averaged components of the velocity field u j .Flow shear stress can be calculated from the relationship: with δ ij the Kronecker delta function.The rest of the parameters are constants of the model.The turbulent kinematic viscosity is estimated as: where α 1 = 0.31, Ω = 2w ij w ij , with w ij are the components of the vorticity vector w and F 2 = tanh (λ 2 ), with: where d is the distance from the field point to the nearest wall.With these elements, the following model parameters were adopted: With this choice, the parameters σ k , σ ω , γ, β, F 1 , F 2 involved in Equations ( 3)-( 7), were easily estimated [32,33].
Modeling domains for the CT and NT were defined in rectangles of L 1 × L 2 = 1150 × 1050 m 2 and 1150 × 950 m 2 , respectively (Figure 3a).The boundaries and bathymetric regions of each domain were differentiated by mesh refinement levels (Figure 3b).Far from the towers, the seafloor and the mean sea level, the cell-size was set to 8 m in all directions.The grid was progressively refined toward the structures to properly include the bathymetric characteristics of the seafloor (Figure 3c).The grid adopted a cylindrical shape with a radius of 60 m just above the center of mass of CT and NT, with a cell-size of 2 m.Piles and its surroundings upto 1 m were modelled with cell-sizes of 0.25 m (Figure 3d).A cell-size of 2 m was also used within a band of 5 m centered vertically at the mean sea level.The choice of the grid and the cells size refinement, particularly near the walls of the numerical domain, was based on the own experience of designers fitting the critical deadlines of the project.Other CFD parameters adopted for these simulations are included in Table 1.Each simulation took around 12 to 18 h with nearly 120 parallel processors.These computations lead to simulation periods ∆T of around 370 s (≈6.2 min) up to 550 s (≈9.2 min), which were enough to observe a completely developed flow around the piles of each tower.The steady regime was reached once no significant variation of the hydraulic parameters of incoming flow (e.g., water depth, velocity) were observed.This usually occurred at the time t * < ∆T, varying between 350 s and 530 s.Once the steady regime was reached, a periodic vortex shedding phenomena was generated downstream each row of piles.The analysis of the most significant hydraulic variables of the flow (e.g., shear stress and velocity components) was undertaken considering this phenomenon into the steady regime (i.e., for t ∈ [t * , ∆T]).Flow velocities u j and flow shear stresses τ ij around the structure were obtained from the software ParaView [34].These fields were calculated for a flood-tide and an ebb-tide scenarios.Velocity and flow depth at NW and SE boundaries of our high resolution models (Figure 3b) were obtained from a large-scale model built in a Mike21, using the non-linear shallow water equations for extreme tidal ranges [13].Table 2 defines the boundary conditions at both sites.The large-scale model was calibrated using four tide gauges, one ADCP and two current meters along Chacao Channel, and then validated with three tide gauges and an ADCP at different sites.The model-to-data comparisons of the large-scale model showed coefficients of determination of R 2 = 96% for tide elevations and R 2 = 73% for depth-averaged speeds at the study site, values which are considered acceptable given the complex hydrodynamic patterns of Chacao Channel.Finally, the boundary conditions of the closure turbulence equations were imposed constant over time and estimated based on the freestream velocity, a turbulence intensity equal to 10%, and a turbulence length scale equivalent to 7% of the depth; typical values used when measurements are not available [35].

Estimation of the Erosion Rate in Earth Materials
To determine the erosion rate, we conducted experimental tests with the Rotational Erosion Test Apparatus (RETA), a device that measures the rate of erosion of different earth materials induced by the application of a shear stress at sample's surface [36][37][38].Five samples were obtained for NT and CT, however, only three were used for the former while no samples at CT were discarded.Figure 4d,e show the state of some samples obtained at both sites sites after tests.Table 2. Boundary conditions at NW and SE sides of the numerical domain of both sites obtained from a large-scale numerical model built in a Mike21, using the non-linear shallow water equations for extreme tidal ranges [29]. (1)A negative value denotes a free surface under the sea level (=0 m.a.s.l.). (2)A negative value denotes a flow direction opposite to flood tide.

Velocity
Free The RETA measures the mass loss of the sample within a given period and a range of shear stresses [38].Its configuration follows a Couette-like geometry, with a cylindrical sample in the center and an outer cylinder rotating at a constant speed ω, which mobilizes the water in between (Figure 4a).Before each test, the samples are trimmed, smoothed and perforated in the center, generating a cylindrical volume (Figure 4b).The dimensions and saturated mass of the sample are measured before immersion.An initial spin lasting up to 12 h, guarantees the elimination of loose particles, thus avoiding unusual high erosion rates.After that, the sample is removed and the water drained (Figure 4c).The dimensions and saturated mass are recorded again, from which its surface area, density, and volume are estimated.In the actual test, the sample is placed in a slightly larger unit filled with water, where an initial torque is applied during an interval that depends on the scour resistance of the sample.Erosion can be determined visually by a trained operator in tests that can last up to 72 h.
After testing, eroded sediment is placed in an oven where water is evaporated.The mass of dry sediment corresponds to the eroded material.The average shear stress τ(ω) is estimated from: where ω is the rotation speed of the cylinder, T the applied torque, R the radius and L the length of the sample.The erosion rate ṙ, measured in mm•s −1 , is estimated as: where ∆r is the average change of the sample radius for each test.This change is calculated assuming that all mass loss ∆m occurs uniformly across the outer surface of the sample within the interval ∆t.The parameter D = 2R denotes the initial diameter of the sample and ρ the dry mass density.To account for the variability of the process, between 3 to 5 tests are performed for each sample.The repeatability of these tests is ensured by operating several RETA units in parallel using the same material.

Erosion and Scour in Central Tower
Soil at CT site is a complex rock material, for which many of existing scour formulas are not directly applicable.In general, determining erodibility and local scour in rocky materials is a difficult task for hydraulic engineering [21,24,25].Erodibility of these materials has been studied to determine bedrock scour depth in dams due to the impact of a high power jet [39][40][41][42], to study the stability of spillway channels [43][44][45] and to determine the power required for excavating earth materials by means of an erosion index [46,47].These findings have been extrapolated to flows around bridges [27], since most of existing methods were thought for non-cohesive granular materials, in which the resistance to erosion of clays, rocks and cemented materials is neglected.According to Arneson et al. [48], two modes of scour and erosion in rock could arise: The first mode, referred to as quarrying and plucking, involves the removal of relatively intact rock blocks in fissured and jointed rock masses.The second mode, more gradual, is caused by bedload abrasion of the rock surface over relatively long periods of time.Here we assumed the first mode is dominant in CT while the second mode is relevant in NT.
Scour prediction at Remolinos rock (CT) was based on the Erodibility Index Method (EIM) [21,22], taking into account its variable geological composition.This method considers the different mechanism involved in rock scour (Figure 5a-c), defining an erosion threshold for the soil when subjected to an erosive current of power P A and introducing a geomechanical dimensionless index known as erodibility index K.According to Keaton et al. [27], scour in rock is essentially controlled by rock type, discontinuity spacing, abrasion resistance and weathering rate.For example, fractured and weathered rocks are more susceptible to scour, while mechanical and chemical weathering can also increase the susceptibility.Then, K quantifies the relative ability of the material to resist erosion.According to Annandale [22], K is defined as: where M s is the mass strength number, K b the block or particle size number, K d the discontinuity shear strength number, and J s the relative ground structure number, all of them dimensionless.M s and J s are estimated based on hardness and consistency of the material.In particular, M s is calculated as follows: where C r = ρ r g/27, 000, ρ r is the density of the rock, and σ ucs is the unconfined compressive strength of the sample.This strength typically falls into the range 1.7 ≤ σ ucs ≤ 212 MPa [24].The parameters K b and K d are calculated as follows: where RQD is the Rock Quality Designation, a standard dimensionless parameter in drill core logging [49], and J n is the joint set number, a dimensionless function of the number of joint sets in a rock mass (Figure 5d,e).A joint set is a family of parallel, evenly spaced joints that can be identified through mapping and analysis of their orientations, spacing, and physical properties [50].When there is not much information about the bed, the parameter RQD can theoretically be estimated as RQD = 115 − 3.3J c , where J c = 3/δ + 3 and δ is the average diameter of rock blocks that are likely to detach from their matrix, measured in m.In turn, δ can be estimated as δ = (J x J y J z ) 0.33 , an expression valid for δ > 0.1 m, where J x , J y , J z represent the mean spacing between contiguous cracks in the rock mass along (x, y, z) directions [24], all of them measured in m.RQD index usually falls in the range 5 < RQD < 100, where 5 and 100 represent rocks of very poor and very good quality, respectively.The index J n falls in the range 1 < J n < 5.The dimensionless parameters J r and J a denote the joint roughness number and joint alteration number, respectively, and can be obtained from [51].
The parameter K d typically varies in the range 1 < K d < 100 and thus, the erodibility index K varies from 1 to 100 for rocky beds.According to Yang [52], the available erosive stream power can be written as P A = γqS f , where γ = ρg is the specific weight, Q the flow rate, q = Q/W the flow rate per unit width, W the width of the channel, the length of the channel in the streamwise direction and s f the friction slope.This last parameter can be calculated as S f = dE/dx ≈ ∆E/ , where ∆E is the energy dissipation.The available stream power P A can be written as [48]: where τ = γR h S f is the shear stress acting on the rocky bed [53], R h is the hydraulic radius and V the mean flow's velocity.Here τ and V are representative scales for velocity and shear stress, so that they can be chosen according to the purpose of the study.For a wide channel, the hydraulic radius satisfies R h ≈ h, where h is the mean flow's depth for a given cross-section.Annandale [25] provided an alternative expression of Equation ( 12), as a function of the shear stress acting in the vicinity of the bedrock.According to this author, it is reasonable to assume that K is a function of the erosive power of the stream P = f (K), traducing into a threshold diagram of erosion represented by the formula: where the fitting parameters are α = 0.48 and β = 0.44 for K < 0.1, and α = 1.0 and β = 0.75 for K > 0.1 [21,22].Consequently if, for a given erosion index K, the available power P A is such that P < P A , there is erosion risk for the material, otherwise not [22].If P A /P < 1 then Z = 0, where Z is the scour depth.However, if P A /P > 1, the scour depth can be estimated from the relation [48]: where Z max is the expected maximum scour depth, ( â, b) = (8.42,0.712) are fitting parameters and B = D is the width of the pile.Annandale [25] proposes a similar formula but with another values of ( â, b), although with significant scattering, revealing a non-unique definition of the scour depth in rocky materials.

Erosion and Scour at North Tower Site
Soil at NT site is constituted by a cohesive cemented soil, highly resistant to hydrodynamic erosion but not quite enough to be considered non-erodible.For this site, an easy way to relate the erosion rate ṙ to the probable scour depth Z(t) is through the expression Z(t) = ṙt, leading to a maximum scour depth (Z s ) max = ṙT, T being bridge's lifetime.This formula assumes, however, that flow conditions are stationary, leading to unrealistic values of scour depth.A more realistic calculation should take into account the soil's nature, varying hydrodynamic conditions and the asymptotic nature of the erosive process [54,55].For such purposes, the SRICOS-EFA method (Scour Rate In COhesive Soil -Erosion Function Apparatus) proposed for bridges by Briaud [56][57][58][59] was used herein.This method has been applied to estimate scour at bridge piers [60,61], in contraction zones [62] or a combination of both, including piers of different cross-sections, with different angles of attack and pier spacing [62].This leads to scour evolution curves Z(t) for a variable velocity V(t), and can be applied to cohesive, non-cohesive, stratified and uniform granular materials, silts, clays and some soft rocks characterized by a threshold shear stress τ c [63,64].According to Briaud [63] once the flow-induced shear stress exceeds τ c , the material erodes following the evolution curve: where Z max = 0.18(Re) 0.635 is a parameter, Re = Vh/ν is the Reynolds number of the upstream flow, V = Q/Wh is the mean flow velocity, h is the mean flow depth, and t is the time such that 0 ≤ t ≤ T.

Flow Field
Figure 6 shows the velocity fields at −6 m LAT (Lowest Astronomical Tide) in the vicinity of the NT for the flood-tide and ebb-tide scenarios in the no-project condition and with the structure.The unperturbed flood-tide (Figure 6a) and ebb-tide flows (Figure 6c) show mild horizontal gradients reaching up to 3.0 m/s and 4.0 m/s, respectively.Once the structure is present, the flood-tide flow (Figure 6b) is drastically perturbed, increasing to localized speeds of up to 3.5 m/s in between piles to nearly stagnant conditions at the leeward side of each pile, as a consequence of individual wakes.For the ebb-tide flow (Figure 6d) speeds locally increase to up to 4.5 m/s, and there is a larger increase of flow speed around the group of piles.In both scenarios, there is also an overall reduction of the speed at the wake of the NT, which is heavily influenced by the pile group effect.
Figure 7 shows side views of velocity fields at vertical planes coinciding with the southern-most line of piles and between lines of piles in the NT for both scenarios in the no-project condition and with the structure.The unperturbed flood-tide flow (Figure 7a,e) shows relatively uniform speeds of nearly 2.5 m/s above a very confined bottom boundary layer.In contrast, for the ebb-tide flow (Figure 7c,g), there is a marked vertical gradient with speeds up to 4.0 m/s immediately above a thin bottom boundary layer which decrease to nearly 2.5 m/s towards the surface.This pattern is caused by the varying topography at the Tower's basis, which tend to contract the upstream streamlines, thus increasing the speed in the lower portion of the water column.Once the structure is present, both the flood-tide (Figure 7b) and ebb-tide flows (Figure 7d) are drastically reduced at the leeward side of each pile.Conversely, the speed increases to nearly 4.0 m/s in the upper portion of the water column for the flood-tide scenario (Figure 7f), and to 4.5 m/s for the ebb-tide scenario (Figure 7h), thus enhancing the speeding effect caused by the shallowing bathymetry.

Shear Stress Around Bridge Piers
Turbulent shear stress around each pile is represented by τ(t), where t is the time falling into the range 0 ≤ t ≤ ∆T, and ∆T is the simulation time introduced in Section 2.1.This stress was computed at every position (x, y, z) of the numerical domain.However, for design purposes, the maximum stress was computed close to the perimeter of each pile within the steady regime interval as follows: where t * is the instant where the incoming flow reaches a steady regime, ≈ 5 cm, r is a radial coordinate measured from the center of each pile, φ a polar angle and z a vertical coordinate directed upwards and centered in the pile (Figure 2c).The parameter τ w (t) was computed with a proprietary MATLAB code at the interception between pile's plane and the seafloor, where the stress is expected to have maximum values.Note that even in the (upstream) steady regime, τ w (t) is not strictly constant because of the vortex shedding downstream each pile.Such dynamics continuously changes the flow patterns emphasizing the need of computing the shear stress at every time of the simulation, as proposed in Equation (18). Figure 8 shows, as example, the iso-contours of shear stress τ for an arbitrary time in the steady regime, around several piles at the NT for both, flood-tide and ebb-tide scenarios.Notice that high intensity stress values concentrate on the sides of the piles with respect to flow's stream-wise direction, where significant abrasion is expected.Figure 9 shows the frequency and cumulative frequency distributions of τ w (t) for an arbitrarily chosen pile at NT, for both flood and ebb-tide scenarios.The frequency here is understood as the number of times that a value of τ w occurs into a given class interval.The highest 5% of each histogram, hereinafter noted as τ λ , was used for further analysis (green shadowed box in Figure 9a,b).Minimum, maximum and average values of this subset were extracted for each pile (τ λmin , τ λmax , τ λave ), the distribution of which is shown in Figure 9c,d for flood and ebb-tide scenarios, respectively.Each set was compared for a range of polar angles θ (Figure 2c).Regardless of the statistical parameter, the highest values of shear stress falls into the range τ λ < 250 Pa, values which were finally used to build the erosion evolution curves with the SRICOS method [63], as explained in Section 3.3.c,d) in the no-project condition and with the structure, respectively.Similar results at a vertical plane (B) between the southern-most and the immediate interior lines of piles for the flood-tide (e,f) and ebb tide-scenarios (g,h) in the no-project condition and with the structure, respectively.Arrows show the local flow direction and black lines delineate the structure's profile.

Erodibility and Scour Depth
Table 3 shows the erosion rate ṙ obtained from RETA tests for NT samples.Results show that this material is always erodible and that ṙ increases with the shear stress applied over the sample.In this case the concept of continuous and homogeneous erosion is coherent with theory described by Briaud [63].In addition, Table 3 reveals that erosion occurs only once a threshold shear stress τ c is exceeded.This threshold is sample-dependent and it is not always possible to clearly identify it.However, these measurements allowed to build an empirical relation ṙ(τ), such that ṙ = 0 and Z(t) = 0 for τ ≤ τ c and, ṙ > 0 and Z(t) > 0 otherwise.With these elements Equation ( 16) is rewritten as: where Z max is the asymptotic erosion depth of the model when t → ∞, measured in meters, and depends on the Reynolds mean flow number estimated with OpenFoam.The erosion depth (Z s ) estimated from this model for t = T is, of course, lower than Z max .In this model t is the real time, such that 0 ≤ t ≤ T where T is bridge's lifetime.As already introduced, a subset of τ w (t) was extracted from each pile to estimate the erosion rates as required by SRICOS method.Such values, denoted as τ λ , correspond to the [95%, 100%] interval from the cumulated frequency distribution shown in Figure 9a,b.By combining the empirical curves ṙ(τ) obtained from RETA tests and the statistical parameters τ λmin , τ λmax , τ λave , we interpolated a minimum, maximum and average rates, ṙmin , ṙmax and ṙave , for each pier.From these parameters, minimum, maximum and average evolution curves of erosion around each pile of the NT were computed (Figure 10a).The distribution of Z max for each pile is depicted in Figure 10b, as a function of the upstream velocity, V = max 0≤z<h V(z) , estimated for ebb and flood-tide scenarios.Figure 2 shows the reference points where such velocity was measured at both sites at a given plane z.The local velocity scale V(z) was calculated as follows: where ( ū, v, w) are the time-averaged velocity components estimated from OpenFoam in (x, y, z) directions, respectively.These averages were, once again, estimated into a steady regime.An advantage of SRICOS approach is the possibility of integrating the erosive action on each pile estimated for flood and ebb-tide scenarios, into a variable approximating flow velocity V(t), and, in turn, the shear stress around each pile τ(t), for t ∈ [0, T].We herein assumed V(t) only takes a square-wave type series with binary values for the ebb and flood-tide scenarios acting during bridge's lifetime, for each pile.The parameter Z max is described by the next fitting function, for both tide scenarios: where (m, n) = (2.10,0.63) are fitting parameters.Figure 10c shows the cumulative frequency distribution of the erosion depth for t = T, Z s (T) obtained when considering ṙ = ṙmin and ṙ = ṙmax .Around 80% of data falls into the range Z s ≤ 1.45 m, with a maximum erosion value of (Z s ) max = 1.85 m.No significant differences are observed between both curves, a consequence of the low differences between the curves in Figure 10a.19) versus mean flow velocity V at each pile ( V is the velocity obtained from OpenFoam).The fitting curve is given by Equation ( 21).(c) Cumulated frequency distribution of Z s for its minimum and maximum range of values.Some referential points were included as a guide to the eye.Figures adapted from [29].

Results for Central Tower 4.1. Flow Field
Figure 11 shows the velocity fields at −4 m LAT in the vicinity of the CT for the flood-tide and ebb-tide scenarios in the no-project condition and with the future assemble of 36 piles.The flood-tide flow is naturally perturbed by Roca Remolinos generating a recirculating pattern with very small speeds on the southern array of 18 piles (Figure 11a).In the northern array of piles, in contrast, speeds reach up to 6 m/s because of the funnelling effect of bathymetry.As in the NT, the ebb-tide in the CT shows larger speeds than the floodtide (Figure 11c), with speeds of up to 10 m/s in some piles and significant depressions of the surface's profile downstream two of the shallows shaded in black.The flood-tide flow is locally reduced by the presence of the northern array of piles (Figure 11b), but the speed reduction is moderate.For the ebb-tide flow (Figure 11d) speeds are locally reduced in the wake of individual piles, being the group effect larger in the northern array of piles.
Figure 12 shows side views of velocity fields at vertical planes coinciding with the southern-most line of piles and between lines of piles in the CT for the flood-tide and ebb tide-scenarios in the no-project condition and with the structure.For this site, mean flow's velocity significantly increase given the small flow's depths in this area.The unperturbed flood-tide flow (Figure 12a,e) shows speeds of nearly 5.0 m/s while for the ebb-tide flow (Figure 12c,g), speeds reach 9.0 m/s under a near critical flow, evidenced by localized hydraulic jumps caused by the varying topography.Once the structure is present, both the flood-tide (Figure 12b,f) and ebb-tide flows (Figure 12d,h) are reduced and spatial gradients become larger due to the presence of piles.Once again, the velocity field is perturbed by the group of piles, but to a lesser degree when compared to the NT.
It is also clear from Figure 13 that the ebb-tide scenario plays a leading role on the erosion process at the CT.Particularly, Figures 13c,d show values of shear stress up to 1500 Pa for an ebb-tide scenario, that is, at least 3 times higher than values observed for NT.In addition, the irregular bathymetry at CT site induces significant changes in the velocity field between the different sections of the group of piles, suggesting that the distribution of shear stress is highly variable throughout the perimeter of the group and around each pile.These observations emphasize the need of locally determining the velocity and shear stress around each pile when dealing with a complex bathymetry, as suggested in [48].

Erodibility and Scour Depth
Table 4 shows the erosion rate ṙ obtained from RETA tests for samples obtained at Remolinos rock.These results significantly differ with those at NT and constitute a clear evidence of the heterogeneous geological nature of Chacao Channel.The samples showed low or susceptibility to hydrodynamic erosion, even days after the RETA tests (samples RR4 and RR5).Only samples RR2 and RR3 showed non-zero erosion rates, a sign that a part of the rock could be eventually eroded by flow conditions.For this part of the rock, the SRICOS method can be applied.However, sample RR4 showed high erosion rates of ṙ > 400 mm/year, values explained by large pieces of rock suddenly removed due to increasing shear stresses, remaining the rest of the sample unperturbed (Figure 5d).For this case, the concept of continuous and homogeneous erosion described by Briaud [63] becomes misleading, forcing to explore alternative methods.The erodibility index K was then used, taking into account the information obtained from geological drilling of the site [6].
Table 5 shows the values used to compute this index, which turned to be K = 11.66.This calculation assumed an almost soft rock whose possible fractures are homogeneously distributed throughout the sample.We also assumed the sample contained cemented fissures, filled with a cohesive-like material or crushed rock elements.Underwater explorations and qualitative analysis conducted on samples provided reasonable support to these hypotheses.
Using this index, the erosion threshold stream power was estimated from Equation (15) as P = K 0.75 = 6.31 kWm −2 .This value can be compared with the available stream power P A , estimated from Equation ( 14) as P A = τV, where τ and V are representative scales for shear stress and velocity, respectively.For this study such velocity scale is given by the upstream flow velocity defined by Equation ( 20), whereas τ corresponds to τ λmin , τ λave and τ λmax .This choice is in agreement with recommendations proposed by Arneson et al. [48] who suggests the use of local variables for a better estimation of P A , thus rock scour around bridges.Using data of Figure 13c,d and the velocities calculated with OpenFoam, the hydrodynamic available stream power P A was estimated for flood-tide and ebb-tide scenarios around each pile.These estimations led to minimum, average and maximum values of P A . Figure 14a compares P A /P with the approximating velocity V of each pile, for both hydrodynamic scenarios.The ratio P A /P >1 represents an estimate of the effective probability of rock erosion, otherwise not.The black line is a fit of average data, represented by the function: where ( m, n) = (0.035, 1.78) are fitting parameters.Note that from the critical condition P A /P = 1, we obtain that V c = (1/ m) 1/ n, a critical velocity above which rock scour could occur around a pile.A simple evaluation gives that V c ≈ 6.6 m/s.To delve into the plot of Figure 14b, the cumulative frequency curve of P A /P was built for the minimum, maximum and average values of P A .Note that only 7% of data at CT site satisfy the condition P A /P > 1, implying the risk of erosion is very low for this site.In addition, this value is grounded on extremely conservative conditions when calculating the threshold erosion power and the strength of the rock.Average joint spacing (along y) 1.00 J z Average joint spacing (along z)  If we focus in the range 1.0 < P A /P < 1.37 and if some scour risk is assumed, from Equation ( 16) and the guidelines by Arneson et al. [48], it is found that Z s /D ≈ 3.43.This ratio would lead to the range 7.48 m < Z s < 8.59 m for D = 2.5 m.However, these values are extremely high and have to be carefully considered for design purposes.In this context, the scour depth for this site was estimated considering the erodible part of the rock according to RETA tests reported in Table 4.For this purpose it was essentially adopted the same calculation procedure as in NT.Estimations of local scour depth led to a more reasonable range 0.01 m < Z s < 0.63 m for this tower.

Discussion
In many engineering projects, the estimation of local scour at bridge piers can be a routine calculation, particularly when dealing with granular beds.However, this is not the case when dealing with cemented cohesive or rocky beds, like those in NT and CT.Most of existing formulas are no longer applicable for these complex soils.To reduce the uncertainty in scour estimates, a combination of 3D computational simulations, experimental tests (RETA) and the adaptation of empirical formulas provide a way to achieve safe designs when dealing with complex materials and environment conditions.In this context, a differentiated analysis is performed at both sites as follows.

About North Tower
RETA tests show that samples at NT are always erodible and the erosion is almost homogeneous (Figure 4e).According to Keaton et al. [27], cavitation is also unlikely for the dominant hydraulic conditions at this site.Consequently, an erosion curve Z(t) was built from the velocity time series V(t) around each pile.These series were assumed as the combination of pure flood-ebb tides acting during bridge's lifetime T = 100 years, from which a maximum long-term scour depth Z max was determined using SRICOS methodology.The combination of various conservative assumptions (e.g., velocity time series defined from the highest tidal ranges on a 19 years cycle and the selection of the highest 5% of shear stresses of the most critical pile) provided relatively high but conservative values of Z s for this type of materials.Another advantage of this approach is the feasibility of building a relation between τ λ and ṙ, which explicitly considers the local conditions of each site.Though the values of ṙ used by Briaud et al. [63] were estimated from piston-type devices (e.g., SERF flume [38]), the hydrodynamic setting in RETA tests is not significantly different, that is, a persistent shear stress τ acting on sample's surface.
We remark that the decreasing slope curve ṙ(τ λ ) proposed by SRICOS methodology was not observed in this study.Instead, a cloud of points with significant scattering was obtained, therefore, the erosion curve Z(t) adopted for design purposes showed an intrinsic uncertainty related to the substrate's nature and the variability of the erosion rate ṙ(τ λ ).
To reduce this uncertainty, minimum, average and maximum erosion curves were built for ṙmin , ṙave and ṙmax , respectively.By following this strategy, we determined a more reasonable maximum scour depth at each pile, falling into the range 0.67 m < Z s < 1.82 m for this site (cf. Figure 10c).Just for comparison, these results were contrasted with the classical HEC18 formulae proposed by Richardson and Davis [65], recently updated by Shan et al. [66], and proven to be applicable to a wide range of conditions and type of bridges by [67].The formulae reads as: where D is the diameter of the pile and Fr = V/ gh the Froude number of the approximating flow.The coefficients K 1 , K 2 , K 3 and K 4 consider corrections for pier nose shape (K 1 = 1.0), the angle of attack of the flow (K 2 = 1.0), bed condition (K 3 = 1.1) and bed material armoring (K 4 = 0.6).In addition, the next parameter was introduced: Figure 15 shows the distribution of ξ versus the Froude number (Figure 15a) and the relative depth h/D (Figure 15b), for flood and ebb-tide scenarios.Notice that most of points fall into the region ξ < 1, meaning that in this region HEC-18 formulae underestimates the scour depth obtained from SRICOS methodology.Just for a few conditions, the parameter ξ is greater than 1.Then, from a conservative point of view, the method proposed by Briaud et al. [58,63] minimizes uncertainties in the design process.Although Equation ( 23) provides a simple way to estimate scour depth for many materials, it must be carefully used according to recommendations proposed by Keaton et al. [27].

About Central Tower
At CT site the scouring mechanism is complex because of the three dimensionality of the flow [68] (e.g., existence of bow waves, horseshoe and wake vortexes), the characteristics of the soil and the structure itself.Then, the estimation of scour requires more elements than those for single piles.The application of the EIM method [25] leads to high erosion threshold powers (P(K)) and a very low probability of erosion at this site.Only 7% of data shows that erosion could arise for the analysed samples, but such probability is a combination of extreme conservative conditions which will be seldom reached during bridge's lifetime.Even more, if some scour is accepted to occur, Equation ( 17) leads to values of 7.48 m < Z s < 8.59 m, which are very high compared with those obtained at the more erodible material found in the NT.This contradicts evidence from underwater explorations around Remolinos rock, which showed few pieces of rock removed by currents.Then, there's a strong reason to consider this theoretical scour range unrealistic.In addition, Annandale [25] proposed another following expression for estimating Z s , which was also used in the present study.This equation led to a similar range of scour depths and, once again, relies on a large amount of samples which were nonexistent in CT.
An important aspect of this study is related with key engineering decisions adopted in the project.The idea that Remolinos rock presents a very-low risk of erosion was difficult to accept by design experts.Thus, scour depth estimations were anyway conducted on those parts of the rock showing some erodibility according to RETA tests.In this context, the procedure used in NT was adopted for this site, leading to a scour depth Z s = 0.63 m, a moderate value that agrees with the order of magnitude reported by Arneson et al. [48].
Several improvements could be achieved to reduce the uncertainty of scour estimates at this site.For example, a rock scanner would provide valuable information of the joints' structure of the rock, the faces prone to erosion and its strength.Unfortunately, budget limitations and the strict deadlines of the project inhibited the acquisition of several samples at the CT site.Additionally, the underwater exploration was restricted to nearly 10 min during slack tide (i.e., the shift from ebb to flood tides when speeds are reduced), which allowed highly qualified divers to recover only few samples from Remolinos rock.This condition imposed a serious restriction for this site, limiting our knowledge about the overall erodibility of the material.

Design Value of Scour Depth
Although the scour depth Z s was estimated for extreme conditions, an overall safety factor f s was considered for each pile.This factor was calculated as f s = f1 f2 , where fk = (1/N) ∑ N i=1 f ki is the mean value of f k , where N is the number of piles at each tower and k = 1, 2. Unfortunately, there is no physically meaningful estimation for these factors and different visions across disciplines involved in their definition persist, as pointed by Briaud et al. [69].In principle, both factors can be estimated as follows: The factor f 1i compares the maximum local scour depth (Z s ) max and the safe scour depth 1 + (Z s ) max of each pile.By contrast, f 2i compares (Z s ) max with the probable erosion obtained from average flow conditions represented by (Z s ) ave , for i = 1, ..., N.Both quantities, (Z s ) max and (Z s ) ave , are obtained from Figure 10a at t = T = 100 years.However, the application of Equation ( 25) leads to extremely high values for f 1 (e.g., in CT reach values of up to 50), which look unrealistic and not practical for design purposes.A statistical treatment was therefore applied by discarding abnormally high values, leading to average values of f1 = 1.15 and 1.35 for CT and NT, respectively.A same procedure was adopted for f 2 , obtaining in average f2 = 1.46 and 1.07 for CT and NT, respectively.With these assumptions, the equivalent safety factor was estimated as f s = 1.67 and 1.44 at both sites.The final scour depth for design Z D calculated as follows: where Z s is the overall maximum scour depth at each tower and δ is an additional depth of foundation chosen according to designers experience.The final scour depth is summarised in Table 6.

Conclusions
The pier scour depth at two towers of the Chacao Channel bridge was estimated considering the non-conventional nature of the substrates, the geometry of piles and local hydrodynamic conditions at both sites.The erodibility properties and scour depth were constrained by the limitations of existing formulas and the lack of data on the geological conditions of the substrates, which altogether forced us to propose an integrated approach to properly account for soil's nature.
At NT site, a simple adaptation of SRICOS method lead to reasonable values of the scour depth, that is, Z s ≈ 1.82 m, which turned to 2.82 m considering safety criteria.This site is characterized by a highly cohesive bed with high resistance to hydrodynamic erosion, but not quite high to be non-erodible.In contrast, the nature of the soil at CT site suggested the use of the EIM erosion index, showing large uncertainties when applied to Remolinos rock.Although erosion risk was estimated to be very low for this substrate, a conservative approach considered a non-zero safety scour depth of Z s ≈ 0.63 m for this site, which turned 1.63 m considering safety criteria.After several discussions with the Ministry of Public Works, we proposed to complement these calculations with a monitoring plan in order to control possible bedrock removals during the bridge's lifetime.This case study reveals that a combination of numerical modeling, experimental tests and long-term scour monitoring plan could provide a reasonable strategy to reduce the uncertainty when estimating the scour depth around bridges in complex and heterogeneous substrates with scarce information about geology.Additionally, a more detailed characterisation of the CT site, and in less degree of the NT site, would have further contributed for a better estimation of scour depth and the susceptibility to erosion.
To conclude, literature still lacks of a general and confident formula for rock scour that considers all the mechanical properties of a rock, the involved erosion mechanisms, the local hydrodynamic patterns, and pile group effects, among other variables.This restriction was a key aspect of this study and is a challenge that deserves to be afforded in future hydraulic engineering studies.

Figure 1 .
Figure 1.(a) Map showing Chacao Channel, Chiloé Island and the Chilean Inland Sea.(b) Aerial view of the narrowest section of Chacao Channel during ebb tide.Remolinos Rock generates a turbulent wake (courtesy of Horacio Parrague).(c) A frontal view of the bridge showing the position of Central and North towers.(d) A close view of Remolinos rock.(e) A current view of bridge piers being mounted over Remolinos rock site.(f) A projected view of the bridge.

Figure 2 .
Figure 2. Top view geometry and dimensions, in m, of (a) North Tower (NT) and (b) Central Tower (CT) of Chacao bridge (adapted from [29]).NT is made up of two groups of piles, each having 9 cylindrical piles with diameters of D = 2.5 m.CT is a larger structure made up of four groups of piles, each formed by 9 piles of the same diameter.NT and CT occupy areas of 50 × 21 m 2 and 59 × 51 m 2 , respectively.The distance between neighbouring piles is s = 7.5 m, measured from their respective vertical axes.Dashed lines show profiles where speed and shear stress are estimated.Points in red (green): locations where the velocity V(t) was estimated for flood (ebb) tide scenario.Points in black: locations for estimating V(t) in both scenarios.(c) Definition of a reference frame for the estimation of the approximating flow velocities V(t) and the shear stress τ w acting around the piles.(d) Front view of CT site.

Figure 3 .
Figure 3. (a) Definition of the modelling domain at the CT.(b) Definition of boundaries of the numerical domain.(c,d) show a views of the refined mesh used for numerical modelling.Figures adapted from [29].

Figure 4 .
Figure 4. (a) A front view of RETA apparatus [38].(b) A view of some of the samples and the involved dimensions.Taps at both extremes were also indicated (arrows in red).(c,d)Some samples of Remolinos rock after an erosion test.(e) A view of a sample from NT after an erosion test (figures extracted from[28]).

Figure 5 .
Figure 5. Consecutive steps of rock scour for a turbulent flow: (a) jacking, (b) dislodgement and (c) displacement.Schematic presentation of joint sets for a rock specimen, (d) a rock with 1 joint plane and (e) 3 joint planes, pointed in red circles.Figures adapted from[22,24].

Figure 6 .
Figure 6.Velocity fields at −6 m LAT in the North Tower site for the flood-tide scenario in the (a) no-project condition and (b), with the future assemble of 18 piles.Similar results for the ebb-tide scenario in the (c) no-project condition and (d) with the piles' structure.Arrows show the local flow direction, black circles the location of each pile and black regions represent seafloor located above −6 m LAT.

Figure 7 .
Figure 7. Side views of velocity fields at a vertical plane (A) coinciding with the southern-most line of piles in the North Tower site for the flood-tide (a,b) and ebb tide-scenarios (c,d) in the no-project condition and with the structure, respectively.Similar results at a vertical plane (B) between the southern-most and the immediate interior lines of piles for the flood-tide (e,f) and ebb tide-scenarios (g,h) in the no-project condition and with the structure, respectively.Arrows show the local flow direction and black lines delineate the structure's profile.

Figure 8 .
Figure 8. Iso-contours plots of shear stress τ (measured in Pa) around the piles of North Tower for (a) flood-tide and (b) ebb-tide scenarios.High intensity shear stress is shown in yellow, and low intensity in blue.For simplicity, only results for the front and rear rows of the tower were shown.Figures adapted from [29].

Figure 9 .
Figure 9. Frequency distribution of τ w (t) for pile A1, obtained for t > t * (upper left pier in Figure 2a) of the NT for (a) the flood-tide (b) and the ebb-tide scenarios.The cumulative frequency curve (in red) is also included.Distribution of τ λ versus the central angle θ (Figure 2c) for (c) the flood-tide (d) and the ebb-tide scenarios.Minimum, maximum and the average of this parameter are shown for reference.Figures adapted from [29].

Figure 10 .
Figure 10.(a) Erosion curves Z(t) obtained by considering different erosion rates, ṙmin , ṙave , ṙmax , for one of the piles of NT (pile A6).The erosion depth Z s was indicated for reference.(b) Parameter Z max involved in Equation (19) versus mean flow velocity V at each pile ( V is the velocity obtained from OpenFoam).The fitting curve is given by Equation (21).(c) Cumulated frequency distribution of Z s for its minimum and maximum range of values.Some referential points were included as a guide to the eye.Figures adapted from[29].

Figure 11 .
Figure 11.Velocity fields at −4 m LAT in the Central Tower site for (a) the flood-tide scenario in the no-project condition and (b), with the future assemble of 36 piles.Similar results for (c) the ebb-tide scenario in the no-project condition and (d) with the piles' structure.Arrows show the local flow direction, black circles the location of each pile, black regions represent seafloor located above −4 m LAT and white regions areas where the water level is below −4 m LAT.Figures adapted from [29].

Figure 12 .
Figure 12.Side views of velocity fields at a vertical plane (A) coinciding with the southern-most line of piles in the CT site for the flood-tide (a,b) and ebb tide-scenarios (c,d) in the no-project condition and with the structure, respectively.Similar results at a vertical plane (B) between the southern-most and the immediate interior lines of piles for the flood-tide (e,f) and ebb tide-scenarios (g,h) in the no-project condition and with the structure, respectively.Arrows show the local flow direction and black lines delineate the structure's profile.Figures adapted from[29].

Figure 13 .
Figure 13.Frequency distribution of τ w (t) for pile A1, for t ∈ [t * , ∆T] for pile A1 (upper left pier in Figure 2b) of the CT site for (a) the flood-tide scenario (b) and the ebb-tide scenario.The cumulative frequency curve (in red) is also included.Distribution of τ λ versus the central angle θ (Figure 2c) for for (c) the flood-tide scenario (d) and the ebb-tide scenario.Minimum, maximum and the average of this parameter are shown for reference.Figures adapted from [29].

Figure 14 .
Figure 14.Erodibility at CT site.(a) Distribution of P A /P versus V obtained for minimum, maximum and average shear stress (τ λmin , τ λave , τ λmax ).The fit corresponds to Equation (21) for averaged data.A critical "scour" velocity V c =6.6 m/s was included.(b) Cumulated frequency distribution curves obtained for the ratio P A /P for the same conditions.The erosion risk region (P A /P ≥ 1) was shown for reference (box in light green).

Figure 15 .
Figure 15.Distribution of parameter ξ for flood-tide and ebb-tides acting at North Tower site, versus (a) the mean Froude (Fr) number and (b), the relative depth h/D, where h is the mean flow depth around each pier and D the diameter.The horizontal dashed line corresponds to ξ =1.

Table 1 .
Simulation parameters adopted in OpenFOAM for both towers.

Table 3 .
Rate erosion ṙ obtained from RETA tests for samples extracted from North Tower (samples: TN1, TN2 and TN3).

Table 5 .
Erodibility index K calculated for Remolinos rock.The parameters included in the table were already described in the text.

Table 6 .
(26)e of scour depth Z s and final scour depth value Z D adopted at NT and CT. ( * ) Estimated considering Equations (25) and(26).