Mixed-Convective Processes Within Seafloor Sediments Arising From Fresh Groundwater Discharge

The dependence of near-shore ecosystems on the freshwater component of submarine groundwater discharge (SFGD) is well recognized. Previous studies of SFGD have typically assumed that SFGD occurs through aquitards that are in direct contact with seawater. These studies provide no guidance on the distribution of freshwater discharge to the seafloor where SFGD occurs through sandy sediments, even though in most situations, seabed sediments are permeable. We find that SFGD may occur in unconfined, seafloor sediments as density-driven flow in the form of fingers, or otherwise, diffusive freshwater discharge is also possible. Unstable, buoyancy-driven flow within seabed sediments follows similar patterns (except inverted) to the downward free convection of unstable (dense over less-dense groundwater) situations. Consequently, the same theoretical controlling factors as those developed for downward mixed-convective flow are expected to apply. Although, there are important differences, in particular the boundary conditions, between subsea freshwater-seawater interactions and previous mixed-convective problems. Simplified numerical experiments in SEAWAT indicate that the behavior of fresh buoyant plumes depends on the aquifer lower boundary, which in turn controls the rate and pattern of SFGD to the seafloor. This article provides an important initial step in the understanding of SFGD behavior in regions of sandy seafloor sediments and analyses for the first time the mixed-convective processes that occur when freshwater rises into an otherwise saline groundwater body.


INTRODUCTION
Subsea fresh groundwater discharge (SFGD) is the release of freshwater from seafloor sediments, and has been identified as an important means to transport dissolved nutrients to the ocean, thereby having a significant influence on marine ecology and benthic organisms (e.g., Johannes, 1980;Moore, 1999). Knowledge of the salinity patterns in subsea sediments is important for several reasons, including for the understanding of benthic ecosystems in regions of SFGD, and for the design of SFGD measurement approaches. The distribution of groundwater discharge to the sea depends on the characteristics of geological formations. Fresh groundwater discharge to intertidal zones is often associated with unconfined aquifers, whereas SFGD, occurring beyond the intertidal zone, typically requires low-permeability confining units that act to limit freshwater-seawater mixing that would otherwise occur near the shoreline (e.g., Jiao et al., 2015;Michael et al., 2016). Predictions of the extent of subsea fresh groundwater typically assume that the lowpermeability layers that preserve freshwater and allow SFGD to occur are in direct contact with seawater (e.g., Kooi and Groen, 2001;Bakker et al., 2017;Solórzano-Rivas and Werner, 2018;Werner and Robinson, 2018;Solórzano-Rivas et al., 2019). However, the seafloor more often comprises high-permeability sediments (>50% of the global shelf sediments have high permeability; Riedl et al., 1972). Yet, this configuration has not been widely considered in studies of SFGD.
The occurrence of freshwater in high-permeability sediments overlain by seawater, as occurs when SFGD passes through permeable seafloor sediments, creates a mixed-convective condition. That is, flow processes are expected to be controlled by both buoyancy-driven flow (i.e., free convection) driven by water density differences and hydraulic-driven flow (i.e., forced convection) arising from groundwater heads in subsea sediments that exceed those of the sea. Thus, solute distributions in subsea aquifers may resemble those of other mixedconvective situations. However, studies of mixed-convective or free-convective processes where buoyancy is created by salinity gradients usually involve descending plumes of higher-density fluid that contaminate underlying lower-density groundwater (e.g., Webster et al., 1996;Smith and Turner, 2001;Stevens et al., 2009;Xie et al., 2011). Conversely, the upward movement of lower-density groundwater (e.g., as expected to arise during SFGD) is rarely explored in the solute transport context, although upward, buoyancy-driven groundwater flow has received significant attention in the field of heat transport and geothermal phenomena. For example, Kurylyk et al. (2018) utilized temperature-based methods to quantify submarine groundwater fluxes in seafloor sediments offshore of eastern Canada. However, their study suggested that the flow patterns inferred from seafloor temperature-depth profiles appear to be influenced by both density-driven and geothermal processes.
Perhaps the most pertinent prior investigation of mixedconvective transport accompanying SFGD is that of Moore and Wilson (2005), who were perhaps the first to recognize that SFGD may be driven by upward buoyancy forces. Based on temperature measurements of subsea groundwater (1.5 m below the seabed), Moore and Wilson (2005) interpreted that a sudden drop in the ocean temperature produced upward motion of warmer (i.e., less dense) groundwater to the ocean. However, previous research offers little guidance on the characteristics of free-convective or mixed-convective flow within seafloor sediments. In particular, whether mixed-convective flow within seafloor sediments can be characterized according to buoyancy theory developed for other unstable, solute transport conditions (e.g., leading to downwardmoving fingers of higher density; e.g., Wooding et al., 1997) is unclear. This includes the application of several dimensionless variables that are used to categorize mixed-convective problems (e.g., Wooding et al., 1997;Simmons et al., 2010). Additionally, to our knowledge, no quantitative evaluation of mixed convection driven by the upward movement of lower-salinity groundwater has been reported in the literature. Rather, upward convection of lower-density groundwater has only been explored where buoyancy is created by temperature gradients (e.g., Irvine et al., 2015).
The aim of this study is to undertake a review of the existing buoyant theory, including non-dimensional numbers used to characterize mixed-convective flow, and explore whether this theory can be applied to SFGD through high-permeability seafloor sediments. A small number of highly idealized numerical simulations are used to provide an initial demonstration of freshwater-seawater mixing accompanying the rise of fresh buoyant plumes through permeable seafloor sediments. We expect that the upward movement of less-dense groundwater will show similar, albeit inverted, characteristics (e.g., buoyancydriven fingers of freshwater) to the downward movement of more dense saltwater that has been comprehensively assessed in numerous prior investigations.

