Similar results have been obtained during experiments of X‐ray irradiation of silica at temperatures of 100 K << 300 K (Griscom,; Revesz,, and references therein). Griscom () reviewed electron spin resonance experiments of silica samples with 1,200 ppm of hydroxyl irradiated by hard X‐rays and subsequently annealed at 77 K and 100 K. These works inferred that the production/destruction of metastable hydroxyl occurred via the following reactions (Griscom,):where (.) and (o) represent a valence electron and electron hole, respectively. The reaction in equation 2 proceeds to the right under irradiation and to the left during annealing. Griscom () showed for such samples that the reaction in equation 3 is limited by diffusion, as the reaction rate occurs more than 4 orders of magnitude faster than the diffusion rate. Various studies examining the mobility of hydrogen in silica have shown that diffusion was best characterized by using a distribution of activation energies to represent the random sizes and of interstitial and trapping sites in irradiated silica (Devine,; Griscom,; Shelby & Keeton,). The mean activation energies derived from (H, H) diffusion in various pristine and amorphous silica samples ranges from 0.18 to 0.73 eV (Fink et al.,, their Figure 6).

Fink et al. ( 1995 ) performed experiments of hydrogen implantation and diffusion in silica annealed in the temperature range of 320 K–550 K. Permanent OH is not removable by degassing at these temperatures (Lee, 1964 ). However, the Fink et al. experiments demonstrated that ion irradiation also enhances diffusion. That is, irradiation at low temperature leads to a dynamic equilibrium between the formation of dangling bonds (immobile atoms with an unsatisfied valence) and diffusing hydrogen (H, H 2 ; Fink et al. 1995 ; Griscom, 1984 ). The dangling bonds become trapping sites for interstitially diffusing hydrogen.

Experiments performed under conditions more applicable for the lunar environment have shown that metastable O‐H is also produced during irradiation of silica at low temperatures (Griscom, 1984 ; Fink et al., 1995 ; Zeller et al., 1966 ). Zeller et al. ( 1966 ) irradiated silica glass with protons and found strong evidence for the enhancement of the OH absorption feature as seen by a monotonic increase of the optical density of the feature with increasing flux. Furthermore, it was observed that with increasing dosage, the optical density approached a constant, thus, indicating that irradiation was both producing and destroying OH. Zeller et al. found that up to 100% of the implanted H could be converted into OH at the onset of bombardment. When a similar experiment was performed with helium and hydrogen‐2 ions, no significant increase in optical density was observed.

It is well known in the silicon microelectronics industry that material defects significantly alter the diffusive properties of silica. Impurities in natural silica alter the spacing between oxygen atoms, which in turn affect the local electron cloud, and the covalent bonds within OH groups. This population of hydroxyls is commonly referred to as metastable (Fink et al., 1995 ; Lee, 1964 ). The existence of metastable OH has been inferred from experiments of loading H 2 gas to silica samples and subsequent annealing (Devine, 1985 ; Lee, 1963 , 1964 ; Shelby & Keeton, 1974 ). Lee ( 1963 ) found the IR absorption feature about 2.75 μm became deeper and more broadened during experiments that loaded H 2 gas to silica laden with impurities at temperatures of 900 K–1,100 K. Upon degassing in vacuum at T ~ 1,073 K, it was found that the IR absorption feature was removed. However, when the same treatment was applied to highly pure samples, it did not produce any significant change in their respective OH feature. Later, Lee ( 1964 ) using experimental data from IR (absorption feature) and mass (H 2 degassing rates) spectrometers showed that the amount of metastable OH could significantly outnumber the amount of permanent OH in samples. For example, they reported a metastable OH content of ~100 ppm compared to <5 ppm of permanent OH in an experiment annealed at 1,100 K. It was concluded that samples with impurities possessed populations of both free hydroxyl characteristic of strong O‐H bonds with a well‐defined potential minimum and metastable O‐H bonds characterized by potential wells of varying depth.

