Typhoon-induced megarips as triggers of turbidity currents offshore tropical river deltas

Tropical cyclones impose stresses on narrow and shallow continental shelves. The interaction of strong wind- and wave-induced currents with the local topography near the shore gives rise to complex flow and sediment transport patterns. Considerable uncertainty remains on the initiation mechanisms of turbidity currents, particularly in coastal oceanic settings subject to extreme weather events. Here we use state-of-the-art numerical models to investigate the implications of tropical cyclone-induced coastal circulation patterns for the generation of turbidity currents. In our simulations tropical cyclones induce megarip currents associated with shoreline curvature and rotation of incoming wave directions. These currents flush water and sediment towards submarine canyons, ultimately triggering turbidity currents into deep waters. Evidence of sediment-laden underflows, which resulted in subsea pipeline displacements, supports our hypothesis that tropical cyclone-induced megarip currents can trigger turbidity currents offshore from tropical river deltas. Intense rip currents caused by tropical cyclones can drive sediment-laden turbidity flows down submarine canyons, according to numerical simulations. Shoreline shape, bathymetry and incoming wave direction are key factors controlling this phenomenon.

T ropical cyclones regularly devastate coastal regions as their landfall induces extreme waves 1 , surges 2 and rainfalls 3 that have resulted in some of the world's most destructive coastal floods. Less well known is that cyclone-induced storms 4 and river floods 5 can trigger sediment gravity flows that move down the continental slope posing a serious threat for seafloor infrastructures 6 such as pipelines and thus a potential risk for natural hazards.
Though it has been recognised that turbidity currents can be initiated by various mechanisms, including submarine sediment failures 7 , rivers floods 8 , downwelling 9 as well as standing edge waves 10 , only recently the nearshore circulation induced by extreme events has been directly linked with turbidity currents 6,11,12 . Particularly, observed pipeline displacements have been related to turbidity currents triggered by two typhoons 13 , disclosing intense alongshore currents directed towards submarine canyons as a novel triggering mechanism for turbidity currents. Indeed, even during moderate storms, many mutually interacting environmental forces, such as buoyant river plumes, waxing and waning tides, storm waves, onshore blowing winds, give rise to complex patterns of currents on inner shelves 14 . These currents are further exacerbated by the nearshore heterogeneity 15 which is invariably determined by the presence of coastal boundaries and the shifting nature of sandy shores underpinned by sediment transport 16 . The impact severity of such processes is amplified during extreme weather events and strongly depends on the target nearshore morphology. Sediment continually evacuated from river deltas to the inner shelf are suspended over the seabed during these events by intense wave stirring and then distributed by coastal circulation patterns to be eventually transported offshore by rip-current cells that can potentially transition into turbidity currents along submarine canyon systems 10 .
Ubiquitous along energetic coastlines 17 , rip currents are funnelled flows of surface water that project out from the shoreline, across the surf-zone to finally broaden beyond the breakers. Rips are driven by alongshore variations in wave breaking as a result of the alongshore variability of the coastal morphology. As waves break they lose energy and thus transfer momentum that forces water in the direction of their propagation. This induces the mean sea level to 'set-up' onshore 18 and water to flow as alongshore currents depending on the breaking wave angle relative to the coastline 19,20 . Alongshore variations of breaking wave heights and angles are due to wave refraction owing to the morphological variability of the coastal region. The resulting patterns of higher and lower pressure gradients 21 , in combination with alongshore currents induced by obliquely-incident waves, dictate the nearshore circulation. Specifically, rips occur wherever flow convergences in the surf-zone are established 22 . During extreme wave conditions, topographically-controlled circulations give rise to large scale rips (megarips 23 ) that intensify as high-energy conditions are enforced. Such megarip circulations, usually characterised by higher velocities than typical rip currents, extend well beyond the breakers (i.e. at a distance from the shore of the order of 1 km 17 ) and are responsible for the erosion of large volumes of sediment 24 which are steered offshore through wide, deep channels. These erosional footprints have been observed adjacent to the headlands of small embayed beaches or in the centre of larger embayments 25,26 .
This study is intended to unravel the conundrum of the initiation mechanism responsible for the triggering of a recently detected turbidity current event associated with the landfall of an intense typhoon. For this purpose, the root cause of the lateral displacement occurred on a submarine pipeline has been analysed. This pipeline is located offshore Philippines along the north-eastern coast of Oriental Mindoro nearby the delta of the Magasawang-Tubig (MT) river. The MT river debouches into a very narrow shelf, only a few hundreds of meters wide, that transitions into a steeper continental slope carved by multiple submarine channels and canyons. Further, this area is prone to experience extreme weather events, such as typhoons. The pipeline lateral displacement was detected through a regular pipeline inspection in early 2016, precisely where it crosses the largest of these submarine features. Subsequent seabed assessments indicated turbidity currents as the most likely cause of this event as opposed to landslides or ground motion due to earthquakes. In particular, deep-water bedforms widespread across the submarine canyon system were interpreted as signatures of down-canyon turbidity currents. The pipeline was not damaged by the turbidity current and so the exact timing of its occurrence is unknown. However, the landfall of typhoon Melor nearby this area in December 2015 suggests it was the likely culprit of the triggering of the turbidity current that caused the pipeline deformation. Nevertheless, the physical processes driving the turbidity current remained elusive. As such, a wide array of data was collected, including high-resolution bathymetric surveys, sub-bottom profiles, sediment cores, rainfall time-series as well as a dedicated tropical cyclone hindcast. These data were then used to inform the coupling of a Reynolds-averaged Navier-Stokes model of turbidity currents with a depth-averaged shallow-water model of the tropical cyclone-induced coastal circulation. Our simulations reveal a typhoon-induced megarip current generated in-between two neighbouring river deltas and extending well beyond the shelf-break as the trigger of a strong turbidity current down the largest canyon that the pipeline crosses and where its displacement was observed.

