Gas terminal velocity from MIRO/Rosetta data using neural network approach

Context. The Microwave Instrument on the Rosetta Orbiter (MIRO) on board the Rosetta spacecraft was designed to investigate the surface and gas activity of the comet 67P/Churyumov-Gerasimenko. The MIRO spectroscopic measurements carry information about the velocity of gas emanating from the nucleus surface. Knowledge of the terminal velocity of the H2O gas is valuable for interpretation of in situ measurements, calibrating 3D coma simulations, and studying the physics of gas acceleration. Aims. Using a neural network technique, we aim to estimate the gas terminal velocity from the entire MIRO dataset of nadir geometry pointings. The velocity of the gas is encoded in the Doppler shift of the measured rotational transitions of o-H16 2 O and o-H 18 2 O even when the spectral lines are optically thick with quasi or fully saturated line cores. Methods. Neural networks are robust nonlinear algorithms that can be trained to recognize underlying input to output functional relationships. A training set of about 2200 non-LTE simulated spectra for the two transitions is computed for known input cometary atmospheres, varying column density, temperature, and expansion velocity profiles. Two four-layer networks are used to learn and then predict the gas terminal velocity from the MIRO nadir measured o-H16 2 O and o-H 18 2 O spectra lines, respectively. We also quantify the accuracy, stability, and uncertainty of the estimated parameter. Results. Once trained, the neural network is very effective in inverting the measured spectra. We process the entire dataset of MIRO measurements from August 2014 to July 2016, and investigate correlations and temporal evolution of terminal velocities derived from the two spectral lines. The highest terminal velocities obtained from H18 2 O are higher than those from H 16 2 O with differences that evolve in time and reach about 150 m s−1 on average around perihelion. A discussion is provided on how to explain this peculiar behavior.


Introduction
One of the objectives of the Rosetta mission was to better characterize and understand the smallest features achievable regarding surface morphology, its properties, and their link to coma of the comet. The Microwave Instrument on the Rosetta Orbiter (MIRO) measured spectral radiance from several preselected rotational transitions of key molecular species (H 2 O, CO, CH 3 OH, NH 3 ) in the 548-579 GHz frequency interval (Gulkis et al. 2007). These measurements contain information on the column density within the MIRO beam and can be used to estimate the respective production rates, Q[H 2 O] [molec s −1 ], which have been measured since the early arrival at the comet 67P/Churyumov-Gerasimenko Lee et al. 2015;Biver et al. 2015;Zhao et al. 2019). Marshall et al. (2017) used a method of pre-calculated interpolation tables to invert H 2 O line areas of the entire MIRO dataset to study spatial Tables with the estimated velocities are available at the Open Science Foundation link: http://dx.doi.org/10.17605/OSF.IO/XMN32 distributions of activity in greater detail. Biver et al. (2019) analyzed limb observations when MIRO controlled the spacecraft pointing to obtain production rates as a function of heliocentric distance for multiple species. Self-consistently retrieved line-ofsight (LOS) profiles of water density, temperature, and velocity from the H 16 2 O and H 18 2 O line shapes for a subset of data (May 2015) are presented in Marschall et al. (2019), which show deviation from the expected characteristics of gas expansion into vacuum. In Lee et al. (2015) terminal expansion velocity was estimated simultaneously with water production rates from the H 16 2 O MIRO data on August 7-9 and 18-19 in 2014, for heliocentric distance >3.5 AU. This ensured optically thin or moderate opacity in the measured line. However, even under such conditions, the optimal estimation algorithm may suffer instabilities and the parameters may not be fully independently determined because of the need for a matrix inversion. Unfortunately, retrieving multiple parameters via the optimal estimation method would not work for observations where the line core is completely saturated as is the case throughout most of the year A&A 648, A21 (2021) 2015. Therefore, as detailed in our paper, we constrain the problem to a single parameter and study the neural network (NN) approach as described below.
A unique property of the MIRO measurements is that they provide a direct measure of the Doppler shift of the spectral lines from which a LOS expansion velocity can be extracted. This information cannot be easily obtained from other instruments on board Rosetta; therefore, this information may be useful for consistent interpretation of their measurements (especially in situ instruments). The directly estimated velocity data also serve to validate 3D gas coma models.
An explicit distinction should be made between estimating the velocity profile along the LOS that may be retrieved under suitable conditions from these data (e.g., Marschall et al. 2019) and the terminal expansion velocity (TEV) of the velocity profile. In the context of near-nucleus operations, when the spacecraft is effectively inside the coma, conditions might arise in which we do not sample final velocity of the neutral gas as is discussed later. An expanding gas reaches the TEV at some altitudes at which collisions with neutrals and/or charged particles are no longer effective and most of the (random) kinetic energy is converted into bulk motion (Combi et al. 2004;Davidsson 2008). Both of these asymptotic values are affected by initial gas and outflow speed. For this reason these values would be expected to exhibit a dependence on the heliocentric distance for near surface sublimation from cometary nuclei (Davidsson et al. 2010). In this paper, although we use the terminology of TEV, the interpretation of data during some periods may need additional consideration due to the changing spacecraft-nucleus distance throughout the mission. Similarly, because of different opacity in the two water transitions, the TEVs are generally sampled from different altitude regions along the LOS (additional discussion and an example is provided in Appendix A).
The goal of this work is to analyze the entire set of MIRO nadir measurements of H 16 2 O and H 18 2 O spectra to obtain TEV estimates and study the patterns of correlations with several key effective geometrical parameters, such as solar zenith angle (SZA), view angle (VA), and local time (LT). This first order comparison may provide valuable indicators (in a statistical sense) of the gas outflow physics and its variations in time and space. Nevertheless, the retrieved TEV may need to be assimilated into 3D kinetic models to learn precisely about the regional outgassing and observed gas coma parameters . In this study, we implemented a neural network for the TEV inversion.
The paper is structured as follows: In Sect. 2, the extracted MIRO data are described as prepared for our analysis. Then we provide a description of the NN structure, training, and performance for the task of extracting terminal expansion velocities in Sect. 3, and we follow that with a description of the key results (Sect. 4). Section 5 provides discussion, summary, and future outlook.

