Liquid water in the Martian mid-crust

Large volumes of liquid water transiently existed on the surface of Mars more than 3 billion years ago. Much of this water is hypothesized to have been sequestered in the subsurface or lost to space. We use rock physics models and Bayesian inversion to identify combinations of lithology, liquid water saturation, porosity, and pore shape consistent with the constrained mid-crust (∼11.5 to 20 km depths) seismic velocities and gravity near the InSight lander. A mid-crust composed of fractured igneous rocks saturated with liquid water best explains the existing data. Our results have implications for understanding Mars’ water cycle, determining the fates of past surface water, searching for past or extant life, and assessing in situ resource utilization for future missions.

Liquid water existed at least episodically on Mars in rivers (1), lakes (1), oceans (2), and aquifers (3) during the Noachian and Hesperian, more than 3 billion years ago.Mars lost its ability to host persistent bodies of liquid water on its surface after the planet lost most of its atmosphere during this time period (4).The ancient surface water may have been incorporated in minerals (5), buried as ice, sequestered as liquid in deep aquifers, or lost to space (4).
Geophysical measurements have the potential to identify water in the deep subsurface.For example, seismic velocities derived from ground motion measured by the InSight (interior exploration using seismic investigations, geodesy, and heat transport) mission and interpreted with rock physics models have been used to constrain water distribution to depths of 20 km beneath the InSight lander, Elysium Planitia.The shear V s and compression V p wave velocities within the upper 300 m beneath InSight are consistent with a dry crust composed of minimally cemented (<2% of the pores) sediments (6).V s in the upper 8 km beneath InSight is lower than expected for an ice-saturated cryosphere (7), though V s may be higher elsewhere (8,9).Kilburn et al. (7) argue that the crust between 8 and 20 km beneath InSight is a) mafic and highly porous or b) felsic and less porous, but with V s alone, could not determine whether the fractures contain liquid water.
We assess whether V s (10)(11)(12)(13), V p (12), and bulk density b (14) data (Table 1) are consistent with liquid water-saturated pores in the mid-crust (11.5 ± 3.1 to 20 ± 5 km) within 50 km of the InSight lander.The mid-crust is one of four robust seismically detectable kilometer-scale layers beneath InSight (10)(11)(12)(13) and may be global (8).V p and layer thickness have been challenging to obtain for other locations on Mars (see ref. 9 and references therein).Temperatures on present-day Mars become warm enough for stable liquid water near the top of mid-crust (15), and pores are expected to have closed at the bottom of the layer (16).We use Bayesian inversion and a Markov chain Monte Carlo (MCMC) algorithm (17) to identify combinations of six lithologic parameters (pore shape aspect ratio , porosity , liquid water saturation w , mineral bulk modulus m , mineral shear modulus m , mineral density m , Table 2) that best reproduce the three observed data points V p , V s , and b (Table 1).Calculations combine the seismic velocity equations, the Berryman self-consistent rock physics model (18), and the Gassmann-Biot equations (19) (Materials and Methods).A mid-crust composed of igneous rock with thin fractures filled with liquid water can best explain the geophysical data.

Results and Discussion
Fig. 1 summarizes inversion results when the MCMC algorithm samples a range of mineral moduli and densities spanning from mafic (14,20) to more evolved igneous rocks (14,21) represented by a range between 100% basalt and 100% plagioclase.Several combinations of parameters produce good fits to the observed V p , V s , and b data within assumed errors (Fig. 1 V -X )., , m , and w are well resolved.A fully liquid water-saturated crust w = 100% is most probable (Fig. 1F ); is estimated as 0.17 ± 0.07 (Fig. 1C) and as 0.19±0.18(Fig. 1A), implying thin fractures.The inversion recovers a nonlinear relationship between and (Fig. 1B).m is not well-constrained by the data (Fig. 1K ).
A mid-crust containing liquid water has implications for the Martian water budget and hydrological cycle.Assuming the InSight location as representative, motivated by similar V p /V s (1.81 to 1.98) and seismically derived (0.1 to 0.17) (8) beneath InSight and areas up to 4,500 km away from the lander, 10 km of crust with porosity of 0.1 to 0.2 translates to 1 to 2 km of water-more than the water volumes proposed to have filled  hypothesized ancient Martian oceans (2).Thus, Mars' crust need not have lost most of its water via atmospheric escape.Liquid water in the pores of the mid-crust also requires high enough permeability and warm enough temperatures in the shallow crust to permit exchange between the surface and greater depths.While available data are best explained by a water-saturated mid-crust, our results highlight the value of geophysical measurements and better constraints on the mineralogy and composition of Mars' crust.

