Merge pull request #244 from VirtualPlanetaryLaboratory/DiffusionLimA…

…tmescChange

Diffusion lim atmesc change

examples/DiffLimWaterEscape/README.rst

@@ -0,0 +1,58 @@
+Accumulation of Atmospheric Oxygen due to Water Photolysis and Hydrogen Escape on the TRAPPIST-1 planets b and e with the LS16 Water loss model
+==================================================
+
+Overview
+--------
+
+Track water loss and oxygen build-up due to hydrodynamic atmospheric escape.
+
+===================   ============
+**Date**              08/21/2023
+**Author**            Megan Gialluca
+**Modules**           AtmEsc
+                      STELLAR
+**Approx. runtime**   less than 1 minute
+===================   ============
+
+M dwarf stars experience particularly long and violent pre-main sequence (PMS) phases.
+At this time, terrestrial planets that may be in the HZ during the main sequence experience
+intense amounts of XUV radiation, typically past the runaway limit. Additionally, planets
+always destined to remain interior to the HZ experience a heightened level of radiation
+across the system's lifetime. In both of these cases, the planets are subjected to strong
+hydrodynamic thermal atmospheric escape; in a pure water atmosphere, this photolyzes water 
+molecules and allows hydrogen to escape, and drag oxygen during energy-limited escape.
+During diffusion-limited escape, the escaping flow of hydrogen is throttled by a build-up of
+oxygen in the upper atmosphere as opposed to limited by the incoming stellar energy. During diffusion,
+oxygen drag cannot continue.
+This example plots the water loss and oxygen build-up over time for the terrestrial planets 
+TRAPPIST-1b (interior to HZ) and e (HZ planet). Additionally it marks the time diffusion-limited 
+escape begins.
+
+This is a recreation of the T1-e curve from Figure 2 of Gialluca et al (2023, in prep), and a 
+recreation of the T1-b curve but with a differing initial water content than in the published
+figure (10 TO instead of 200 TO initial for T1-b). The other change to the published version is
+T1-b and e use red and blue coloring, respectively, for better contrast.
+
+To run this example
+-------------------
+
+.. code-block:: bash
+
+    python makeplot.py <pdf | png>
+
+
+Expected output
+---------------
+
+.. figure:: T1e-b_EscapeThroughTime.png
+   :width: 600px
+   :align: center
+
+   Water loss [TO] (top) and oxygen produced/retained [Bar] (bottom) over time [Gyr] for TRAPPIST-1b (red) and e (blue)
+   for an initial water content of 10 TO [Terrestrial Oceans]. In this particular case, hydrodynamic escape is 
+   halted on both planets within 1 Gyr; for T1-e escape ends due to HZ entrance at ~370 Myr and for T1-b escape ends
+   due to desiccation at ~900 Myr. Both planets experience diffusion-limited escape during this simulation, which begins
+   at the time indicated by the vertical dashed lines; this occurs when the mixing ratio of free oxygen in the atmosphere
+   is equal to that of water in the LS16 water loss model. 
+
+

examples/DiffLimWaterEscape/b.in

@@ -0,0 +1,23 @@
+# Planet b parameters
+sName                     b              # Body's name
+saModules                 atmesc         # Modules to apply, exact spelling required
+
+# Physical Properties (from Agol et al. [2021])
+dMass                     -1.377           # Mass, negative -> Earth masses
+dRadius                   -1.117           # Radius, negative -> Earth radii
+
+# AtmEsc Properties
+dXFrac                    1.0
+dSurfWaterMass            -10
+dEnvelopeMass             0.0
+sWaterLossModel           ls16
+bInstantO2Sink            0
+sAtmXAbsEffH2OModel       bolmont16
+dJeansTime                -10
+
+# Orbital Properties
+dEcc                       0.00401        # Eccentricity
+dOrbPeriod                 -1.510826
+
+
+saOutputOrder Time -SurfWaterMass -OxygenMass WaterEscapeRegime

examples/DiffLimWaterEscape/e.in

@@ -0,0 +1,23 @@
+# Planet b parameters
+sName                     e              # Body's name
+saModules                 atmesc         # Modules to apply, exact spelling required
+
+# Physical Properties (from Agol et al. [2021])
+dMass                     -0.693           # Mass, negative -> Earth masses
+dRadius                   -0.922           # Radius, negative -> Earth radii
+
+# AtmEsc Properties
+dXFrac                    1.0
+dSurfWaterMass            -10
+dEnvelopeMass             0.0
+sWaterLossModel           ls16
+bInstantO2Sink            0
+sAtmXAbsEffH2OModel       bolmont16
+dJeansTime                -10
+
+# Orbital Properties
+dEcc                       0.00637         # Eccentricity
+dOrbPeriod                 -6.101013
+
+
+saOutputOrder Time -SurfWaterMass -OxygenMass WaterEscapeRegime