From equation 1 an approximate timescale for diffusive loss can be defined as τ 1 = h 2 / D ( U, T ; e.g., Starukhina, 2006 ). As described in Starukhina ( 2006 ), the degassing rate decreases significantly in cooler regions due to the exponential dependence on temperature. Therefore, at high latitudes where the surface has a high U / T ratio more hydrogen is accumulated even though the solar wind flux is lower. In contrast, at low latitudes hydrogen is more easily degassed even though the solar wind flux is higher because the U / T ratio is lower and D 0 exp(− U / T ) is larger. Starukhina ( 2006 ) emphasized that the lunar surface would have a broad range of activation energies characteristic of trapping by a variety of physical and chemical defect sites. Estimates of the diffusional lifetime for atoms trapped physically and chemically are given in Table 1 . The lifetimes were evaluated using the τ 1 with D 0 = 10 −6 m 2 /sand h = 100 nm from Farrell et al. ( 2015 ). Within crystalline SiO 2 , the energy barrier to vacancy diffusion and H bound in OH groups is ~1–3 eV. This barrier to diffusion of H bound in OH groups is high compared to the thermal energy of surface bound atoms ~0.03 eV, (Starukhina, 2006 ). However, diffusion via interstitial hopping and along grain boundaries is more favorable (energy barrier ≲0.1 eV) due to increased interatomic spacing and lower atomic binding energies. Starukhina ( 2006 ) estimated that activation energies of U ≳1.3 eV and U ≳0.8 eV for the lunar surface at low and high latitudes, respectively, are required to trap hydrogen concentrations on the order of ≲160 ppm. Therefore, it was suggested that vacancy and chemical trapping would be required to build up the H concentration at low latitudes, whereas even physical trapping could be important at high latitudes.

Starukhina () predicted that solar wind elements are favorably retained at high latitudes by investigating the diffusional timescales of implanted hydrogen. The steady state concentration of atoms was balanced by incident solar wind flux using Fick's first law= −(d/d), whereis the diffusive flux,exp(−) is the diffusion coefficient (m/s),is a preexponential factor,is the activation energy,is temperature in eV, and d/dis the concentration gradient. From Fick's first law the diffusive flux in terms of the total density,, within the implantation layer of depth,, is approximated by equation 1

In warmer regions the loosely bound surficial hydrogen attached to a metastable OH (Fink et al., 1995 ) is possibly outgassed as H 2 (Starukhina, 2006 ). Stern et al. ( 2013 ) reported on a substantial H 2 exosphere with densities between 1,200 and 9,000 cm −3 . Both LAMP and the Apollo 17 orbiter both performed UV observations of the exosphere constraining the near‐surface atmospheric density to <9,000 cm −3 (Feldman & Morrison, 1991 ). Cook et al. ( 2013 ) analyzing LAMP spectra noted a dusk to dawn asymmetry of 1,700 ± 400 and 2,100 ± 300 cm −3 , respectively. Partial pressures measured by the Chandraayan‐I Chandra's Altitudinal Composition Explorer (CHACE) mass analyzer were used to infer near‐surface densities of 500–800 cm −3 (Thampi et al., 2015 ). In that study, the estimate of the neutral densities depended on the assumption of the surface temperature. Indirectly, the Neutral Mass Spectrometer (NMS) on the Lunar Atmosphere and Dust Environment Explorer inferred consistent neutral densities by detecting H 2 ions produced by photoionization (Halekas et al., 2015 ).

Hydroxyl groups have been identified on the Moon by observation of infrared absorption features in the wavelength range of near ~2.8 μm (Clark, 2009 ; Pieters et al., 2009 ; Sunshine et al., 2009 ). Sunshine et al. ( 2009 ) reported on observations by the Deep Impact (DI) High‐Resolution Instrument‐infrared spectrometer of absorption features at low latitudes 20–60°N that indicated a diurnal dependence of water content with the strongest absorptions occurring near the terminators. From morning to noon, the band depth decreased by ~70%. At high latitudes >70°N the strongest absorption features were observed but with less diurnal variation. Li and Milliken ( 2017 ) constructed the first quantitative global maps of OH content from M 3 observations that were consistent with the DI findings. They found that the water content varied by 200 ppm over the lunar day at latitudes between ~30° and 70°. However, at latitudes above 70° and below 30° a diurnal variation was not identifiable. The maximum content of ~500–750 ppm was observed in the Northern Hemisphere at latitudes above 70° with values consistently larger than that in the Southern Hemisphere, whereas below 30° the surficial OH content was constrained to values <~100 ppm (Li & Milliken, 2017 ).

One pathway for thermal hydrogen to leave the surface is as molecular hydrogen through recombinative desorption (Starukhina, 2006 ). In this process, the diffusing H atoms escape the surface binding potential by combining with another H via the sharing of valence electrons. The resultant H 2 molecule has enough kinetic energy to escape the 1–2‐eV surface binding energy (Starukhina, 2006 ). Observations of the H 2 lunar exosphere obtained with LAMP are consistent with this hypothesis (Hurley et al., 2016 ; Stern et al., 2013 ). Additional observations of energetic neutral H (Wieser et al., 2009 ) and methane (Hodges, 2016 ) suggest that H atoms can find a variety of pathways to escape the surface potential, but molecular hydrogen appears to be the dominant pathway.

