The multi-scale modelling framework developed for this study utilizes a nested construct of the Regional Ocean Modelling System (ROMS). ROMS is a free-surface, terrain-following primitive equation model widely used to study ocean dynamics in diverse environments, from the sea surface to deep ocean ridges (e.g., Shchepetkin and McWilliams, 2005; Warner et al., 2010; Vic et al., 2018). ROMS implements efficient nesting schemes to transfer data between different grids in a multi-resolution simulation. The nesting can either be one-way, where a coarse grid provides the lateral boundary conditions for the embedded finer grid(s), or two-way, where the fine grid(s) also provide(s) feedback to update the coarse-grid solutions. For this study, we adopt two-way nesting to ensure a smooth transition of model fields and energy cascades between nested grids.

The multi-resolution domain of our modelling framework ( Figure 2 ) consists of two layers of nested grids having spatial resolutions that increase from 360 m (fine) in the outer layer ( Figure 2B ) to 72 m (ultra-fine) in the inner layer ( Figure 2C ). Both grids are rotated to align with the ridge axis (20° from true north). The computational cells in each grid are uniform in the horizontal direction and each grid consists of 32 vertically stretched layers designed to enhance the resolution near the ocean surface and the seafloor, yielding a vertical resolution of approximately 15 m at the crest of the Endeavour Segment. The dynamics within the model domain are driven at the lateral boundaries of the fine grid and at the surface of both grids. The model relies upon lateral boundary conditions that are derived from the daily averaged output from a stand-alone pilot simulation, which covers a larger domain ( Figure 2A ) and has a lower horizontal resolution of approximately 1800 m (coarse). The lateral boundary forcing for the pilot simulation is constructed from the daily averaged product of the Copernicus Marine Environment Monitoring Service (CMEMS) global ocean eddy-resolving (1/12°) reanalysis (Fernandez and Lellouche, 2021). The surface forcing used in the pilot simulation is derived from the 3-hourly winds, atmospheric pressure, and fluxes from the fifth generation European Centre for Medium-Range Weather Forecasts atmospheric reanalysis of the global climate (ERA5, Copernicus Climate Change Service (C3S), 2017).

The pilot simulation serves as a transition between the global reanalysis and our regional simulations, maintaining a grid refinement factor (i.e., the ratio of spatial resolutions between two consecutive grids) of approximately five. To account for tidal forcing, which is absent in the global reanalysis and the pilot simulation, we have added tidal elevation and barotropic currents acquired from the TPXO global tidal solution (Egbert and Erofeeva, 2002) to the lateral boundary forcing for the fine domain. We derive the surface forcings for both fine and ultra-fine domains from the same ERA5 3-hourly dataset used in the pilot simulation. The initial model fields used for the pilot simulation are constructed from the CMEMS global reanalysis data. The initial fields for the nested simulations are interpolated from the pilot simulation output onto the fine and ultra-fine grids. The computational time increments are 150 sec, 30 sec, and 6 sec on the coarse (pilot), fine, and ultra-fine grids, respectively.

The seafloor topography in our model is constructed from the Global Multi-Resolution Topography (GMRT) Synthesis of the Marine Geoscience Data System (Ryan et al., 2009). The original GMRT bathymetry has a spatial resolution of 250 m. Before interpolating on our model grids, we first apply a spatial median filter to smooth out features with length scales smaller than twice the size of a computational cell (i.e., the smallest resolvable wavelength). To further enhance numerical stability, we then apply the smoothing techniques described in Sikirić et al. (2009) to reduce the maximum stiffness ratio (Shchepetkin and McWilliams, 2005).