CONCEPTUAL MODEL
Two hypothetical conceptual models of SFGD are considered, both of which involve freshwater-seawater mixing within highpermeability, seabed sediments, as shown in Figure 1. Model A represents the situation where the seafloor sediments are underlain by a higher-permeability freshwater source, whereas Model B involves underlying sediments of lower permeability (i.e., an aquitard). The former is intended to reflect the situation of SFGD where an underlying aquifer sub-crops to the seafloor and is overlain by sand, while the latter represents SFGD passing through a leaky aquitard into more permeable seafloor sediments.
An important unknown for the conceptual models considered here is the distribution of freshwater within the seafloor sediments. While Figure 1 shows buoyant freshwater plumes that take the shape of fingers, similar to classic unstable flow patterns for downward free convection, freshwater-seawater mixing may alternatively be dominated by dispersive process (at least theoretically), leading to different salinity patterns to those represented in Figure 1. Most field studies of SFGD adopt methods that assume that SFGD occurs in a distributed form (i.e., presuming dispersive freshwater-seawater mixing in seafloor sediments), with devices located on the seafloor (i.e., seepage meters) in regular grid patterns (e.g., Taniguchi et al., 2003;Michael et al., 2005;Kurylyk et al., 2018). However, the buoyant freshwater plumes that are hypothesized here (Figure 1) may require alternative measurement strategies. For example, the spacing between, and size of, freshwater plumes, should they occur, may influence the deployment of seepage meters.
Although freshwater plumes within seafloor sediments are expected to behave similarly to the well-studied downwardmoving saltwater plumes of prior studies (e.g., Xie et al., 2011), there are important differences between the conceptual model of Figure 1 and previous analyses of downward-moving plumes. Consider for example, the Elder problem (e.g., Elder et al., 2017), which was transformed to free-convective solute motion by Voss and Souza (1987), and is a common benchmarking problem for assessing numerical models. Various modifications of the Elder problem have occurred to explore other aspects of free convection problems. For example, Xie et al. (2010) introduced mechanical FIGURE 1 | Two conceptual models of a subsea unconfined aquifer subject to SFGD. Blue is freshwater, and purple is seawater. Model A (left) involves sandy sediments overlying a subcropping, higher-permeability aquifer, and Model B represents an aquitard overlain by sand.
dispersion and changed the lower boundary condition to represent more realistic conditions, similar to the approach of Post and Kooi (2003). The conceptual model applied to SFGD (Figure 1) is, conceptually at least, an "inverse" of the solute-based Elder problem, where the aquifer is initially seawater-filled and a source of constant freshwater is introduced at the bottom boundary. The initial seawater conditions in the aquifer allow for the transience of salinization to be observed. However, the boundary conditions differ to those of the Elder problem to improve the representation of the expected SFGD process. That is, specified-head boundaries are imposed at the top (seawater) and the bottom (freshwater), whereas the Elder problem uses no-flow boundary conditions at all edges and specified hydraulic heads at the top corners (e.g., Guo and Langevin, 2002). Thus, the behavior of SFGD within seafloor sediments, under the conditions presented in Figure 1, cannot be inferred from previous studies. Further details of numerical models are provided in section "Numerical Simulation of SFGD."