Schematic of proton implantation and the degassed hydrogen corona. Incident SW protons (p+) implanted in the regolith excite atoms leading to ionization, excited electrons, high‐energy neutrals producing various types of physical and chemical defects. The diffusion of H atoms is slowed when they form metastable bonds with oxygen O‐H. Direct diffusion of hydrogen and recombinative desorption leading to the formation of H 2 are pathways for hydrogen to degas into the exosphere.

Experiments have shown that when silica is irradiated, an infrared absorption feature, near 2.8 μm, grows as the result of hydroxyl (OH) production (Fink et al., 1995 ; Griscom, 1984 ; Zeller et al., 1966 ). The absorption feature has been attributed to a combination of permanent and metastable hydroxyls present in the silica (Fink et al., 1995 ; Lee, 1963 , 1964 ). It appears that this process is occurring on the Moon due to solar wind‐implanted hydrogen as observed by three independent spacecraft missions (Clark, 2009 ; Pieters et al., 2009 ; Sunshine et al., 2009 ). Within regolith grains implanted hydrogen is present in atomic form and combined with various molecular species. As shown in Figure 1 , implanted protons diffuse and can chemically combine with other regolith atoms, like oxygen, or become trapped in physical defects. This process is driven by the energy deposited along the path of penetrating protons that energizes atoms leaving the regolith oxides with dangling bonds (Griscom, 1984 ; Zeller et al., 1966 ). Subsequently, the interstitial diffusion of hydrogen is hindered through forming permanent and metastable bonds with regolith oxides. For example, in pristine silica (damage free), hydrogen has a derived diffusion coefficient with a preexponential factor of D 0 ~ 10 −9 m 2 /s (Lee, 1963 ); however, experiments of diffusion in irradiated silica indicate that the mobility of hydrogen is significantly reduced in damaged SiO 2 , for example, D 0 ~ 10 −12 m 2 /s (in Fink et al., D ( U = 0.52 eV, T = 350 K) = D 0 exp(− U / T ) ~ 10 −19 m 2 /s).

Lunar soils are predominantly composed of olivine (Mg, Fe) 2 SiO 4 , pyroxene (Ca, Mg, Fe)SiO 3 and plagioclase feldspars (Ca, Na)AL 2 Si 2 O 8 (Williams & Jadwick, 1980 ). Down to 60 cm below the surface the bulk density of the regolith is estimated to be in the range of 1,500–1,800 kg/cm 3 . Oxygen and silicon are the two most abundant elements accounting for an average atomic weight percent of ~60% and ~16%, respectively (Turkevich, 1973 ). Typically, they are in the form of silica (SiO 2 ) with wt. % ~45 (McKay et al., 1991 ).

The solar wind is the tenuous plasma emitted from the solar corona. For nominal conditions, the typical flux at Earth is ≈2.0 × 10 8 cm 2 /s, with the composition consisting of primarily protons and electrons, combined with a few percent of He ions and fractional percentages of heavier ions. The protons have energies in the keV range, and the nominal value is typically taken to be 1 keV. Because the Moon does not have a significant global magnetic field, the exposed surface is under constant bombardment by the solar wind particles, except for ~5 days when the Moon is in the geomagnetic tail. When the protons penetrate the surface, they are neutralized after undergoing several interatomic collisions while penetrating regolith grains. The penetration or implantation depth depends on the energy of the incident ion, the angle of incidence, and the composition of the target surface and is typically about 20 nm for a 1 keV proton (Farrell et al., 2017 ). During this process, SW ions energize atoms and molecules via momentum transfer events and excite electrons along its path (Johnson, 1990 ). On the other hand, the concomitant radiation damage creates physical and chemical trapping sites that also hinder the mobility of hydrogen (Fink et al., 1995 ; Starukhina, 2006 , 2012 ; Zeller et al., 1966 ).