Materials and Methods
Constraining the Mid-Crustal Bulk Density.The bulk density of the midcrust has not been directly constrained by the gravity, seismic velocity, and mineralogical data used to derive the average bulk density and thickness of the crust beneath InSight (14).We can, however, infer the bulk density of the mid-crust using three constraints.First, the average bulk density within the upper 1.2 km is 1,600 ± 360 kg/m 3 and 2,300 ± 130 kg/m 3 between 1.2 and 11.5 km.These numbers are based on the estimated average bulk densities within the upper few hundred meters below the surface (26) and ∼5 km below the surface (27) of the adjacent Gale Crater on Mars.Second, the bulk density of the crust increases with depth (22).Third, the bulk density of the layer beneath 20 km ± 5 km is the same as its mineral density due to pore-closure (16).An average bulk density of the mid-crust can be obtained by solving a constrained problem to reproduce the average bulk density of the crust, 2,580 ± 209 kg/m 3 (14).
Rock Physics Models.Seismic velocities V p and V s depend on bulk density b and effective shear e and bulk e moduli: Berryman's rock physics model (18) provides dry-frame shear d and bulk d moduli of fractured rocks [see ref. 7 for a list of Berryman's equations (18)].The model uses a self-consistent approach and long-wavelength scattering theory that allows inclusions to overlap (18).Model inputs are , m , m , m , and .e = d (19).We use Gassmann-Biot fluid substitution theory (19) to estimate e from d , , m , and the bulk moduli of the fluid in a dry ( f 1 = 0 kPa for gas) versus partially to fully liquid-saturated ( f2 ) rock, With constraints on e and e from Berryman and Gassmann-Biot equations (18,19), we then estimate V s and V p via Eq. 1.
Bayesian Inversion.We perform a Bayesian inversion, which requires that we specify a prior p 0 (x) and a likelihood p l (y|x), where x are the six unknown parameters that we invert for ( , , w , m , m , and m , which control e , e , and b ) and y = (4.1 km/s, 2.5 km/s, 2,589 kg/m 3 ), [3] are the three data (V p , V s , and b ) we seek to explain.The prior is a uniform distribution over the parameter bounds in Table 2, combined with the constraint that V p > V s .The likelihood follows from assuming Gaussian errors in the data p(y|x) ∝ exp −0.5 W(y − m(x)) 2 2 , [4] where m(x) is the rock physics model (i.e., the forward model) and where W is a diagonal matrix whose diagonal elements are the reciprocals of the standard deviations of the data ( V p = 0.2 km/s, V s = 0.3 km/s, b = 157 kg/m 3 , derived from Table 1 to render all reported data points as likely).Jointly, the prior and likelihood define a Bayesian posterior distribution, p(x|y) ∝ p 0 (x)p l (y|x), which we sample via an affine invariant MCMC ensemble sampler (17).Sensitivity analyses confirm that water saturation does not significantly influence V s (19) and most strongly influences the V p , followed by b (18).

Fig. 1 .
Fig. 1.Summary of inversion results.Panels (A-U): Histograms of marginal posterior distributions of model parameters, computed from 5 × 10 5 iterations of the MCMC (17).The area under each histogram is equal to one.In the 2D histograms, cold colors (blues) indicate low posterior probability, and warm colors (yellows and whites) indicate regions of high posterior probability.In the 1D histograms, black stair plots show results for our default parameters bounds (Table2).The light gray stair plots in panels (C) and (F ) illustrate results obtained with widened bounds on mineralogical parameters (Results and Discussion).Water content is nearly uniformly distributed (F ) under these assumptions, but the porosity takes on unreasonably large values (0.29 ± 0.07).Panels (A) and (C) show that and are tightly constrained by the data.Panel (B) reveals a nonlinear relationship between and .Panel (F ) indicates that a high water saturation is likely in view of the data.Panel (J) shows that m is not constrained by the data.Panels (V -X ): Data fits.Histograms show model responses (V p , V s , and b ) for each of the parameters in panels (A-U), normalized so that the area under the graph is one.The orange error bars (horizontal) illustrate the mean of the data (filled dot) and expected errors (two SD).

Table 1 . Geophysical data for the mid-crust beneath the InSight lander
See Materials and Methods for b calculations, which assume crustal mineralogies ranging between 100% plagioclase and 100% basalt.