REVIEW OF MIXED-DENSITY NON-DIMENSIONAL NUMBERS APPLICABLE TO SFGD THROUGH SEAFLOOR SANDY SEDIMENTS
The characterization of buoyancy-induced flow associated with unstable systems has been the subject of extensive research. One the most common non-dimensional numbers applied for this purpose is the Rayleigh number (Ra) (e.g., Simmons et al., 2010). Ra is the ratio between buoyancy forces (inducing instabilities in the form of solute fingers) and dispersion forces (tending to resist the formation of unstable solute fingers). When freshwater and seawater interact, Ra can be defined as (e.g., Post and Kooi, 2003): where K is the aquifer hydraulic conductivity (L T −1 ), ρ is the difference between the seawater (ρ s ) and freshwater (ρ f ) densities (M L −3 ), H is the aquifer thickness (L), and D is hydrodynamic dispersion (L 2 T −1 ). D is equal to the sum of molecular diffusion [D 0 , (L 2 T −1 )] and mechanical dispersion, which is often simplified to α L v L , where α L is the longitudinal dispersivity (L) and v L the advective flow velocity (L T −1 ) (i.e., transverse dispersivity, α T , is often neglected).
There are several alternative Ra formulations that have been reported in the literature. For example, some studies have presumed that D consists only of D 0 , with mechanical dispersion neglected (e.g., Post and Kooi, 2003;Stevens et al., 2009;Simmons et al., 2010;Xie et al., 2011). These cases involve only free convection, i.e., hydraulic-driven convection is not considered. Alternative expressions for Ra have been developed for mixed-convection processes. For example, in mixed-convective problems examined by Simmons and Narayan (1997) and Smith and Turner (2001), both mechanical dispersion and molecular diffusion were included in their definitions of D. However, Simmons and Narayan (1997) adopted α T and Smith and Turner (2001) adopted α L in defining D. This study adopts the Smith and Turner (2001) The value of Ra has been used to predict the occurrence of unstable solute motion, typically in the form of solute fingers (e.g., Wooding et al., 1997). The critical Ra is defined as the threshold value at which buoyancy forces overcome dispersion forces that inhibit finger formation. This threshold value differs depending on conditions in which free convection occurs. For example, in the case of free convection in hydrogeological settings, Stevens et al. (2009) and Moore and Wilson (2005) refer to the critical Ra of 4π 2 (i.e., for the appearance of solute fingers) reported by Lapwood (1948) for temperature gradients in porous media. However, alternative critical Ra values are offered by other authors for free convection problems. According to van Reeuwijk et al. (2009), the historical incongruency in critical Ra estimates is related to such aspects as the numerical approach, governing equations, and slight differences in the initial conditions. They suggest that the critical Ra for free convection involving solutes is zero (in other words, solute fingers form without the need for a "dispersive boundary layer"; see below).
Applications of Ra theory to mixed-convective processes include the study of Simmons and Narayan (1997), who imposed linearly varying hydraulic heads across the top boundary, to induce lateral flow (i.e., flow parallel to the solute source boundary). The transverse dispersion created by this lateral flow was incorporated into the stabilizing dispersive force in assessing the critical Ra. By incorporating α T v L in D, Simmons and Narayan (1997) include forced convection (i.e., v L ) in the definition of Ra. They found that the critical Ra ranged between 300 and 500 under these conditions. Smith and Turner (2001) considered a somewhat similar situation to that of Simmons and Narayan (1997), except fresh groundwater discharged to the upper solute boundary condition (i.e., the estuary) in Smith and Turner's (2001) analysis of estuaryaquifer interaction. They obtained a critical Ra value of 5 for the onset of finger development. Solórzano-Rivas and Werner (2018) adopted the Smith and Turner (2001) formulation for the analysis of freshwater discharge through subsea aquifers, and found a critical Ra of about 2 for the salinization of submarine aquitards overlying fresh offshore aquifers. The application of Ra to mixed-convection problems introduces challenges in establishing general values for the critical Ra, given differences in the representation of forced convection in defining Ra across the abovementioned studies. Wooding et al. (1997) analyzed the development of boundary layers within mixed-convective flow systems (i.e., where freshwater flows upwards towards a saltwater boundary). They derived an expression for the steady-state boundary layer thickness, δ (L), defined as the thickness of solute formed (by dispersion) in the presence of upward freshwater flow, q z (L T −1 ), as: Wooding et al. (1997) combined Equations 1 and 2, by equating H to δ (i.e., conceptually, this is the case where the boundary layer encompasses the entire aquifer thickness), to define a boundarylayer Rayleigh number, Ra δ , namely: The thickness of δ that leads to the onset of unstable solute fingers defines the critical Ra δ . Wooding et al. (1997) undertook laboratory experiments using a Hele-Shaw cell to determine the critical Ra δ , including different rates of lateral inflow at the right boundary and vertical outflow to the top, and using various inclination angles of the cell. They also performed numerical experiments that reproduced the two-dimensional flow and solute transport behavior observed in the Hele-Shaw experiment. From the laboratory and numerical experiments, they found Ra δ values between 8.9 and 9.8, concluding that a good estimate for the critical Ra δ is approximately 10. Another widely used non-dimensional parameter for the characterization of mixed-convective processes is the mixed convection ratio (M). M describes the relationship between buoyancy-driven forces and hydraulic-driven forces as (e.g., Simmons et al., 2010): where h l (-) is the hydraulic gradient over a distance l. An alternative expression for M, through substitution of Darcy's Law, is given by Smith (2004), as: By combining Equations 2 and 3, it is apparent that Ra δ and the formulation for M given in Equation 5 are the same, i.e., M = Ra δ . Surprisingly, this has not been reported previously to the authors' knowledge. The implications of this are that critical values for Ra δ also apply to M. However, there is no evidence that the equilibrium value of M = 1 suggested by Simmons et al. (2010) (whereby if M is larger than 1 the problem is said to be free-convection dominated, and forced convection is thought to be the dominant process if M is less than 1) has application in terms of Ra δ . That is, Ra δ = 1 has not been reported as a significant value previously. It follows that the onset of instabilities is not predictable through comparison of free and forced convective forces (i.e., through the use of M), as expected given the role that dispersion plays in instability initiation. Nevertheless, Stevens et al. (2009)