MIRO data preparation
The MIRO spectroscopic measurements are nominally provided with 30 s integration time as a result of three cycles of the frequency switch mode 1 . At full resolution of 44 kHz or about 24 m s −1 at 557 GHz, the random noise level is about 10 K after 30 s integration time. Although accurate line shapes are not the key requirement in this work, we have to apply additional time averaging to decrease the random noise for reliable estimations. Selecting the appropriate time averaging window is important in order to preserve a good spatial resolution. Therefore, the averaging time window is set to minimum 5 min and maximum 15 min (resulting in a typical random noise of about 4 to 1.5 K). The variable averaging time window allows us to obtain more data points that a constant time window might have to reject.
To automate nadir spectral line identification for the averaging, we use a combination of three guidelines: (i) the measured continuum temperature has to be over 30 K; (ii) the line center (line core) cannot be above the continuum temperature to reject pure limb emission lines; and (iii) the entire MIRO beam, taken as twice the half-power beam width (HPBW), must be within the nucleus shape model. This calculation uses the SPICE toolkit (Acton 1996) and relies on their accuracy. Despite this effort, there are traces of spectra in our database that visually appear as if a portion of the beam is seeing limb emissions (half-nadir half-limb). This occur mostly when the spacecraft is further than 200 km from the surface of the comet. Nevertheless, these traces of spectra do not pose a difficulty for the NN to properly identify the nadir component of the line and extract the correct TEV. With

Neural network modeling
A NN approach has two main advantages as a functional approximator: first its ability to naturally deal with nonlinear mapping between inputs and outputs, and second its speed in evaluating new examples after the training phase is done. A NN model consists of nodes (neurons) interconnected through weights (Gardner & Dorling 1998). Each neuron receives an input data that is transformed into an output value according to an activation function. Neurons are grouped in two or more layers and the network is said to be fully connected if every neuron in a layer is linked with every neuron in the successive layer through a weight value: the larger the weight, the larger the contribution given from the corresponding connection to the output of the layer. A graphical representation of a NN is sketched in Fig. 1.
Our first step is to determine the NN structure (Sect. 3.1), then train it to learn how to infer TEV from measured spectra A21, page 2 of 12 L. Rezac et al.: MIRO terminal velocity x 3 x 2 Example of activation function computation for a single neuron. The value of σ(z +1 i ) is computed from Eq. (1). The symbols have the meaning specified in Sect. 3.1. In the example shown, n = 3 (number of neurons in layer ), i = 1 (neuron considered in layer + 1), k = 1, 2, 3 (neuron identifiers for layer ). The value of σ(z +1 i ) becomes x +1 j , which is used to compute the activation function value for neurons j in the successive layer + 2.