Before the ground‐breaking observations of the lunar OH veneer in 2009 (Clark, 2009 ; Pieters et al., 2009 ; Sunshine et al., 2009 ) pioneering theoretical investigations and experiments proposed that the solar wind (SW) is a dynamic source for delivering hydrogen to the surfaces of rocky bodies without atmospheres (Mattern et al., 1976 ; Starukhina, 2006 ; Zeller et al., 1966 ). In this study, we examine the role of the solar wind as a widespread source of hydrogen and hydroxyl to the lunar surface by quantitatively estimating the SW‐derived H (i.e., OH) surface content and the H 2 exosphere densities. In our approach we use a Monte Carlo model to calculate the dynamic buildup of OH over time due to populations of molecules with high activation energy diffusional pathways. To this end, we apply a statistical mechanics approach developed in Farrell et al. ( 2015 , 2017 ) to estimate the diffusional lifetime of H atoms in the silica surface hindered both by physical defects and by forming temporary bonds with oxygen atoms. The statistical distributions are used within a global Monte Carlo simulation to track the temporal OH surface content and the subsequently degassed H 2 exosphere. Li and Milliken ( 2017 ) quantified the global distribution of OH surface content using refined observations of infrared absorptions near 2.8 μm in the surface reflectance by the Moon Mineralogy Mapper (M 3 ) instrument on the Chandrayaan‐1 mission. The observations indicate that OH content increases with latitude up to 500–750 parts per million (ppm) at high latitudes, and at low latitudes there appears to be a diurnal variation of up to 200 ppm over the lunar day linked to warmer surface temperatures. The model results are in good agreement both with the global maps of OH surficial content produced from M 3 observations discussed in Li and Milliken ( 2017 ) and observations of H 2 in the exosphere by the Lyman Alpha Mapping Project (LAMP) on the Lunar Reconnaissance Orbiter, (Cook et al., 2013 ; Stern et al., 2013 ).

Guided by previous experimental work as discussed in section 1.3 , our approach examined the effect of diffusion on the exospheric hydrogen budget. Therefore, to estimate an upper limit of the molecular hydrogen density, we assumed that all particles escaped as H 2 by reducing the degassed particle weight by half. This is equivalent to removing two hydrogen atoms from the surface inventory. The degassed molecules were assigned a random velocity determined from the local surface temperature. Therefore, we emitted molecules into the exosphere using a thermal speed randomly selected from the Maxwell‐Boltzmann flux‐speed distribution (Brinkmann, 1970 ), and the velocity vector was directed with a cosine distribution about the surface normal. The upper radial bound of the grid was specified to be 5–10 scale heights above the surface. If a molecule traversed the upper boundary with a speed larger than the escape speed, it was removed from the simulation. Otherwise, the molecule was tracked until it returned ballistically back into the averaging domain. We also tracked the shadow of the Moon projected onto the exosphere to consider photo‐induced losses. The lifetime of H 2 against photodestruction is taken to be τ photo ~ 6.8 × 10 6 s (2.5 months) from Huebner and Mukherjee ( 2015 ). Even though there was a probability for molecules to be lost from photochemical processes, r > exp(− dt /τ photo,ADS ), the primary loss mechanism was thermal escape. We did not consider the lifetime for H 2 adsorbed on the surface because above ~15 K its sticking probability approaches 0 (Acharyya, 2014 ; Penteado et al., 2017 ).

For the solar wind we applied typical conditions of density n sw = 5 × 10 6 H/m 3 and speed v s w = 400 km/s. The local flux incident to each surface element is estimated using its zenith angle with n sw v sw cos( Z ) assuming that the Sun is in the z ‐plane. The implanted atoms and exospheric molecules are represented statistically by using a particle weight. The particle weight was derived from the SW flux so that each computational particle represented 10 24 implanted hydrogen atoms. Leveraging the studies of Farrell et al. ( 2015 , 2017 ) each implanted particle was given an activation energy, U , selected from a Gaussian energy distribution (equation 5 ) using a Monte Carlo choice. Likewise, each atom was also prescribed a random depth, h , chosen from a Gaussian distribution based on TRIM profile. Then at each integration step, the diffusion lifetime was determined using by evaluating the probability for an implanted particle to degas as r > exp(− dt/ τ num ), where r is a uniformly distributed random number ranging from 0 to 1. We considered that all surface loitering hydrogen bond (permanently or metastable) with oxygen in order to obtain an upper limit estimate of the hydroxyl surface concentration. We did not directly consider the effect of the reaction rate of H in regolith oxides, which is another limiting factor to the production of OH. However, based on the hydroxyl observations of 10–1,000 ppm a lower bound estimate of the reaction times is on the order of > ~50 s and > ~1 day for the subsolar point and terminator regions, respectively (Farrell et al., 2015 ).

