Hydrodynamic outflows and the lunar volatile element depletion
Abstract
These three worksheets written for Maplesoft's Maple 2020 software package were used to produce FIgures 1, 2, and 3 (C-O-H-adiabats.mw) and Figure 4 (C-O-H-Na-adiabats.mw, adiabatic-solutions.mw) of the accompanying publication. PDF versions are also included for those readers without access to the Maple 2020 software package or those who wish to reproduce the presented results in other software packages or platforms.
Full text
> > > > > > restart : # start with a clean slate with linalg : # load the linear algebra package with plots : # load the plots package with CurveFitting : # load curve fitting package R0.008314 : # kJ mol.K kB 1.38 10 23 : # J K mp 1.67 10 27 : # kg # Thermodynamic Data: II. Vapor speciation, 0 3000 K fits # Equation a - Iron-wustite (IW) Buffer # 2FeO(s) = 2Fe(s) O2(g) dSa 0.127 : # kJ mol.K dHa 541.0 : # kJ mol # Equation 1 - H2O dissociation equilibrium # H2O(g) = H2(g) 1/2 O2(g) dS1 0.0561 : # kJ mol K dH1 247.2 : # kJ mol # Equation 2 - OH dissociation equilibrium # OH(g) = H(g) 1/2 O2(g) dS2 0.0438 : # kJ mol K dH2 184.2 : # kJ mol # Equation 3 - H2 dissociation equilibrium # H2(g) = 2H(g)
> > > > > > dS3 0.1161 : # kJ mol.K dH3 444.4 : # kJ mol # Equation 4 - O2 dissociation equilibrium # O2(g) = 2O(g) dS4 0.1312 : # kJ mol.K dH4 504.7 : # kJ mol # Equation 5 - CO2 dissociation equilibrium # CO2(g) = CO(g) 1/2 O2(g) # dS5 0.0853: # kJ mol.K # dH5 281.5: # kJ mol # Entropies of gas species, kJ/mol.K sH T0.02083 ln T0.00407 : # 1-H sO T0.0210 ln T0.0419 : # 2-O sH2 T0.0314 ln T0.0499 : # 3-H27 sO2 T0.0353 ln T0.000766 : # 4 O2 sOH T0.03212 ln T0.00131 : # 5-OH sH2O T 0.0346 ln T0.00843 : # 6-H2O sCO2 T0.0547 ln T0.10592 : # 7-CO2 sCO T0.03375 ln T0.0025 : # 8-CO sNa T0.020834 ln T0.0349 : # 9-Na # Parameter choices # NatoH 0.5: # the sodium-to-hydrogen ratio BY NUMBER ( critical number according to H2018 and MS1995) # CtoH 1/12 : # The carbon-to-hydrogen ratio BY NUMBER! (=12 according
> > (2)(2) > > > > > > (1)(1) to Hirschmann and Dasgupta, 2009) sT 5 107: # surface density of magma, kg/m2= 5 x107 mks for 2 lunar mass disk within 5 Earth radii, 1010for TMO sHT 100 10 6sT : # size of H inventory in kg/m2= 5000 mks or106 mks for 100 ppm H in lunar disk and terrestrial MO sCT 1 / 12 sHT 12 : # size of C inventory in kg/m2nominally equal to the H inventory by mass sNaT 0.35 2700 10 6sT : # size of the Na inventory in kg/m2 nominally 20 times the H inventory by mass S0.175 : # entropy per volatile element atom (H and C and Na), kJ/mol.K 3.14 10 4: # Keplerian frequency, 2.4 10 4 at Roche radius, s1 n5 : # number of atmospheric layers # initialize matrices # xCO2 vector n : # mole fraction, CO2 xNa vector n : # mole fraction, Na xCO vector n : # mole fraction, CO xH2O vector n : # mole fraction, H2O xH2 vector n : # mole fraction, H2 xO2 vector n : # mole fraction, O2 xOH vector n : # mole fraction, OH xH vector n : # mole fraction, H xO vector n : # mole fraction, O xLH vector n : # mass fraction H T vector n : # Temperature, K p vector n : # Pressure, bars u vector n : # mean molecular weight, amu eval sNaT ;eval sCT ;sHT; 47250.00 5000 5000 # do the buffered calculation of last equilibration to determine the atomic C-O-HNa abundances of the basal gas eq1 p 1 2xo2 1 2xh2 xh2o = exp dS1 Rexp dH1 R T : # H2O dissociation equilibrium eq2 p 1 2xo2 1 2xh xoh = exp dS2 Rexp dH2 R T :
> > > > (2)(2) > > # OH dissociation equilibrium eq3 p xh2 xh2 = exp dS3 Rexp dH3 R T : # H2 dissociation equilibrium eq4 p xo2 xo2 = exp dS4 Rexp dH4 R T : # O2 dissociation equilibrium # eq5 p 1 2xo2 1 2xco xco2 = exp dS5 Rexp dH5 R T : # CO2 dissociation equilibrium eq5 p xo2 = exp dSa Rexp dHa R T 10 2: # specify fO2 relative to IW; eq6 xh xo xoh xh2 xh2o xo2 xco xna = 1 : # all species explicitly accounted for eq7 S = xh sH T xo sO T xoh sOH T xh2 sH2 T xh2o sH2O T xo2 sO2 T xco sCO T xna sNa T R xh ln xh xo ln xo xoh ln xoh xh2 ln xh2 xh2o ln xh2o xo2 ln xo2 xco ln xco xna ln xna ln p 1 / xh xoh 2 xh2 xh2o xco xna : # this is the entropy per hydrogen atom, so that an adiabat in which the # of particles changes has the same entropy (per H atom, which is conserved along an adiabat) eq8 u = xh 1 xh2 2 xo 16 xoh 17 xh2o 18 xo2 32 xco 28 xna 23 : # mean molecular weight in amu, ignore the electrons, they're weightless eq9 shv = 2 / 18 sh2o 1 / 17 soh sh2 sh : eq10 sCT = 12 / 28 sco : # single species diatomic vapor eq11 sNaT = sna : # one-species atomic vapor eq12 xh2o = sh2o geff p 105 u 18 : eq13 xh2 = sh2 geff p 105 u 2: eq14 xoh = soh geff p 105 u 17 : eq15 xh = sh geff p 105 u 1: # eq16 xco2= sco2 geff p 105 u 44 : eq16 xco = sco geff p 105 u 28 :
> > > > (2)(2) > > eq17 xna = sna geff p 105 u 23 : eq18 geff = 2 3.14 kB T ump 1 / 2 : # Effective gravity is Pv/ v # eq18 geff=10: # Effective gravity for the Earth is prescribed eq19 sHT = shv xlh sT : # statement of mass balance eq20 p xh2o = 2.9 104xlh 0.01 2 : # solubility law, Sossi et al. 2023 sol fsolve eq1, eq2, eq3, eq4, eq5, eq6, eq7, eq8, eq9, eq10, eq11, eq12, eq13, eq14, eq15, eq16, eq17, eq18, eq19, eq20 , xh, xo, xoh, xh2, xo2, xh2o, xco, xna, xlh, p, u, sna, sco, sh2o, soh, sh2, sh, shv, T, geff , xh = 0.00000000000001 ..1, xo = 0.0000000000000000000000000001 ..1, xoh = 0.0000000000000000000001 ..1, xh2 = 0.00000000001 ..1, xo2 = 0.00000000000000000000000000000001 ..1, xh2o = 0.0000000001 ..1, xco = 0.00000000001 ..1, xna = 0.000000000001 ..1, xlh = 0.0000000001 ..0.01, p = 0.001 ..1000, u = 1 ..44, sna = 0.000000000025 ..250000, sco = 0.000000025 ..250000000, sh2o = 0.0025 ..250000, soh = 0.0025 ..250000, sh2 = 0.0025 ..250000, sh = 0.0025 ..250000, shv = 0.000025 ..2500000, T = 300 ..12000, geff = 0.001 ..1000 ; results eval xh, xo, xoh, xh2, xo2, xh2o, xco, xna, xlh, T, p, u , sol : # FdissH results 1 results 3 results 1 results 3 2 results 4 results 6 ; xH 1 results 1 : xO 1 results 2 : xOH 1 results 3 : xH2 1 results 4 : xO2 1 results 5 : xH2O 1 results 6 : xCO 1 results 7 : xNa 1 results 8 : xLH 1 results 9 : T 1 results 10 : p 1 results 11 : u1results 12 ; OtoH xH2O 1 xOH 1 xO 1 2 xO2 1 xH 1 xOH 1 2 xH2O 1 xH2 1 ; CtoH xCO 1 xH 1 xOH 1 2 xH2O 1 xH2 1 ;
> > (2)(2) > > > > (3)(3) > > NatoH xNa 1 xH 1 xOH 1 2 xH2O 1 xH2 1 ; sol T = 3002.603965, geff = 0.1794759856, p= 0.1206874420, sco = 11666.66667, sh = 1263.507259, sh2 = 3216.169130, sh2o = 2771.914532, shv = 4830.253176, sna = 47250.00000, soh = 723.9668128, u= 12.09337801, xco = 0.07493429933, xh = 0.2272320748, xh2 = 0.2892016564, xh2o = 0.02769486311, xlh = 3.394936486 10 6, xna = 0.3694586759, xo = 0.003681398237, xo2 = 0.0001382162741, xoh = 0.007658815887 u112.09337801 OtoH 0.04525410042 CtoH 0.08626186898 NatoH 0.4253085193 # calculate the thermal and speciation structure of an adiabatic column for i from 1 to n1 do p i 1p i / 10 : eq1 p i 1 1 2xo2 1 2xh2 xh2o = exp dS1 Rexp dH1 R T : # H2O dissociation equilibrium eq2 p i 1 1 2xo2 1 2xh xoh = exp dS2 Rexp dH2 R T : # OH dissociation equilibrium eq3 p i 1xh2 xh2 = exp dS3 Rexp dH3 R T : # H2 dissociation equilibrium eq4 p i 1xo2 xo2 = exp dS4 Rexp dH4 R T : # O2 dissociation equilibrium eq5 OtoH =xh2o xoh xo 2xo2 xh xoh 2xh2o xh2 : # isochemical atmospheric column eq6 CtoH =xco xh xoh 2xh2o xh2 : # isochemical atmospheric column eq7 NatoH =xna xh xoh 2xh2o xh2 : # isochemical atmospheric column
> > (2)(2) > > > > (3)(3) > > eq8 xh xo xoh xh2 xo2 xh2o xco xna =1 : # all species explicitly accounted for eq9 S = xh sH T xo sO T xoh sOH T xh2 sH2 T xo2 sO2 T xh2o sH2O T xco sCO T xna sNa T Rxh ln xh xo ln xo xoh ln xoh xh2 ln xh2 xo2 ln xo2 xh2o ln xh2o xco ln xco xna ln xna ln p i 1 1 / xh xoh 2 xh2 xh2o xco xna : # this is the entropy per volatile atom (C-H-Na), so that an adiabat in which the # of particles changes has the same entropy per atom sol fsolve eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8,eq9 ,xh,xo,xoh, xh2,xo2,xh2o,xco,xna,T,xh = 0.00000000000000000000001 ..1, xo = 0.0000000000000000000000000000000000001 ..1, xoh = 0.000000000000000000000000000000000001 ..1, xh2 = 0.0000000000000000001 ..1, xo2 = 0.00000000000000000000000000000000000000000000001 ..1, xh2o = 0.000000000000000000000001 ..1, xco = 0.00000000000001 ..1, xna = 0.00000000001 ..1, T= 30 ..10000 ; results eval xh,xo,xoh,xh2,xo2,xh2o,xco,xna,T,sol ; xH i1results 1 ; xO i1results 2 ; xOH i1results 3 ; xH2 i1results 4 ; xO2 i1results 5 ; xH2O i 1results 6 ; xCO i1results 7 ; xNa i1results 8 ; T i 1results 9 ; end: FdissH results 1results 3 results 1results 3 2 results 4results 6; Error, invalid input: eval received fsolve({0.4525410042e-1 = (xh2o+xoh+xo+2*xo2)/(xh+xoh+2*xh2o+2*xh2), 0.8626186898e-1 = xco/ (xh+xoh+2*xh2o+2*xh2), .175 = (xh*(0.2083e-1*ln(T)-0.407e-2)+xo* (0.210e-1*ln(T)+0.419e-1)+xoh*(0.3212e-1*ln(T)-0.131e-2)+xh2* (0.314e-1*ln(T)-0.499e-1)+xo2*(0.353e-1*ln(T)+0.766e-3)+xh2o* (0.346e-1*ln(T)-0.843e-2)+xco*(0.3375e-1*ln(T)+0.25e-2)+xna* (0.20834e-1*ln(T)+0.349e-1)-0.8314e-2*xh*ln(xh)-0.8314e-2*xo*ln (xo)-0.8314e-2*xoh*ln(xoh)-0.8314e-2*xh2*ln(xh2)-0.8314e-2*xo2*ln (xo2)-0.8314e-2*xh2o*ln(xh2o)-0.8314e-2*xco*ln(xco)-0.8314e-2*xna* ln(xna)+0.7501145606e-1)/(xh+xoh+2*xh2+2*xh2 ... e-10 .. 1, xo = 0.1e-36 .. 1, xo2 = 0.1e-46 .. 1, xoh = 0.1e-35 .. 1}), which is
> > (5)(5) > > > > (2)(2) > > > > > > > > (3)(3) > > (4)(4) not valid for its 2nd argument, eqns FdissH 0.005136186371 eval xH ; eval xO ; eval xOH ; eval xH2 ; eval xO2 ; eval xH2O ; eval xCO ; eval xNa ; eval T ; eval p ; 0.2272320748 0.09812843695 0.005035117600 `?`4`?`5 0.003681398237 0.0003412021273 2.545143104 10 7`?`4`?`5 0.007658815887 0.001999215444 0.00002946193989 `?`4`?`5 0.2892016564 0.3791848316 0.4459034485 `?`4`?`5 0.0001382162741 0.00001486529293 1.155560605 10 8`?`4`?`5 0.02769486311 0.04011072514 0.04459344475 `?`4`?`5 0.07493429933 0.08097563518 0.08505923757 `?`4`?`5 0.3694586759 0.3992450883 0.4193790236 `?`4`?`5 3002.603965 2423.604369 1754.877979 `?`4`?`5 0.1206874420 0.01206874420 0.001206874420 0.0001206874420 `?`5 eval sNaT ; 54000.0