The near-nucleus gas coma of comet 67P/Churyumov-Gerasimenko prior to the descent of the surface lander PHILAE

Context. The European Space Agency (ESA) Rosetta mission was the most comprehensive study of a comet ever performed. In particular, the Rosetta orbiter, which carried many instruments for monitoring the evolution of the dusty gas emitted by the cometary nucleus, returned an enormous volume of observational data collected from the close vicinity of the nucleus of comet 67P/Churyumov-Gerasimenko.
Aims. Such data are expected to yield unique information on the physical processes of gas and dust emission, using current physical model fits to the data. We present such a model (the RZC model) and our procedure of adjustment of this model to the data.
Methods. The RZC model consists of two components: (1) a numerical three-dimensional time-dependent code solving the Eulerian/Navier-Stokes equations governing the gas outflow, and a direct simulation Monte Carlo (DSMC) gaskinetic code with the same objective; and (2) an iterative procedure to adjust the assumed model parameters to best-fit the observational data at all times.
Results. We demonstrate that our model is able to reproduce the overall features of the local neutral number density and composition measurements of Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) Comet Pressure Sensor (COPS) and Double Focusing Mass Spectrometer (DFMS) instruments in the period August 1–November 30, 2014. The results of numerical simulations show that illumination conditions on the nucleus are the main driver for the gas activity of the comet. We present the distribution of surface inhomogeneity best-fitted to the ROSINA COPS and DFMS in situ measurements.