We have adapted the Monte Carlo model applied in Tucker et al. ( 2015 ) to the Moon in order to track the surface inventory of atomic hydrogen and subsequently degassed molecular hydrogen to the exosphere. The computational domain is centered on the Moon and discretized in spherical coordinates in the azimuth and zenith angles using 8,000–20,000 surface elements. The corresponding surface areas ranged from 45 to 6,600 km 2 from the poles to the equator. In the radial direction we used a nonlinear grid composed of spherical cells. The exosphere grid was discretized into 60,000 cells with volumes of 40,000–3 × 10 6 km 3 with the lower volumes corresponding to cells near the surface and the larger for cells near outer boundary of ~6,000 km. Each cell was given a radial width chosen to be a fraction of the scale height corresponding to the lowest surface temperature of 100 K. We used an inertial frame of reference centered on the Moon, and the surface temperature was calculated using T (Z) = 280 K*cos 1/4 ( Z ) + 100 K (Butler, 1997 ; Crider & Vondrak, 2000 ). At temperatures of 100 K–380 K the mean thermal speed of H 2 ranges from ~1,000–2,000 m/s. Based on the ratio of the lowest radial cell width divided by the thermal speed, we obtained a time step on the order of tens of seconds. The integration period for averaging macroscopic properties during the simulation was limited to time intervals <10 min, so that molecules traveled distances less than a scale height per integration period. More detailed models could include the effect of the lunar rotation in the Earth‐Sun system. For simplicity, we have used an inertial frame of reference and considered only the near‐surface exosphere. The Moon's rotation in a three‐body system results in 1–2° variation of subsolar point (e.g. Tenishev et al., 2013 ). However, including such effects would not alter the conclusions of this study.

After diffusing to the surface, many hydrogen atoms leave in a molecular form to escape the surface potential (Starukhina, 2006 ). This is supported by Apollo era measurements that place an upper limit of ~10 cm −3 of thermal atomic hydrogen in the exosphere. Ground‐based and remote observations indicate that the dominant pathway of hydrogen loss from the surface is in the form of molecular hydrogen (Cook et al., 2013 ; Hodges, 1973 ; Stern et al., 2013 ). Helium measurements, which have been well characterized, strongly suggest that nonthermal sources produced by micrometeoroid impacts and sputtering are not significant (Hodges, 1973 ). Therefore, the source hydrogen to the exosphere depends on the diffusion of hydrogen and its conversion into H 2 within the surface. To this end, we have examined diffusion‐limited degassing as the limiting factor to the exosphere. The advantage to this Monte Carlo approach is that both the population of hindered H atoms can be tracked and those that are lost in the H diffusion and H 2 escape process (H‐H 2 pathway).

Subsequent modeling of timescales longer than a lunation resulted in a buildup of surface concentration because not all of the implantations contribute to an immediate dynamic equilibrium over the lunar day. Given the distribution of activation energy with U 0 = 0.5 eV and U W = 0.1 eV, those few atoms with energies at the high‐energy tail of the distribution can remain trapped over periods longer than a lunation. These populations of trapped H atoms accumulate with time, and the accumulated effect of this small but ever‐growing population is shown in Figure 2 c. Hence, while a dynamic equilibrium condition applies in warm regions to most of the H atoms undergoing diffusion, in cooler regions there is a population of high U valued implantations that accumulate and dynamic equilibrium does not apply.

We used the Monte Carlo approach to track the diffusion lifetimes of hydrogen atoms in the surface. In Figure 2 a, we show a contour map of the surface concentration calculated with equation 7 using the parameters U 0 = 0.5 eV and U W = 0.1 eV. As discussed, the concentration is largest in cooler regions along the terminator. The analytical solution assumes dynamic equilibrium. At solar zenith angle of 0° and 89° the volumetric surface concentration is 6 × 10 21 and 2 × 10 23 H/m 3 , respectively, where the subsolar point is defined at 0° latitude and local time 6 hr. As discussed below in section 2.1 , each implanted atom was assigned a random implantation depth and activation energy. These values are then used to calculate the diffusion lifetime on an implantation‐by‐implantation basis. This lifetime is used to make a Monte Carlo choice to determine if an H 2 molecule escapes the surface. In Figure 2 b we plot the bulk average H concentration using an energy‐averaged diffusion time based on the Monte Carlo implantations in a given region. The activation energy‐averaged numerical Monte Carlo model and the activation energy‐integrated analytical result are shown to be consistent for timescales much shorter than a lunation, providing a check that the model was performing as expected.