Network structure
We built two separate NN models: one trained for H 16 2 O spectra and another for H 18 2 O. We experimented with different NN layouts in order to obtain the simplest structure (i.e., the smallest number of neurons and number of hidden layers), which ensures at least 95% accuracy in estimating the TEV from synthetic spectra, for which we know the correct TEV value by definition. These requirements were satisfied with NN having two hidden layers, each containing 28 neurons. The networks are the same for both spectral transitions.
The activation function σ(z) defines how the output of each neuron is computed from its input (shown in Fig. 2). We adopted the rectified linear unit (ReLU) activation function, defined for neuron i of layer as where w ik is the weight value from neuron k of layer to neuron i of the successive layer + 1; x k is the input value from neuron k of layer ; b i is the bias of neuron i in layer , that is, a constant value used to offset the result of the activation function; and n is the number of neurons in layer . The quantity σ(z +1 i ) becomes the input value x +1 i of the considered neuron i of layer + 1 to be used to compute σ(z j, +2 ) of neuron j in the successive layer + 2. Compared to other activation functions, the ReLU ensures faster computational time and convergence (Glorot et al. 2011). Initial values for weights w ik and biases b i , which are changed during the learning process known as backpropagation, were initialized with a normal distribution.

Training set
The performance of the NN strongly depends on how well the training set represents the real data to which the network will be applied. During the training phase, the network is fed with synthetic (calculated) spectra, for which true TEV is known, and learns by iteratively adjusting the values of weights and biases in order to match its output to the true TEV values. The synthetic spectra database is generated by applying our a 1D spherical symmetric general non-LTE forward model based on the accelerated lambda iteration (ALI), which has previously applied to the MIRO measurements (Marshall et al. 2017;Marschall et al. 2019) and also validated for icy-moon atmospheres (Yamada et al. 2018). The 1D model radial profiles of density follow the Haser distribution, whereas velocity and temperature profiles are parameterized through two variables hyperbolic tangent function and a power law, similarly to Lee et al. (2015) and Zhao et al. (2019). The only additional step is to vary the free variables to create a database that samples the necessary range of densities and velocities as described below.
The training set consists of synthetic nadir spectra for each water isotopolog transition, sampling a variety of input coma conditions using a 1D Haser model. The generated atmospheric inputs have a range of TEV from 0.1 to 1.2 km s −1 and water column density range from 10 11 to 10 19 cm −2 to capture optically thin and thick conditions. The NN is trained with a 80/20 approach: 80% of randomly selected inputs are used for training and the other 20% for validation.
About one-half of the training set (1100 cases) of coma inputs was generated by the 1D Haser model, and the other half is supplemented by scaling sampled atmospheres from direct simulation Monte Carlo (DSMC) simulations of 67P/CG Marschall et al. 2016). The Haser model for density profile already depends on assumed velocity and does not provide any guidance on temperature altitude dependence, which in our case is also an empirical modification of a power law similar to Lee et al. (2015) (Eq. (1)-(3)). We find that by including the DSMC samples we get an increase in NN validation accuracy of about 15%. The DSMC model itself is set up for the illumination conditions of the comet on 2015-07-10T01:32:39; subsolar point is −30.48 • , such that the southern hemisphere of the cometary surface is illuminated. The boundary conditions for the sublimation of H 2 O ice from the nucleus is driven by the insolation at the surface using a simple thermal model from Liao et al. (2016) and Marschall et al. (2016). The H 2 O ice is homogeneously distributed all over the nucleus with an effective active fraction of 10%, meaning that 10% of every unit area is active and responds to solar irradiation in the same way. The DSMC simulation domain used extends up to 10 km from the center of the nucleus. To get the gas properties at the spacecraft position, we extracted information from the closest point within the simulation domain and perform a linear interpolation to obtain values for the number density, temperature, and velocity along the LOS as done in Marschall et al. (2019). Two DSMC outputs for production rates of 10 26 and 10 28 [molec s −1 ] were used, sampling profiles around the nucleus at 10 • steps over the plane spanned by the vector pointing to the Sun and that pointing to the spacecraft position in the comet frame for the chosen date. The extracted profiles from the DSMC provide samples for day, night and terminator conditions. By scaling the density and shifting the velocity of these profiles, we increased their number to about a half of the total number of inputs.
Improvements in the NN performance were expected because the DSMC provides radial profiles of temperature and density near the nucleus that leads to a better agreement in spectral slope of line wings just outside the saturated line cores in the case of the H 16 2 O transition. Such inputs cannot be generated by the adopted Haser 1D model with the assumed scaling laws A21, page 3 of 12 of temperature and velocity. In addition, by sampling the different illumination conditions from the model output, as described below, a realistic variability of inputs is obtained and provides a further improvement. The benefit of including the DSMC inputs, it should be made clear, does not rest on the expectation that the 3D model provides accurate absolute values of the inputs. First, a physical inconsistency is guaranteed by shifting the extracted velocity profiles and scaling the water densities. In this work this fact has no effect on the NN performance, the output accuracy, or the latter interpretation. Density and temperature profiles do not influence the center frequency Doppler shift, which is the parameter of interest. Second, the line center frequency of the water transitions, already for moderately optically thick conditions, forms well above the actual computational domain of the DSMC (∼10 km). Therefore, the benefit lies in the fact that by sampling DSMC outputs (including modifications we described) we obtain a set of synthetic line shapes that the 1D Haser parametrization cannot capture. It is in principle possible to create a more complex 1D parametrizations that would be able to provide better and more complete approximations observed in the MIRO lines. At the same time, we would expect that it is possible to create a comprehensive training database using only DSMC samples, if more runs were available. Nevertheless, the hybrid approach we describe gave us a perfect trade-off managing computational speed, flexibility in range, sampling density of inputs, and reducing implementation complexity. In the case of future work, when a NN is tasked with a more difficult problem, such as entire vertical profile retrievals, the hybrid approach may have to be re-evaluated.

