Probabilistic Landslide Tsunami Estimation in the Makassar Strait, Indonesia, Using Statistical Emulation

This paper presents a significant advancement in the understanding of tsunamigenic landslide hazard across the length of the Makassar Strait in Indonesia. We use statistical emulation across the length of the continental slope to conduct a probabilistic assessment of tsunami hazard on a regional scale, across 14 virtual coastal gauges. Focusing on the potential maximum wave amplitudes (distance between the wave crest and the still‐water level) from possible tsunamigenic landslide events, we generate predictions from Gaussian Process emulators fitted to input‐outputs from 50 training scenarios. We show that the most probable maximum wave amplitudes in the majority of gauges are between 1 and 5 m, with the maximum predicted amplitudes reaching values of up to 10 m on the eastern coast, and up to 50 m on the western coast. We also explore the potential use of Gaussian multivariate copulas to sample emulator prediction input values to create a more realistic distribution of volumes along the continental slope. The novel use of statistical emulation across a whole slope enables the probabilistic assessment of tsunami hazard due to landslides on a regional scale. This area is of key interest to Indonesia since the new capital will be established in the East Kalimantan region on the western side of the Makassar Strait.

as an area of potential risk from submarine landslide-induced tsunamis. Notably, this is an area of particular risk due to its proximity to the announced location of the de facto capital city. This is further exacerbated by a number of additional factors all contributing to the oversteepening and destabilization of the continental slope: sediment erosion and deposition generated by the Makassar Throughflow currents, sediment influx from the Mahakam River, carbonate growth, and frequent proximal seismic events (Brackenridge et al., 2020;Gumbira et al., 2021;Hutchings & Mooney, 2021).
Research aiming to simulate the tsunamigenic landslide events in the Makassar Strait has generally followed scenario-based methodologies (e.g., Gumbira et al., 2021;Pranantyo et al., 2021). To develop the aforementioned scenario-based deterministic methodologies into a probabilistic tsunami hazard analysis within the Makassar Strait region requires the production of a catalog of synthetic (but realistic) tsunami events, involving the probabilistic sampling of key parameters.
Tsunami modeling packages generally require significant time and computational resources to run, a limitation only more prevalent when needing to simulate a large number of synthetic tsunami events (Abrahams et al., 2023;Ai et al., 2021;Davies et al., 2022). One novel, yet increasingly popular method to overcome this issue is the utilization of statistical emulation now breaking new frontiers in terms of the number of scenarios, locations, or scales used (e.g., Ehara et al., 2022;Gopinathan et al., 2021;Salmanidou et al., 2021). Statistical emulation provides a surrogate model which is fitted to a simulator's input-outputs overcoming the computational complexity and expense of carrying out simulations.
In this work, we make use of emulation to produce a set of probabilistic tsunami events within the Makassar Strait. This paper builds upon previous work focusing on single landslide event emulation (Salmanidou et al., 2017), or for specific slide locations (Løvholt et al., 2020;Salmanidou et al., 2019;Zengaffinen-Morris et al., 2022). To our knowledge, this is the first attempt to integrate statistical emulation into landslide-driven probabilistic tsunami hazard assessments for an entire region. We fit an array of Gaussian Process (GP) emulators to the relationship between the input parameters of a small set of simulated training scenarios and the maximum wave amplitude (distance between the wave crest and the still-water level; Sorensen (2005)) observed for a series offshore virtual gauges. We then use these fitted emulators to predict maximum wave amplitude at each virtual coastal gauge (Figure 1), given a large set of sampled input parameters. This enables us to observe the probabilistic distribution of maximum wave amplitudes at each gauge, with particular attention being paid to those surrounding the location of the new capital city.
We present the various steps that form this methodology in Section 2, and the respective outputs in Section 3. Later, we discuss the impact of the outputs of this work and touch upon its limitations and proposed future work in Sections 4 and 5.

Methods and Data
The workflow of this study consists of the following chronological steps: (a) A series of deterministic simulations are conducted for a subset of representative synthetic scenarios; (b) Quantities of interest, specifically maximum wave amplitude, are extracted from each simulation output; (c) Emulators are fitted to the relationship between the quantities of interest and input parameters; (d) Predictions of quantities of interest are generated from a larger sample of possible landslide scenarios; (e) The predicted outputs are interpreted to form the basis of a probabilistic assessment of tsunamigenic landslide hazard in the Makassar Strait.