Results
Turbidity current signatures in the MT canyon system. Since its installation in 2001 the pipeline experienced only one displacement off the MT delta ( Fig. 1) detected in 2016. Geophysical seabed surveys were carried out in 2017 and indicated turbidity currents as the most likely cause of this incident. Three main features characterise the continental slope off this delta, namely a sea-floor headwall, several submarine canyons and diffuse seafloor bedforms. These deep-water dunes found throughout the canyon system were interpreted as a signature of turbulent underflows, specifically turbidity currents. The very existence of the canyons strengthened such an interpretation. The MT canyon system develops off the homonymous delta whose two outlets debouch onto a narrow shelf that links the river directly to the continental slope, leaving little room for fluvial sediment to accumulate in shallow waters. Starting at~50 meters water depth some channels are interrupted by the headwall, which marks the continental slope at~120 meters water depth and has been interpreted as a landslide scar. Sub-bottom profiles have also shown that the submarine channels have eroded this headwall and are thus clearly younger than the associated landslide. Only one of these channelised features, the longest south-east canyon (LSEC), extends beyond 250 meters water depth and bounds the system to the south-east, where the pipeline displacement was recorded ( Figs. 1 and 2).
Bedforms display wavelengths ranging from 50 to 100 m and heights of the order of 5 m. Their orientation on the left side of the LSEC indicates that turbidity currents converge into the canyon, leading to these flows become confined and thus accelerating. The stronger the flow, the larger the bedforms, whose crests invariably grow perpendicular to the forcing currents. These bedforms show asymmetric profiles with steep stoss sides and longer lee sides with gentler gradients and thus were interpreted as antidunes. Further confirmation of such distinctive nature was provided by visual analyses of their internal architectures which revealed these bedforms were migrating up the slope counter the forcing flow direction ( Fig. 2a and b). As antidunes are generated by supercritical flows, their presence suggested the occurrence of intense turbidity currents.
On the basis of these findings, it was hypothesized that sediment-driven gravity flows developed across the upper-slope off the MT delta and reached the pipeline as violent turbidity currents due to flow convergence within the LSEC and consecutive in-channel self-acceleration. Although no indication of possible leading mechanisms for the triggering of gravity flows were found, seabed assessments point to turbidity currents as the root cause of the pipeline displacement.
Typhoon melor triggering turbidity currents off the MT delta. Some of the planet's most powerful typhoons follow generally straight westward trajectories, crossing the Philippines at their greatest intensity. Here, due to the increasing roughness of the terrestrial surface, their strength diminishes. Nevertheless, Mindoro Island experienced violent typhoons as it is only~300 km from the easternmost islands of the Philippines Archipelago.
Since 2001 there have been eleven typhoons which have passed nearby the North of Mindoro. These were all included in a dedicated hindcast which provided meteocean conditions for Tayabas Bay, with an increasing resolution towards the coast of Mindoro. Besides, 17 years of rainfall observations (2000-2016) from a weather station in Calapan (North-East Mindoro) were available. Data indicated that the typhoons' proximity was the main factor determining the severity of meteocean and hydrological forcing nearshore the MT delta. The closer the typhoons' path, the higher the resulting waves and the river discharges. Typhoons' intensity, defined as their maximum sustained wind speed, was the second most important factor. Waves exceeding four meters heights were predicted for each typhoon ranked 'category 3' or higher and passing within 50 km from the trigger area 13 . Furthermore, the landfall of such typhoons resulted in nearshore currents stronger than tidal currents and rainfalls comparable to those induced by monsoons. Three most relevant typhoons were identified, namely Durian (2006), Melor (2015) and Nock-Ten (2016), whose eyes passed within 25 km from the MT delta with an intensity ranked 'category 3' (Fig. 1).
A key differentiator among these typhoons was their approach latitude, Melor passing south of the trigger area as opposed to Durian and Nock-Ten passing north (Fig. 1b). Such difference probably affected the typhoon-induced rainfalls as the larger portion of rain usually occurs in tropical cyclones' right-front quadrants 27 . This might explain why Melor resulted in river floods larger than those associated with Durian and Nock-Ten. Furthermore, tropical cyclones' approach latitude with respect to the target coastline dictates the rotation of the incoming wave direction (Fig. 3), which is primarily due to the counter-clockwise wind patterns of tropical cyclones in the Northern Hemisphere.
This analysis revealed a variety of factors which back the capability of these three typhoons to trigger turbidity currents off the MT delta by inducing exceptional river floods and storms. Accordingly, the likelihood of the occurrence of two processes have been further assessed: indirect triggering due to flood discharges from the MT river and direct triggering due to nearshore circulations off the MT delta.
Likelihood of hyperpycnal flows from the MT river. The very first condition for a river plume to plunge into the sea as a hyperpycnal flow is to have a concentration of carried sediments larger than a certain threshold value 28 , which depends on salinity and temperature of the coastal waters. In order to assess the potential occurrence of hyperpycnal flows from the MT delta during the lifespan of the pipeline, the rating curve of the MT river, which relates water and sediment discharges at the catchment outlet, was investigated. A rainfall-runoff model was used to estimate peak flood discharges in the river basin. In addition, the outcomes of this hydrological analysis were employed as inputs for a sediment transport model providing estimates of the suspended sediment concentrations at the delta. The flood induced by typhoon Melor was characterised by a large return period of 35 years, as opposed to Durian and Nock-Ten whose floods were predicted to be lower. This notwithstanding, the sediment concentrations associated with the flood peak induced by Melor resulted in values well below the critical threshold for the inception of hyperpycnal plumes 29 .
Megarip circulation nearshore the Naujan embayment. The interplay between tropical cyclone-generated river outflows and meteocean conditions at the MT delta was addressed by a highresolution, depth-averaged application of the hydrodynamic model Delft3D with superimposed wave-current interactions obtained by its coupling with the spectral wave model SWAN. Three simulations were performed, including rising, peak and falling stages of each storms as well as the peak river floods that followed the passage of typhoons with a time-lag of~12 h. Time-series of wind and wave conditions are provided in Fig. 3 for comparative purposes. Wave roses describing the three storms are also presented with arrows designed to highlight the evolution of the incoming wave direction, which represents a key aspect in regulating coastal circulations. Snapshots of the peak currents predicted for each storm are depicted in Fig. 4. During typhoon Durian, oblique wave breaking generated a persistent southeastward alongshore current that dominated the nearshore circulation ( Fig. 4a and d). This current reached its maximum at the storm peak when the incoming wave direction was~0 ∘ N with velocities exceeding 2 ms −1 offshore the MT delta, leading to the southeastward deflection of both river outlets. Coastal and fluvial sediments were thus transported alongshore towards the Naujan embayment. Despite the wave angle of attack rotated counter-clockwise a short time-lag after the storm peak (Supplementary Movie 1), the southeastward flow persisted even during the tail of the storm with the decreased winds and waves now acting from the west inducing a reduction of the nearshore circulation, but no current reversal.
The clockwise rotation of the direction of winds and waves associated with Melor resulted in a completely different scenario ( Fig. 4b and e). During the storm waxing stage, northerly incoming waves generated a southeastward current similar to that seen dominating Durian. However, as Melor's eye approached the coast, winds and waves turned to the east-northeast to attain an overall shore-normal direction at peak conditions (~45 ∘ N). This gave rise to a northwestward alongshore current on the southeast side of the Naujan embayment that acted in opposition to the southeastward, antecedent current still prevailing on the northwest side of the embayment. The convergence of these currents initiated a megarip circulation in the centre of the embayment, where a strong jet-like flow propagated offshore towards the LSEC. This intense seaward transport of water and sediment is thus interpreted as being the leading mechanism for the triggering of the turbidity current. The waning stage of Melor was then dominated by a general northwestward flow owing to the ongoing clockwise rotation of winds and waves (Supplementary Movie 2). The remainder of the megarip circulation was deflected northward and deviated offshore the MT delta as a boundary-controlled rip current 30 .
Lastly, looking at the outcome of the 'Nock-Ten' modelled scenario (Fig. 4c), winds and waves did not originate any megarip circulation but, rather, the predicted circulation closely resembled that induced by Durian.
Typhoon-induced megarip as trigger of turbidity currents. Predicted velocity components and sediment concentrations were extracted along the 10 m water-depth contour to inform the modelling of turbidity currents, which was performed by the  computational fluid dynamics code TCsolver. In line with the timing of the pipeline shifting, numerical simulations predicted the triggering of a severe turbidity current within the LSEC only by typhoon Melor as opposed to Durian and Nock-Ten. The former induces mass gravity flows within the canyons offshore the MT delta, while the latter did not generate any intense gravity flow. Figure 5b shows an instance of the modelled turbidity current triggered by Melor (Supplementary Movie 3). Despite the inflow conditions were forced along the entire southwestern boundary, the flow was triggered with strong velocities only nearshore the centre of the Naujan embayment, as only there orientation and strength of the offshore-directed megarip provides enough momentum and suspended sediment concentration to trigger intense turbidity currents. Conversely, along the rest of the southwestern boundary, where the nearshore current was mainly alongshore-directed, powerful turbidity currents were absent. Offshore the shelf-break, the predicted turbidity current was mainly confined within the LSEC. The confined current impacted almost perpendicularly the pipeline. Contrarily, the results of the 'Durian' and 'Nock-Ten' modelled scenarios revealed the absence of any intense current pattern within the LSEC.
These outcomes matched precisely the available evidence of the observed pipeline displacement, giving confidence that the physical processes involved in the turbidity current ignition were qualitatively and quantitatively well captured. In particular, the role played by the tropical cyclone-induced coastal circulation was found to dictate where the turbidity currents were strongly locally triggered.