The model incorporates the buoyancy generated by hydrothermal venting as bottom heat flux originating from the five major vent fields located within the axial valley. Specifically, each of the vent fields is represented as a 216 m by 216 m (3 by 3 cells) heat source in the ultra-fine domain. The corresponding heat flux is obtained from ship-based hydrographic surveys conducted in September 2004 (Kellogg and McDuff, 2010; Kellogg, 2011) and are prescribed as follows: 135 MW for Mothra, 306 MW for MEF, 387 MW for HR, 67.5 MW for Salty Dawg, and 4.5 MW for Sasquatch. In the model, the sources are treated as having no explicit volume flux. Consequently, the hydrothermal discharge is represented as thermally driven plumes. This simplified treatment of the plume sources is reasonable, considering that hydrothermal plumes are typically characterized as ‘lazy’ plumes (Hunt and Kaye, 2005). In such cases, the volume flux of the source has negligible influence on the plume’s behaviour and dynamics downstream. The vertical eddy viscosity and diffusivity in the model are determined using the Mellor-Yamada 2.5 level turbulent closure scheme (Mellor and Yamada, 1982), which has been shown in previous studies to enable simulation of vertical motions of plume fluid and the associated entrainment of ambient seawater into the plume, despite the use of the hydrostatic approximation (Thomson et al., 2005; Thomson et al., 2009).

In this study, we have conducted simulations covering two calendar-year periods: 2016 and 2021. We selected these years primarily due to the availability of concurrent observational data, which can be used for model validation. Those observations include the time series of flow velocity measured by current meters on cabled moorings deployed within the axial valley of Endeavour as part of the NEPTUNE seafloor observatory infrastructure by Ocean Networks Canada (Kelley et al., 2012) ( Figure 1B ). Additionally, comparing the model outputs between these two one-year periods enables us to investigate the interannual variability in ocean circulation and plume dispersion near the ridge segment.

To investigate how ocean tides and buoyancy from seafloor venting impact plume dispersion and ocean circulation near the ridge segment, we have conducted three runs for the 2016 simulation. The primary simulation output for model validation and interpretation was obtained from the first run, referred to as ‘baseline’ hereafter, which utilized the full set of forcings described in Sections 3.1 and 3.2. The remaining two runs employed a similar model configuration, but with a reduced set of forcings. Specifically, the second run, referred to as ‘no tide’ hereafter, included buoyancy from venting but not ocean tides, while the third run, referred to as ‘no vent’, included ocean tides but not buoyancy from venting. Additionally, we have conducted one run of the 2021 simulation using the full set of forcings (2021 baseline).

Our evaluation of the accuracy of simulated flow velocities utilizes concurrent, collocated comparisons with measurements obtained from the three current-meter moorings (RCM-NE, RCM-SE, RCM-SW) deployed within the axial valley of the Endeavour Segment ( Figure 1B ). Specifically, we utilized the current-meter time series from all three moorings to assess model skill during the 2016 simulation period. On the other hand, for the 2021 simulation period, we only used data from the two southern moorings (RCM-SW and RCM-SE) since no data was available from RCM-NE for most of the simulation period. We quantify the statistical agreement between the model and the observations using the Willmott Skill Score ( W S S ), which measures the normalized mean squared error in the forms of

and

for the zonal ( U ) and meridional ( V ) flow velocities, respectively. In equations (1) and (2), the subscript ‘m’ denotes model output, while ‘d’ indicates observed data. The angled brackets (e.g., ( e . g . ,   〈 | U m − U d | 2 〉 ) . ) represent time averaging of the quantity enclosed within the brackets. The score is a measure of the similarity between the model output and the observed data. The score ranges from 0 to 1, with higher values indicating better model performance (Willmott, 1981). The W S S listed in Table 1 suggest that the model performs better at simulating the meridional velocity ( V ) than the zonal velocity ( U ).

