The first problem with the acquired data is the unusually high measurements towards the poles that are considered as incremental noise scattered over the entire longitudinal axis. The second problem is incomplete coverage along the latitude and hence the poles themselves. To estimate the missing information along the latitudinal region, the following strategy is applied. If ilmenite abundance is correlated with the classification of the highlands/Mare regions and if the pole region geology features highland characteristics, then the expected values of the known highland regions serve as estimates of the ilmenite content at the poles.

As shown in Figure 2, the distribution characteristics of these two regions deviate considerably, where the average abundances also vary from 3.38 wt% in Mare to 1.1 wt% in the highland regions. Therefore, the ilmenite content correlation is given, and the estimates over the missing latitudinal areas of the highlands are set to 1 wt%; this matches with the original assumption of the WAC for values under the detection ratio. To remove the incremental noise at the extreme latitudes, a mask is created from the Mare boundaries as per the method of Nelson et al. (2014) and merged with a constant separation at ϕ = ±56°. The replacement values for the mask are set equally to 1 wt%. After applying both estimates, a low-noise global ilmenite map is obtained as shown in Figure 3.

The global ilmenite abundance map of Figure 3 is now used as the input to Eq. (2), which results in the location-dependent ISRU hardware mass shown in Figure 4.

The initial radial distance to the Moon’s centre of mass r(ϕ, λ) influences the ideal Δv demand directly, as shown in Eq. (3), for ascent into a circular orbit at r _orbit with standard gravity g ₀ and standard gravitational parameter μ. Δ v ideal ϕ , λ = v orbit 2 + v ascent 2 = μ r orbit + 2 ⋅ g 0 ⋅ r ϕ , λ − r ϕ , λ 2 r orbit (3)

The global ground elevation data are now used in the form of a displacement map (Wright and Petro, 2019), which originates from the Lunar Orbiter Laser Altimeter (LOLA) measurements (Smith, 2015). The elevation ranges from −9.115 km to 10.757 km with regard to the reference radius r _ref of 1737.4 km and therefore defines r(ϕ, λ) globally. The displacement map is shown in Figure 5.

Evaluating the extreme values on a low lunar orbit (LLO) at 100 km altitude (r _orbit = 1837.4 km) using Eq. (3) yields Δ v min = Δ v ideal max r ϕ , λ = 1725.187 m s Δ v max = Δ v ideal min r ϕ , λ = 1725.204 m s .

The influence of ground elevation on Δv is therefore of the order of 0.001%, which is extremely low.

The second celestial influence, namely, the initial surface velocity v ₀, is either an additional Δv demand or a Δv reduction depending on the shared velocity components in the launch direction. Together with the sidereal rotation period and assumption of a spherical lunar surface of r _ref , the surface velocity can be derived as a function of the latitude ϕ, as displayed in Eq. (4). v 0 ϕ = 2 π 27.322 d a y s ⋅ cos ϕ ⋅ r ref (4)

Evaluating the extreme points of the polar and equatorial locations on Eq. (4) yields Δ v min = Δ v 0 ϕ = ± 90 ° = 0.000 m s Δ v max = Δ v 0 ϕ = 0 ° = ± 4.624 m s .

Comparing this range |Δv _max| − Δv _min with the ascent from the reference radius to a circular LLO of 100 km as Δ v ideal ( r ref ) = 1725.196 m s gives the Δv influence of the surface velocity to be of the order of 0.27%, which is significantly more than the influence of elevation but still considerably low.

Explicit transfer from any lunar geodetic point to an NRHO and vice versa entail a high-fidelity problem that is usually solved non-analytically, as in Trofimov et al. (2020). Additionally, there are multiple transfer strategies that can be deployed for different optimisation goals. Between optimisation of the required Δv and transfer time, two transfer options were analysed for the chosen scenario.

First, a long-duration transfer that features a very low required Δv of only 664.9 m s to an LLO at an altitude of 100 km, which is very close to the theoretical limit of 654.8 m s , requires minimum energy change (Whitley et al., 2018). Moreover, it features an almost complete independency with the surface location that is achieved by something similar to a three-impulse transfer, where the lunar sphere of influence is left to circle once around the Earth before reinsertion. This allows removal of any inclination restrictions but at a cost of a long transfer time of 100.1 days. If this transfer option is chosen, the influence of the transfer efficiency is extremely low and marginal to the ISRU dependencies derived earlier. In this case, the transfer dependencies can be neglected and location selection can be simplified based on only ISRU efficiency.

