Characterizing Liquid Water in Deep Martian Aquifers: A Seismo‐Electric Approach

Deep Martian aquifers harboring liquid water could hold vital insights for current and past habitability. We show that with seismo‐electric interface responses (IRs) we can quantitatively characterize subsurface water on Mars. Full‐waveform simulations and sensitivity analyses across diverse Martian aquifer scenarios demonstrate the technique's effectiveness. In contrast to how seismo‐electric signals often appear on Earth, Mars' desiccated surface naturally removes co‐seismic fields and exposes useful IRs that allow us to characterize several aquifer properties. Changing the aquifer depth, thickness, or quantity changes the IR arrival times or shape: aquifer depth is a strong control on evanescent IRs, thickness affects the relative timing of IRs, and increasing the number of aquifers introduces more dipole sources to the waveform. Other factors, such as aquifer saturation, chemistry, and salinity, strongly affect IR amplitude but have minimal or no effect on waveform shape. Notably, for a deep low‐porosity aquifer, the salinity and brine chemistry (perchlorate vs. chloride) are the strongest controls on signal amplitude. Analyzing the effects of epicentral distance shows that radiating and evanescent IRs separate at large source‐receiver offset, allowing analyses of both signals and accurate event distance derivation. From this numerical investigation of the sensitivity of IRs to deep Martian aquifers, we anticipate future analyses of electromagnetic data from the InSight lander or future missions to Mars and other planets.