Validation
The validation step allows us to perform a statistical analysis of a trained NN outputs with the true values, assessing the performance of the network. As noted above, 20% of the entire training data are set aside for validation. In the validation phase, as well as in subsequent application to real measured MIRO spectra, we do not include absorption lines whose line area is smaller than 4.5 K km s −1 , which cannot be reliably distinguished from noise. Figure 3 shows  predictions is set within ±0.03 km s −1 with respect to the true TEV. The 30 m s −1 corresponds roughly to the uncertainty in the MIRO calibration of the frequency grid (i.e., the width of a single spectral channel; 25 m s −1 ). The NNs used in this work over several runs perform with an accuracy of 95% for H 16 2 O and 96% for H 18 2 O. Figure 4 further illustrates the NN performance, showing nine randomly selected validation spectra, the true TEV, and the predicted TEV. The uncertainty in NN outputs is estimated by performing 50 complete runs with fixed training and validation dataset. The estimated 1σ is a measure of the uncertainty inherent in the NN itself due to its nonlinear nature and design (e.g., shuffling the order of inputs in which they are presented to the NN during training, random initialization of weights and biases).
As shown in Fig. 4, the trained NN is able to very accurately perform the inverse mapping from spectra to TEV. In fact it is visually not possible to distinguish the vertical lines as they overlap each other, thereby confirming the excellent agreement between the H 16 2 O and H 18 2 O TEV estimates and the true value. We want to emphasize that even in cases of completely saturated line cores prevalent for H 16 2 O lines for high water column densities, the NN is able to extract TEV with great accuracy and statistical stability. We find no inherent biases in the trained NN for predicting TEV from the calculated H 16 2 O and H 18 2 O spectra. In other words, TEVs predicted from either spectral line are always within the uncertainty estimate to be inherent in the NN itself (1σ error), which is always smaller than the single channel uncertainty in velocity calibration of MIRO spectra.

Application to measurements
As in the validation phase, we estimate TEVs from the measured MIRO data by averaging 50 runs of the NN, each automatically differing in the selection of training profiles and weights initialization. We opt for this approach to capture and quantify the uncertainty in TEV estimation due to NN setup itself. The 1σ error is computed as the standard deviation of the 50 runs output. Additionally, we exclude from the prediction step any spectra with a line area smaller than 4.5 K km s −1 , which cannot be reliably distinguished from random noise. that occur for early and late mission phases, and for night side observations. All estimates of TEVs and their uncertainties are made available online 2 . Below, we present the main results of correlations and time evolution analysis of the retrieved TEVs.  Fig. 5 plotted as a 2D joint density distribution. The color bar in the plot is qualitative with bright (yellow) color indicating higher frequency of occurrence. We can see that the cloud of points aligns as expected for a positive correlation between the two TEV estimates. However, it is also evident that TEVs from H 18 2 O spectra show a considerable number of data points with larger absolute values, especially when TEV(H 16 2 O) is above 650 m s −1 . The estimated uncertainty on the individual data points can be considered small, ranging from 5 to 30 m s −1 in most cases. There are a few statistical outliers with larger uncertainty, mainly in the case of the weaker H 18 2 O isotope with lower signal-to-noise ratio (S/N). In addition, we notice that the range of estimated TEVs is between 0.5 and 0.9 km s −1 , as shown in Fig. 5. As the low TEV values (<500-600 m s −1 ) generally occur for conditions of lower number densities (e.g., shadowed regions and/or large heliocentric distances), the sampling is worse by the weaker H 18 2 O measurements than for main isotopolog H 16 2 O. This is also apparent in subsequent figures.