A one-to-one comparison between the simulated flow velocities and measurements for all three moorings demonstrates an encouraging level of agreement. Considering flow at RCM-SW in 2021 for illustration ( Figure 3 ), the observed currents predominantly flow in directions that are roughly aligned with the orientation of the ridge-axis at RCM-SW. This alignment is likely caused by the steering of flow within the topographic confinement of the axial rift valley. At a depth of 2162 m, just 1 m above the seafloor ( Figure 3G ), the direction of the observed currents oscillated between south-southwest (SSW) and north-northeast (NNE). This included an episode of intensified southwestward flow around Oct 27th, 2021. During that episode, the magnitude of the 2-day averaged flow velocity reached a maximum of 7.4 cm/s, compared with a median magnitude of 0.8 cm/s over the measurement period. At a shallower level, at 2117 m depth ( Figure 3E ), the direction of the observed currents was primarily SSW and only occasionally switched to north-northwest (NNW). Higher in the water column, at 2042 m and 1967 m depths which are shallower than the confining depths of the axial rift valley ( Figure 3A, C ), the westward component of the observed flow velocity was more prominent and the predominant flow direction was rotated toward the west-southwest (WSW), likely indicating a weakening of topographic steering effects at these higher altitudes. Encouragingly, the magnitude and direction of the simulated velocities ( Figure 3B, D, F, H ) are largely consistent with the observations at all depths at RCM-SW. There are two notable exceptions to this general case, 1) at 2042 m depth, the change of flow direction from SW to NW that occurred around Oct 16, 2021 is much less pronounced in the model velocity ( Figure 3C vs 3D ); 2) at 2162 m depth, the episode of intensified currents around Oct 27th, 2021 is absent in the simulated velocity ( Figure 3G vs 3H ).

At mooring RCM-SE, ( Supplementary Figure 1 ) the observed currents at 2093 m and 2167 m depths are not as well aligned with the ridge segment as they are at similar depths at RCM-SW, despite the proximity of the two moorings ( Figure 1B , Supplementary Figure 1 ). In particular, the observed currents at 2167 m depth predominantly flow in the NNW direction ( Supplementary Figure 1C ), differing from the currents observed at approximately the same depth at RCM-SW that oscillated between SSW and NNE ( Figure 3G ). In comparison, at 2093 m depth, the observed flow direction at RCM-SE is primarily SSW, occasionally shifting to NNW ( Supplementary Figure 1A ). This vertically sheared flow structure between 2167 m and 2093 m depths at RCM-SE is not observed at RCM-SW. However, the simulated flow velocity for RCM-SE does faithfully reproduce a similar vertical sheared structure, with currents predominantly flowing in the NNW and SSW directions at 2167 m and 2093 m depths, respectively. One notable discrepancy between measured and simulated flows at the RCM-SE locale is that the simulated SSW flow at 2093 m depth is more uniformly consistent than the observations ( Supplementary Figure 1a vs 1b ). The differences in observed currents between RCM-SW and RCM-SE reveal evidence for a complex flow pattern that exhibits both horizontal and vertical shear along the ridge segment, especially within the axial valley. Overall, these features are reproduced, at least qualitatively, by the model.

For the 2016 simulation period, the flow velocities obtained from the simulation generally match the measurements taken at the RCM-NE ( Supplementary Figure 2 ), RCM-SE ( Supplementary Figure 3 ), and RCM-SW ( Supplementary Figure 4 ) moorings. However, one notable distinction is evident for the current meter situated at a depth of 1956 m on RCM-NE, where the observed currents ( Supplementary Figure 2A ) predominantly flow southward throughout most of the simulation period. In comparison, the simulated flow velocity ( Supplementary Figure 2B ) primarily exhibits a westward direction, transitioning from NNW to SSW during the initial few months of the simulation period, and subsequently shifting towards a southward direction for the remainder of the simulation period. At the RCM-SE mooring for the same period, both the observed and simulated currents at 2027m predominantly flow southward, but the observed flow velocity exhibits a larger magnitude compared to the simulated flow velocity, which also has a stronger westward component ( Supplementary Figure 3A,B ).

To assess the model’s proficiency in reproducing the oscillations in flow velocity at various frequencies, we have computed the power spectra for both observed and simulated velocity time series at RCM-SW ( Figure 4 ) for the 2021 simulation period and at RMC-NE ( Supplementary Figure 5 ) for the 2016 simulation period. The spectrum of the measured flow velocity exhibits prominent peaks at semi-diurnal (S, ~2 cycles/day) and diurnal (D, ~1 cycle/day) tidal frequencies, as well as near the local inertial frequency (f, ~1.48 cycles/day) at the mooring sites. Moreover, the flow velocity measured at all but the lowest current meter (5 m above bottom) on RCM-SW in 2021 ( Figure 4 ) displays significant oscillations in the observational data within the ‘weather’ band previously noted, with a peak near a period of 4 days in the clockwise component of the spectrum (Cannon and Thomson, 1996; Thomson et al., 2003).