Tsunami Modeling
A synthetic database of 50 landslide events was created to act as a training set ( Figure 2b). The Tsunami Open and Progressive Initial Conditions System (TOPICS), developed by Watts et al. (2003), was implemented to initialize the free surface deformation, which was used as input in the propagation model. Given the physical parameters of a submarine landslide, TOPICS uses the curve fits obtained from a nonlinear model to provide an approximation of the initialized water displacement field (e.g., Pranantyo et al., 2021;Watts et al., 2003). For each of the training events, location points (x and y) were sampled along a mid-continental slope depth contour of 1,500 m to generate a cross-slope profile. A second sampling process was used to sample a point along this profile. This two-stage process ensures that we can sample a point along the irregular shape of the continental slope sampling space, as defined by depth contours. Values for landslide volume (v) were sampled with minimum and maximum values of 50 and 600 km 3 , respectively. We used a Latin Hypercube design to sample the three variables in a way that efficiently explores the input space for building a GP emulator. For each training scenario, we also calculate the Great-circle distance between the landslide location and gauge as an additional emulator input (dist).
The sampled landslide volume and location constrain the physical parameters relating to the initialized water displacement. At each site, GEBCO Bathymetric data (GEBCO Compilation Group, 2022) was used to obtain the average slope angle (θ) and water depth (d) at the landslide point of origin, and slide direction (ϕ) perpendicular to the depth contour. The dimensions of the slide (i.e., its thickness, length and width) are maintained at the default ratio set within TOPICS (to create a set of initial surface conditions to be used by the wave propagation model), with the additional limiting of the maximum thickness (th max ) depending on the location of the slide along the shelf. This is done to account for the varying possible slide thicknesses, in kilometers, as determined by the study of seismic sections (Brackenridge et al., 2020), and are limited simply according to latitude as: For the generation, propagation, and runup of the tsunami training events, we use the Basilisk framework, an open-source Computational Fluid Dynamics (CFD) library specializing in spatially adaptive solutions. This library has been developed and used in a range of CFD applications, but specifically, the free-surface solvers (including relatives such as Saint-Venant and Boussinesq-type schemes) have previously been used and validated for tsunami applications Hayward, Whittaker, Lane, Power, Popinet, & White, 2022;Lane et al., 2017;Lee et al., 2021;Popinet, 2015).
From this framework, we use the non-hydrostatic multilayer finite-volume numerical scheme to simulate the generation and propagation of tsunamis. The discrepancies between non-hydrostatic and hydrostatic numerical schemes become more evident in long distance wave propagation, where over-steepening of the leading wave may occur in hydrostatic solvers (Popinet, 2020). The model vertically incorporates incompressible Euler equations with a free surface and gravity using n layers, which are vertically discrete and horizontally gridded. This layer integration can be represented by the following equations: where, in the x-z reference frame, k is the layer index, h k layer thickness, ϕ k the non-hydrostatic pressure, u k , w k the horizontal and vertical velocity components, g gravitational acceleration, η the free-surface height (sum of layer thicknesses and bathymetry height z b ), and the height of layer interfaces.
The numerical scheme implemented is described in further detail by Popinet (2020), and has been validated for wave propagation and runup Hayward, Whittaker, Lane, Power, Popinet, & White, 2022) and against various benchmarks, from simplified evolution of tsunami-like perturbations and periodic wave propagation to real cases such as the modeling of the dispersive features of the 2011 Tohoku earthquake tsunami (Popinet, 2020).
The Basilisk framework provides options for parallelization and a grid refinement algorithm using an adaptive tree walk protocol. This allows for the computationally and memory-efficient simulation of nonlinear runup and inundation with complete control over grid resolution, balancing the use of high resolution where required and computational efficiency where it is not (Popinet, 2013).
A digital elevation model (DEM) is constructed of the region comprised of a combination of SRTM and GEBCO datasets, reprojected and aligned to fit in the numerical model. We run each of the 50 training scenarios on systems containing 16 logical CPU cores. To initialize each simulation, the domain is filled with a bathymetry of the DEM and fluid layers are emplaced to represent sea-level at rest. The water surface displacement specified from TOPICS for the given landslide scenario is subtracted from the fluid layers, which results in the initial condition of the simulation. Boundary conditions are set to open flow to prevent the impact of wave reflections from the domain edges. Each scenario is compiled and run separately for a simulated time of 3 hr, where TOPICS is implemented at the initial time step, and the maximum spatial resolution is 33.75 arc-second.