TEVs and line area variations
We also analyzed whether a relationship exists between the extracted TEVs from each transition and their corresponding spectrum line area, but we found no strong correlations. Part of the reason for the weak correlation stems from the high opacity of the H 16 2 O line for most of the observational period. Therefore, the growth of line area is nonlinear with increasing production rate. we also did not find a strong correlation. Therefore, from a global dataset perspective the TEVs are not connected to the line area by a simple linear relationship.

Continuum temperature variation
The MIRO instrument also provides direct measurements of continuum temperatures in two channels. Full details on the relation between the MIRO antenna and the physical temperature of a target is available in the user manual (Frerking et al. 2018). 3 Since we explore only trends and correlations we do not convert the continuum antenna temperatures to the physical temperature. Center frequencies of the continuum receivers for submillimeter (sub-MM) and millimeter (MM) wavelengths are 552 (∼0.53 mm) and 188 GHz (∼1.6 mm), respectively (Gulkis et al. 2007). Since the instrument operated mostly in dual mode, we have simultaneous spectra and continuum measurements, which are certainly connected by physical relationship. For instance, the (sub)surface temperature at which a molecule undergoes its last collision before it escapes from the surface determines its initial velocity. Figure 6 shows correlation plots of TEV with both sub-MM and MM continuum measurements for the two transitions. The color bar represents a density function, with brighter colors indicating a higher frequency of occurrence. In this overall picture we do not find a notable relationship with the sub-MM or MM continuum temperature for either of the water transitions. The low TEVs (H 16 2 O) (<0.6 km s −1 ) tend to occur predominantly when sub-MM antenna temperatures are below 175 K; however, these would correspond to times of earlier phases of the mission. Similarly, we can see that for TEVs (H 18 2 O) the data form a point cloud without any strong correlation. To highlight this point we present more detailed analysis of density plots split into several time segments in Appendix B. Nevertheless, even then it is difficult to extract a general conclusion from these plots. This may be not surprising since the different illumination and observation conditions are not separated as they evolve with time. A more elaborate analysis may be required to extract more information on the dependence of TEV and continuum temperatures.

Solar zenith angle correlations
Next, we consider the TEV variation with the effective SZA, LT and the VA 4 . The effective SZAs, LTs, and VAs represent the mean values of the respective variables determined as an average for all facets that lie within two times the HPBW of MIRO's beam projected on the surface using the 125k SHAP7 shape model (Preusker et al. 2017). The projected size of MIRO's beam varies with time as the spacecraft distance changes during the mission phases. Despite the complexity of considering effective mean variables, it is natural to look for correlations from a point of view of the entire dataset, to spot whether there are prevalent statistical relations among the variables. For example, it would be of interest to find an analog of spatial and temporal variations of gas TEV with SZA, as has been established for dust jets (Lai et al. 2019). Figure 7 summarizes the TEV(H 16 2 O) variation with LT (left), SZA (middle), and VA (right). The blue nuances correspond to the joint probability distribution between the variables: the darker the color, the higher the joint probability. The 1D distribution over a given parameter from the entire dataset is shown on top of each panel (LT, SZA, VA) and for TEV on the vertical axis of the rightmost plot.
The TEV 1D probability distribution curve reflects changes in the physical conditions of the coma (e.g., heliocentric distance) and the observation geometry. It is clearly not Gaussian shaped but has a general resemblance of a bimodal distribution. As already noted, the low TEV values would predominantly correspond to the early and late mission phases (months of 2014 and 2016), while the high velocities occurred in 2015, that is, near perihelion. In the 1D distribution over the LT we clearly get a bimodal structure with long tails on either side. This is due to a combination of nadir geometry (of this dataset) and a preference to operate the spacecraft close to the terminator orbit whenever possible. This is also reflected in the effective SZA probability distribution (middle panel), which peaks near 75 • . The VA distribution has a single, but broad, peak around 50 • . The shape of the joint probability is not symmetric around noon (LT 12) and shows a slight tendency for higher TEV values before noon, but the lower values (<0.5 km s −1 ) have non-negligible probability shifted toward later afternoon LT. The TEVs in the range 0.7 ± 0.05 km s −1 seem to occur with the broadest range of LT (from 5 to 20 h). Similarly, this range of TEVs shows a wide non-negligible probability range over SZA and VA distribution. The most frequent TEV over the entire dataset is about 0.8 ± 0.05 km s −1 , occurring mostly at around 9-10 h in the morning or 15 h afternoon. The TEV and VA joint probabilities appear more uniform, dominated by the single peak in the VA distribution. Nonetheless, TEV values above 0.7 km s −1 can also occur even at VAs larger than 70 • , unlike low TEV values, which tend to cluster between 20 and 40 • .
The TEV(H 18 2 O) joint probability distributions with LT, SZA, and VA are shown in Fig. 8. The panels, axis description, and color scale are the same as in Fig. 7. We note that minor isotopolog line strength is about 500 times weaker than the main line strength. Therefore, most of the TEV values in the early and late mission phases (2014 and 2016) are not sampled, giving the TEV 1D distribution a much narrower form, although still appearing multimodal. The missing data from early and later phases also present modified 1D distributions of LT, SZA, and VA; however, in a broad sense they retain the same characteristics as those for H 16 2 O. We note that H 18 2 O TEV distribution between 0.8-0.9 km s −1 displays a double-peak structure at ∼9 and ∼15 h LT. We can also see non-negligible occurrence frequency of 0.9 km s −1 velocities for LT around 5-7 h. On the other hand, TEVs of around 0.7 km s −1 shows a flatter distribution over a narrower range of LT (9-16 h). Compared to H 16 2 O data shown in Fig. 7, the joint distributions appear rather symmetric around the region of maximum probabilities for LT, SZA and VA parameters. This is because the H 18 2 O dataset does not sample the early and late mission phases properly owing to generally lower water vapor production rates during these times.    Figure 9 shows the time evolution of TEVs as a function of days from perihelion. The left panel indicates daily mean values for H 16 2 O (blue) and H 18 2 O (red) along with corresponding standard deviations. The weaker isotope transition is not detectable for low production rates in the early and late mission phases as pointed out previously. The range of TEVs is also smaller for the minor isotope since measurements of low production rate conditions (night/shadow) cannot be made with reliable sensitivity. Importantly, we can clearly observe the aforementioned difference in TEVs magnitude from the two isotopic lines. The right panel in Fig. 9 presents daily averages of the H 16 2 O and H 18 2 O TEV difference with their standard deviation shown as a thin vertical line. The differences are possible only when we have the same data point in time for both measurements, which explains the gaps in the data for the early and late mission phases. The data point clusters of the same color represent a particular month of the year (with perihelion passage occurring in August 2015).