In comparison, the spectra derived from simulated flow velocity time series demonstrate overall agreement with the observations, with a few noticeable discrepancies. Firstly, the simulated flow velocity lacks the 4-day ‘weather’ band oscillations present in the observations. Secondly, the model appears to overestimate the near-inertial oscillations (f) at the current meters located close to the seafloor, specifically at altitudes of 50 m and 5 m ( Figures 4C, D and Supplementary Figure 5C ). As indicated by Thomson et al. (2003), the near-inertial oscillations experience significant attenuation within the rift valley of the Endeavour Segment. This previous finding is consistent with our current observations, which demonstrate a decrease in the magnitude of the near-inertial spectral peak as depth increases ( Figure 4 and Supplementary Figure 5 ). On the other hand, the simulated flow velocity shows less prominent attenuation of near-inertial oscillations within the axial valley. This discrepancy can likely be attributed to the smoothing of the bathymetry used in the simulations, which is done to enhance computational stability and reduce numerical errors. The bathymetric smoothing process within the model potentially suppresses the ruggedness of the axial valley, thus reducing its damping effect on near-inertial oscillations. The reason behind the absence of 4-day oscillations in simulated flow velocity remains unclear. Cannon and Thomson (1996) attributed the occurrence of these oscillations observed along the Juan de Fuca Ridge to local storms. Therefore, the lack of such oscillations in our simulations may indicate an inadequate representation of local storm events in the surface forcing used by the model (Section 3.1).

To evaluate the model skill in simulating plume dispersion, we compare the characteristics of the non-buoyant cap of the plume generated by the model with observations from previous studies (Baker and Massoth, 1987; Thomson et al., 1992; Kellogg and McDuff, 2010). Specifically, we focus on three key properties of the plume’s non-buoyant cap: height, core potential temperature anomaly, and core density. Our analysis shows that in the two instances demonstrated in Figures 5 , 6 , the buoyant plumes rising from the five major vent fields coalesce and reach neutral buoyancy at a height of 100 – 200 m above the floor of the axial valley. In the interior of the ultra-fine domain, within tens of kilometers of the ridge segment ( Figures 5B , 6B ), the non-buoyant plume cap disperses as a cohesive band in the NW and SW direction on Oct 2, 2016 and Aug 16, 2021, respectively. In comparison, as the plume cap spreads farther away from the ridge segment, it becomes more diffuse and displays meandering patterns ( Figures 5A , 6A ). In both cases, the core of the plume’s non-buoyant cap is situated at depths between 2000 m and 2050 m, between the 27.7 and 27.71 k g / m 3 isopycnic surfaces ( Figures 5C , 6C ). We calculate the isohaline potential temperature anomaly ( Δ θ S ) using the method described in Kellogg and McDuff (2010) as the deviation of absolute potential temperature ( θ ) from a reference for a given salinity ( S ). To obtain the reference temperature, we fit a second-order polynomial to the model θ / S data between 1700 m and 2300 m depths at three background stations ( Figures 5B , 6B ). These stations are located east of the MEF at (47.95°N, 128.94°W), (47.95°N, 128.85°W), and (47.95°N, 128.81°W), respectively. The resulting Δ θ S at the core of the non-buoyant plume exhibited a maximum, at a depth of 2050 m, of approximately 0.08°C and 0.10°C on Oct 2, 2016 and Aug 16, 2021, respectively ( Figures 5D , 6D ). In comparison, a hydrographic survey conducted along the axial valley in June 2004 showed hydrothermally influenced water with Δ θ S ≈ 0.10°C at the same depth (Kellog and McDuff, 2010).