Introduction
Several works have undertaken the derivation of the distribution of the surface activity of comet 67P/Churyumov-Gerasimenko (hereafter 67P) from the in situ measurements made near to the nucleus by the Rosetta Orbiter Spectrometer for Ion and Neutral Analysis (ROSINA) consisting of two mass spectrometers and a pressure sensor (Balsiger et al. 2007). This requires (1) the use of a coma gas model relating the gas flow variables at any point to a set of adjustable parameters meant to represent the nucleus surface gas production (i.e. defining the so-called flow "surface boundary conditions"), and (2) an iterative procedure to derive the best-fit values of these parameters. The coma gas model has two requirements: (1) compliance with the well-established rules of physical gasdynamics and (2) plausibility of the surface boundary conditions including the complex shape of the nucleus and the associated intricate illumination conditions. The inhomogeneous gas flow in the coma is characterized by the juxtaposition of regions with different flow regimes from free-molecular to fluid, the presence of multiple shocks, and non-equilibrium. The universal approach to handle this situation is the direct simulation Monte Carlo (DSMC) kinetic method (Bird 1994), but this may become very expensive in terms of computational time and required memory space. However, as has been demonstrated in many previous works, inviscid and/or viscous fluid methods (e.g. Euler equations, EE, or Navier-Stokes equations, NSE) can describe the flow in the coma with sufficient precision while preserving a high computational efficiency (e.g. Lukyanov et al. 2005).
The second aspect (the relevance of boundary conditions) consists in the underlying physical model and the way parameters for this model were derived. Complex, precise physical model with multiple parameters may have no benefit owing to the large uncertainty in the parameters magnitudes due to the lack of precise observational information. The other limit -an oversimplified model -is not worth the effort since it has no physical meaning.
In Bieler et al. (2015), both kinetic and fluid approaches were used (we do not consider the third, purely geometrical model) to compute the gas flow, assuming a pure water production from a homogeneous surface. At each surface point the gas flux and the surface temperature varied between imposed minimum and maximum values according to the cosine of the local solar zenith angle (i.e. gas production and surface temperature are decoupled). This model was able to reproduce the overall features of the local neutral number density measurements for the time period between early August 2014 and January 1, 2015, although some details in the measurements were not reproduced. Also, it should be noted that additional latitude-dependent correction factors were applied in a post-processing manner to the model outputs to improve the correlation to the observations. In Fougere et al. (2016a,b), the gas coma was simulated by the DSMC method and took into account the most abundant observed components: H 2 O and CO 2 or H 2 O, CO 2 , CO, and O 2 , respectively. The model of activity was similar to Bieler et al. (2015) but with the addition of surface inhomogeneity factors (for each species) to each surface element. The heterogeneity pattern of the surface outgassing activity was constrained by in situ coma density and composition measurements performed from August 2014 to the end of February 2016. A formal numerical data inversion method based on the extension of the geometrical model from Bieler et al. (2015) was used for the derivation of the surface fluxes. Agreement of the adjusted model and measurements was quantified by the root mean square deviation normalized by the mean measured value (for H 2 O and CO 2 , these are 1.14 and 1.08, respectively). It should be noted that (1) during a large part of this time period the Rosetta spacecraft was at distances greater than 100 km from the nucleus (where many flow structures are washed away) and (2) derivation of the surface properties based on long period measurements implicitly suggests that surface properties are invariable on this timescale (which is unlikely for more than one year).
In Marschall et al. (2016) the DSMC method was used assuming a pure H 2 O coma (though there is a high relative abundance of carbon dioxide in the ROSINA data set). The model of gas production assumed surface ice sublimation from an inhomogeneous surface. The distribution of surface inhomogeneity was postulated by a set of patches. The geometry of patches was inspired by the observed surface morphology (i.e. it is not connected explicitly with surface activity). Marschall et al. (2016) also investigated activity from cliffs versus plains, setting the difference between the two at a certain angle of the surface with respect to the local gravity field. The total number of defined regions on the surface of the nucleus is 26 and each is of a considerable size, therefore adjustment by patches is rather crude. The outgassing activity of patches at the surface was constrained by in situ coma density measurements performed from the end of August to the end of September 2014. It is worth noting that in this work simulation was done for the 10 km size region only. Beyond 10 km the simulated data were extrapolated assuming radial outflow with constant velocity.
In Kramer et al. (2017), a "simplified analytical model" (a kind of a collisionless expansion from point sources placed at shape model facets) was used for the definition of the gas coma parameters distribution. It is clear, however, that the gas coma was anything but collisionless. Furthermore, the model considered the time-averaged (over several weeks) surface emission rate and did not consider oscillations caused by momentary illumination conditions. The gas model of gas production was constructed by fitting tens of thousands of measurements to thousands of potential gas sources distributed across the entire nucleus surface. The relative difference with in situ density measurements was 11-18%.
It should be noted that the period of fitting was 2015-2016 and at this time Rosetta was at more than 100 km from the nucleus and therefore the numerous gas structures that are formed by collisions in the nucleus vicinity become in large part smoothed out. Therefore, the good agreement between an artificial collisionless model with quite homogeneous surface properties and smoothed-out structures are a natural consequence. However, the correlation to the footpoints of dust outbursts identified by the Optical, Spectroscopic, and Infrared Remote Imaging System (OSIRIS), even one year after the occurrence found in this work, is rather remarkable.
This short overview of currently published results on the distribution of surface activity shows that (1) the gas production mechanism is formulated in different ways; (2) the methods used to constrain the models are different; (3) the resulting patterns of surface activity are quite different and no one model fits perfectly the observational data for the close distances over a long period. Therefore, the derivation of surface activity still deserves attention.
In the present work we present in detail our own model and the results of its adjustment to the in situ measurements. We restrict ourselves to the so-called "prelanding data" collected in August-November 2014, which we had to interpret in real time to predict the aerodynamic force exerted on the lander during its descent (Jurado et al. 2016). We use both fluid and kinetic methods alternately in order to reduce computational expenses but keep the physical accuracy. The method based on fluid approach uses Euler equations and was described in Rodionov et al. (2002). The kinetic approach was implemented on the basis of the DSMC method. The same methods are used in the fitting procedure.

RZC model description
The observational nucleus data for the model were (a) a nucleus rotation model, and (b) a nucleus surface shape model. To these we added our coma gas flow model(s) and our surface production model.

Nucleus surface model
Many variants of the nucleus shape model were derived during the mission. In the present work we adopted the so-called "RMOC shape 3" (see Fig. 1), despite the fact that there are more precise shape models of 67P at the present time. This was done for the following reasons: (1) we started our work in 2014 and at that time it was the most precise model; (2) we study the period August-November 2014 when the Sun was in the northern hemisphere and the northern part of the shape was mapped with sufficient precision; (3) the more precise model would have needed to be degraded anyway in order to reduce the number of surface elements to match the requirements of the gas outflow solvers.