Oftentimes, a transfer time of 100 days is simply too long for certain applications as it may, for example, induce general system lag times and therefore poor dynamics in propellant delivery adjustments for the target missions. For this reason, a second transfer option is analysed as a direct transfer trajectory between the NRHO and surface as per Trofimov et al. (2020), featuring the shortest transfer time of only hours but at the cost of a higher Δv and greater location dependency. The direct transfer, illustrated in Figure 6, is the subject of the analyses hereafter and serves as a worst-case scenario for an NRHO transfer in terms of the location dependency.

In the previous work by Trofimov et al. (2020), a set of possible direct descent trajectories and their associated landing points and Δv demand were identified. The resulting map of scatter points for the southern 9:2 NRHO was taken as the starting point to derive a global Δv map. As this result does not contain the solutions of the cheapest trajectory for each location but rather all the solutions for direct descent, the data points can have both low- and high-cost solutions for the same location. Since the lowest cost option of a location is chosen during mission planning, a minimum estimation is performed by splitting the map into 20° square tiles, where constancy is assumed and the lowest value is set for the entire tile. This tiling on an equirectangular projected map gives higher resolutions towards the poles causing the problem; accordingly, the solution coverage is so low in the southern polar region that only high-cost solutions are present in a tile even if the neighbouring tiles may feature low-cost solutions. To mitigate this, the data were removed from particular high-cost trajectories of Δ v > 2985.65 m s , leaving a few non-defined tiles on the southern pole. If such data removal is not performed, up to 3300 m s transfer options would be carried over to the final map, which are clearly high-cost solutions, and would not be considered in a real mission. This process is visualised in Figure 7.

The non-defined tiles are estimated from their longitudinal neighbouring tiles via linear interpolation, which is only necessary in tiles at the southern polar region where there is already a higher geodetic resolution. This results in the final Δv map depicted in Figure 8.

Even though the data in Figure 8 are computed for the descent only, they also serve as estimates for the ascent, which is biased because these problems are not entirely symmetric. Additionally, it should be mentioned that even though a 2414 m s transfer is very viable, a transfer of around 2900 m s may be used in a real-world scenario based on a different transfer strategy through an LLO with a waiting time to reduce Δv. However, in the present analysis, Δv(ϕ, λ) is globally defined in Figure 8 with Δ v min = 2414.35 m s and Δ v max = 2985.65 m s .

To derive the associated transport mass costs, the previously determined Δv has to be applied to a specific launcher. The Argonaut, formerly known as the European Large Logistics Lander (EL3), is chosen as a starting point for this scenario. Its initial configuration is based on the published information from the European Space Agency (2023) of wet mass of 10,000 kg, dry mass of 1,600 kg, and payload of 2,100 kg as of the time of writing this article. Additionally, a hydrogen and oxygen propulsion system with an oxidiser fuel ratio of 6 and a specific impulse of 400 s is assumed. When this original configuration is applied the proposed mission with Δv _max, the fuel is depleted before the round-trip can be completed. Therefore, the launcher configuration has to be altered according to the needs of the scenario.

Launcher upscaling is performed by maintaining the dry mass constant at 1,600 kg but adding fuel until the mission can be completed. The minimal viable system features an empty H₂ tank upon arrival at the gateway and an empty O₂ tank at the lunar surface. The iteration scheme for upscaling is visualised in Figure 9, which can converge to a minimal viable launcher for any given payload. The iteration commences with an undersized launcher starting from the moon with “initial H₂” as the amount of hydrogen in the tank at the lunar surface, “wet mass” that together with a constant dry mass indirectly represents the amount of oxygen in the tank on the lunar surface, and “refill H₂” as the amount of hydrogen refuelled upon arrival at the lunar gateway. When Δv is applied to this simulated round-trip, the undersized launcher experiences one of the defined failure cases, which indicates a particular missing propellant at some point of the mission. “O₂ empty” and “H₂ empty” describe the depletion of the oxygen and hydrogen tanks during flight, which results in an increase of the desired propellant for the next iteration, at which point “O₂ leftover” and “H₂ leftover” trigger the analogue opposite. “H₂ insufficient” is a less obvious failure case in which the launcher returns to the lunar surface but does not have enough hydrogen to perform the next round-trip run as it can only be refuelled at the lunar gateway and not on the ground. Owing to the fact that the launcher used in the iteration is undersized, the value of “refill H₂” may be considered unidirectionally while still reaching convergence with a small enough increment of the propellant.