Discussion
Turbidity currents in submarine canyons initiate as sedimentladen underflows controlled by the suspended-sediment supply at the canyons' heads and thus the triggering of the flow is the key for their development. Despite the ongoing advancing of turbidity current modelling, considerable uncertainty still remains regarding the formulation of realistic boundary conditions for their field-scale numerical simulation, particularly in coastal oceanic settings. Down-canyon density currents have been reported to be caused by complex combinations of air-sea-land interactions, including the alongshore pile up of water and its motion owing to strong onshore blowing winds 31 , high incident surface waves 32-34 , tidal currents 35 , dense sea water cascades 9 (downwelling) 36 , internal waves 37,38 and excitation of standing infra-gravity waves resulting in canyon-oriented rip-current cells 10 . In line with these previously disclosed mechanisms, our results reveal a new initiation mechanism for turbidity currents in submarine canyons in which the passage of a cyclone generates a sequence of opposing alongshore currents in an embayed zone between two protruding river deltas (Fig. 6). The resulting transitory megarip-current cell produces substantial sediment flux directed towards an offshore canyon triggering an ignitive turbidity current that is inferred to have been responsible for the shifting of a seafloor pipeline (Fig. 5). As rips are ubiquitous along energetic, rugged coastlines, this new initiation mechanism for turbidity currents may have general application to worldwide submarine canyon systems and thus extensive implications for both the risk assessment of seafloor facilities and the current understanding of sediment transport processes across continental margins and their associated deep-sea deposits.
In light of this newly identified initiation mechanism, we reanalyse the nature of two previously detected typhoon-induced turbidity currents that impacted upon the same pipeline in the slope located just offshore the northern coast of Mindoro in the Verde Island Passage 13 . There, the Malaylay and Baco canyons are fed by the homonymous rivers that debouch straight into the canyons' heads in an area with also a very narrow shelf. Multiple canyons on the upper continental slope coalesce further downstream into a larger main canyon where the pipeline is located. Analogously to the MT canyon system, turbidity current signatures in the form of seafloor bedforms characterise the Malaylay and Baco canyons. The pipeline was displaced by a turbidity current event in 2006 and in 2016 another event severely damaged the rock berm that had been installed to protect the pipeline after 2006. The first displacement was detected after the passage of Durian where the Malaylay canyon merges into the main canyon. The second was recorded by a pipe survey in early 2017 after the passage of Nock-Ten. Similarly to the MT river, the Malaylay and Baco rivers have been shown to experience floods that are unlikely to generate hyperpycnal plumes and thus the triggering of turbidity currents near the Malaylay river mouth has been related to typhoon-induced waves and currents. In particular, Sequeiros et al. 13 predicted turbidity current triggering associated with Durian and Nock-Ten and the lack of any intense gravity flow associated with Melor in relation to the time evolution of typhoon-induced alongshore currents, which strongly depends on the typhoons' approach latitude with respect to this trigger area. Durian and Nock-Ten passing north generated peak currents directed towards the Malaylay canyon, while Melor passing south induced a current reversal at the storm peak that prevented intense sediment transport towards this canyon. The Delft3D application used to investigate the coastal circulation off the MT delta was applied to the Malaylay and Baco delta. Simulated river plumes were found to be deflected by windand wave-induced surges and currents. All the modelled scenarios revealed the occurrence of megarip circulations across the beach compartment on the west of the Malaylay and Baco delta. Figure 7 shows snapshots of the circulation patterns at peak conditions. Opposite alongshore currents formed on both sides of the western embayment and converged at its centre where megarip currents flushed the entire beach compartment. Depending on their incoming direction, winds and waves drove the dominance of one of the feeder currents and thus the deflection of the offshore-flowing rips. Specifically, Durian and Nock-Ten peaks were dominated by north-west directed forcing and an associated eastward flow that resulted in megarips deflected towards the Malaylay canyon. These are thus interpreted as an important contributing mechanism for the transport of sediment towards the canyons and the related triggering of the aforementioned turbidity currents. In contrast, the peak of Melor was dominated by a westward flow due to winds and waves coming from north-east that deflected the reminder of the megarip away from the Malaylay and Baco canyon system. Therefore, the role of the coastal circulation was found to be critical in differentiating between the most relevant typhoons also in this case. It is worth noting that, in terms of location, the MT and Malaylay-Baco are opposite systems. In the former, the main canyon is in the middle of the embayment, and the rip current must be offshore-directed to trigger a turbidity current. On the other hand, in the latter, the main canyons are at the edge of the embayment, just in front of the delta, and the rip current must be re-directed to the point of resembling an alongshore current to feed the canyons.
Generally, the triggering of turbidity currents has been related to a particular combination of wind and wave conditions appropriate to generate megarip currents directed towards the heads of submarine canyons. The clear matching with field evidence in the form of pipeline displacements associated with three different tropical cyclones in two morphologically different settings substantiates the importance of megarip circulation patterns in driving the cross-shelf transport of sediment during major storms and the related local triggering of previously unforeseen turbidity currents offshore tropical river deltas. Although in the cases presently investigated these deltas play a critical role not only in supplying year-long sediment to the shelf but also in dictating the nearby shoreline concavity, this newly identified mechanism may similarly apply to the triggering of turbidity currents offshore embayed/concave shorelines, prominent headlands, pocket beaches, even in the absence of any river outflow.