Nucleus rotation model
From the orbiter camera images, the nucleus rotation was accurately determined as a single-axis rotation with constant period (12.40 h) and invariant axis (OZ in the present figures) inclined in the period under discussion from 45 • to 55 • from the Sun direction.

Structure of the gas coma computational grid
This model is given as an unstructured triangular grid with a resolution of about 30 metres. In our computations we use a block-structured grid consisting of six blocks joined in a "football" manner (see Fig. 2 of this paper and Rodionov & Crifo (2006) for more details) and this grid covers the whole surface. Since the shape model is non-stellar-like, the grid generation procedure used in Rodionov & Crifo (2006) was modified. In the present version the initial point for radial directions is floating along the X-axis between −1 and 2 in order to ensure that the radial direction from the initial point crosses the surface only once. In addition, the resolution of the shape model is excessively fine for gas outflow simulations. Therefore, for computations we generated three versions S i of an "effective surface model" to be used by the gas code differing by the characteristic angular resolution of cells (1 o × 1 o , 2.5 o × 2.5 o and 5 o × 5 o ). This procedure evidently erases all small-scale details (see images in Fig. 1) but also creates a surface slightly different from the "real" one, which has, for example, a different shadow pattern under oblique solar illumination. Figure 1 shows a comparison of the initial shape model with the "effective surface model" of different resolution. The volume grid is constructed having the effective surface model as the inner boundary (ground level) and a sphere as the outer boundary (top level, see Fig. 2). The radial spacing of levels is non-uniform to ensure that the radial steps near the surface and near the outer boundary are comparable with inner and outer surface cell resolution, respectively. The model allows for the rotation by a succession of steady-state solutions separated by a constant time interval δt such that (1) during δt the computed surface temperatures and gas fluxes vary negligibly, and (2) the typical gas transit time across the computational domain is smaller than δt.