The analysis of the 2016 baseline solution ( Figure 7A ) shows that the flow field, when averaged over a five-month period from July to November 2016, exhibits an anti-cyclonic circulation (clockwise in the northern hemisphere) that is centered around the high point on the western side of the ridge segment. Specifically, there is a poleward along-ridge flow over the western flank of the ridge, while an equatorward along-ridge flow is present over the eastern flank ( Figure 8A ). At the core of these ridge flank jets, the velocity reaches a maximum of approximately 0.04 m/s at a depth about 2050 m. A similar anti-cyclonic circulation has been observed around the North Cleft Segment of the southern JDFR and was attributed to two potential generation mechanisms: hydrothermal forcing and tidal rectification (Cannon and Pashinski, 1997). More specially, the entrainment of ambient seawater into buoyant hydrothermal discharge can drive an anti-cyclonic circulation at the spreading level of the non-buoyant plume cap (Speer, 1989). Additionally, tidal rectification has been shown to generate anti-cyclonic residual currents around pronounced topographic features such as Georges Bank (Chen and Beardsley, 1995) and Axial Seamount (Xu and Lavelle, 2017).

A comparison between the baseline solution and the results of the no-tide run indicates that the anti-cyclonically sheared mean flow across the ridge segment substantially weakens after the removal of the tidal forcing ( Figure 7A vs 7B , Figure 8A vs 8B ). Furthermore, a comparison between the baseline solution and the results of the no-vent run also shows a reduction in the strength of the anti-cyclonic flank jets along the ridge-axis, albeit to a lesser extent ( Figure 7A vs 7C , Figure 8A vs 8C ). Taken together, these comparisons suggest that tidal forcing is likely the primary driver of the anti-cyclonic circulation around the ridge segment, with hydrothermal venting playing a secondary role. This finding is consistent with Cannon and Pashinski (1997)’s conclusion that tidal forcing is the dominant generation mechanism for the anti-cyclonic circulation around the North Cleft Segment.

The buoyant hydrothermal plumes that rise from the axial valley draw in surrounding water and drive an inflow, which is most pronounced around the HR vent field with which exhibit the highest heat outputs ( Figures 9A, C ). This inflow alters the near-bottom mean flow within the axial valley, which is primarily directed along the ridge axis in the SSW direction in the absence of venting ( Figure 9B ). In comparison, the venting-induced local circulation causes the mean flow to become cyclonically sheared within the valley, with an equatorward current on the western side and poleward current on the eastern side ( Figures 8A , 9A, B ). This cyclonically sheared flow pattern is also evident in current-meter measurements from the RCM-SW and RCM-SE moorings in 2016. At approximately 50 m above the seafloor, the current at RCM-SW on the western side of the valley primarily flowed in an equatorward direction along the ridge axis, while the current measured at RCM-SE on the eastern side of the valley was primarily poleward along the ridge axis ( Supplementary Figure 3E vs Supplementary Figure 4C ). The model currents sampled at the same locations as these two moorings show reasonable agreement with the measurements ( Supplementary Figures 3, 4 ). Previously, current velocity measurements within the axial valley had been interpreted in such a way as to suggest that venting drives a steady near-bottom flow into the valley from its northern and southern ends (Thomson et al., 2003). However, the long-term current-meter measurements and concurrent simulation results presented here reveal a more complex venting-induced flow structure that is cyclonically sheared across the valley. Our findings are consistent with those from a previous numerical study, which showed a similar cyclonically sheared near-bottom flow within the axial valley, using realistic ridge topography, a near-steady cross-ridge flow and hydrothermal venting as forcing (Thomson et al., 2009). Based on previous laboratory and numerical studies (e.g., Speer, 1989; Helfrich and Battisti, 1991; Speer and Marshall, 1995; Fabregat Tomàs et al., 2016), it has been observed that the influx of surrounding water into a turbulent buoyant plume drives a local circulation that is anti-cyclonic near the spreading level of the non-buoyant plume cap and cyclonic around the plume stem below. In our simulations that account for realistic ridge topography, the deeper forced circulation, driven by the convergent flow towards the plume within the confinement of the axial valley, leads to the formation of the cyclonically sheared mean flow within the valley ( Figures 8A , 9A, C ). Meanwhile, the divergent flow out of the non-buoyant plume cap contributes to the formation of the anti-cyclonically sheared mean flow at depths overlying the axial valley floor ( Figure 7A ).