Heliocentric evolution of TEV
The observed differences in TEVs is the most puzzling aspect of the two datasets. It is unambiguous that the differences grow larger approaching the perihelion. There, these differences reach about 0.15-0.2 km s −1 , which is more than five spectral channels of MIRO. The MIRO frequency calibration uncertainty is better than one channel, which indicates that these deviations are significant. For another view of the TEV differences, we plot nine randomly selected MIRO measured spectra in Fig. 10, one for each month in 2015 as shown in labels. The shallow absorption line (orange) is for the H 18 2 O and the deeper and the wider spectral line is for the simultaneously measured H 16 2 O isotopolog. The vertical lines in each panel correspond to the estimated TEVs by the NN. Even for individually randomly selected spectra, the differences are easily perceived as the separation between the vertical lines is growing when approaching August-September 2015. We scrutinized the NN inversion for potential biases and can confidently rule this out as a source of the discrepancy. spectral line. The H 16 2 O spectra with its fully saturated line core, makes the visual identification of correct Doppler shift more difficult. However, as shown in our simulations (see Fig. 3 and also Fig. A.1 with contribution functions), the correct Doppler offset qualitatively falls in the middle of the saturated (flat) portion of the line. The NN also accurately recovers the TEVs for cases in between optically thin and complete line core saturation. As shown in the NN validation (Sect. 3.3), the NN is able to obtain accurate (i.e., within one channel width) TEV from the H 16 2 O spectra independent of optical thickness or line shape complexity.
Neither frequency calibration nor NN inversion for TEV seem to be the reason for the observed differences. Hence, we need to consider the TEV differences estimated from the two different spectral lines arising through a physical mechanism. The heavier isotopolog (H 18 2 O) spectra yields a faster expansion velocity. Also, because of its smaller opacity, these spectra sample regions at around 10-20 km from the nucleus surface even for large water column density (see Appendix A). The line center of H 16 2 O spectra typically forms at higher altitudes, even above 100 km altitude for larger production rates. Therefore, if the H 2 O gas velocity begins to decrease for some reason above around 10-20 km, the TEV difference would be manifested directly in MIRO measurements of the radiance of the main water isotopolog. The cometo-centric distance of Rosetta varies with time and happens to reach consistently above 150 km altitude around and through the perihelion passage. Therefore, owing to the higher water production rates the MIRO measurement of the optically thick H 16 2 O are formed at higher altitudes (see Appendix A) than the optically weak H 18 2 O. Hence, if there is an extended source of water gas, this could very well fit the observed behavior of TEV. The extended sources and its effects (including velocity radial profiles) were already modeled to explain the observed slow decay of water rotational temperatures and densities for the comet 73P/Schwassmann-Wachmann 3B in Fougere et al. (2012).
Therefore, our first hypothesis builds on the extended source of gas occurring around the perihelion passage. As the comet approaches perihelion, the gas production rate is increased such that icy aggregates are lifted and start outgassing as they become embedded in the faster outflow (Keller et al. 2017). The velocity of released gas, while depending on the temperature of these larger aggregates, would be lower compared to the velocity of gas that has accelerated from the comet surface up to a given height. In addition, if we suppose they are not oriented on average such that the released gas in the direction of the bulk flow, it would also result in a slower average outflow velocity within the MIRO sampled volume. This hypothesis would also suggest that kinetic temperature of the gas should be higher, which has been already noted in several MIRO observations in May 2015 (Marschall et al. 2019).
A second plausible hypothesis to explain the implied deceleration of natural gas above the line center formation of H 18 2 O line is the ion-netural gas drag at the boundary of diamagnetic cavity (Gombosi 2015). The diamagnetic cavity region in coma of the comet is manifested when the incoming (solar) plasma is strongly decelerated as to become stagnant at some altitude range with an associated build-up of a magnetic barrier around it. This typically happens in the close nucleus region (approximately below 100 km depending on Q[H 2 O]) because the ion-neutral collisions tend to slow the in-flowing plasma. There are already reliable hybrid model simulations of this effects in the context of 67P/CG in Koenders et al. (2015). The observational evidence for the presence of the diamagnetic cavity, or a structure similar to it, in the coma of 67P/CG is provided by the Rosetta Plasma Consortium suite of instruments monitoring the proton A21, page 8 of 12 and alpha particle bulk properties as a function of time (Behar et al. 2017). It was found that from around mid-April 2015 until mid-December 2015 there is complete lack of solar wind particles at the location of Rosetta spacecraft because these are deflected by the magnetic barrier associated with the diamagnetic cavity. Nevertheless, there is still little knowledge about the coupling effects of neutral gas deceleration by the stagnant plasma, especially the magnitude, altitude extend of this mechanism, and whether the kinetic temperature increase would also be affected.
A global 3D DSMC hybrid models could be used to weigh in as to which scenario is a better explanation of the observed TEV differences obtained from the two isotopolog water transitions. However, a model capable of handling the coupling of gas-dust and ion-gas to obtain self-consistently the neutral gas velocity and temperature is not yet available to our knowledge. The desired model domain for this problem extends at least several 100 km with adaptive grid to allow at least 2 km resolution below 20-30 km. It is desirable to perform such simulations in the future to disentangle the TEV evolution of measurements from MIRO.