NUMERICAL SIMULATION OF SFGD
Here, the conceptual models illustrated in Figure 1 are adopted, except the situation is assessed whereby the sandy seafloor sediments are presumed to be originally filled with seawater. Numerical modeling, using SEAWAT (governing equations are found in Guo and Langevin, 2002; not repeated here for brevity), was used to assess two simple situations of SFGD through sandy seafloor sediments. The two cases (Models A and B; see Figure 2) involve contrasting mixed-convective forces arising from the inclusion or omission of an underlying aquitard. This is reflected in the corresponding values of Ra and Ra δ , which are reported for each case and compared to previously reported critical values (Wooding et al., 1997;Smith and Turner, 2001) in section "Review of Mixed-Density Nondimensional Numbers Applicable to SFGD Through Seafloor Sandy Sediments." The calculation of Ra and M for Model B does not include the aquitard thickness (i.e., only H is regarded as per the Equation 1 definition). Figure 2 is a simplified section showing two possible configurations of a subsea aquifer subjected to SFGD, which is homogeneous, isotropic, and of length L (L) and thickness H (L). The upper boundary represents seawater immediately above the seafloor, simulated by a high hydraulic conductivity (i.e., 100,000 m/d) and a specified-head condition equal to z s , assigned to the top row of the model. SEAWAT converts z s to a constant pressure (Langevin et al., 2008) reflecting the column of overlying seawater, of density ρ s . The solute concentration condition in this top layer depends on the flow direction, FIGURE 2 | Conceptual models for the numerical implementation of SFGD through an unconfined subsea aquifer, where the left side shows sandy sediments overlying a subcropping, higher-permeability aquifer, and Model B represents an aquitard overlain by sand. Red and blue represent seawater and freshwater boundary conditions, respectively. The equivalent freshwater head at the base of the aquifer, when the aquifer is seawater-filled (i.e., the initial condition), is shown as h f . whereby seawater concentration (i.e., solute concentration = 1) was assigned to any inflowing water, whereas discharge to the sea occurs at the ambient groundwater concentration. This type of solute concentration condition avoids salt accumulation at the boundary in an unrealistic manner (i.e., upstream dispersion; Irvine et al., 2021) and is consistent with the approach of Solórzano-Rivas and Werner (2018). The lower freshwater boundary condition differs between Models A and B. In Model A, the lower boundary reflects a situation where an underlying layer of higher permeability occurs, containing freshwater. That is, there is no restriction to the entry of freshwater through the lower boundary of Model A, represented by a specified-head condition equal to h f , which was converted to a constant pressure (by SEAWAT) based on a water density of ρ f , as shown in Figure 2. In Model B, the lower boundary is composed of lower-permeability sediments (i.e., reflecting aquitard material) of thickness H L , with the specified head h f imposed on the bottom row of the model (i.e., the base of the aquitard). The solute concentration of the lower boundary for both Models A and B was set to freshwater (i.e., solute concentration = 0), although this depends on the flow direction in a similar manner to the upper boundary.
No-flow boundaries on the left and right edges of the model domain reflect mostly vertical flow processes. Post et al. (2007) and Langevin et al. (2008) offer the following formulation for vertical flow, q z , under mixed-density conditions: where ρ a is the average density between ρ s and ρ f (M L −3 ) and h f l is the hydraulic gradient (in terms of equivalent freshwater heads) over a distance l.
The aquifer is assumed to be initially full of seawater but is underlain by a fresh groundwater source. The lower boundary head of both models is chosen so that the initial condition is hydraulically stable, primarily to reflect the freeconvective situation (i.e., neglecting forced convection) that occurs in the Elder problem and its many variants, thereby providing opportunities to compare upward buoyancy-driven flow to the downward movement of solute fingers observed by others (e.g., Xie et al., 2011). That is, the hydraulic head of the fresh groundwater source is equal to the equivalent freshwater head, h f , at the bottom of aquifer (at the beginning of each scenario). Considering an initial seawater hydrostatic condition, h f is given by ρ s z a (where z s and z a are defined in Figure 2). Consequently, the only forces driving flow (initially) are buoyancy forces, because the hydraulic gradient is initially zero. Thus, at the beginning of the scenarios, Equation 6 reduces to: Equation 7 is usually referred to as the convective velocity (e.g., Simmons et al., 2010). Vertical flow velocities in the numerical model, from the initial time-step, were found to be consistent with Equation 7 (results not shown for brevity).
Although the initial condition is a free convective condition, density changes due to the introduction of freshwater through the lower boundary create head gradients that induce boundary inflows, and the system becomes mixed convective. This differs to the free convection Elder problem, which is bounded by noflow conditions, including the upper salt source boundary. Thus, although our situation is initially free convective, the situation becomes mixed convective during the course of simulations.
The model is a two-dimensional domain represented by a finite-difference grid of uniform discretization, both vertically and horizontally, of z = x = 0.1 m. Three cases are evaluated for each model (i.e., a total of six numerical simulations) to briefly explore the role of dispersion in the mixed-convective processes associated with SFGD through the seafloor sediments, namely a Base Case with α L = 0.1 m, Case 1 with α L = 0.5 m, and Case 2 with α L = 1 m. These values fall within the range of "moderate" and "high" reliability values recommended by Zech et al. (2015). Other parameters used in all three cases for Model A are: {H, L, z a , z s , h f , K, ρ f , ρ s , D 0 , α L /α T , ε} = {50 m, 100 m, 80 m, 160 m, 162 m, 2.5 m/d, 1000 kg/m 3 , 1025 kg/m 3 , 8.64 × 10 −5 m 2 /d, 10, 0.2}, where ε (-) is effective porosity. The three cases of Model B differ to those of Model A in that K a is 0.001 m/d and H L is 1 m. The parameters adopted here correspond to typical, field-scale values used in other SFGD-related studies (e.g., Knight et al., 2018;Solórzano-Rivas and Werner, 2018). Timescales for numerical models to reach steady state or quasisteady state conditions differed between Models A and B, as shown in section "Results and Discussion." Consequently, Model A was simulated for 5 years and Model B for 100 years. Timescales for freshwater to reach the seafloor depend on the velocities of rise of buoyant freshwater fingers. We adopted the tip of the highest buoyant finger (HBF) to characterize mixed convective velocities (assuming time-constant velocities) in a similar way to the approach of Xie et al. (2011). The tip was defined by the 0.85 isochlor (i.e., 85% of seawater concentration).