The 2016 baseline solution shows that, within the ultra-fine domain, the non-buoyant plume cap disperses primarily towards the north within a few tens of kilometres from the Endeavour Segment ( Figures 10A, C, E ). As the plume travels further away ( Figures 11A, C, E ), it splits into two branches near the Endeavour Seamount (ESM) and the southern end of the West Valley Segment (WVS) of the JDFR, which lies immediately to the north of the Endeavour Segment. One of those branches continues its northward movement, to the east of Endeavour Segment, while the other turns westward flowing south of the Heck Seamount chain ( Figures 11A, C, E ). In comparison, the no-tide run shows a similar northward dispersion within the vicinity of the Endeavour Segment ( Figures 10B, D, F ). Farther afield, however, instead of bifurcating into northward and westward moving branches as in the baseline solution, the plume veers to the northeast past Endeavour Seamount ( Figures 11B, D, F ). This disparity in plume dispersion patterns can be attributed to the existence of tidally rectified anti-cyclonic circulation in the baseline solution around Endeavour Seamount, the West Valley Segment, and the Heck Seamount (HSM) chain ( Supplementary Figure 6A ). In particular, the baseline solution shows a strong northwestward current flowing along the south-facing flanks of these topographic structures, which drives the westward moving branch of the plume fluid. This current is significantly weaker in the no-tide run, explaining the observed deviation in plume dispersion ( Supplementary Figure 6B ).

To investigate the interannual variability of plume dispersion patterns, we have calculated the plume occupancy rates for both the 2016 baseline and 2021 baseline simulation results. These rates represent the proportion of time that a given location is occupied by a hydrothermal anomaly at the depth of the non-buoyant plume cap ( Figure 12 ). The results reveal a shift in dispersion direction between 2016 and 2021. Specifically, the 2016 simulation results show that the non-buoyant plume cap primarily dispersed towards the north of the Endeavour Segment before splitting into two branches that continued northward and westward, respectively. In contrast, the 2021 simulation results demonstrate that the plume cap primarily dispersed towards the southwest of the Endeavour Segment then southward farther away from the ridge segment. This southwest dispersal direction is more consistent with previous field observations of the non-buoyant plume cap by Baker and Massoth (1987). Meanwhile, the north-to-south change of dispersion direction between 2016 and 2021 suggests that the trajectory of the plume can vary significantly on interannual time scales.

What might be the cause of such pronounced inter-annual variability? Large-scale abyssal flow plays a critical role in driving the dispersion of hydrothermal plumes, particularly in off-axis regions. Consequently, a plausible explanation for the change in dispersion direction between 2016 and 2021 could be a shift in the direction of abyssal flow near the depth of the non-buoyant plume cap, around 2100 m ( Figures 5C , 6C ). To investigate this proposed explanation, we have examined the flow field at a depth of 2100 m from the CMEMS global reanalysis (Fernandez and Lellouche, 2021) for July to September of 2016 and 2021, respectively ( Figure 13 ). The results exhibit a noticeable alteration in the flow direction within the regions surrounding the Endeavour Segment. Specifically, in 2016, currents flowed northward across the ultra-fine domain ( Figure 13A ), whereas in 2021, the predominant flow was southwestward over and to the east of the ridge axis, and westward over the western flank ( Figure 13B ). These flow directions are largely consistent with the plume dispersion directions observed in our simulations for both 2016 and 2021.