Since this method of upscaling effectively increases the mass ratio ( r m = m wet m dry ) of the launcher, this assumption becomes increasingly unrealistic. Additionally, from the perspective of the fuel depot, oxygen is delivered but hydrogen is also removed, thereby effectively trading their masses according to the exchange ratio ( r e x = m payload m H 2 ,refill ) . To decide upon the range of launchers that must be compared to derive the cost differences, Figure 10 presents the parameter space between the exchange ratio r _ex and mass ratio of the launcher r _m .

Here, two sections are greyed out, where r _ex < 1 for economic reasonableness and r _m > 10 as the soft border of the mass ratio for a realistic single-stage launcher. The chosen launcher frames were selected through a sequence of movements in the parameter space, starting with economically reasonable exchange ratios of 1.5 and 2.0 (solid black lines, also called milestones) that are then projected to their required mass ratios (dashed black lines). This currently holds a set of constant exchange ratios over the Δv range; however, to obtain comparable results, the mass ratio r _m has to be constant over one set as it represents the efficiency of the launcher. Therefore, the maximum value of the mass ratio (at Δv _max) is maintained constant over the Δv range to yield the chosen frame for the mass ratio (solid yellow line). When this set is then back-projected onto the exchange ratio (yellow dashed line), a span of exchange ratios can be achieved over the Δv range.

This determines the two chosen frames (yellow lines) as r m = 8.555 w i t h 1.5 ≤ r e x ≤ 3.303 r m = 10.688 w i t h 2.0 ≤ r e x ≤ 3.889 .

Analysing the problem on two frames with different mass ratios provides insights into the sensitivity towards more efficient launchers in general and their influences on location selection.

Using another iteration scheme targeting the chosen mass ratio r _m, a launcher can be converged for any Δv. The expended fuel is directly drawn from the simulated round-trip and normalised by the payload size, which can be combined into a direct mapping from the required Δv to spent fuel k _Flight in kg per kg of the payload. This dependency can be seen in Figure 11 for both mass ratios r _m = 8.6 and r _m = 10.7. In this comparison, the higher mass ratio r _m = 10.7 features a smaller absolute and relative growth, indicating that the differences in the spent fuel decrease with increasing launcher efficiencies.

To meet the additional demand for oxygen that the launcher requires for transport every year, the ISRU facility is upscaled linearly by its ISRU costs per kilogram oxygen k _ISRU for each location and corresponding fuel requirements. The scaling factor originates from the base configuration cost in Figure 4 and the annual base production of m _base = 23.9 t, which gives k ISRU ( ϕ , λ ) = K base ( ϕ , λ ) m base . The additional oxygen demand every year is derived from the yearly payload of m _pl,y = 8 t that is set to minimise scaling on the base configuration, oxidiser–fuel ratio of r _of, and spent fuel k _Flight depending on the two selected mass ratios r _m ∈ {8.6, 10.7}. Therefore, the fix costs representing the mass supplied from Earth towards the construction of the ISRU facility are determined through Eq. 5 K Fix ϕ , λ , r m = K base ϕ , λ + k ISRU ϕ , λ m pl,y ⋅ k Flight ϕ , λ , r m ⋅ r of r of + 1 + m year − m base . (5)

The expended hydrogen is fully considered as a mass cost since it is retrieved from the lunar gateway depot and assumed to be delivered from Earth. Here, the cost levels of the lunar gateway and lunar surface are simplified to be equal to those from the Earth for comparability, which rather overestimates the cost of hydrogen compared to costs on the lunar surface when supplied from Earth. Therefore, the dynamic costs representing the mass of hydrogen supplied from the Earth to the lunar gateway every year t are determined through Eq. 6 K Dynamic ϕ , λ , t , r m = t ⋅ m pl,y ⋅ k Flight ϕ , λ , r m ⋅ 1 r of + 1 . (6)

Combining both K _Fix and K _Dynamic, the final total costs of the mission over the location and time are K Total ϕ , λ , t , r m = K Fix ϕ , λ , r m + K Dynamic ϕ , λ , t , r m . (7)

Applying Eq. (7) to the earlier location-dependent results gives the total cost maps for the mission time of 20 years for r _m = 8.6 (Figure 12) and r _m = 10.7 (Figure 13). In a direct comparison between Figures 12 and 13, the flight cost influence is visibly less pronounced for the higher mass ratio. Moreover, a reduced variation in the total mission cost is observed.