RESULTS AND DISCUSSION
The Influence of the Lower Boundary right column) at different times. The bottom row in all three figures represents the steady-state conditions, which in the case of Model B is a dynamic equilibrium, also referred to as quasi-steady sate. That is, the Model B distribution of SFGD within seafloor sediments is temporally unstable, involving finger patterns that change in time.
Figures 3-5 highlight substantial differences in mixedconvective processes caused by the existence (or not) of a leaky aquitard below the subsea aquifer. The most important difference between Models A and B, in the context of SFGD characteristics, is that unstable solute patterns persist under quasi-steady state in Model B, whereas in Model A, the flow instabilities are temporary, and the system reaches a steady state condition in which the aquifer is completely fresh. The occurrence of the leaky aquitard creates other important differences in mixed convective processes. For example, freshwater fingers produced from Model A are fresher, for fingers of a comparable height, to those produced by Model B (e. g., Figures 3c,g). This leads to sharper freshwater-saltwater interfaces in Model A results. The timescales for the rise of buoyant fingers also differ between Models A and B. For example, the time to obtain freshwater fingers of a roughly similar height differs substantially, as evident in the timing of 292.3 days for Figure 4c (Model A) and 2191.9 days for Figure 4g (Model B). Thus, the lower aquitard in Model B significantly restricts finger speeds, which is an intuitive outcome. HBF in Model A for the Base Case, Case 1, and Case 2 reaches elevations of 119.15, 110.95, and 110.75 m in Figures 3c, 4c,  5c, respectively. These represent average finger speeds of 48.9, 38.7, and 38.4 m/year. Lower velocities are produced in all three cases of Model B. That is, HBF heights were 129.85, 127.15, and 120.95 m in Figures 3g, 4g, 5g, leading to velocities of 8.31, 7.86, and 6.82 m/year, respectively. There is evidence of boundary effects in Model A, associated with the vertical noflow boundaries, that was also observed by Xie et al. (2011). Greater velocities of fingers adjacent to no-flow boundaries are attributed to the absence of the host fluid (i.e., seawater) moving in the opposite direction to the buoyant fingers, which is otherwise expected to retard the finger upwelling speed. For that reason, fingers adjacent to vertical boundaries were neglected in assessing HBF speeds.
According to Xie et al. (2011), the theoretical velocity of fingers is 13.1 m/year, based on ( ρK/ρ f ε) × f , where f (0.115) is the Xie et al. (2011) corrective factor to the theoretical convective velocity (i.e., the ratio between Equation 7 and ε). Interestingly, the Xie et al. (2011) corrected convective velocity, which neglects the leaky layer in Model B, is closer to the results in Model B (i.e., 8.31, 7.86, and 6.82 m/year), compared to the finger velocities in Model A (i.e., 48.9, 38.7, and 38.4 m/year). However, finger velocities in Model A could be better predicted applying values of f that fall within the range of corrective factors proposed by Post and Kooi (2003) and Wooding (1969) (i.e., 0.22 and 0.446, respectively). Figure 6 illustrates the temporal variability in total solute mass within model domains, for the six scenarios considered in this study. The fluctuations in total mass in Model B are caused by mixed-convective instabilities, which persist as quasi-steady state conditions, as discussed above. The time for Model A to reach steady-state conditions (approximately 2 years) is much less than the time required for Model B to reach quasi-steady state conditions (approximately 10 years). This is consistent with the difference between buoyant fingering speeds found in Models A and B, as discussed above. The larger velocities in Model A (e.g., 48.9 m/year, 38.7 m/year, and 38.4 m/d in Figures 3c, 4c, 5c) led shorter timescales to reach steady state compared to Model B, in which velocities were 8.31, 7.86, and 6.82 m/year (Figures 3g, 4g, 5g). Figure 6 shows that variations in dispersion had no noticeable impact on the time required to reach steady and quasisteady state conditions. This is consistent with the relatively   small change between cases (i.e., from Base Case to Case 2) in buoyant fingering speeds, as observed also by Xie et al. (2011). The little influence of dispersion on buoyant fingering speeds is also in accordance with the definition of the theoretical convective velocity proposed by Xie et al. (2011), which neglects dispersion. While dispersion has an unimportant role in the speed of buoyant finger rise, dispersion does appear to influence the quasi-steady state buoyant fingering patterns associated with mixed-convective processes (Figures 3-5). Figure 7 compares the three cases (i.e., Base Case, Case 1, and Case 2) of Model B, showing quasi-steady state salinity distributions in profile, and the temporal variation in SFGD solute concentration at the seafloor using a timeframe of 85 years. Subplots Figures 7a-c correspond to Figures 3h, 4h,  5h, respectively. Figure 7 shows that the number of fingers reaching the seafloor decreases with increasing dispersion due to the widening of fingers ; a similar observation to those of Xie et al. (2011) for dense, downward-moving fingers.