The hydrothermal anomaly, which is quantified by the concentration of the passive tracer released from the source vent sites, exhibits higher values in the non-buoyant plume cap in our no-tide simulation results when compared to the baseline solution results ( Figures 10 , 11 ). Generally, a larger hydrothermal anomaly indicates a lower degree of plume fluid dilution and thus less entrainment of ambient seawater. Previous laboratory and numerical studies have demonstrated that the presence of a cross-flow enhances the entrainment of ambient water into a buoyancy-driven plume, yielding increased dilution factors (e.g., Xu and Di Iorio, 2012). To quantify the mixing of plume fluid with ambient seawater, we calculate the volume integrated variance of the concentration of the passive tracer ( C ) released from the source vents as

In equation (3), C ′ = C − C m , C m being the tracer concentration averaged over the integration volume ( V ), which encloses the interior of the final domain, bounded by the seafloor and the sea-surface. The loss of variance of a conservative tracer is a well-established approach for defining mixing in studies of ocean turbulence (e.g., Burchard and Rennau, 2008; Wang and Geyer, 2018). More recently, destruction of salinity variance has been used to quantify mixing in estuaries (e.g., MacCready et al., 2018; Broatch and MacCready, 2022).

In our simulations, during the time between when the tracer is released and when it reaches the lateral boundaries of the integration volume, tracer variance is produced at the source vents on the seafloor and destroyed within the volume of the plume due to mixing of plume fluid with ambient seawater. A comparison among the four simulation runs demonstrates that: V a r ( C ) computed from the 2016 no-vent simulation results has the highest overall value over a 70-day period starting from the onset of tracer release; V a r ( C ) computed from the 2016 no-tide run has the second highest overall value; V a r ( C ) computed from the 2016 baseline solution and the 2021 baseline solution have the lowest overall values ( Figure 14 ). This suggests that the buoyancy flux from venting and tidal forcing both enhance mixing of hydrothermal plumes with ambient seawater, but that the buoyancy flux effect has the larger impact.

The preceding simulation results demonstrate the model’s proficiency in capturing the dispersion of hydrothermal plumes within a complex flow field that arises from the dynamic interaction between abyssal flow, ocean tides, and the topography of the ridge segment. Nonetheless, the comparison of the simulated flow velocity with direct measurements has revealed certain discrepancies, which may affect the accuracy of simulated plume distributions and trajectories with implications for future research efforts. The most apparent discrepancy is the absence of the 4-day ‘weather’ band oscillations in the simulated flow velocity ( Figure 4 ). Previous observations by Cannon and Thomson (1996) have indicated that these oscillations could have a significant impact on the dispersion of hydrothermal plumes near the ridge axis. Therefore, the absence or misrepresentation of these oscillations in our simulations has the potential to undermine their overall accuracy, particularly concerning the dispersion of hydrothermal plumes in close proximity to the ridge segment.

As mentioned in Section 4.2, the lack of 4-day oscillations in the simulated flow velocity could potentially be attributed to inaccuracies in the surface atmospheric forcing, such as wind stress, employed in the model. While a definitive understanding of the generation mechanism for the 4-day oscillations remains elusive, previous studies have indicated that these oscillations are associated with local storms (Cannon et al., 1991; Cannon and Thomson, 1996) and could be characterized as ridge-trapped subinertial waves (Allen and Thomson, 1993). In order to gain further insights into the generation of the 4-day oscillations along the Endeavour Segment and other regions of the JDFR, a potentially effective approach would involve diagnosing the factors contributing to the absence of these oscillations in our simulations. By making adjustments to the surface forcing and other model parameters, we can strive for an improved representation of these oscillations in future numerical studies. This approach holds promise for enhancing our understanding of the mechanisms underlying the 4-day oscillations and their behavior along the JDFR, making an important avenue for future research.

Additionally, we can achieve a more robust and comprehensive validation of our modeling framework by leveraging the continuously advancing capabilities of deep-diving autonomous underwater vehicles (AUVs). A carefully designed field program, employing an AUV equipped with dedicated in-situ instruments, would be able to provide direct observations of plume-specific tracers across a wide spatial area. Such future AUV-based plume surveys would be able to complement the use of the data collected by stationary instruments, such as the current-meter mooring data employed here, enhancing our assessment of model performance by enabling us to compare simulated plume distributions and trajectories with field observations.