The Gaussian distribution has the form ofwhereis the energy of the peak of the distribution andis the width about the peak. Over timescales much less than the diffusive timescale, a quasi‐steady state density of implanted H atoms can be obtained using the continuity equation:whereandare the solar wind density and velocity respectively andis the solar zenith angle (Farrell et al.,). In equation 6 the diffusive flux is defined in terms of a mean value because of the use of the Gaussian distribution. Upon integration over all energies, an effective diffusion coefficient is obtained. Similar approaches have been used in the studies of hydrogen diffusion in silica glass (Devine,; Shelby & Keeton,). Using equation 6 , Farrell et al. () define an analytical expression of the implanted density as a function of solar zenith angle assuming dynamic equilibrium on timescales much less than a lunation:whereis defined as the diffusive velocity, or in terms of a diffusive lifetime. A comparison of expressions ofandorandshows that the diffusion rate is decreased by the factor ofwhen considering a distribution of activation energies. This is due to the fraction of the distribution that possesses activation energies less than the mean value of, which is characterized by the width of the distribution,, as shown in Figures 2 a and 2 b. On the other hand, 50% of the distribution also has, and therefore, over long time periods the net result is an increase in the surface concentration compared to the monoenergetic case, as shown in Figure 2 c. That is, over time atoms and molecules with high activation energies will dominate the surface concentration.

Farrell et al. ( 2015 ) showed that using a distribution of activation energies to characterize hydrogen retention in the regolith could simulate the observed diurnal variation of the IR absorption feature. As discussed therein, it is expected that implanted H atoms would have a range of activation energies in space‐weathered grains due to variable interatomic spacing. Each atom can be considered to have activation energy dependent on its individual diffusion path but characteristic of a mean value. Without knowledge of this distribution a priori, a Gaussian distribution of energies was used to characterize both populations of atoms that are easily outgassed and those that are trapped for long time periods.

Here we model the experimental and theoretical approaches of hydrogen diffusion discussed in Fink et al. ( 1995 ) and Starukhina ( 2006 ). That is, the implanted hydrogen diffuses by hopping from O to O in the bound SiO 2 ‐rich regolith. The O then hinders the H diffusion (Fink et al., 1995 ), to form metastable OH. Once on the surface, the mobile H atoms find each other and leave the surface via recombinative desorption H + H = H 2 . Therefore, a dynamic equilibrium exists between the diffusion of atoms and molecules slowed by trapping in defects and the subsequently degassed molecular hydrogen to the exosphere.

3 Results