Introduction
Water is essential to life as we know it.By understanding the location and distribution of liquid water on Mars, we can assess the potential for planet to both harbor past or present life and provide valuable resources for future human exploration.
The study of deep groundwater on Mars continues to be primarily theoretical, largely due to limited field observations.Past studies have sought to answer a brought range of general questions: what is the volume of deep subsurface water, both locally in aquifers (Head et al., 2003) and globally (Carr & Head, 2015)?What is the latitudinal extent and thickness of the cryosphere which could trap subsurface aquifers (Clifford et al., 2010)?How might the water be stored within the rock matrix (Kilburn et al., 2022)?What are the possible chemical and electrical properties of the groundwater (Burt & Knauth, 2003;Michalski et al., 2013;Rivera-Valentin et al., 2020)?How do all of these factors change through time in a Martian hydrologic system (Carr & Head, 2015;Clifford, 1993)?
Very few quantitative methods have been developed for deep Martian groundwater detection and characterization; those methods that have been proposed generally utilize electromagnetic (EM) or seismic sounding.EM sounding theory, sensitivity, resolution, and instrumentation have been thoroughly assessed for groundwater on Mars (Grimm, 2002;Grimm et al., 2009).EM sounding is highly sensitive to the large resistivity contrasts created by Martian groundwater, and the method can generally be used to derive subsurface conductivity and can inform on other EM parameters and aquifer porosity.For EM sounding using natural EM sources, the exact nature and of the source is important in determining the level of success it could provide in passive EM surveys (Grimm, 2002).Grimm et al. (2009) developed a model for an active-source time-domain EM sounder instrument that could detect Martian saline groundwater aquifers up to 5 km deep.Seismic data from the Interior Exploration using Seismic Investigations, Geodesy and Heat Transport (InSight) mission is actively being used to investigate the existence of groundwater and the cryosphere beneath the lander (Banerdt et al., 2020;Kilburn et al., 2022;Manga & Wright, 2021;Wright et al., 2022).Cryosphere thicknesses beneath the InSight lander may vary from 0 to 9 km, below which liquid water may be present (Clifford et al., 2010).Seismic S-(shear) wave velocities suggest no detectable shallow ice or fluid, but pore gas or liquid may exist 8-20 km deep (Kilburn et al., 2022).Observed Swave anisotropy could be created by a dry-or liquid-filled fracture system or igneous protrusions ∼8 km below InSight (Li et al., 2022).Seismic methods are useful for constraining mechanical differences in the subsurface, but these S-wave analyses alone cannot distinguish between pore-water and pore-gas (Kilburn et al., 2022;Li et al., 2022).
We present the seismo-electric (SE) phenomenon for deep Martian aquifer detection and characterization.As seismic waves traverse saturated porous media or interact with interfaces in saturated layers, they induce EM waves through pore-scale electrokinetic coupling (Frenkel, 2005;Pride, 1994).The resulting EM wavefields are termed SE fields.SE wave amplitudes are highly responsive to variations in poro-elastic properties (e.g., density, bulk modulus, pore size, and rock permeability) and electrochemical properties (e.g., electric permittivity contrasts between rock and water, salinity, electrolyte chemistry) (Dzieran et al., 2019;Gao & Hu, 2010;Jufar et al., 2017;Ma et al., 2022;Zyserman et al., 2017).SE fields are characterized by two main types of EM fields: co-seismic waves and interface responses (IRs) (Butler et al., 1996(Butler et al., , 2018;;Dzieran et al., 2019;Ren et al., 2016;Thompson & Gist, 1993;Zhu et al., 2000).Figure 1 is a cartoon of the SE effect on Mars.The co-seismic field, shown in Figure 1i, is an EM signal coupled into the seismic wavefield traveling at the same velocity and with amplitudes proportional to the ground acceleration (Garambois & Dietrich, 2001).Interface responses (IRs) are generated when a seismic wave creates a charge imbalance across a saturated interface (Haartsen & Pride, 1997).There are two main types of IRs: radiating IRs ("RIRs") radiate from the interface independently at EM velocities, regardless of fluid content in the second medium (Gao et al., 2017).These can be easily seen in Figure 1e: these RIRs are created as the deep marsquake impinges on the bottom and top of an aquifer.They can arrive at a receiver before any seismic waves and are characterized by a zero-moveout signal across large distances (Wang et al., 2021).The second type of IR is the evanescent IR ("EIR") first described by Ren et al. (2016).EIRs are generated when a seismic wave impinges on a saturated interface at an angle equal to or greater than the critical angle: where V seis , V EM are the seismic and EM velocities, respectively.EIRs diffuse into the second medium and decay exponentially in amplitude.Both IRs are created only in the presence of mobile fluids but do not require a saturated layer to propagate or diffuse, making them useful to find deep groundwater on Mars (Olhoeft, 2003).Further, the frequency-dependent attenuation of RIRs is low compared to their seismic counterparts, allowing them to travel through kilometers of rock without significant energy loss (Haartsen & Pride, 1997;Thompson & Gist, 1993).We introduce the naming convention of a "P-" or "S-" prefix for IRs produced by P-and S-waves, respectively.
SE methods have been used to study groundwater on Earth.Electric and magnetic components of SE fields from earthquakes have been recorded, and IRs from local seismicity have been proposed to monitor aquifers in Turkey (Dzieran et al., 2019(Dzieran et al., , 2021;;Gao et al., 2016;Kasdi et al., 2022;Romano et al., 2019).SE signals were used to map thin subglacial water strata at high resolution (Kulessa et al., 2006).Lab experiments by Zhu et al. (2008) and Zhu and Toksöz (2012) measured key SE coupling coefficients.Later, Peng et al. (2018) revealed that these SE coupling coefficients are approximately constant in the seismic frequency band.Reservoir-scale SE surveying with high-power electrical sources has imaged reservoirs to kilometer-scale depths (Fourie et al., 2000;Thompson et al., 2006).
On Mars, natural seismic sources primarily include marsquakes and surface meteor impacts (Clinton et al., 2021;Daubar et al., 2018;Dzieran et al., 2021;Giardini et al., 2020;Jacob et al., 2022;Kasdi et al., 2022;Kedar et al., 2021).Earthquakes have been shown to naturally excite SE waves (Dzieran et al., 2019(Dzieran et al., , 2020(Dzieran et al., , 2021;;Gao et al., 2016), so as an analogy marsquakes and impacts could act as natural seismic sources for the generation of SE signals on Mars.One of the biggest problems for SE studies on Earth is removing the co-seismic signal, seen clearly in the Earth-analog model of Figure 1, panels g, h, and i, where the co-seismic signal overpowers IRs by several orders of magnitude.On Mars, we would not expect such a co-seismic wave to be observed at the surface (seen in Figure 1f).The generation of the co-seismic wave is inherent only to seismic waves traveling through saturated media (like the near-surface on Earth), whereas the Martian near-surface media are desiccated and do not support the generation of co-seismic waves.
The remaining sections of this paper demonstrate the power of the SE method to investigate deep groundwater beneath the surface of Mars.This method has the potential to: 1. Unambiguously identify the existence of kilometers-deep liquid water aquifers, 2. Constrain aquifer volume through pore saturation and aquifer thickness, 3.And characterize aquifer depth, poroelastic parameters, pore-fluid chemistry, and the state of chemical equilibrium and number of aquifers.
We carry out simulations of various aquifer saturations, depths, thicknesses, chemistries, geometries, and source distances.For all of these models, we track relative change in the amplitude and phase of various components of the SE wavefields, providing an analysis of the wavefield sensitivity.This sensitivity analysis helps quantify the contribution of each parameter to the overall SE amplitude.

Mars Geologic Model
The simulations in this paper calculate wavefields in a combined solid-porous layer model, illustrated in Figure 2.
The surface of the model is simplified from recent interpretations from satellite studies and the InSight lander.Shallow lithology was investigated and interpreted via crater wall studies and satellite imagery (Golombek et al., 2018), retro-rocket upheaval of the upper regolith (Golombek et al., 2020), compliance measurements (Kenda et al., 2020), and ambient noise (Hobiger et al., 2021;Wright et al., 2022).In the top few hundred meters, highly fractured crater ejecta overlays shallow layers of Amazonian and late-Hesperian basalts which themselves overlay older Noachian sedimentary layers.We simplify this near-surface geology to three layers: regolith, a basalt flowpath, and upper fractured sediments.Beyond these layers, lithology is speculative and based on velocity measurements and interpretations from marsquake-induced surface waves (Kim et al., 2022) and body waves (Knapmeyer-Endrun et al., 2021).We represent these results as increasingly competent Noachian sediments down to 8 km.
Below this, what Kilburn et al. (2022) terms the "deep crust," S-wave velocity measurements are consistent with consolidated basalts or more fractured feldspar-rich rocks with porespace filled with gas or fluid (Kilburn et al., 2022;Kim et al., 2022).Our model represents this as the bottom bedrock halfspace.Thermodynamic models most favor the stability of liquid water within this deep crustal layer (Clifford et al., 2010).A deep crustal aquifer could extend in depth down to around 20 km, where Mars' gravity likely closes all porespace (Kilburn et al., 2022).As a general approximation, the porosity at the surface is 20% and decreases with depth with a scale height of 2.82 km (Clifford, 1993).

Seismo-Electric Governing Equations and Numerical Implementation
Seismic waves moving through saturated porous media are commonly described by the Biot equations (Biot, 1962).These equations treat the combined solid-fluid medium as a complex system with properties that depend on the pore-filling fluid.Pride (1994) demonstrated that when these equations are coupled with Maxwell's equations of electromagnetism, they can explain the pore-scale SE phenomenon and numerically solve for SE wavefields in horizontally layered media (Gao & Hu, 2009, 2010;Haartsen & Pride, 1997;Pride & Haartsen, 1996;Ren et al., 2012Ren et al., , 2015)).
This work utilizes the method originally presented in Haartsen and Pride (1997) to simulate the seismic and EM wavefields generated by marsquakes in the presence of deep aquifers.Assuming an e ωt time dependence, the governing equations can be stated: These first four equations are related to Biot's formulations to determine the seismic wavefield, where τ is the bulk stress tensor, ω is radial frequency, and ρ, ρ f are the bulk and fluid densities, respectively.The field variable u is the average solid displacement and is counterpart to w, the relative solid-fluid phase displacement.F, f, are applied force densities acting on the bulk material and fluid phases, respectively.Gassmann's bulk modulus, K G , is joined by C, M, which are coefficients related to the bulk moduli of the fluid and solid phases.The stiffness G is a shear modulus related to rock matrix grains.I is the identity matrix, P is pore-fluid pressure, k is the wavenumber, η is the pore-fluid viscosity, L is the electrokinetic coupling coefficient which heavily modulates SE converted amplitude (Jaafar et al., 2011;Pride, 1994;Revil, 1999), and E is the electric field.
Equations 5-9 compute the EM wavefield from Maxwell's equations, and the transport Equations 4 and 5 govern the coupling between the seismic and EM fields.In these equations, H is the magnetic field, B is the magnetic flux density, J is the electric current density, D is the dielectric displacement, σ is the bulk conductivity, ϵ is the electric permittivity, and μ is the magnetic permeability.Two source terms also appear: M and C are applied magnetic-current and electric current-density source terms, respectively.
The physical importance of pore-fluid on the results of these equations is clear by the presence of the numerous fluid-related parameters.This fundamentally tells us that a pore-filling electrolyte fluid is necessary for this phenomenon to occur.Thus, the observation of an unambiguous SE signal is diagnostic of mobile fluids.Although the SE phenomenon has never yet been studied in relation to Mars, current understanding of Martian geology give us no reason to believe the deep pore structure of the Martian crust would behave any differently to mechanical stresses than the porous media here on Earth.We move forward under the assumption that these methods and the underlying theory are just as applicable on Mars as they are on Earth.
The treatment of the uppermost desiccated layers of Mars requires an extension of SE theory to account for partially saturated and un-saturated media.Saturation here refers to the porespace fraction occupied by mobile fluid and is notated as S w (Bordes et al., 2014).We adopt the linear model of Guichet et al. (2003) that introduces a saturation functional factor, f(S w ) = S w , to the electrokinetic coupling coefficient L. The static saturationdependent electrokinetic coupling coefficient can be written: where ϕ is the porosity, α ∞ is the tortuosity, σ w is the fluid conductivity, ϵ 0 is the permittivity of free space, ϵ f is the fluid permittivity, ζ is the zeta potential (which is dependent on fluid molarity and chemistry), Λ is a pore geometry factor, and the length scale d is less than or equal to the Debye length.The second exponent of Archie's Law, n, is kept constant at 2 for all models, as is appropriate for generalized or unknown grain structures (Glover, 2017).
The saturation fraction affects the effective fluid modulus, the poroelastic moduli B, C, M and H, the bulk density and effective fluid density, the effective fluid viscosity, the effective electrical conductivity of the medium, the effective electrical permittivity, the frequency dependent dynamic permeability, the complex density, the electrokinetic coupling coefficient, and the P-and S-wave slownesses.For the wavefield modeling results presented below, the most noticeable difference that saturation causes is a change in P-wave velocity.The full equations of the listed parameters' saturation dependence can be found in Bordes et al. (2014), Appendix A.

Results
We simulate Martian aquifers of various saturations, depths, thicknesses, and chemistries under the excitation of marsquakes.We also analyze some model results at epicentral distances up to 100 km and one model implementing two separate confined aquifers.In scenarios where the placement or extent of the aquifer changes, the saturation-dependent parameters mentioned earlier also change in the newly saturated layers.The default model parameters and scenario-dependent changes that are not shown here are all listed in Supporting Information S1, Tables 1-4.Unless otherwise specified, the default configuration of the model are all follows: • Saturation of all non-aquifer layers are set to near-zero (0.01%-0.0001%) • Marsquake source is 20 km deep at x = y = 0 m.
• Receiver line is set along the coordinates {x = 0 km …20 km, y = 0} with a total of 200 receivers at a spacing of 100 m.Radial and transverse seismic and EM wavefields are measured.The large array allows easy visualization of SE signals, though the bulk of our analyses focus on results from one single receiver.• The marsquake is represented by an x-z double-couple rupture source with a moment of 1e14 N m.This seismic moment correlates to an earthquake with a moment magnitude of 3.3 (Hanks & Kanamori, 1979).Its source function is a ricker wavelet with a central frequency of 1 Hz and a central time delay of 1 s.Double-couple sources are frequently used in earthquake modeling to represent fault-slip behavior (Gao & Hu, 2010).
Whenever the maximum IR amplitude is referenced, this is calculated as the maximum of the absolute value of the trace in question.There is no specific time window used, though the maximum value reliably comes from the S-RIR or P-EIR generated at the aquifer interface nearest to the surface.

Saturation
The aquifer saturation was modeled at 100%, 80%, 50%, 10%, and 1%. Figure 3 compares selected traces for aquifers saturated 100%, 50%, and 10%.Before we compare the results of each saturation model, it is worthwhile to describe the visible arrivals in each component.In Figure 3a, the vertical component of the particle velocity, V z , for each saturation model is shown.Due to the x-z orientation of the double-couple source, the vertical component is dominated by P-wave energy.Figure 3b shows the radial particle velocity, V r , which has strong SHwaves.Figures 3c and 3d show the radial electric (E r ) and transverse magnetic (H t ) components, respectively.The first signal observed in both components is the RIR created by the marsquake as its SH-wave impinges on the bottom of the finite aquifer.The marsquake SH-wave (velocity 3,000 m s 1 ) would be expected to reach the bottom of the aquifer 2 km above by 0.66 s.As the seismic wave continued upwards, it would impinge on the top of the aquifer at 4 s.This aligns with the strongest RIR.P-RIRs are present just preceding the S-RIRs, though they are much smaller and difficult to see.
We call IRs that arrive after the largest RIR as "coda IRs," which begin around 7 s in Figures 3c and 3d.Most large coda IRs can be traced back through travel time analysis to downward-propagating SH and P-S seismic reflections that interact with the deep aquifer.S-RIRs are much larger than the P-RIRs over all of our tested scenarios.Thus we conclude that the SE coupling for a double-couple source is dominated by the SHTE (SH and horizontally polarized electromagnetic) mode.
The percent change in the maximum IR amplitudes of each scenario relative to a 100% saturated aquifer are shown in Figure 3e.The amplitudes depend strongly on the values reported in Table 1.Changing the pore fluid saturation fraction changes primarily the effective fluid viscosity and bulk conductivity, both of which also change the electrokinetic coupling coefficient.The electrical permittivity is also changed to a smaller degree.

Aquifer Chemistry
Pore fluid chemistry on Mars is not well understood.We test three different solutes that have been theorized to exist on Mars: calcium perchlorate (Ca (ClO 4 ) 2 ) (Primm et al., 2017), magnesium perchlorate (Mg(ClO 4 ) 2 ) (Rivera-Valentin et al., 2020), and sodium chloride (NaCl) (Thomas et al., 2019).Each of these brines are modeled at 5% concentration, and NaCl and Ca (ClO 4 ) 2 are also modeled at 1% concentration to test salinity sensitivity.Finally, an NaCl gradient is tested.Directly beneath the cryosphere, liquid water might form a salinity gradient as water freezes out of the brine solution into the base of the cryosphere.Denser hyper-saline brine would descend  deeper into the aquifer, replaced by comparatively dilute brine at the top (Clifford, 1993).The choice of the gradient severity in our model is arbitrary, as we know of no existing theories on these characteristics.Beginning at the 8 km deep aquifer top, salinity increases in 85 equal steps.Each step is a discrete ∼118 m layer.The top of the aquifer starts at 1% salinity and increases linearly to 20%.These layers have no elastic differences from the default aquifer.We initially tested gradients with between 2 and 85 discrete layers and found that the higher end, 85 layers, was necessary to more accurately represented a physical gradient.Brief results from the model tests varying gradient granularity are shown in Figure S1 in Supporting Information S1.
Pore fluid composition plays a large role in the SE conversion through the electrokinetic coupling coefficient calculated with Equation 11, which is a function of brine chemistry through several factors.First, the length scale d ≤ d, where the Debye length d is: where N i is the ionic concentration in ions per meter cubed of ion type i, z i is the electrolyte valence of ion type i, e is the electron charge equal to 1.60 × 10 19 C, and kT is Boltzmann's constant k (1.38 × 10 23 J C 1 ) multiplied by the temperature T. For a z:z binary symmetric electrolyte (NaCl) the summation in the denominator simplifies to 2Ne 2 z 2 , as used in Haartsen and Pride (1997).However, perchlorate electrolytes are not symmetric and require the separate consideration of their cation and anion species in the summation.
The fluid conductivity (neglecting surface conduction) σ 0 and pure electrolyte conductivity σ f are where b ± = 3 × 10 1 1 m s 1 N 1 are the ionic mobilities of the cations and anions.
The zeta potential ζ in this equation has been formulated in the literature in two broad domains.The first is an approximation often used for pore fluids with a cation molarity of less than 0.5 M: where C is the electrolyte molarity.
The second formulation of the zeta potential is far more chemically involved, but allows the modeling of pore fluid concentration greater than 0.5 M.Many of the empirical chemical parameters for this formulation are not known for the perchlorates used here, so this is only used in our models that use a NaCl concentration greater than 0.5 M. The entirety of the second formulation is described in Glover et al. (2012), though in brief the zeta potential is approximated by: where ϕ d is the Stern plane potential: where N a is Avogadro's constant, K Na = log(7.5)is the binding constant for sodium on quartz, K ( ) = -log(7) is the dissociation constant for dehydrogenization of silonol surface sites, Γ o s =10 nm 3 is surface site density, C f is the variable fluid molarity (in mol L 1 ), and the pH = 7.These values are appropriate for sodium chloride brine chemistry and are used by Glover et al. (2012) and Gao et al. (2014).The shear plane distance χ ζ is typically equal to 2.4 × 10 10 m.The basic numerical parameters dependent on the fluid chemistry and salinity are summarized in Table 2 for each of the chemical scenarios.For the non-symmetric brines in Table 2, anion molarity is double cation molarity and anion valence is 1.Brine chemistry does not affect the rock elastic properties in our simulations, as seen in Equations 12-17, thus we show only E r and H t in Figure 4. EM traces in panels a-b show that the NaCl gradient is the only trace with changes to the RIR wavetrain shape.The gradient traces lack S-RIRs generated at the base of the aquifer, likely due to the higher saline concentrations there effectively dissipating the bottom S-RIR as it moves through the gradient.All other models varied only in amplitude.
Figure 4c shows the maximum amplitude is very sensitive to chemistry and particularly the existence of a chemical gradient.The EM-relevant parameters are shown in Table 3. Changing the pore fluid chemistry affects the zeta potential, effective conductivity, electrokinetic coupling coefficient, and fluid conductivity.Each of these parameters are larger for the 5% NaCl case compared to the default 5% Ca(ClO 4 ) 2 case.This leads to a larger SE signal generated at the interface, but the higher conductivity also means that the RIR generated at the lower aquifer interface attenuates faster as it passes upwards through the aquifer.
Between the three 5% salinity models, the sodium chloride showed the largest overall amplitude.This is due to the increase in the aquifer electrokinetic coupling coefficient, reported in Table 3, where 5% NaCl has the highest value of all three brines tested at 5% concentration.The effect of decreasing salinity can be seen by comparing the 5% and 1% calcium perchlorate cases and the 5% and 1% NaCl cases in Figure 4e.Consistent with prior understanding, the amplitude of the SE response varies inversely with increasing salinity (Gao & Hu, 2010).The difference between 5% magnesium perchlorate and 5% calcium perchlorate is very small due to their similar molarities and densities, which leads to negligible differences in EM coefficients.There is a much larger amplitude contrast between the perchlorates and the sodium chloride.
The amplitude of the gradient concentration aquifer should most aptly be compared to the 1% constant concentration NaCl aquifer, as the gradient is topped with 1% NaCl brine.Figure 4e shows a slight dampening in the maximum amplitude of the gradient model, suggesting that a linear saline gradient created by overlying cryospheric ice may lead to a less-detectable aquifer than one with constant concentration.

Depth
Three aquifer depths are exhibited here: 8, 5, and 2 km, visualized in Figure 5a.These depths were chosen as reasonable estimates based on cryosphere thicknesses across various Martian latitudes estimated by Clifford et al. (2010).In each case, the aquifer saturation extends from the noted depth downward to 18 km depth.Figure 5 gives an overview of the velocity, electric, and magnetic wavefields compared between the three depth cases.The electric field (Figures 5c, 5e, 5g) has a stronger P-EIR signal than the magnetic field.The EIRs grow significantly stronger at shallower depth, seen in panel g, dominating E r in the shallowest case.Single receiver traces from directly above the marsquake source are shown in Figure 6.
The SE traces in Figures 6b and 6c show increasing time delay as the top of the aquifer is farther from the marsquake source.Only the 8 km depth case has a visible S-RIR from the aquifer bottom, as the thicker aquifers of the other models heavily attenuate this RIR.The coda IRs are more compressed in shallower aquifer models due to the shorter time required for downward seismic reflections to reach the aquifer top and generate coda IRs. Figure 6d  shows that depth is a significant control on IR max amplitude.The electric and magnetic fields increase with different trends as the aquifer gets shallower: we attribute this to the stronger growth of EIRs in the electric field.

Aquifer Thickness
We compare five thinner aquifer scenarios to the default model of 10 km thickness to investigate IR sensitivity: 5 km, 2 km, 1 km, 750 m, and 100 m.For all thickness models, the aquifer top remains at 8 km depth and the bottom aquifer interface is moved upwards to produce the desired thickness.Figure 7 shows timing changes in the IR wavetrains for these thickness models.The first thinner aquifers in Figure 7c shows the S-RIR induced at the aquifer bottom appearing closer to the S-RIR generated at the top.This trend continues in Figures 7d and 7e.For these thickness models, the RIR produced from the bottom interface is generated only fractions of a second before the RIR from the top, thus there is not two distinct IR arrivals for our central frequency of 1 Hz.Later RIRs are also affected, as seen most easily in Figures 7b and 7c between 6 and 9 s.
Figure 8 shows that the maximum IR amplitudes of the models only change when the aquifer becomes thin enough to cause some overlap between the RIRs produced at its two interfaces.Thus, we can attribute the amplitude changes to constructive and destructive interference between the multiple SE sources.

Source Distance
While all our previous models have simulated a receiver array of 200 receivers over the span of 20 km, we have not focused on the effects of epicentral distance on SE waveform and amplitude.In this section, we compare waveforms and amplitudes from receivers at distances of 0, 25, 50, 75, and 100 km.In this scenario, we use 200 receivers set 500 m apart over 100 km Figure 9 shows that the EM components decrease quickly away from the source, with the RIRs in E r decreasing slightly faster than those in H t .
As can be seen in Figures 9b and 9c, EIRs play a larger role in the electric field wavetrain.Locally generated EIRs should increase in prevalence beyond the critical radius as defined by Ren et al. (2016), r c = h ⋅ sin(Θ crit ), where r c is the radius of the circular concentration area on the interface where RIRs would be generated, h is the source depth measured from the generating interface, and Θ crit is the SE critical angle.Because the EM velocity is frequency dependent, the SE critical angle and thus the critical radius varies with frequency.These values also differ between incident P-and S-waves, as they have different velocities.Table 4 records the ranges of these values for our default model.4 show that, for the 500 m receiver spacing used in this model, EIRs should be present primarily beyond one or two receivers.The RIRs are generated only when the incident seismic wave is essentially normal to the top of the aquifer.

The critical radii in Table
Figure 10 compares receiver traces from the sampled distances along the array.Even at 100 km distance, the RIRs generated directly above the distant source are the most powerful SE signal.These same RIRs are present in all traces.EIRs, on the other hand, are generated locally beneath each receiver by the traveling seismic wave.They decay exponentially away from the top of the aquifer, with higher frequencies decaying faster.Comparing the Stransforms between 0 and 1 Hz of the receiver directly above the source and the receiver at 100 km in Figure 11 reveals the presence of low-frequency EIRs at the distant receiver between 20 and 60 s.
Because both RIRs and EIRs are generated at different locations, detecting both IRs with receivers at large epicentral distances could allow for the twofold inversion of the local subsurface as well as the subsurface directly above the seismic source.The trend in the IR amplitudes as a function of epicentral distance, shown in Figure 10e, appears as expected, with a general geometric spreading decay seen in both EM components.

Two Aquifers
The choice of a single aquifer used in the rest of this paper has been arbitrary, allowing us to focus on fundamental effects of various parameters.In reality,  Journal of Geophysical Research: Planets 10.1029/2024JE008292 multiple layers of impermeable rock acting as aquifer seals can exist in a single setting, creating several stacked confined aquifers.So, we also briefly demonstrate the effects of introducing a second aquifer layer to the simulation.Taking the default model, we replace the top two km of the aquifer (depths of 8-10 km) with a 1 km thick aquifer 1 km above the top of the existing aquifer (placed between 10 and 18 km depth), shown in Figure 12a.Figures 12b-12d summarizes the results.
As seen in Figures 12b, 12d, 12e, adding a second aquifer introduces additional SE sources that overlap, effectively widening the RIR bunch that appears between 3 and 6 s.Similar to the thickness models, the position and thickness of the second aquifer layer can be expected to have a frequencydependent effect on the observed signal caused by the interference of multiple RIRs radiating from the aquifer interfaces.These additional SE sources end up creating a unique signal that can be used to distinguish configurations of multiple-aquifers from single aquifers.E r experienced a 330% increase in its maximum IR amplitude, while H t had a 97% increase compared to the default, singe aquifer model.These increases are certainly due to the addition of more SE sources generated at the bottom and top of the new aquifer.

Discussion
For an x-z double-couple source and receivers along the x-axis, the SHTE mode of the SE conversions is much stronger than the PSVTM mode.This carries the implication that, while S-waves alone cannot uniquely characterize pore fluid with a traditional seismic approach, they are very useful for SE generation.The results also showed that saturation does not severely modulate the amplitude of the signal for a low-porosity deep aquifer unless the saturation fraction is very low, and it does not affect the SE waveform shape or arrival times from a double-couple source.The SE amplitude peaked for an 80% saturated aquifer, with an increase of less than 15% in E r and 8% in H t relative to the fully saturated aquifer.This is a small sensitivity compared to the large effects of other parameters, most likely due to the low porosity of our aquifers.This insinuates that the ability to characterize aquifer saturation may be difficult, but when characterizing many variables, such as depth, thickness, and chemistry, it could be safe to assume a fully saturated aquifer.The low-saturation case can also inform on our understanding of how these results might be affected by thin films of liquid water coating mineral surfaces in generally ice-saturated cryospheres that might overly some of the groundwater on Mars (Clifford et al., 2010;Grimm, 2002).While our models did not simulate an ice-filled cryosphere, our results of ultra-low aquifer saturations show that the electrokinetic coupling coefficients of such layers are orders of magnitude beneath even marginally more saturated media.So, theorized semi-liquid layers at the base of icy cryospheres would likely not contribute meaningfully to an SE signal.
Brine composition is a vital factor in determining the aquifer's astrobiological significance and habitability.Our models showed that both the electrolyte species and the brine concentration modulate the SE signal amplitude, and not the waveform shape.Thus, the SE method is not suitable for distinguishing between these two aspects of the pore fluid composition, but it can constrain estimates of the combinations of these two variables.The large amplitude difference between perchlorate and chloride brines in Figure 4e is promising for differentiate between the two electrolyte structures (symmetric vs. asymmetric), a prominent question in the study of water on Mars (Primm et al., 2017).The ability to discriminate between solutes with the same electrolyte structures seems unlikely.Decreasing the calcium perchlorate salinity from 5% to 1% led to an amplitude shift in H t of the same magnitude as moving an 8 km deep aquifer toward the surface by 3 km.This shows that salinity is likely one of the strongest controls on the amplitude of deep Martian aquifers.Our arbitrary NaCl saline gradient slightly dampened IR amplitude when compared to the aquifer with a constant composition of 1% NaCl.The gradient aquifer also showed a distinct lack of S-RIRs generated at the bottom of the aquifer, which demonstrates that the SE method can distinguish this configuration from an aquifer of constant salinity.Wavefields in all panels are normalized to the same color scale, so the colorbar is unitless.
Altering the depth of the aquifer top changes both the IR arrival times and their amplitude.The arrival times of the IRs are dependent on the relative distance between the aquifer top and any other geologic layers that reflect seismic rays, conceivably allowing the inversion of whole subsurface geology through a method analogous to seismic travel-time inversion.Even results from a single receiver could be used to infer a 1D geologic model.A network of receivers, as shown in Wang et al. (2021), can allow 3D tomographic modeling of non-flat interfaces.Depending on the depth of the aquifer and the properties of the subsurface, EIRs can play a prominent role in the SE wavefield (shown in Figure 5g).EIRs are believed to be the main contributor to IR wavefields on Earth because of the prevalence of shallow fluids, and schemes using them for subsurface inversion have already been developed (Dzieran et al., 2019(Dzieran et al., , 2020)).In contrast, our models suggest RIRs will be the dominant SE signal for deep Martian aquifers.
Along with aquifer saturation, analyzing aquifer thickness can tell us about its fluid volume.Altering the thickness of the aquifer subtly changes the relative arrival times of the bottom and top RIRs and can cause them to overlap.This is showcased in the H t amplitude progression in Figure 8: the bottom RIR is first separate, then constructively interferes for a 1 km thick aquifer, then destructively interferes for a 100 m thick aquifer.The amplitude's thickness dependency should also be wavelength dependent, so the exact sensitivity quantifications given here should only be taken for our source characteristics.
As was experienced by the InSight lander, marsquakes can be detected from great distances (Banerdt et al., 2020).We showed that increased epicentral distance separates the RIR generated at the source and the EIR generated beneath the receiver.This raises the possibility of both local subsurface and near-source subsurface investigation, as the RIRs carry information about the subsurface stratigraphy and aquifer at the source location and EIRs carry information only from beneath the receiver.
The addition of a second, shallower confined aquifer added several RIRs to the SE wavetrain.Because of our layer geometry, these RIRs constructively interfered and increased the maximum amplitude of the SE signal in both components.The visual arrival of the aquifer top RIRs between 3 and 6 s of the simulation (shown in Figure 12a) is elongated due to the addition of multiple new sources, differentiating these results from other models.Thus, we show that stratified or multi-layer aquifers create visually distinct signatures.Learning about the specific aquifer geometry can hint at the region's geologic history and tell us more about the aquifer's hydrodynamics.
The RIR generated at the bottom of our finite aquifers can shed light into how the SE waves are modified by their passage through saturated aquifer layer.RIRs, like any EM field, attenuate in conducting media in a manner dependent on layer conductivity, permittivity, permeability, and frequency.So, while the RIR at the bottom of the aquifer is closer to the seismic source than the RIR at the top of the finite aquifer, the 10 km of conducting saturated rock attenuates the bottom RIR significantly.By the time both RIRs reach the surface, the bottom RIR is much weaker than its upper counterpart which was not heavily attenuated by the desiccated crust above the aquifer.The variations in conductivity and other important EM parameters for the saturation scenarios (Table 1) and chemistry scenarios (Table 3) give quantitative insight into how different aquifer configurations attenuate the bottom-interface RIR that propagates through them.
On Earth, practical use of the SE method is limited by the ability to separate IRs and co-seismic waves.On Mars, where the near-surface is certainly desiccated, no such separation is needed.The co-seismic response is entirely absent, as shown in all the simulations of this paper, allowing the IR to be investigated alone.While co-seismic responses are often orders of magnitude larger than RIRs (Haartsen & Pride, 1997;Kulessa et al., 2006;Ren et al., 2012), Ren et al. (2016) showed that EIRs can themselves be order of magnitude larger than co-seismic fields depending on measurement proximity to the generating interface.The only model in this paper that showed EIRs stronger than RIRs was the 2 km deep aquifer top case.On Earth, SE signals from earthquakes have been inverted for fluid parameters (Guan et al., 2013;Jardani et al., 2010;Ma et al., 2022) and used in imaging kilometers-deep reservoirs (Fourie et al., 2000;Thompson et al., 2006).Inversion of earthquake-based SE data has been done with various algorithms that use data from either a single receiver or an array as input (Dzieran et al., 2020;Guan et al., 2013;Jardani et al., 2010;Ma et al., 2022).These studies inverted for: fluid salinity (Dzieran et al., 2020;Jardani et al., 2010;Ma et al., 2022) and subsurface permeability (Guan et al., 2013;Jardani et al., 2010), as well as porosity, conductivity, and material moduli (Jardani et al., 2010).
A full analysis of the detectability of these signals is not undertaken here, as such a study should account for instrument sensitivities, receiver locations, specific source locations and mechanisms, and natural EM noise sources.The amplitudes reported here can be varied significantly by parameters that are poorly constrained for the subsurface of Mars, for example, porosity.Kilburn et al. (2022) estimated that the porosity of the structure around 8 km below InSight could range from 1% to over 20%.Many of our parameter choices, for example, the rock permeability, were chosen to show a conservative amplitude estimate.The permeability of the Martian subsurface is not well understood (Clifford & Parker, 2001), and our used value of 1e 13 m 2 is likely high for the deep subsurface layers.Because permeability varies inversely with the magnitude of the SE transfer functions reported in the work of Bordes et al. (2014), our use of a high permeability gives a conservative estimate for SE amplitudes in these regards.The detection of SE signals from earthquakes is traditionally done by large electric dipole equipment as part of magnetotelluric stations, with dipole lengths on the order of tens or hundreds of meters (Dzieran et al., 2021;Gao et al., 2016).On Earth, magnetotelluric equipment can detect electric signals on the order of 1 μV m 1 and magnetic signals on the order of 1 pT.The InSight lander's IFG magnetometer has a similar magnetic field noise floor of ∼30 pT at 2 Hz.In our models, E r and H t fell between the orders of 10 6 10 4 μV m 1 and 100 1,000 pT at receiver zero.Thus, for this work, the transverse magnetic component is most reasonably detectable with signal-to-noise ratios well over 10 for the IFG magnetometer.
A large magnetotelluric station has never yet been deployed on the surface of another body, but may be possible to create in the context of a distributed seismic network such as the Cerberus Fossae network proposed by Stähler et al. (2021).Cerberus Fossae is the site of geologically recent near-surface flooding possibly caused by deep subsurface aquifers cracked by intrusive dike formations (Burr et al., 2002;Head et al., 2003).Cerberus Fossae has been shown by the InSight lander to be extremely seismically active, which would provide bountiful seismic sources for SE generation below or within deep local aquifers (Banerdt et al., 2020;Clinton et al., 2021;Jacob et al., 2022;Kim et al., 2022;Stähler et al., 2021).Other intriguing theories on Martian seismicity could be tested with the SE method, such as one recently posed by Manga et al. (2019).They investigate a relationship between tidally increased pore pressure in deep aquifers and a prevalence of local seismicity.This could lead aquifers to create repeatable seismic sources.
The SE phenomenon may have applications beyond Mars as planetary research directives increasingly look toward the outer solar system.Probing the chemistry of icy world subsurface oceans or crustal brine pockets (Chivers et al., 2021) are possibly within the capabilities of this method.

Figure 1 .
Figure 1.Cartoon of the SE effect produced by the excitation of a deep marsquake in three aquifer configurations.(a-c) the seismic ground velocity and electric field measurements that might be expected for a completely desiccated Mars where there is no deep aquifer.(d-f) the result in the case of a deep aquifer.A P-EIR can be seen originating at the same time as the P-RIR and increasing in moveout at large epicentral distance.The S-EIR is much weaker and is not visible here.(g-i) a fully saturated model, analogous to what we observe on Earth.The P-wave arrival times are variable between these examples due to changes in layer saturation.

Figure 3 .
Figure 3. (a, b) V z and V r for the saturation cases.(c, d) E x and H t , with coda regions and visual amplifications labeled.All other EM components are much smaller or zero.The initial S-RIRs are labeled in the insets.(e) Relative percent change in the maximum IR amplitude of the SE wavefield components of receiver zero for the saturation models relative to the default scenario of 100% saturation.Right axes track the actual amplitudes of E r and H t .

Figure 4 .
Figure 4. (a, b) E r and H t field results.Examples of S-RIRs and P-RIRs are given and any visual amplifications are noted.The earlier S-RIR is produced at the aquifer base.(c) Relative percent change in the maximum IR amplitude of the SE wavefield components of receiver zero for the varied brine chemistry models, relative to the default scenario of 5% calcium perchlorate brine.Right axes track the actual amplitudes of E r and H t .

Figure 5 .
Figure 5. (a) Stratigraphic column demonstrating the locations of the aquifer tops in the various depth models tested.The depths are based on the range of theoretical aquifer depths theorized by Clifford et al. (2010).(b-d) Results for the 8 km deep aquifer.(b) V r shows strong S-waves and preceding P-waves.(c, d) E r and H t , respectively.(e, f) EM results from the 5 km case.Changes in travel times are due to the modified aquifer geometry.Stronger EIRs can be seen, even appearing in the magnetic field.(g, h) Results for the 2 km deep case.(g) E r with powerful EIRs.(h) H t .While the EIRs show significantly more, the RIRs still dominate.Wavefields are normalized to the same color scale, so the colorbar is unitless.

Figure 6 .
Figure 6.(a) V r for the depth cases.(b, c) E r and H t , respectively, with notable arrivals and coda marked.The S-RIRs from the aquifer top are circled.(d) Percent change in the maximum IR amplitude (the circled S-RIR) of the SE wavefield components of receiver zero for the varied depth models, relative to the default scenario at 8 km deep.Right axes track the actual amplitudes of E r and H t .

Figure 7 .
Figure 7. (a) V r in the 10 km thickness model.(b-e) H t for 10, 5, 1, and 100 m thick models.The early RIR from the bottom aquifer interface can be seen moving closer to the RIR generated at the top of the aquifer interface, until the two RIRs overlap in the 1 and 100 m thick cases.

Figure 8 .
Figure 8. Relative percent change in the maximum IR amplitude of the SE wavefield components of receiver zero for the varied thickness models, relative to the default scenario of a 10 km thick aquifer.Right axes track the actual amplitudes of E r and H t .

Figure 9 .
Figure 9. Results for 200 receivers spread over a 100 km line.(a) V r (b) E r .RIRs appear to maintain their amplitude better over larger distances compared to EIRs.(c) H t , showing similar patterns to the electric field.Wavefields in all panels are normalized to the same color scale, so the colorbar is unitless.

Figure 10 .
Figure 10.(a) V r at several distances.Solid and dotted lines trace the P-wave and S-wave arrivals, respectively.These traces are normalized to focus on waveform shape (b) H t , also normalized.Initial RIRs are the highest amplitude portion of the signal.EIRs are difficult to identify visually in these traces.Light pink dotted lines show the P-EIR generation times.The inset for (b) shows the first seconds of the normalized H t .The waveform shapes are almost identical.(c) Relative percent change in the maximum IR amplitude of the SE wavefield components of receiver zero, relative to the default scenario at x = 0 km away from the source.Right axes track the actual amplitudes of E r and H t .

Figure 11 .
Figure 11.(a) S-transform of E r at x = 0 km.Little low-frequency content is present.(b) S-transform of E r at x = 100 km.Here, low-frequency energy is present at the expected arrival time of the P and SH waves at the top of the 8 km deep aquifer below the receiver.This is diagnostic of the locally generated EIR.The expected arrival times for the P-EIR and S-EIR are marked, respectively, by pink solid and dashed lines.

Figure 12 .
Figure 12.(a) Stratigraphic column demonstrating the location of the second aquifer implemented in this section.H t across entire receiver array.The wide RIR signal arrival between 3 and 6 s is comprised of RIRs from the top of the lower aquifer and both interfaces of the shallower, newly placed aquifer.(c) V r from receiver at x = 0. (d, e) E r and H t , respectively.The addition of the new SE sources from the second aquifer is evident in the more complex RIRs.
SE measurements give us a way to detect and image Martian groundwater kilometers below the surface.We discussed the general characteristics and uses of seismo-electricity on Mars, as well as providing amplitude sensitivity analyses.The desiccation of the upper crust of Mars removes the co-seismic wavefield from the expected SE signals, which leads to a completely different SE wavefield than what we see on Earth.The SE effect Journal of Geophysical Research: Planets 10.1029/2024JE008292 could be a singular experimental tool capable of answering many of the questions driving research into deep Martian aquifers.Beyond unambiguous detection of mobile fluids in deep aquifers, the SE method offers expansive characterization of many unknowns: aquifer depth and geometry, water volume, source distance, and brine chemistry and gradient.If a long-lasting SE station is deployed on the surface of Mars, repeated SE signals over the course of months or years could also allow the study of how each of these factors change under the effects of a broader Martian hydrologic cycle.As SE exploration become more widespread on Earth, this study represents the first foray of the method to other worlds.

Table 1 (
Columns Left to Right) Effective Relative Permittivity, Viscosity, Zeta Potential, Effective Conductivity, Static Electrokinetic Coupling Coefficient, and Fluid Conductivity for the Three Different Modified Saturation Aquifers • Aquifer thickness is 10 km, from 8 km depth to 18 km depth.

Table 3
Effective Relative Permittivity, Viscosity, Zeta Potential, Effective  Conductivity, Electrokinetic Coupling Coefficient, and Fluid Conductivity,  Respectively,for the Five Different Modified Chemistry Aquifers Two higher-salinity NaCl gradient layers are showcased at the bottom.