examples/DiffLimWaterEscape/makeplot.py

@@ -0,0 +1,75 @@
+"""
+Reproduces the T1-e curve from Figure 2 of Gialluca et al (2023, in prep)
+and the T1-b curve from that Figure except with 10 TO initial instead of 200 TO
+
+Megan Gialluca, UW, 2023
+
+"""
+
+import sys
+import os
+import pathlib
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+import vplanet
+import vplot
+
+# Get Path
+#path = os.getcwd()
+path = pathlib.Path(__file__).parents[0].absolute()
+sys.path.insert(1, str(path.parents[0]))
+from get_args import get_args
+
+model = vplanet.run(path / 'vpl.in', quiet=True)
+
+# Define Time Array
+Time = model.b.Time
+
+# Find time of diffusion-lim escape for both bodies
+bdifftime = Time[np.where(model.b.WaterEscapeRegime == 4)[0][0]].value
+edifftime = Time[np.where(model.e.WaterEscapeRegime == 4)[0][0]].value
+
+# Plot figure
+fig = plt.figure(figsize=(8,10))
+wl = fig.add_subplot(2,1,1) # Water Lost
+ob = fig.add_subplot(2,1,2) # Oxygen Buildup
+
+wid=2
+
+wl_b = [10-x.value for x in model.b.SurfWaterMass]
+wl_e = [10-x.value for x in model.e.SurfWaterMass]
+
+bcolor = vplot.colors.red
+ecolor = vplot.colors.dark_blue
+
+wl.plot(Time*1e-9, wl_b, color=bcolor, lw=wid, label='TRAPPIST-1b')
+wl.plot(Time*1e-9, wl_e, color=ecolor, lw=wid, label='TRAPPIST-1e')
+# Doesn't matter what this is, just to get the label in the legend:
+wl.plot([-8,-9], [0,0], color='xkcd:grey', lw=wid, linestyle='--', label='Onset of Diffusion-Lim Escape')
+
+ob.plot(Time*1e-9, model.b.OxygenMass, color=bcolor, lw=wid)
+ob.plot(Time*1e-9, model.e.OxygenMass, color=ecolor, lw=wid)
+
+for ax in [wl, ob]:
+    ax.tick_params(length=4, width=2, labelsize=18)
+    ax.tick_params(which='minor', width=1, length=2)
+    ax.spines['top'].set_linewidth(2)
+    ax.spines['bottom'].set_linewidth(2)
+    ax.spines['left'].set_linewidth(2)
+    ax.spines['right'].set_linewidth(2)
+    ax.axvline(bdifftime*1e-9, color=bcolor, linestyle='--', lw=wid)
+    ax.axvline(edifftime*1e-9, color=ecolor, linestyle='--', lw=wid)
+    ax.set_xlim([0,1])
+
+wl.legend(fontsize=16)
+wl.set_xticklabels([])
+wl.set_ylabel('Water Lost [TO]', size=18)
+ob.set_ylabel('Oxygen Produced [Bar]', size=18)
+ob.set_xlabel('Time [Gyr]', size=18)
+
+plt.tight_layout()
+
+ext = get_args().ext
+fig.savefig(path / f"T1e-b_EscapeThroughTime.{ext}", dpi=200)

examples/DiffLimWaterEscape/star.in

@@ -0,0 +1,9 @@
+sName	          star
+saModules	      stellar
+dMass           0.09
+dAge            1e6
+sStellarModel   baraffe
+dSatXUVFrac     1.e-3.03
+dSatXUVTime     -3.146
+dXUVBeta        1.172
+saOutputOrder 	Time -Luminosity -LXUVTot Temperature HZLimRecVenus HZLimRunaway HZLimMaxGreenhouse HZLimEarlyMars

examples/DiffLimWaterEscape/vpl.in

@@ -0,0 +1,17 @@
+
+sSystemName               trappist1
+iVerbose                  0
+bOverwrite                1
+saBodyFiles               star.in b.in e.in
+sUnitMass                 solar
+sUnitLength               AU
+sUnitTime                 YEARS
+sUnitAngle                d
+bDoLog                    1
+iDigits                   6
+dMinValue                 1e-10
+bDoForward                1
+bVarDt                    1
+dEta                      0.1
+dStopTime                 0.95e9
+dOutputTime               1e6