3.1 Surface Concentration Three case studies were simulated using a Gaussian activation energy distribution each defined with a different mean energy of U 0 = 0.3, 0.7, and 0.5 eV, respectively, for comparison to Farrell et al. (2015, 2017). In the Farrell et al. studies a range of mean activation energies was considered to characterize surfaces defined with U 0 < 0.3 eV as emissive, U 0 > 0.9 eV as retentive, and 0.3 eV < U 0 < 0.9 eV as diurnal varying. Emissive surfaces degassed their hydrogen content in less than a lunation, whereas retentive surfaces possessed more high activation energy trapping sites retaining atoms longer than a lunation period. For the intermediate energy range Farrell et al. (2015) found that the surface concentration could vary diurnally. However, these previous studies considered the dynamic buildup of the hydrogen surface concentration on timescales much smaller than a lunation. Here we simulate the dynamic buildup of hydrogen over several lunations. The adjustable parameters are D 0 , U 0 , and U W for the activation energy distribution and h 0 and h W for the implantation depth. To this end, we used parameters derived/inferred from experiments and theoretical studies: U W = 0.1 eV, D 0 = 10−12 m2/s, h 0 = 20 nm, and width in the implantation distribution of h W = 13 nm (Farrell et al., 2017; Fink et al., 1995; Griscom, 1984; Shelby & Keeton, 1974; Starukhina, 2006; Zeller et al., 1966). The simulations were started without any implanted H atoms, and we simulated the progression to a quasi‐steady state. In order to obtain an upper limit estimate of the hydroxyl content, it was assumed that all implanted hydrogen paired with oxygen to form a metastable O‐H pair like that defined in Fink et al. (1995). We calculated distributions of the surface content as a function of latitude and local lunar time. The results of the three simulation cases are presented below in parts per million (ppm) calculated using a regolith density of 1,700 kg/m3. Case Study 1, f(U 0 = 0.3 eV, U W = 0.1 eV): Result for the activation energy distribution centered at 0.3 eV with width 0.1 eV is shown Figure 3. For this distribution Farrell et al. (2015) obtained a result for which no atoms were retained during 8:00 am to 12:00 pm (lunar local time) at the equator, and at 80° latitude only 0.2% of atoms were retained. We note that Farrell et al. (2015) used a much larger D 0 = 10−6 m2/s derived from diffusion experiments with pristine silica. For comparison we show our result obtained with the same energy distribution in Figure 3a but using the smaller preexponential factor of the diffusion coefficient derived from the Fink et al. (1995) experiments of diffusion in irradiated silica glass. We obtained a quasi‐steady state concentration balanced by the SW flux after two lunations. At latitudes below ~85°, the surface concentration was balanced by the degassing of H 2 molecules and implanted solar protons during the lunar day. On the illuminated surface the concentrations below this latitude were nearly uniform with values < ~ 1 ppm. Above this latitude, independent of the local lunar time, there was a net accumulation on the order of ~19 ppm (OH) after 18 lunations, Figure 3a. However, it is expected that eventually, the surface concentration at high latitudes will reach dynamic equilibrium as well (discussed below in section 3.3). Figure 3a shows a diurnal variation that is largest near midlatitudes ~45°. While a time of day dependence is also obtained at low latitudes, the variation is not very significant, Figure 3b. The surface concentration was largest at the morning terminator because of the addition of freshly implanted atoms added to the concentration of atoms retained over the lunar night. Because most atoms degassed during the day (14 Earth days), the concentrations were lower at the evening terminator compared to morning. The day/night asymmetry obtained in our study was solely due to implantation and not the migration of hydrogen above the surface. Figure 3 Open in figure viewer PowerPoint (a) f(U 0 = 0.3 eV, U W = 0.1 eV) snapshot contour of dynamic surface concentration over the lunar surface as a function of lunar local time. The dashed lines identify the location of the morning and evening terminators and local noon. (b) Diurnal surface concentration extracted from (a) along 0° latitude. Case Study 2, f(U 0 = 0.7 eV, U W = 0.1 eV): Result for the activation energy distribution centered at 0.7 eV with width 0.1 eV is shown in Figure 4. This distribution results in retention times that exceed several lunations for a significant fraction of the implanted atoms. Therefore, the very large surface concentration (much larger than reported in Li & Milliken, 2017) is mostly reflective of the solar wind source. As shown in Figure 3a, the concentration increases from the poles with decreasing latitude to ~40°. However, at lower latitudes the accumulation is partially slowed due to the increased degassing occurring from the warmer surface regions. Although the residence time at local noon is less than a lunation for implanted atoms with U = 0.7 eV, atoms with higher activation energies stick around longer and the cumulative effect over several lunations masks any diurnal variations, for example, Figure 4b. Figure 4 Open in figure viewer PowerPoint (a) f(U 0 = 0.7 eV, U W = 0.1 eV) snapshot contour of dynamic surface concentration over the lunar surface as a function of lunar local time. The dashed lines identify the location of the morning and evening terminators and local noon. (b) Diurnal surface concentration extracted from (a) along 0° latitude. Case Study 3, f(U 0 = 0.5 eV, U W = 0.1 eV): Result for the energy distribution centered at 0.5 eV with width 0.1 eV is shown in Figure 5. The snapshot of the OH concentration after 18 lunations shown in Figure 5a is consistent with the Li and Milliken (2017) globally averaged map constructed from the M3 observations, for example, their Figure 1a. Farrell et al. (2017) used this same distribution and the preexponential factor to the diffusion coefficient of D 0 = 10−12 m2/s with equation 7 but estimated a dynamic mass fraction that was less than 0.3 ppm at latitudes above 85° for timescales much less than a lunation. Therefore, Farrell et al. showed that if 90% of the diffusing H atoms had an activation energy of 0.5 eV and 10% had an activation energy of 0.7 eV, the amount retained would be >30 ppm. In contrast, using the Monte Carlo approach with a single activation energy distribution (U 0 = 0.5 eV), we obtained a much larger OH mass fraction of ~700 ppm at high latitudes and ~45 ppm near the subsolar point after 18 lunations. The difference is due to the accumulation of those very slowly diffusing H atoms with higher activation energy that are not in strict dynamic equilibrium during a lunation period. Consistent with Li and Milliken (2017), we obtained concentrations < ~100 ppm at latitudes up to 30°, and above latitudes of ~70° the concentrations were constant with values of ~750 ppm. However, we obtained a diurnal effect even along the subsolar longitude. Li and Milliken noted in their analyses that a clear enhancement between dawn and dusk was not identifiable because the observations of local evening occurred closer to noon than the morning observations. In our study, we obtained a relative enhancement in concentration at dawn compared to dusk due to the summation of atoms retained throughout the night and newly implanted atoms added in the early morning hours. This population represents a fraction of newly implanted atoms with low barriers to diffusion characteristic of physical trapping sites such that they degas during the lunar day. For a surficial point rotating near the equator its concentration increases by ~15 ppm from 6 to 7 am, followed by a decrease of 65 ppm as it rotates to noon, and finally, there is another increase by 40 ppm as the point rotates to the night. On the cold night side with the source turned off the surface maintains a constant concentration of ~95 ppm. We did not obtain a variation in surface concentration of ~200 ppm at midlatitudes as shown for two spectra in Li and Milliken (2017), their Figures 6c and 6d. Nevertheless, our reported variation is within the range of the standard deviation of the observation. Figure 5 Open in figure viewer PowerPoint (a) f(U 0 = 0.5 eV, U W = 0.1 eV) snapshot contour of dynamic surface concentration over the lunar surface as a function of lunar local time. The dashed lines identify the location of the morning and evening terminators and local noon. (b) Diurnal surface concentration extracted from (a) along 0° latitude. From the three case studies, we find that the H concentrations match the IR observations for D 0 = 10−12 m2/s, U 0 ~ 0.5 eV, and U w ~ 0.1 eV. Activation energies below U 0 ~ 0.5 eV tended to release the hydrogen too fast, resulting in low surface concentrations. Values far above 0.5 eV tended to retain the hydrogen, resulting in high surface concentrations and a loss of the latitude and diurnal variations.