Numerical Gauges
An array of 14 virtual numerical gauges were created across the study area ( Figure 1). Aside from the five central gauges located in the deeper, central waters of the Makassar Strait, the remaining gauges are focused offshore from concentrations of human populations and infrastructure. Particular attention is placed on the use of numerical gauges around the entrance to Balikpapan Bay (labeled as Balikpapan and Balikpapan (2) in Figure 1). This is to best record, and later probabilistically quantify, the hazard posed to the recently announced new Indonesian capital city of Nusantara, East Kalimantan. Each of the virtual gauges were configured to output water surface elevation at each computed time step (Figure 3), and we compute the maximum wave amplitude at each gauge for each simulation, to form the output of interest for emulator fitting. In some cases, the initial grid refinement of the areas where numerical gauges are located results in non-zero surface elevations at the first timestep (e.g., Sangatte as shown in Figure 3). In these cases, there is at least one training scenario that generates an initial displacement responsible for altering the grid refinement at the point of the gauge. In cases where the largest computed wave amplitude is at the initial time step, we discount this value to ensure we are exclusively training on the maximum amplitude of the simulated wave.

Statistical Emulation
For each of the training scenarios detailed in Subsection 2.1, maximum wave amplitude values were simulated for each numerical gauge location. Input-outputs consisting of scenario-specific landslide location (x and y), distance between gauge and landslide (dist) and landslide volume (v), and simulated maximum wave amplitude outputs were used to fit an emulator at each gauge location. To fit the emulators to the inputs-outputs of the 50 training scenarios, we use the Python package, Multi-Output Gaussian Process Emulator (MOGP) (Daub et al., 2022).
Once a fitted emulator was created for each numerical gauge location, a set 2,000 probabilistic scenarios were used to generate a set of predicted maximum wave amplitudes based on landslide location and volume. An important consideration when generating outputs from emulators is that inputs must fall within the same range as those used to fit them. Previous work has shown that the prediction uncertainty grows as inputs deviate from the emulator training points, and so to use inputs beyond the range of the training points would result in increasing prediction error moving further from the minimum and maximum training parameter values (Mohammadi et al., 2019). With this in mind, a spatial sampling space was generated along the continental slope, between 113 and 1,657 m depth contours (Figure 2a), within the depth ranges of the 50 training scenarios. This sampling space was then discretized into cell nodes 100 m apart (in both directions), generating 2,000 landslide locations from a possible approximately 1 million potential locations, to then be sampled from randomly (

Parameter Sampling With Gaussian Multivariate Copulas
Beyond the use of the independently-sampled emulator inputs outlined above, we explore the use of alternative emulator input parameter sampling techniques to improve the output probabilistic wave amplitude distributions. Submarine mass failures are notoriously difficult to catalog, and as such, most knowledge of the slope failures within the Makassar Strait comes from the aforementioned analysis of mass transport deposits from surveys of the seafloor. We use the mapped mass transport deposits published in Brackenridge et al. (2020) to construct a correlation between latitude and landslide volume, from which we calculate a new respective volume for each of the 50 training scenarios.
Using these latitudes and calculated volumes, we fit a Gaussian Multivariate Copula to the data using the Copulas library (datacebo, 2023). Copulas is a Python library used for modeling and sampling from multivariate distributions using Copula functions. This fitted copula was used to generate 2,000 synthetic latitude-volume parameter pairs. To ensure that the latitude was matched with a longitude value to place the landslide along the continental slope, we use the same spatial sampling grid as outlined in Section 2.3, and randomly sample along the longitudinal axis within the bounds of the sampling space grid contours at the given latitude ( Figure 5). It is worth noting that as this copula was fitted to the 50 uniformly sampled training landslides, the longitude and latitude sampling still remains relatively uniform, and the significant change is the sampling of volume. As it is crucial that all input parameters fall within the ranges of those used to train the emulators, we filter out any volume values that are less than 50 km 3 , or larger than 600 km 3 . After the landslide-gauge distance was computed for each landslide, these new copula-derived synthetic landslide parameters were used as alternative inputs for the fitted emulators. The result of this sampling approach is that we place a slight bias on larger landslides taking place at lower latitudes, as is reflected in previous work (Brackenridge et al., 2020).

Emulator Outputs
For each virtual gauge, we use the array of fitted emulators to produce a set of maximum wave amplitude predictions based on 2,000 synthetic scenarios along the Makassar Strait continental slope. On a gauge-by-gauge basis, we can observe the relationships between the input parameters (x and y location, distance between gauge and landslide, and landslide volume) and the output quantity of interest, maximum wave amplitude, albeit with variability due to the other three parameters varying (not being held constant) so this is not a clear sensitivity analysis. We show examples of these emulator outputs, for gauges surrounding Balikpapan, in Figure 6.
When observing the correlations between emulator input parameters and predicted maximum wave amplitudes, we focus on the inputs of volume and landslide-gauge distance. As shown in Figure 6, the latitude and longitude show a relationship with predicted wave amplitude that suggests higher uncertainties when the tsunami is required to interact with headlands and shallower waters. The relationships with volume appear more complex, with deeper-water gauges demonstrating larger uncertainty in prediction ability, yet some degree of positive correlation between volume and wave amplitude  is shown across fitted emulators, agreeing with previous work (Murty, 2003). We also note that closer landslides show greater uncertainty when observing the correlation between wave amplitude and Euclidean landslide-gauge distance. As identified in Dignan et al. (2020), the directionality of the sliding mass has significant potential to influence the directionality of the largest tsunami waves, and this could explain these results. For the majority of landslide scenarios, the easterly sliding mass causes the resulting tsunami to need to reflect off the eastern coastline before reaching the Balikpapan gauges (which can also be seen in the delay in gauge arrival of scenarios included in Figure 3), introducing wave amplitude stochasticity, and therefore uncertainty, in the process. The lack of linearity in the distance-wave amplitude plots in Figure 6 can also be attributed to the process of the tsunami interacting with nearby headlands and local bathymetry to reach the Balikpapan gauges, which are relatively well protected from tsunamigenic landslides generated at higher latitudes. In future works, empirical relationships for each parameter correlating with maximum wave amplitudes could be developed, similar to previous work, for example, in Sabeti and Heidarzadeh (2022).

Emulator Validation
We evaluate the accuracy of the emulator predictions by using a leave-one-out cross validation (LOOCV) method (as in Sarri et al. (2012)) with the 50 training scenarios. LOOCV is a thorough, yet computationally intensive, technique for model evaluation, wherein each data point is iteratively used as a single test sample, while the remaining data is used for training. This provides an unbiased estimate of a model's performance and addresses potential issues of over-fitting and bias.
We fit the emulator 50 times for each gauge, each time leaving out one scenario from the fitting data. The removed scenario is then used to have the emulator predict the output, resulting in 50 predicted values that we compare to the corresponding simulated values. The proximity of the emulator-predicted values to the actual simulated values is used as a measure of its accuracy. We show this comparison of predicted and simulated maximum wave amplitudes for each gauge in Figure 7. Whilst there were some, albeit negligible, differences between actual and predicted maximum wave amplitudes for the Bontang gauge showed the closest match (root mean square error, RMSE = 0.001), followed by Balikpapan (RMSE = 0.02). The gauges with the highest RMSE, and therefore poorest prediction accuracy, were Koealabadak (RMSE = 2.17, Mean Error = 19.63%) and Sangatte (RMSE = 1.50). We propose that a possible source of the fluctuation in prediction accuracy links to the linearity of the fluid flow required to reach the respective gauge. The Koealabadak gauge is located in relatively shallow water, and is also somewhat protected by the Mahakam Delta from tsunamis generated on the southern end of the continental slope, the processes of shoaling and refraction introduce non-linear processes to tsunamis reaching the gauge, introducing difficult-to-predict stochasticity.
When we put the RMSE values from the LOOCV analysis into context, the prediction ability of all fitted emulators was generally good. Excluding Koealabadak, the largest RMSE values ranged between 1.50 and 2.17, broadly meaning that the least accurate gauge had a prediction error of 0.71 m in maximum wave amplitude. The Koealabadak gauge had an RMSE of 2.17, which, as can be seen in Figure 7, is due to a minority of erroneous data points; the vast majority of data points within the LOOCV for Koealabadak show good prediction ability. As RMSE can be translated into the same units as the target value, in this case, centimeters, we can state that the majority of the gauges showed an RMSE of less than 71 cm, an error that is acceptable, particularly for the purposes of this probabilistic assessment. Balikpapan and Balikpapan (2) gauges also showed very good prediction accuracy, particularly the former, which is important for assessing landslide-based tsunami risk to the new capital. We must, however, consider this error when we analyze and extract results from the probabilistic wave amplitude distributions.

Wave Amplitude Distributions
Using the predicted maximum wave amplitude outputs from the fitted emulators for the independently-generated 2,000 landslide scenarios (as described in 2.3), we are able to produce a probabilistic histogram for each gauge location (Figure 8). From a probability perspective, the amplitude of each bar in the histogram can be interpreted as an estimate of the probability density of the corresponding maximum tsunami wave amplitude range.
When the emulator input parameters of the landslide scenarios are sampled using a fitted copula instead of independently, the resulting predicted maximum wave amplitude distributions show some differences. In all cases along the eastern coast of the Makassar Strait, we see that using copula-derived parameter sampling shifts the distribution toward smaller wave amplitudes and leaves independent sampling of parameters producing distributions with more frequent larger wave amplitudes. This is replicated in the southern center-channel gauges. The comparison is either negligible or mixed along the western coastal gauges. Broadly speaking, there is a more noticeable difference in gauges located in shallower waters, an effect that can be seen around the two gauges offshore Balikpapan (Figure 9).
The resulting probabilistic wave amplitude histograms (Figure 8) vary spatially. For all gauges, the most likely maximum wave amplitude was around 0-2 m, but the wave amplitude distributions vary from gauge-to-gauge. The Strait (Center) gauge shows significant probability of larger maximum wave amplitudes, but this is attributable to the proximity of the gauge to the landslide locations relative to the other deep water gauges. On the eastern coast, the shape of the distribution was broadly similar with maximum predicted wave amplitudes being 6.5-10 m. Along the western coast, the wave amplitude distributions have more prominent peaks around 0-2 m, however, exhibit much higher maximum predicted wave amplitude values (up to 50 m). It is widely accepted that there is a strong relationship between water depth and maximum wave amplitude as the process of shoaling begins (Li & Raichlen, 2002;Synolakis, 1987), and so it is probable that these larger (yet relatively unlikely) tsunami events occur along this coast due to the shallower gauge depth and shoaling across the continental shelf. Exceptionally high maximum wave amplitude values, such as for Bontang, can generally be attributed to the maximum recorded wave amplitude recorded from the training scenarios ( Figure 3). In the case of the Bontang gauge, a combination of large maximum amplitude training scenarios (attributed either to the largest volume landslides or to the gauge proximity to the source) and the shallow water depth at the gauge location are responsible for the small number of very large maximum amplitude values. Out of the proximal cluster of training events close to the Bontang gauge, the maximum landslide volume was just over 500 km 3 , which generated a maximum wave amplitude of approximately 21 m. When the increased number of synthetic scenarios fill the sampling space more completely, we see larger maximum wave heights where these proximal events have volume values closer to 600 km 3 . Similarly high maximum wave amplitude values can be seen along this western coastline caused by the same impact of scaling up the number of events and the full use of the volume parameter range. The impact of this is particularly prominent in Bontang (and less so in neighboring virtual gauges in Sangatte and Koealabadak) due to the Bontang gauge being located directly on the top edge of the continental slope and is, therefore, able to capture the maximum amplitude as shoaling occurs and before shallow water wave energy dissipation occurs.
Across the gauges, the distribution is typical of tsunami maximum wave amplitudes computed in other works, with decreasing frequency as maximum wave amplitudes increase (Choi et al., 2002;Fukutani et al., 2015;Zhang et al., 2020). As previously mentioned, we must consider the less than ±0.71 m influence that the prediction error has on these distributions when interpreting the most probable maximum wave amplitudes. For Koealabadak, the gauge with the poorest prediction performance in LOOCV (RMSE = 2.17), this error represents a 12% error in maximum probable wave height. Comparatively, the Balikpapan gauge LOOCV error (RMSE = 0.02) represents a 0.001% error in maximum probable wave height.

Discussion
Experimenting with the use of fitted copulas to sample emulator input parameters provides an interesting source of discussion; when conducting a probabilistic assessment, we must, to some extent, make some assumptions regarding both the ranges and sampling of the model input parameters. In this work, we assumed a uniform random sampling to select the landslide locations. Although previous studies on the landslide deposits indicate that there is a prevalence of more voluminous landslides in the South of the Strait, landslide locations are spread preconditioning slope instability along the Strait axis (Brackenridge et al., 2020). Future slope stability assessments could shed light on potential landslide locations and could be therefore used to inform hazard assessment in the area. As the sampling of latitude and longitude remained random, the use of a copula only affected the sampling of volume and doing so reduced the frequency of larger events when compared to independent sampling. Whilst the copula-derived inputs are arguably more realistic, using these inputs with smaller high-volume landslide frequency naturally influences the output maximum wave height histograms. As we also concentrate the high-volume landslide events around lower latitudes, we also see the differences between copula-derived inputs and uniformly-sampled inputs more clearly in higher latitudes where there is a more significant reduction in high-volume events. Obviously, the use of the most accurate sampling distributions is critical to ensure the validity of the probabilistic assessment. However, when it comes to submarine landslides, this can often be problematic as we are restricted by the availability of data and previous work done to map the volumes and extents of historic slides within the study area. Even still, submarine landslide catalogs carry significant error in identifying past slides and the presence of a past landslide does not equal the presence of another in the future of equal size. A balance must be found between making use of existing knowledge and geological records to create an informed sampling distribution, and using independent sampling to ensure we are not introducing excessive bias and error into the model output and, ultimately, the tsunami hazard assessment. In this paper, we argue that, due to the incomplete data available to fully identify the parameters of a larger sample size of landslides across the continental slope, an independent sampling strategy should be used to create the final maximum wave height frequency histograms.
The tsunami risk to the future capital city of Balikpapan is difficult to define from this analysis alone, however, we are able to observe the possibility of wave heights in some scenarios in excess of 15 m. Without simulating final stage runup and inundation, it is difficult to estimate the land area likely impacted by a tsunami of this height, yet the gauge's proximity to the coastline would suggest that there is a significant risk of landslide-generated tsunami inundation in the Balikpapan region from credible event scenarios.

Conclusions
In this work, we develop the use of statistical emulation for wave amplitude prediction for tsunamigenic submarine slope failures. We also evaluate the effectiveness and accuracy of emulators for the purpose of expediting the process of tsunami simulations in the Makassar Strait, with particular attention given to the sampling of emulator input parameters. The results of this paper represent a first step in conducting a probabilistic tsunami hazard assessment from tsunamigenic landslides in the region, providing maximum wave amplitude distributions for an array of gauges. We pay particular attention to gauges in the area offshore Balikpapan, in proximity to the new Indonesian capital, and the implications of our methodology when understanding the tsunami risk posed to the future city.
From predictions generated from a GP emulator fitted to input-outputs from the 50 training scenarios, we observe that the most probable maximum wave amplitude in the majority of gauges is between 0 and 2 m. On the eastern coast, the maximum predicted amplitudes were up to 10 m, and up to 50 m on the western coast. We show that the emulators are able to predict maximum wave amplitude with a good degree of accuracy, with the majority of gauges having an RMSE between simulated and predicted values of less than 0.7 m.
This paper also explores the potential use of Gaussian multivariate copulas to sample emulator prediction input values to create a more realistic distribution of volumes along the continental slope. We conclude that, despite the importance of sampling distributions used for probabilistic analysis being representative of the hazard itself, using a sampling distribution on incomplete information is something that should be avoided. On this basis, our analysis used independent sampling of individual distributions of the input parameters (not necessarily uniform), and in doing so, provides somewhat of a "relatively agnostic" probabilistic analysis.
Further work will advance this methodology by including the emulation of late-stage runup and inundation, allowing us to conduct end-to-end emulation of tsunami impact on land from input parameters. We also demonstrate the benefit of using statistical emulation for tsunami simulators of increasing complexity, in an attempt to exploit their increased efficiency to maximum effect.

Data Availability Statement
The training simulation scenario data (including gauge time series outputs, wave propagation animation, and files required to run simulations) used to fit the emulators, as described in this work, are freely available on Zenodo in Dignan et al. (2023)

Table B1
Latitude and Longitude Location of All Virtual Numerical Gauges