The Influence of Dispersion
The temporal and spatial variability in the salinity distribution across the sea floor have important implications for the direct measurement of SFGD. For example, the deployment and size of seepage meters (e.g., Burnett et al., 2006) in areas where mixed-convective processes occur require consideration of the spatial variability in freshwater discharge, given the irregular distributions of SFGD in Figures 3-7. Mixed-convective flow leads to regions of freshwater upwelling and seawater downwelling, and therefore, SFGD measurement methods that can detect seafloor fluxes in both directions (inflow and outflow) may assist in detecting mixed-convective phenomena. Additionally, the duration of seepage meter placement should consider the temporal variation of SFGD created by mixed convective processes. Freshwater upwelling appeared to be relatively stable over durations of days-to-weeks, but may vary substantially over longer timeframes, at least for the cases illustrated in Figure 7.
The non-dimensional parameters Ra and M were determined, using Equations 1 and 5, respectively, from the quasi-steady state results of Model B. The steady-state solution of Model A has no density gradients within the aquifer (i.e., the aquifer is freshwater-filled), and therefore, both Ra and M (and Ra δ ) are zero. Given that the system is forced-convection dominated, SFGD to the seafloor under the Model A conditions will be diffusive or uniformly distributed.
For the application of Equation 5 to estimate M for Model B, an average q z (across the bottom aquifer) over about 85 years of quasi-steady state conditions was used, based on the numerical results. The need to know q z in applying Equation 5 to estimate M (or Ra δ ), as a predictor of buoyancy-driven flow, is potentially problematic for the purposes of designing SFGD monitoring systems (e.g., deployment of seepage meters), because rates of SFGD are typically not known prior to seepage meter deployment. For example, Stevens , which requires the vertical hydraulic gradient (based on localized measurements of hydraulic heads and salinities from bore data) rather than q z . Nevertheless, where mixed convective processes occur in seafloor sediments, the vertical hydraulic gradient is also difficult to ascertain considering the changes in location and temporal variations of SFGD caused by buoyancy forces. Alternative approaches to estimating q z , such as the application of geochemical tracers (e.g., Taniguchi et al., 2019), may overcome difficulties in measuring hydrogeological variables within seafloor sediments.
The steady-state value of q z is 2.9 × 10 −3 m/d in the Base Case (Figure 7a), 2.9 × 10 −3 m/d in Case 1 (Figure 7b), and 2.3 × 10 −3 m/d in Case 2 (Figure 7c). While the values of M for the three cases of Model B vary within the same order of magnitude (i.e., 22-27), the range of values of Ra is wider (i.e., 2031, 431, and 264 for the Base Case, Case1, and Case 2, respectively). These values indicate that the flow system in Model B is free-convection dominated. That is, M > 1, and the critical values of Ra δ (or M) and Ra (i.e., 10 and 5, respectively) are also exceeded, indicating that unstable solute motion is likely to occur. Therefore, the occurrence of unstable flow in the current cases are consistent with critical values proposed in the literature for downwards salinization (e.g., Wooding et al., 1997;Smith and Turner, 2001). However, further investigation is warranted to ascertain if those critical values can be generally applied to SFGD through sandy seafloor sediments.