3.2 Results: Exosphere Densities Previous models of the Moon's H 2 exosphere have been presented by Hodges (1973), Hartle and Thomas (1974), Crider and Vondrak (2002), Hurley et al. (2016), and references therein. Hartle and Thomas (1974) used a Monte Carlo model to predict upper limits of H 2 densities to account for degassed hydrogen from solar wind implantation. The model of Hurley et al. (2016) examined solar wind and micrometeoroid sources to demonstrate that chemical sputtering is the most likely source mechanism supplying H 2 in the exosphere. Chemical sputtering refers to solar wind irradiation of the surface that leads to the formation of H 2 , which then degases at the local surface temperature. The release of atoms via chemical sputtering has a lower energy threshold than physical sputtering, being a thermal release. Hurley et al. (2016) reported that source rates in the range of 2.2–17.7 g/s were required to produce near‐surface exospheric densities in the range of 1,200–9,000 cm−3. A key difference of our approach is inclusion of the implantation model. That is, the previous studies derived exospheric source rates of H 2 by constraining models using observations of the exosphere. Here we follow the H to H 2 pathway from solar wind implantation, surface accumulation, and exospheric emission. Following the suggestion of Hurley et al. (2016), we have examined the effect of varying diffusion on exospheric content as a function of temperature. As expected for an exosphere in equilibrium with the surface temperature, the H 2 densities are more extended on the dayside, and largest densities occur near the surface on the nightside (Hodges & Johnson, 1968). On the dayside the exospheric densities peak near the subsolar point where degassing occurs most efficiently. On the nightside the densities peak after local midnight. This is due to the hopping distance of H 2 as molecules near the dawn terminator hop back to the nightside once the surface temperature warms. We obtained similar exospheric densities as described by Hodges and Johnson for the activation energy distributions centered at U 0 = 0.3 eV and 0.5 eV because in both cases a significant amount of hydrogen is degassed. However, the density distribution was significantly reduced in the Case Study 2 using the activation energy distribution centered at U 0 = 0.7 eV. Figure 6 shows the change in exospheric H 2 content in the case of slow diffusion in a highly reactive regolith. The exosphere is in balance between its surface source and escape. The density of the H 2 exosphere is primarily limited by thermal escape, for which H 2 has a lifetime of 4,200 and 9,600 s at the subsolar point and at the terminators (Z = 89°), respectively (Johnson, 1971). The diffusion lifetimes range from 4,200 s to 4 days for U = 0.5 eV. However, for U = 0.7 eV the diffusion lifetime can be much larger than thermal escape, for example, 1 day to thousands of years. Figure 6 Open in figure viewer PowerPoint Equatorial slice of exospheric density distribution after 18 lunations from simulation with (left) U 0 = 0.5 eV and (right) U 0 = 0.7 eV, respectively. The blue arrows indicate the direction of the solar wind. We can estimate the activation energy above which the exosphere is increasingly limited by diffusion by equating the thermal escape and diffusion lifetimes. For thermal escape , equating τ esc = τ 1 , we obtain , where v m = (2πkT/m)1/2, g = 1.9 m/s2 is the lunar surface gravity, and is the Jeans escape parameter (Johnson, 1971). For activation energy distributions with mean energies above 0.52 eV the abundance of the exosphere would be limited by diffusion out of the surface. For example, we obtain steady state source rates from the surface to the exosphere of 27 and 13 g/s for the activation energy distribution centered at U 0 = 0.5 eV and 0.7 eV, respectively. These loss rates correspond to thermal escape rates of ~8–9 × 1024 H 2 /s and ~3–4 × 1024 H 2 /s, respectively.