Nucleus gas production model
From Earth-based observations of comets, yielding (1) an approximate global production rate Q i of the most abundantly produced molecules and (2) a coarse image of the resulting dusty atmosphere -the so-called "coma" -it has long since been considered that the less volatile molecules (e.g. H 2 O) are produced from the surface of the nucleus by sublimation of (sub)surface ice inclusions and that the more volatile molecules are produced by diffusion through the surface after having been sublimated at variable depths inside the nucleus. It follows that the H 2 O emission must peak at sunlit points, while that of the other molecules need not be limited to such points and could even be nearly uniformly produced over the surface if deep production is involved. Thus, due to the nucleus rotation, the coma is expected to have two time-dependent components: one confined to the sunward hemisphere and the other co-rotating with the nucleus.
Considering the fact that inside each surface grid element there may be a complex topography and therefore temperature inhomogeneity, the distribution of ice may assume any pattern, and the CO and CO 2 production may also have stochastic variations, it was considered simply impossible to perform a gas-dynamically exact modelling of the near-surface nonequilibrium outflow region. Instead, our treatment was a trivial extension of that used in Rodionov et al. (2002) where stochastic fluctuations are absent: (1) The initial gas parameters T 0 , P 0 were related to the initial equilibrium flow Mach number M 0 by the plane-parallel homogeneous vapour sublimation relations derived by Cercignani (1981).
(2) Whether condensation or sublimation of H 2 O occurred was selected by the Riemann problem solution of the fluid EE or NSE, which also provided M 0 .
In the absence of reliable information on the interior of the comet nucleus surface, the number of parameters to introduce for a proper representation of the gas production is potentially very large. However, based on our in-depth analysis of the 1986 Halley comet fly-by data , we used a minimal number of parameters. For H 2 O we assume that each elementary surface area contains some fraction f H 2 O of exposed ice (≤1 and it should also be noted that since no large patches of ice were observed on the surface, this factor might be quite small). The upward H 2 O flux is therefore equal to f H 2 O F HK (T n ), where F HK is the ice Hertz-Knudsen sublimation rate, and T n is the surface ice temperature, computed from a slightly modified classical sunlit ice energy budget equation: where F H 2 O is the net sublimating gas flux, M 0 is the initial gas Mach number, L S is the latent heat of sublimation, ε is the surface infrared emissivity, σ is the Boltzmann constant, E in is incoming energy at instant t, w is a proportion of instantaneous and rotation averaged energy input, and T r is the period of rotation. The rotation-averaged term allows us to also consider the time lag of the heat transfer. The case in which w = 0 corresponds to an instantaneous adjustment of the surface temperature to the illumination (i.e. high thermal conductivity). The case in which w = 1 corresponds to the case when the thermal state is defined by the mean illumination (low thermal conductivity), where A V is the dirty ice visible albedo, z is the solar zenith angle, c is the solar flux, and κ is an heuristic dimensionless model parameter representing a heat input from the nucleus interior. Figure 3 shows an example of the input energy distribution over the surface in two limiting cases: purely instantaneous and rotation-averaged. The presence of shadowing of one element by the nucleus orography is taken into account by a ray-tracing method (without regard for reradiation For CO and CO 2 we adopt total production rates Q CO and Q CO 2 and assume that the fluxes F CO and F CO 2 are distributed over the surface according to laws of the kind F = a + b max[cos z , 0] where a and b are free parameters: where A ext is the total surface area, A sun is the sunlit crosssection, a 0 is the fraction of gas production distributed uniformly over the surface, and f sun is a shadow step function equal to 0 in shadow and to 1 otherwise. For CO and CO 2 , we also attribute to each elementary surface area two factors, f CO and f CO 2 , which affect the upward CO and CO 2 fluxes in the following way: f CO F CO and f CO 2 F CO 2 . The physical nature of these factors is different from the fraction of exposed ice f H 2 O and they could vary from 0 to greater than 1. For example, these factors could be interpreted as an effect of the insulating dust layer thickness over the surface. When f CO > 1 or f CO 2 > 1, it means that the corresponding total production rate Q CO or Q CO 2 was underestimated. In order to have f CO ≤ 1 and f CO 2 ≤ 1, the imposed total production rates should be multiplied by the peak values f max CO and f max CO 2 and at the same time the inhomogeneity factors of surface elements should be divided by these peak values. In the following, we will refer to f H 2 O , f CO 2 and f CO as "the factors of surface inhomogeneity". Figure 4 shows an example of the computed surface temperature for an icy and nonicy surface (for homogeneous surface, i.e. with constant f H 2 O , f CO , and f CO 2 ), for one orientation of the nucleus.

Gas outflow model
The approach we used is similar to that described in detail in Rodionov et al. (2002). However, special efforts had to be made to meet three essential constraints: (1) provide temporal and spatial resolutions comparable to those of the observations; (2) handle the coexistence in the computational domain of widely different rarefaction regimes; (3) minimize the central processing unit (CPU) time to enable iterative adjustment to the observational data.
To meet the mentioned constraints, it was necessary to develop three independent codes: one based on Euler equations (EE), one based on NSE, and one DSMC. The three codes assumed emission of the molecules H 2 O, CO, and CO 2 . All codes have a three-dimensional time-dependent (3D + t) capability, but could be solved under a quasi-steady approximation: at outflow velocities of the order of 100 m s −1 (extremely low estimate), the time for the gas to reach 5 km is of the order of 1 min, during which the nucleus rotates by 0.5 • , inducing sufficiently small surface flux and temperature changes.
The transport coefficients entering the NSE were derived from the Physics Handbook temperature-dependent heat conductivities and viscosities of the three molecules, using classical mixture laws. In EE and DSMC, it was assumed that all molecules are vibrationally relaxed, in view of the range of derived temperatures (<200 K). The Euler equations were solved by a shock-capturing second-order Godunov-type method in the whole computational domain, with Courant number equal to 0.4. The right-hand sides of the NSE were approximated by an explicit central difference algorithm. In the DSMC code, variable hard-sphere intermolecular collision potentials were used, adjusted to the viscosities adopted for the NSE. To describe translational-rotational energy exchanges, the Larsen-Borgnakke model was used.

Model adjustment to observational data
To adjust the model parameters, we have at our disposal the total production rate measurements by the Ultraviolet Imaging Spectrometer (ALICE) and the in situ measurements by the nude gauge (NG) of Comet Pressure Sensor (COPS) and a Double Focusing Mass Spectrometer (DFMS) of ROSINA. The COPS provides data on the overall density, and DFMS provides data on the relative abundance of molecular species.
The adjustment of the model is separated into two consecutive stages. In the first stage, we adjust the integral parameters: the total production rates and the proportion of the instantaneous and averaged energy inputs. In the second stage, we adjust the model to the in situ observational data. Figure 5 shows the trajectory of Rosetta in 1-15 of October 2014.
For a selected time period from the positions of the measurements (i.e. positions of the orbiter), we trace back the gas flowlines down the surface and determine from which surface cell these flowlines originated. For the surface cell we collect statistics on the ratios  where n mod and n NG are the overall gas densities at the position of the orbiter computed and measured by the nude gauge of COPS, c mod H 2 O , c mod CO , c mod CO 2 are the relative abundances of H 2 O, CO, and CO 2 from the simulation, and c DFMS H 2 O , c DFMS CO , c DFMS CO 2 are the relative abundances measured by DFMS.
Since the flux along a flowline is F = nv (where n and v are local gas density and velocity) and taking into account that v tends to const. with increasing nucleocentric distance, the flux at the positions of the orbiter is F ∝ n. Therefore, the mean ratio n mod /n NG shows if the flux from the element was in general under-or overestimated. The current (at iteration i) flux from the surface element is The new (i.e. for the next iteration) inhomogeneity factors f i+1 H 2 O , f i+1 CO , f i+1 CO 2 for each surface cell are computed from the equations where i is the index of iteration and the bar inr means averaging over flowlines originated from a given surface cell. As variations in the flux affect the flow in general it is necessary to repeat iteratively simulations for the whole rotational period with the new flux distribution.
Since observational data are limited, some of the surface elements may have no data (i.e. are empty cells). Therefore, it is necessary to expand the data from non-empty cells in order to cover the whole surface. This procedure consists in consecutive iterations in which, in the loop over empty cells, we assign for the cell the averaged data from the neighbouring non-empty cells. The iterations are repeated until all cells contain data. An example of this process is illustrated in Fig. 6.

Results
For most of our computations we used our coma grid with the resolution 2.5 • × 2.5 • , which gives 7776 surface cells.
We have fitted the data acquired in August-November 2014 by the ROSINA instrument (pressure gauge COPS and DFMS), which are described in detail in Le Roy et al. (2015). In the course of this period, the heliocentric distance changed from 3.62 to 2.87 AU. In order to take into account accurately the variation of gas production, we divided this period into ten-day segments, which corresponds to ∼2% variation of r h . For each time segment, we performed simulations of the gas outflow (with corresponding r h ) for one rotational period split into 72 instants (i.e. 10-min intervals).
During the period August-November 2014, the COPS made about 125 000 measurements (cleaned from the firing events) and the DFMS made about 6000 measurements. For adjustment we used time instances corresponding to the COPS measurements and the DFMS data were linearly interpolated (this implies that the composition of the gas mixture does not vary significantly on shorter timescales). In one time segment, we had typically about 10 000 measurements. Since for a grid with resolution 2.5 • × 2.5 • the number of surface elements is 7776, for one 10-day time segment the number of samples per surface element is very poor and a large part of the surface may not be covered at all. To improve the statistics, we usually combined several successive time segments. This implicitly assumes that inhomogeneity factors do not depend on the variation of r h . We started from the assumption of a homogeneous surface (i.e. inhomogeneity factors were the same all over the surface) with two types of energy input: instantaneous (w = 0) and completely rotation-averaged (w = 1). The H 2 O production was assumed to be set by an icy fraction f H 2 O = 0.033, and the total CO and CO 2 production rates were, respectively, 2.0 × 10 25 and 3.0 × 10 24 molec s −1 (with f CO = 1 and f CO 2 = 1). Figure A.1 shows the comparison of the simulated gas density and velocity distribution. In the case in which w = 1, owing to the more uniform income energy distribution over the surface, the bottom of the "neck" region is relatively more active and the flow field is more symmetric with respect to the axis of rotation Z. Figure A.2 shows a comparison of the gas density in the position of the Rosetta orbiter predicted by models and COPS measurements for the period August-November 2014, with zoom on the period October 8-12. The model with averaged energy input shows less amplitude of density variation during the rotation period and it corresponds better to the measured data when the orbiter passes above the southern hemisphere (sub-spacecraft latitude, Z < 0).
To numerically quantify the agreement between the model and the observations, we computed the mean (δ x ) and the dispersion (σ x ) of relative difference (∆ i x ) in the simulated and measured data x: where N is a number of points in the segment of comparison, and x mod i , x obs i are comparable parameters (n, c H 2 O , c CO 2 , c CO ) in our simulations and in the measurements. The results for the models with homogeneous nucleus when w = 0 and w = 1 are given in Table 1 (cases f033c1w00 and f033c1w10, respectively). According to these criteria, the model with w = 1 practically fits the measured densities two times better than the model in which w = 0. Figure A.3 shows an example of the sample volume distribution for surface elements for the period August-November 2014. As expected in the southern (Z < 0) hemisphere, statistics are very poor. In the northern (Z > 0) hemisphere, the number of samples per surface element reaches 100.
The distribution of the f H 2 O , f CO , and f CO 2 after two iterations of adaptation (with gas outflow simulation based on Euler equations) is shown in Fig. A.4. The EE method was used for the initial iterations of adjustment due to its high computational efficiency and acceptable validity in the prediction of the spatial gas density distribution (see Lukyanov et al. 2005;Crifo et al. 2016). Figure A.5 shows a comparison of the measured and simulated gas densities (after adaptation of the surface inhomogeneity factors). The mean and dispersion of relative difference in simulated and measured gas densities are given in Table 1 (cases f002a1w00 and f002a1w10). After adaptation of the surface inhomogeneity factors, the model in which w = 0 improves the agreement with observations on gas density nearly twofold. The model in which w = 1 produced the same agreement in the overall density but very much improved agreement for the relative abundance of CO.
For further adjustment, we used only the model in which w = 0 and for the gas outflow simulations we applied the DSMC method. For the last iterations we used the DSMC method, which is less computationally efficient than the EE but more physically adequate. In order to reduce computational expenses, for the adjustment we restricted the period of study to September 15-October 30, 2014. The distribution of the f H 2 O , f CO , and f CO 2 after two additional iterations of adaptation is shown in Fig. A.6. Additionally, for this distribution of inhomogeneity factors, the gas outflow was computed by the EE method as well. Figure A.7 shows the distributions of gas density and velocity simulated by the EE and the DSMC methods. The results of the density simulation by gas-dynamic and gas-kinetic codes are in a reasonable agreement but the velocity distribution near the nucleus differs noticeably. Figure A.8 shows a comparison of the gas density in the position of the Rosetta orbiter predicted by models based on the Euler equations and the DSMC approach with instantaneous energy input and COPS measurements for the period September 15-October 30 and the period October 8-12 after adaptation of the surface inhomogeneity factors. The mean and dispersion of the relative difference in simulated and measured gas density are given in Table 1 (cases f004a3w00 for the EE method and m004a3w00 for the DSMC). From these criteria, it follows that for the used inhomogeneity pattern, both methods (EE and DSMC) are almost equally close to the measured density but the DSMC gives better agreement with relative abundances. The density variation predicted by the EE method has more sharp gradients than the DSMC. In both predictions there is an extra or split peak at about 3 h on October 10. Figure A.9 shows a comparison of the relative abundances of H 2 O, CO, and CO 2 in the position of the Rosetta orbiter predicted by the DSMC model after four iterations of the surface inhomogeneity factors adjustment for the period from September 15, 2014to October 30, 2014and during October 8-12, 2014. In the period September 15-October 30, the model initially underestimates the relative abundance of CO but then overestimates it at the approach to the end of the time segment. This is a consequence of the constant total production rate of CO postulated in the model. Figure A.10 shows the variation of total gas production during one rotation period at r h = 3.21 AU. Dashed lines show gas production for the homogeneous nucleus, solid lines show gas production after adjustment of the surface inhomogeneity factors. We stop the adjustment procedure after four iterations when the mean relative difference between simulated and measured parameters is about 30%. This was caused by the following reasons: (1) each iteration of the adjustment procedure needs a considerable amount of computations (13 time segments × 72 time instants); (2) since we started our work, the shape model used in our simulations has improved (especially in the southern hemisphere) and this may affect the gas distribution; (3) the A71, page 7 of 16 A&A 618, A71 (2018) Notes. Here method denotes the method of gas flow description (Euler equations or DSMC ), w is the proportion of instantaneous and rotationaveraged energy input (see Eq. (1)), it is the iteration number of inhomogeneity factors adaptation, δ n is the mean, σ n is the dispersion of relative differences in the simulated and measured overall gas densities, and δ H 2 O , δ CO 2 , δ CO are the means and σ n , σ H 2 O , σ CO 2 , and σ CO are the dispersions of relative differences in the simulated and measured relative abundances of H 2 O, CO, and CO 2 at the position of the orbiter (see Eq. (8)).
measurements have error bars of the order of 10-20%. Therefore, we decided that further adjustment is not worth the effort.

Conclusion
In the case of nuclei with homogeneous surface (i.e. inhomogeneity factors are the same all over the surface), the model with energy input averaged over period of rotation (w = 1) better fits the ROSINA data than the model with instantaneous energy input (w = 0): the mean relative difference in total gas density is about 50% versus about 100%. The amplitude of density variation is less in the model in which w = 1 since it depends only on the instrument location in the coma (the gas production remains practically constant during the rotation). As follows from a distribution of the surface temperature (Fig. 4), in the model in which w = 1 more active surface (due to a higher temperature) concentrates on top of the shape, as the part which is exposed for the longest time (this behaviour is consistent with the results of a further adjustment of surface inhomogeneity factors). However, the adjustment of surface inhomogeneity in the case of w = 1 does not improve the agreement with measured data and this means that gas production should be essentially dependent on instantaneous illumination conditions. Adaptation of surface inhomogeneity strongly improves the fit of the model in which w = 0. Final (after four iterations of adaptation) distribution of the surface inhomogeneity factors shows that for all components (H 2 O, CO, and CO 2 ) the most productive surface is located around the northern (Z > 0) peaks and crests of the shape. We emphasize the difference between "more active" and "more productive" surfaces. The more active surface has a greater production rate per square metre in comparison with other surfaces (possibly having other surface properties or/and having other illumination conditions).
The more productive surface has greater production rate per square metre in comparison with other surface at the same illumination conditions.
The final distribution of the surface inhomogeneity (Fig. A.6) is quite different from the results presented in Fougere et al. (2016a,b) and Marschall et al. (2016). These studies came to the conclusion that the most productive surface is located in the neck region of the shape. According to our results, the neck is also relatively highly productive (only for H 2 O) but the most productive areas (for all species) are located around the northern (Z > 0) peaks and crests of the shape. The reason for this discrepancy could be that the flow from the bottom of the neck region and the flow formed by the interaction of fluxes from the lobes (see Fig. A.7) may result in a similar density in the points of measurements, that is, it is possible to fit the ROSINA data by several models of gas production. This naturally brings us to the necessity of fitting to the observational data from other Rosetta instruments as well (and this work is currently in progress).
The results of our analysis show that in the considered time span the outgassing is defined mostly by instantaneous illumination conditions (which are strongly dependent on the nucleus shape) and local properties (inhomogeneity) of the surface. At the same time, owing to the rather rarefied conditions of gas outflow (and therefore fast dissipation of flow structures), smallscale topographical features do not affect the flow noticeably for the present in situ observations (made at altitudes of 10-100 km). In the considered conditions, "small-scale" means comparable with the mean free path of the molecules near the surface, that is, on the order of 100 m. However, due to the lack of information about the near-nucleus gas coma of other comets, it is not without problems to generalize these results to all comets.