Conclusions
We implemented and analyzed the performance of a NN for retrieving gas terminal velocity from (H 16 2 O and H 18 2 O) rotational transition measured by the MIRO instrument in nadir geometry. The high spectral resolution of MIRO allows us to infer this information from the Doppler shift of the two spectral lines independently. Although the Doppler effect is measured along the MIRO beam LOS, it is expected that the TEVs represent an actual radial outflow values as long as the line center weighting function peak around about 10 km altitude or higher . As demonstrated in Sect. 3.3 the NN can be trained on pre-calculated examples to reach a very high accuracy (95%) extracting TEV even from saturated line cores in the H 16 2 O absorption line. Importantly, the NN method for TEV retrieval is free from systematic biases. We stress that characteristics of the training set, such as the size and coverage of the expected domain, plays a key role in the NN performance. In our case, we found that using even modest samples from 3D DSMC simulations had substantial effects (∼15% improved accuracy), which is mainly due to more physical vertical profiles of kinetic temperature and velocity than achievable with power law approximations typically used in 1D coma models.
We performed 50 independent iterations of training and TEV estimations for each measured spectra to estimate uncertainties originating from the NN method itself. The same approach was also used in the validation phase for the NN, for which we know the true values of velocities. We find that TEV uncertainties from the NN itself are smaller than one channel (30 m s −1 ) and are not systematically biased. The calibration uncertainties of MIRO frequency are also smaller than a single channel. The chirp transform spectrometer of MIRO has an absolute frequency stability of 4 KHz and is unlikely to get frequency inaccuracies of more than 10 kHz (0.005 km s −1 ) from the USO even in the worse case assumptions. The frequency stability of the CTS versus temperatures are documented 5 , and are also smaller than the channel width (44 kHz). Therefore, the derived TEV differences from the two isotopic lines of more than (0.15 km s −1 ) are highly unlikely a consequence of the NN method or the measurement uncertainties themselves.
We apply the NN to all suitable MIRO spectra (see Sect. 2) for both water transitions, based on a minimum line area criterion. We retrieve about 28 000 TEVs from H 16 2 O and about 12 000 from H 18 2 O spectra; the tables are available online 6 . The key points from our global TEV dataset analysis can be summarized as follows: -From the perspective of the aggregate dataset, we found no correlation between TEV and MM or sub-MM continuum temperatures measured by MIRO that would suggest a direct physical link with the surface. However, a more detailed study of this relationship is warranted in the future that would, for instance, focus predominantly on the data points taken under the same illumination and observation conditions. -Joint probability plots of TEV and LT reveal that the highest velocities are almost equally likely to occur in the range between 9-15 h LT. Nevertheless, for H 16 2 O, the most frequent TEV range of 0.8 ± 0.05 km s −1 has a bimodal distribution in the LT range of 9 and 15 h. This distribution is asymmetric and large velocities in excess of 0.7 km s −1 may appear over a broad range of LTs (5-20 h). The TEVs estimated from H 18 2 O spectra follow similar bimodal distribution characteristics of the main isotope. We add a word of caution to this local correlation analysis, since it has been shown that a surface area intercepted by MIRO beam may not be responsible for the bulk of molecules in the upper regions of coma where the spectral lines form. This statement has a general validity. These correlation-based conclusions should be viewed in a statistical sense and may not apply to specific data points over a particular time. More elaborate analyses, considering similar illumination conditions over a particular region, should be used to draw more solid and definite conclusions. -The TEVs from H 18 2 O do not always agree with those from H 16 2 O transition. There is a time dependence to these differences, starting at nearly perfect agreement from August 2014 until beginning of 2015 , but growing in May 2015 onward and peaking in August-September 2015 (∼150 m s −1 ). Subsequently, the differences decrease with time again. We ruled out the possibility that this discrepancy is due to either the NN method itself or frequency biases in MIRO calibration.
-Understanding that the MIRO measured H 16 2 O spectral line core forms at higher altitudes than the H 18 2 O (Appendix A), we considered two plausible mechanisms that could potentially explain the difference in their TEVs. First, we discussed an extended source of water molecules in the coma as the comet approaches perihelion, supplying slower and warmer gas above 20-50 km altitude. This effect was studied and proposed to explain measured decay profiles in rotational temperatures for a comet 73P-B/Schwassmann-Wachmann 3B Fougere et al. (2012). A second hypothesis we discussed suggests that an ion drag might be responsible for the deceleration of the neutral gas in the region around the diamagnetic cavity, which occurs in the coma for higher production rates. This type of structure was observed in 67P/CG (Koenders et al. 2015;Behar et al. 2017). However, the magnitude of this effect and whether it would also provide an increase in neutral gas kinetic temperatures, as was shown to occur in the upper coma regions in Marschall et al. (2019), remains to be modeled.
-In relation to the previous point, we also observed that the maximum TEV values from the two transitions do not necessarily coincide in time or space, especially for nearperihelion data. This might be expected if the information on the TEV from the Doppler shift of each isotopolog line comes from very different altitudes as discussed in the text. We hope that future model studies capable of combining dust-gas and gas-ion coupling may bring new knowledge about the source of these differences. -The NN approach for TEV retrieval could be also applied to other measurements of MIRO (H 17 2 O, CO, CH 3 OH, or NH 3 ). However, because of the low optical thickness in these spectral lines they also form in the lower altitude in the coma and the estimated TEV should follow closely the estimated TEV from H 18 2 O. However, it is still desirable to obtain TEVs from other molecular transitions whether or not they agree with our prediction, perhaps identifying conditions where the TEVs imply only a weak coupling with the water gas. This could shed light on other aspects of cometary activity and gas release from the nucleus. As a general remark we find the NN method surprisingly effective in handling the problem of estimating TEVs from the MIRO spectra, especially in the case of saturated lines. We emphasize the need to carefully consider and build the training datasets covering the expected domain to which the NN will be applied. The NN is a promising tool in combination with the 3D forward models in many domains of science application, for instance to approximate radiative transfer calculations, retrieve parameters from nonlinear relationships to measurements, or simply provide a quick first estimate to constrain the initial conditions of more computationally intensive retrievals.