CONCLUSION
This research highlights that the occurrence of SFGD in permeable seafloor sediments potentially involves unstable flow processes, with important implications for SFGD measurement and the understanding of seafloor ecosystems. This study provides insight into SFGD measurement approaches, since we have demonstrated that for the cases considered here, SFGD measurements through the deployment of seepage meters will depend on the placement and measurement duration of seepage meters. Predicting whether unstable flow processes occur is theoretically plausible based on our overview of the most common non-dimensional numbers (i.e., Ra δ , Ra, and M) used previously to characterize mixed-convection processes in groundwater. Simplified numerical models that represent submarine aquifers settings show that SFGD can occur in the form of upward freshwater fingers, analogous to the inverse downward movement of dense solute fingers. We found that the critical values of Ra δ and Ra proposed in the literature for solute convection apply to the two cases analyzed in this study, although, further investigation is needed to generalize the use of those critical values to SFGD through sandy seafloor sediments. That is, further research into SFGD through high-permeability sediments is warranted to constrain the application of non-dimensional numbers for predicting unstable conditions, given the substantial differences between stable and unstable salinity patterns and distributions of SFGD to the seafloor. The results of this study show that the pattern of SFGD through high-permeability sediments containing seawater is controlled by the lower boundary, intended to represent either an underlying aquifer or aquitard. The numerical results also showed quasi-steady temporal fluctuations in the SFGD pattern behavior, for the case involving a low-permeability layer beneath the seafloor aquifer.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
SS-R: numerical modeling, investigation, analysis, visualization, and writing -original draft and editing. AW: conceptualization, supervision, visualization, and writing -review and editing. DI: supervision, visualization, and writing -review and editing. All authors contributed to the article and approved the submitted version.