Methods
Data. In the following, the data that have been used and analysed in this study are presented.
Bathymetry and topography. Two historical bathymetric surveys, a geophysical seabed survey and information from satellite imagery were available. Sub-bottom profiler data for shallow reflection seismic profiling were collected during the 2017 geophysical survey. The survey extended from a minimum water depth of roughly 60 m down to a maximum depth of approximately 220 m along the northeastern edge of the surveyed area. Main survey SE-NW lines were acquired with a spacing of 160 m with only two perpendicular SW-NE lines. As shown in Supplementary  Fig. 1, the geophysical survey area lies between the two historical surveys and was designed to infill the gap between these two antecedent data sources. Tidal corrections were performed with a tidal model of the region. Offsets related to this correction were suitably accounted for to obtain a more precise representation of the seafloor. The so-created merged bathymetry gridded on 2 × 2 m regular mesh was first extended till the coastline using as anchoring points the EOMAP shallowwater satellite-derived data. Further, a digital terrain model (DTM, 30 m resolution) was used to improve the representation of the coastline and the MT river banks as well as to estimate the spatial extent of the MT catchment basin and then evaluate the shape of the river network based on information about the spatial distribution of contributing areas within the landscape.
Sediment cores. Several push cores and piston cores were collected during the geophysical survey. Some of the key observations from the studies of the cores are summarised below. First, short push cores (<0.6 m long) were collected within and outside the submarine canyons. Cores within these canyons contain one or two beds of mostly fine sand, ranging in thickness from 10 to 30 cm. Cores outside the canyons instead contain somewhat thinner sand layers or no sand. In particular, one single core in the more proximal part of a canyon contains a relatively stiff mud and no sand layers. The stiff mud represents older (radiocarbon ages >2000 years), more compacted layers that have been exposed by erosion. As such, the sand beds were interpreted to be deposited from turbidity currents, whereas the stiff mud containing abundant planktonic foraminifera were interpreted as the result of hemipelagic deposition at relatively low sedimentation rates. Furthermore, radiocarbon dates from planktonic foraminifera in the mud indicate sedimentation rate of the hemipelagic mud outside of the canyons of the order of 0.3 cm/year, while erosion in the canyons is evident from the much slower mud accumulation rates (about 0.02 cm/year). Finally, the recurrence frequency of turbidite sand bed emplacement in the area is estimated to be as high as 1 event every 43 years (3 events over the last 130 years), although there is some uncertainty due to erosion and poor preservation of the core tops. required accounting for the spatial variability of the rainfalls. However, this was not possible because of data limitation in Mindoro Island and thus the 3-h rainfall signal at Capalan was assumed to be spatially uniform in the catchment. Of note, the rainfall signal was cautiously increased (doubled) in correspondence of typhoon events. This was done to compensate for potential underestimation of rain measurements as a consequence of high wind speeds.
Hyperpycnal assessment. A state-of-the-art physics-based rainfall-runoff model was used to investigate the main physical drivers of high-discharge events in the MT catchment basin. The outputs of such hydrological modelling were then coupled with a sediment transport model to assess whether the MT river could generate downstream sediment concentrations high enough to trigger hyperpycnal plumes across the inner shelf.
Rainfall-runoff modelling. In the absence of direct measurements of river water and sediment discharge in Mindoro Island, the rainfall-runoff model was used to explore the main physical drivers of high-discharge events in the MT catchment.
The model combines information available on the landscape morphology, the land cover/land use and the underlying climatic forcing to represent, in a modular and spatially explicit framework, soil moisture dynamics in a spatially distributed setting 40 , incorporating information on the geomorphic complexity of the river network. The concept of "unit hydrograph" allows a proper spatial representation of hydrological processes, especially for applications in poorly gauged areas 41 , such as the MT system, where no calibration is feasible and model parameters must be specified based on the limited available information. Supplementary Fig. 2 shows the time-series of the modelled discharges spanning the entire simulation period that characterised the MT catchment basin as a highly dynamical system, with sudden and pronounced increments of its river stream in response to each considered rain events. In particular, heavy floods were predicted in conjunction with rainfalls when the antecedent moisture of the catchment was relatively high. As a result, the flooding potential of the river turned out to be quite high and not only restricted to typhoon-induced events (see the 2008 monsoon flood).
Sediment transport modelling. A sediment transport model was developed to provide estimates of the instantaneous values of suspended sediment concentration at the MT river outlets. In the absence of a river bed survey, the transport capacity of the river was estimated under the assumption of uniform flow conditions. The starting point for the formulation of the sediment transport predictor relation was a generalised form of the Engelund-Hansen sediment transport theory 42 , modified to properly model the dynamics of fine-grained sediments 43 . Besides the instantaneous water discharge, the river bankfull hydraulic geometry and the flow conductance needed to be specified. A combination of available satellite images, topographic data and a comprehensive river database 44 were used to estimate the following characteristics of the MT river upstream its delta: a river width of 25 m, a river depth of 1.8 m and a river bed slope of 2‰. In the absence of direct sediment data in the river bed and in the shallow shelf, core samples collected from the offshore areas were used to evaluate a reference value for the mean fluvial grainsize. As offshore core samples taken from within the submarine canyon system showed a mean sediment grain-size typically ranging in 100-200 μm, an overall average sand diameter of 150 μm was selected as a proxy for the representative grain-size of the river bed. Given that the quantification of the fraction of mud within fluvial sediment transport is an uncertain task, it was explored a range of sediment sizes from 20 μm to 150 μm. a range of sediment sizes from 20 μm to 150 μm was explored. Based on the above parameters, the formulation proposed by van Rijn 45 provided for the total flow conductance (comprising both the skin and form drag). Before debouching into the ocean, the MT river equally split in two branches. The instantaneous value of the depth-flux averaged suspended sediment concentration at each outlet was then computed as the ratio between the sediment load provided by the sediment transport model and the water discharge (divided in half) provided by the hydrological model.
The predicted flood induced by the landfall of typhoon Melor was characterised by the largest return period. This notwithstanding, the suspended sediment concentrations, evaluated at the catchment outlet during the peak of the flood associated with Melor, resulted in values well below the critical threshold for the inception of hyperpycnal plumes. Importantly, the critical concentration was not exceeded by the predicted sediment fluxes even when a 100%-mud scenario was considered, this suggests that the largest sediment discharges from the MT river experienced since the installation of the pipeline were likely not sufficiently high to generate hyperpycnal flows (Supplementary Fig. 3).
Coastal circulation modelling. High-resolution numerical modelling of the coastal circulation were carried out by means of the Delft3D software, an open-source state-of-the-art process-based model capable of simulating hydrodynamic and morphodynamic processes due to rivers, tides, waves, winds and atmospheric pressures in complex geomorphological settings. In particular, the model has been demonstrated to accurately reproduce along-46 and across-shore 47 currents and sediment transport 48,49 , embayment circulation patterns 50 and related headland sediment bypassing on both micro-51 and macro-tidal embayed beaches 52 .
Although validated numerical studies of energetic coastal circulations have been performed with other modelling platforms 53,54 , Delft3D was chosen as it allows the inclusion of river plumes and their contribution to the nearshore hydrodynamics.
In this study, a 2DH (two-dimensional over the horizontal directions) application of Delft3D was setup by coupling 'online' Delft3D-flow and Delft3Dwave modules. Every 15 min (coupling time-step employed) a depth-averaged flow field was supplied from the flow module to the wave module that, in turn, solved the wave action equation to transform input spectral wave conditions across the Delft3D-wave domain while accounting for several physical processes involved in the propagation of wave spectra as sinks or sources of wave energy. Information on the wave conditions were then provided to the flow module that quantify the waveinduced source term of momentum by computing the variation of the radiation stresses across the smaller Delft3D-flow domain. Here, the depth-averaged shallow water equations were then solved to calculate a new instantaneous flow field each 10 s (computational time-step employed). Water levels and mean currents were then extended to the Delft3D-wave grid within the bounds of the Delft3D-flow grid closing the computational loop.
The wave computation uses the spectral, phase-averaged, third-generation wave model SWAN 55 that is based on the discrete spectral wave action equation. The JONSWAP model 56 was assumed for the parameterized wave spectra as its variation of the wave spectral density with the frequency has been shown to be rather universal. Wave dissipation due to bottom friction was computed accordingly, assuming the suggested value for the dimensional bottom friction coefficient equal to 0.067 m 2 s −3 for fully developed wind sea in shallow waters 57 . Wave breaking dissipation was computed by the Battjes and Janssen formulation 58 with the default value for the wave height-to-depth ratio for irregular wave breaking equal to 0.55. White-capping and nonlinear wave-wave interactions were included. The wave computation was carried out in stationary mode, as preliminary computations performed in both stationary and non-stationary mode did not show important differences.
The flow computation as implemented in this study solves the depth-averaged nonlinear shallow water equations. Bottom friction due to depth-averaged currents was computed by the Chézy formulation with a constant bottom roughness coefficient equal to 65 m 1/2 s −1 . The superimposed effect of the wave-current interaction was included using the model by Fredsøe 59 . The horizontal eddy viscosity was set equal to the default value of 1 m 2 s −1 as a sensitivity analysis conducted assuming different order of magnitude of the eddy viscosity estimated from maximum and minimum expected flow velocities and employed spatial discretizations revealed its limited effect on the hydrodynamics, in line with the shallow water theory.
A series of simulations were setup for a domain encompassing Naujan embayment, including the MT delta ( Supplementary Fig. 4). Two nested computational grids were employed. The wave computation was performed by an outer curvilinear grid with a resolution increasing landward from 100 m × 50 m offshore to 50 m× 25 m at the coast. Significant wave heights, peak wave periods and peak wave directions were linearly interpolated along the offshore boundary of the Delft3D-wave grid from the ADCIRC + SWAN model produced by Oceanweather Inc. (OWI) whose resolution around Naujan embayment was~1 km × 1 km. A higher resolution grid with an increasing refinement in correspondence of the river outlets (minimum grid size 15 m × 15 m) was nested in the outer grid and then employed for the flow computation. Its offshore boundary was driven by water levels interpolated from the ADCIRC + SWAN model, while the lateral boundaries (shore-normal) were set as water level gradients (Neumann boundaries). Wind and pressure forcing were included. Pressure forcing obtained from the National Oceanic and Atmospheric Administration (NOAA) Climate Forecast System Reanalysis (CFSR) at 0.5°resolution were linearly interpolated and applied across the outer and inner grids as spatially and temporally varying inputs. Similarly, mean wind velocities and directions were extracted from the ADCIRC + SWAN hindcast nodes shown in Fig. 1 and applied temporally varying but spatially uniform across the flow computational domain. Furthermore, the total discharge of the MT river associated with the passage of each typhoon was enforced at an open boundary corresponding to a cross-section just upstream the bifurcation featuring the river's delta. The imposed time-series of the total discharge were reconstructed from the rainfall-ranoff model of the MT catchment basin, which however lacked the details necessary to provide an evaluation of the flow partitioning at the bifurcation upstream of the river delta. The geometry of the two downstream river forks was assigned on the basis of an estimate of the widths of these two river branches which was done by visual analysis of satellite images of past floods.
Preliminary simulations were repeated varying computational time windows, time step, outer and inner grid sizes proving mesh independence of the model output for the range of tested numerical discretizations. These numerical tests also indicated that a warm-up of~24 h was sufficient for initialising the ultimate simulations.
Additional knowledge on the transition of storm-induced rips into turbidity currents could be obtained by an extension of the depth-averaged Delft3D application that accounts for the vertical structure of the flow field. However, there exist laboratory measurements 60,61 , field observations 62 and numerical simulations 63 indicating that the vertical structure of rip currents is nearly depth-uniform inside the channels they erode where the waves are breaking (i.e. the surf-zone). Of note, larger offshore-directed volume fluxes induced by higher incoming waves are known to increase near-bottom velocities, leading to weaker vertical variations of rip current profiles 60 . Remarkably, in-situ measurements of storm-induced rips revealed large, near surface, vertical shear due to onshore flowing Stokes drift that can potentially result in depth-uniform seaward currents farther offshore the breakers 62 . Therefore, it seems reasonable to assume the depth-averaged shallow-water approximation when informing turbidity current modelling with inlet conditions reproducing rip currents as long as the inlet boundary lies within the surf-zone, even in the case of typhooninduced megarips. Nevertheless, the straightforward, vertical extension of the proposed Delft3D application would allow for a better description of the structure of both shallow-water rip currents and related deep-water sedimentladen underflows and the full understanding of the transition between these two intriguing phenomena. After a proper validation, such 3D model would allow to overcome existing limitations of the current state-of-the-art-modelling of turbidity currents, dynamically representing their initiation mechanisms in coastal oceanic settings.
Turbidity current modelling. Detailed three-dimensional numerical simulations were performed by means of the computational fluid dynamics code TCsolver developed at the University of South Carolina by Jasim Imran and his research group. The mathematical formulation along with an extensive validation of this code is described in several scientific publications [64][65][66] . Noteworthy, this model has been fully validated against a large number of precisely measured flows in the laboratory under various conditions (2D-and 3D-flumes with both channelconfined and unconfined flows), and against full-field scale turbidity currents using evidence of flow conditions from detailed bathymetry and sediment grain size, including canyons 67,68 .
This model employs a multiphase flow approach for modelling poorly sorted turbidity currents. It is governed by the mass and conservation equations for the incompressible bulk fluid, by the mass conservation equation of each sediment grain-size class and by the transport equation of mean turbulent kinetic energy and its dissipation rate. The solution of the governing equations requires the specification of physical boundary conditions at the computational domain edges. Thus, at the inlet, velocity, sediment concentration for each sediment class and the turbidity current thickness were specified.
Inflow conditions were applied across the southwestern, shallow-water boundary of the computational domain located roughly at the ten meters depth contour where the flow originated and propagated offshore in deeper waters, along the continental slope, through the submarine canyon system. They were specified on the basis of the results of the high-resolution numerical modelling of the coastal circulation and the sediment transport modelling. Seabed sediment was modelled using two sizes, a finer class representing the mud component and a coarser class representing the sand component of the seabed, consistently with the seabed surveys and the sediment transport model employed for the hyperpycnal assessment. The initial seabed composition was assumed to be predominantly coarser across the shelf, and well-mixed along the slope, being composed by the same two sediment classes assigned at the inflow. Sediment entrainment in the submarine canyons is crucial in determining the spatial development of the turbidity current within submarine canyon systems. Numerical simulations showed that the García and Parker equation 69 was a realistic choice. Conversely, the mud to sand ratio range (<5%) assumed in the canyons did not have an important effect on flow evolution. Finally, the turbulence closure in the numerical simulations has only a moderate effect on the results when the flow is not particularly strong, whereas it is more pronounced in stronger events.

Data availability
The data analysed and generated during this study belong to Shell Global and restrictions apply to the availability of these data, which were used under license for the current study, and so are not publicly available. Data are however available from the authors upon reasonable request and with the permission of Shell Global.

Code availability
The application used in this study to generate results on typhoon-induced coastal circulations has been carried out by means of the open-source software Delft3D developed by Deltares in close collaboration with Delft University of Technology. The Delft3D 4 suite used in this study corresponds to the latest tag version available (https:// svn.oss.deltares.nl/repos/delft3d/tags/7545). The suite is open-source upon registration on the Deltares' website and after requiring an individual graphical user interface license.