Keywords: asteroids, shapes, spins

We introduce techniques for characterizing convex shape models of asteroids with a small number of parameters, and apply these techniques to a set of 87 models from convex inversion. We present three different approaches for determining the overall dimensions of an asteroid. With the first technique, we measured the dimensions of the shapes in the direction of the rotation axis and in the equatorial plane and with the two other techniques, we derived the best-fit ellipsoid. We also computed the inertia matrix of the model shape to test how well it represents the target asteroid, i.e., to find indications of possible non-convex features or albedo variegation, which the convex shape model cannot reproduce. We used shape models for 87 asteroids to perform statistical analyses and to study dependencies between shape and rotation period, size, and taxonomic type. We detected correlations, but more data are required, especially on small and large objects, as well as slow and fast rotators, to reach a more thorough understanding about the dependencies. Results show, e.g., that convex models of asteroids are not that far from ellipsoids in root-mean-square sense, even though clearly irregular features are present. We also present new spin and shape solutions for asteroids (31) Euphrosyne, (54) Alexandra, (79) Eurynome, (93) Minerva, (130) Elektra, (376) Geometria, (471) Papagena, and (776) Berbericia. We used a socalled semi-statistical approach to obtain a set of possible spin state solutions. The number of solutions depends on the abundancy of the data, which for Eurynome, Elektra, and Geometria was extensive enough for determining an unambiguous spin and shape solution. Data of Euphrosyne, on the other hand, provided a wide distribution of possible spin solutions, whereas the rest of the targets have two or three possible solutions.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 1 Introduction

Traditionally, asteroid shapes have been modelled using triaxial ellipsoids, which are completely defined by only two parameters, namely, the ratios of the lengths of the principal axes.

Lately, complex shape models, such as arbitrary polyhedra and even non-convex shapes, have become common, raising a question about how to describe the characteristic features of the models. One of the lightcurve inversion methods developed in the 20 th century is the convex inversion method (Kaasalainen andTorppa 2001 andKaasalainen et al. 2001, hereafter KT01 andK01, respectively), which offers an arbitrary polyhedron shape model. Also, radar and even spacecraft-based non-convex models have become more and more common (e.g., [START_REF] Magri | Radar observations and a physical model of Asteroid 1580 Betulia[END_REF]Demura et al. 2006, respectively). Although these models are naturally much more informative than ellipsoids, their weakness is that there is no way to express the exact shape solution in a simple numerical form due to the large number of parameters. We can, however, characterize the shape using various quantities, and grasp the most important shape features with a small number of parameters, which is required, e.g., when carrying out statistical analyses of the shape characteristics.

Previous studies considering spin states and shapes of asteroids in statistical sense have used triaxial ellipsoid shape models. Ellipsoid dimensions of 22 main-belt asteroids were determined by [START_REF] Magnusson | Spin vectors of 22 large asteroids[END_REF] who conclude, based on the triaxial models, that it is not likely that centrifugal forces have affected the shapes. On the other hand, [START_REF] Drummond | Photometric geodesy of main-belt asteroids. IV -an updated analysis of lightcurves for poles, periods, and shapes[END_REF] studied the axis distributions of 26 asteroids concluding that shapes of some of the objects may be jacobi ellipsoids, i.e., the stretching of the ellipsoids is due to the rotation of the body. [START_REF] Hartmann | Trojan and Hilda asteroid lightcurves. I -Anomalously elongated shapes among Trojans (and Hildas?)[END_REF] compared the axis ratios of a sample of 26 Trojan and Hilda asteroids to those of main-belt asteroids. They suggest, that Trojans and Hildas are more elongated that main-belt asteroids of comparable size.

Shape characterization is an important part in many branches of science, e.g., in small-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

particle analysis for which commercial image analyzers have already been developed. The image of the particle is often projected to a plane and, hence, the parameters describe the two-dimensional projection of the shape. These parameters usually include the perimeter, area, and maximum and minimum dimensions which are compared to those of a circle with the same area or circumference, or to the convex hull of the particle to define its roundedness, convexity, or dimensions. One example of small particle shape determination is the shape analysis of lunar dust particles as described in [START_REF] Liu | Shape analysis of lunar dust particles for astronaut toxicological studies[END_REF], who used 2-D images of the sample particles to define their areas and perimeters, and the circular-equivalent diameter from the measured areas. They also fit ellipses to the measured areas in the least-squares sense. Fractal analysis for particle shape characterization has been applied by [START_REF] Carey | Using fractal analysis to quantitatively characterize the shapes of volcanic particles[END_REF] when analyzing the shapes of volcanic particles. [START_REF] Muinonen | Inversion of small-particle silhouettes for Gaussian-sphere parameters[END_REF] presents a method to derive statistical particle shape paramters from a set of 2-D particle silhouettes using a Gaussian shape model consisting of two parameters. A method to characterize macroscopic objects is presented, e.g., by Torppa et al. (2006), who determined the characteristic shape parameters of potato tubers. They used 3-D image data to find the ellipsoid that best fits the tuber shape. They also used the residual radii of the ellipsoid fits in a statistical sense to determine the variance of the Gaussian distribution describing the local irregular surface features. Gaussian statistics have been applied to asteroid shape characterization also by [START_REF] Muinonen | Inversion of shape statistics for small solar system bodies[END_REF], who determined the parameters for the Gaussian random sphere that best mimics asteroids and other small bodies.

In this study, we introduce and compare techniques to express the overall dimensions and irregularity of an arbitrary convex polyhedron and apply the techniques to a set of 87 shape models for studying the statistics of basic shape features of asteroids. We also present new spin and shape results for eight asteroids of which new lightcurve observations have been carried out.

In Secs. 2 and 3, we describe the methods used in lightcurve inversion and characteri-

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
zation of arbitrary shapes. Section 4 is dedicated to results and Sec. 5 to conclusions and future prospects.

Convex inversion for spins and shapes

Spin and shape analysis for selected targets with new photometric observations was carried out using the convex inversion method. Theory of the method is described in detail in KT01 and K01, and is not repeated here. In the following, however, we consider a few developments.

One fact that has not been discussed in the earlier papers is the number of parameters (i.e., facets) in the Minkowski minimization, which is used for determining the polyhedron vertices from the set of facet normals. It has been mentioned in previous papers (e.g., KT01) that, when searching for the best-fit facet areas, the number of facets in the polyhedron shape model has to be of order 800 to guarantee that the shape solution does not depend on the choice of the facet normal directions. In the next step of finding the vertices of the shape from the facet information, however, the number of facets can be significantly reduced.

Without affecting the final shape, this recombination significantly speeds up the Minkowski minimization which previously was the most time consuming part of the procedure. We found it most convenient to combine the facets around each vertex determined from the standard triangulation. Normal directions of the combined facets are calculated by taking a weighted mean of the original normal directions. Due to the way of discretization and combination, some facets remain uncombined, but this does not affect the final shape solution.

Another change to the inversion procedure was that, from the beginning, we aimed at a distribution of spin solutions, not at a single unambiguous solution, since the data sets were not as extensive as in the previous analyses (e.g., [START_REF] Torppa | Shapes and rotational properties of thirty asteroids from photometric data[END_REF]hereafter T03). We call this technique semi-statistical, since we obtain the distribution of possible solutions, but

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
not the probability density. As before, we carried out the initial spin state determination by scanning through a set of spin axis directions and rotation periods, while fitting a low-order functional series shape model. The time step of the period scanning was ΔP/4, where

ΔP = P 2 2T , (1)
P being the period, and T the total time span of the observations. This equation is based on the fact that, if P is changed by ΔP , the lightcurve is phase-shifted by π during T .

Thus, if we have model and data lightcurves with two equal maxima and minima (typical for an ellipsoid) of which the first and the last coincide for period P , then they coincide also for period P + ΔP . Although the data lightcurves usually differ significantly from those produced by an ellipsoid, this equation gives a rough estimate about how densely the χ 2 minima are located in the χ 2 , P -space, and how densely P should be sampled. Spin axis was sampled in eight uniformly distributed directions. From the initial spin state search, we accepted all the solutions with rms < min(rms) × 1.5, which is an empirical criterion for finding the final solution or solutions from the set of initial ones. Starting from the set of initial solutions, we refined each by adding more shape parameters (a higher-order functional series) and setting the spin axis direction and rotation period as free parameters. For the improved solutions, we finally applied the standard Gaussian distribution empirical rule for 3-sigma deviation as the criterion of acceptance. Due to the effects of the uncertainties in the shape and scattering models, an error of 5 • should be added to each spin longitude and 15 • to each latitude. This uncertainty is an empirical value, based on this study and earlier applications (e.g., T03).

The last aspect to consider was the sampling density of the data. Since some of the lightcurves are over-densely sampled, while others are sparse, the weight of each data point could be set inversely proportional to the lightcurve sampling density in order to prevent

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
dense lightcurves from dominating in the χ 2 computation. The function to be minimized would then be

χ 2 = i (L d,i -L m,i ) 2 σ 2 i , ( 2
)
where L d,i and L m,i are the corresponding data and model brightnesses.

σ 2 i is proportional
to the number of data points in and inversely proportional to the overall duration of the i th lightcurve.

Shape characterization

The definition of dimensions of convex inversion models is not straightforward, since they are irregular, and not necessarily in the principal-axis rotation state. In the previous papers,

where convex inversion has beed applied to asteroid lightcurve data (e.g., Kaasalainen et al., 2002a and T03), descriptive axis ratios for the model shapes are given as the mean of two measures: the axis ratios computed from the largest dimensions in the rotation axis direction and in the equatorial plane (hereafter the OD technique, section 3.1 of this paper), and the axis ratios of the best-fit ellipsoid as defined in section 3.2 of this paper (hereafter the E1 technique). There are, however, also many other ways to define the dimensions. In this paper, we compare four approaches to characterize the shape of an arbitrary convex polyhedron. Before presenting the techniques, we discuss the rotational state of the shape models, since the models are rarely in the exact principal-axis rotation state, although the target asteroids with only one lightcurve period should be. When finding the characteristic shape features, it is not clear whether the body should be rotated to the principal-axis state.

Here we include both approaches: methods that fix one axis to be aligned with the rotation axis and ones that allow the coordinate system to be rotated arbitrarily. In the first two methods, i.e., the one finding the overall dimensions of the shape (the OD technique) and

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
the one defining the best-fit ellipsoid in the sense of facet distances (the E1 technique), the model shape is not rotated to the principal-axis state. The third technique, that gives the best-fit ellipsoid in the sense of radii lengths (hereafter the E2 technique, section 3.3) the shortest axis is not forced to be aligned with the rotation axis of the model. Thus, only the shape is fitted, disregarding the rotational state. The fourth technique, computation of the principal moments of inertia (hereafter the PMI technique, section 3.4), also considers only the shape without fixing the original rotation state.

Overall dimensions (OD)

The overall dimensions of the shape were computed in three stages (Fig. 1). First, the largest extent in z-axis (spin axis) direction was calculated from

2c = max(z i ) -min(z i ), (3)
where z i are the z coordinates of all the polyhedron vertices. Second, the largest extent of projection to xy-plane was found from 2a = max(x i cos φ ax + y i sin φ ax )min(x i cos φ ax + y i sin φ ax ), (4) where x i and y i are the values of the x and y coordinates of all the polyhedron vertices.

φ ax is the azimuthal angle between x-axis and the trial direction in which the extent of the projection was calculated. φ ax was sampled with one degree step over the range [0

• , 180 • ].
Last, the largest extent, 2b, of projection to xy-plane, perpendicular to direction of a, was computed analogously with 2a. F i g u r e 1

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Ellipsoid I (E1)

With the first of the ellipsoid-based techniques, we compare the support functions of the model and the ellipsoid. The support function is defined as

ρ(θ, φ) = n(θ, φ) • r(θ, φ), ( 5
)
where θ and φ are the spherical coordinates, n is the surface normal, and r is the radius vector. Thus, ρ(θ, φ) is the distance of the tangent plane at point (θ, φ) from the origin. As in [START_REF] Kaasalainen | Interpretation of lightcurves of atmosphere less bodies. I. General theory and new inversion schemes[END_REF], we calculated the difference

d f = i |ρ m,i -ρ e,i | ρ e,i g i , ( 6
)
where g i is the area of the corresponding facet element on a unit sphere, ρ m,i is the distance of the model shape's facet i from the origin, and ρ e,i is the distance of corresponding ellipsoid facet (with the same surface normal) from the origin.

We chose the orientation (origin and z-axis) of the ellipsoid to coincide with that of the model. The ellipsoid axis ratios were searched by setting the shortest (polar) axis length c=1, and sampling the azimuthal axis lengths a, b within the range [1,3], with a step δ = 0.1.

Axes were then scaled so that the ellipsoid and the model were of the same volume. The azimuthal rotation φ ax of the longest ellipsoid a-axis with respect to the x-axis was sampled over the range [0 • , 90 • ] with a step of δ = 10 • .

Ellipsoid II (E2)

Another way of finding the characteristic ellipsoid dimensions is to compare the radii of the model shape and the ellipsoid. We did not force ellipsoid z-axis to coincide with model's rotation axis or ellipsoid origin with that of the model, and minimized function χ 2 e of nine

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
free parameters: semi-axis of the ellipsoid (a, b, c), transformation of the origin of the model shape (r c ) to coincide with the origin of the ellipsoid, and rotation of the model shape (α, β, and γ) for its axes to coincide with the principal axes of the ellipsoid:

χ 2 e = i (r m,i -r e,i ) 2 r2 m,i , ( 7
)
where r m,i is the model shape's i th radius in a system whose origin translated by r c and which is rotated by Euler angles α, β, and γ, and rm,i is the mean of the model radii. r e,i are the ellipsoid radii calculated for the same vertex directions as the model radii r m,i . The minimization was carried out using Nelder and Mead's simplex method described in Press et al. (1994).

Principal moments of inertia (PMI)

The inertia tensor of an arbitrary shape can be calculated by summing up the moments and products of inertia for different parts of the shape. In this case, we found it most practical to divide each shape into polygonal cones and furthermore into tetrahedra so that the total inertia matrix is obtained from

I jk = i t jk,i , (8)
where t jk,i is the moment (j = k) or product (j = k) of inertia of a single tetrahedron.

j, k = [1, 3] correspond to the x, y, and z coordinates and i stands for the i th tetrahedron.

For calculating t i , we followed the formulas given by [START_REF] Tonon | Explicit exact formulas for the 3-D tetrahedron inertia tensor in terms of its vertex coordinates[END_REF]. The principal moments and principal axes were obtained by diagonalizing the inertia tensor I using Jacobi transformations (Press et al., 1994).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT 4 Results

New spin and shape solutions

New spin and shape solutions are presented for eight asteroids, for which new lightcurves have been observed during the past few years. The new data are summarized in Table 1 and lightcurves are shown in Fig. 2. The data are stored in SAPC (Standard Asteroid Photometric Catalog, http://www.astro.helsinki.fi/SAPC) which is the updated version of the Uppsala Asteroid Photometric Catalogue (UAPC, [START_REF] Lagerkvist | Asteroid Photometric Catalogue, 5 th update[END_REF]. We also used previously published data that is stored in UAPC and SAPC. We considered all the lightcurves as relative photometry, since this is adequate for spin and shape determination, and the number of absolute observations was too small for studying phase-curve behaviour. Table 1 Figure 2 For (79) Eurynome, (130) Elektra, and (376) Geometria, we obtained an unambiguous spin and shape solution, for (54) Alexandra, (93) Minerva, (471) Papagena, and (776) Berbericia a few possible spin solutions and, for (31) Euphrosyne, a constrained distribution of spin solutions. As mentioned above, errors of 5 • (λ) and 15 • (β) apply to each spin solution due to inaccuracies in the model shape and scattering law. The uncertainty of the period depends on the period itself and on the total time span of the observations, as explained in connection to Eq 1. The number of digits in the rotation periods corresponds to their accuracy in Table 2, where all the spin solutions are summarized. Some of the targets show double solutions, implying two solutions with the same rotation period and pole latitude, and pole longitudes roughly 180 • apart.

Table 2 The accuracy of spin solution is usually considered to depend on the number of apparitions but at least in this case, where the number of apparitions was five or more in every case, it did not correlate with the period accuracy. Neither did the total time span of the observations which was of the same order for all the targets. We, however, found a correlation between the number of lightcurves and the accuracy of the period solution. In general, the period

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
uncertainty for each spin solution is ΔP from Eq. 1. For all eight targets of this paper, ΔP for a single spin solution was about 0.00001. However, for a majority of the targets, there were more than one spin solution, and thus the real period uncertainty was larger than ΔP . The plot in Fig. 3 shows how the period uncertainty depends on the number of lightcurves. It seems that, when the number of lightcurves is around 10-13, the period uncertainty decreases sharply, then reaching more slowly the level of ΔP , which is reached at around 30 lightcurves. This, of course, is highly sensitive to the quality of the lightcurves, but for a set of lightcurves covering close to one rotation, this kind of behaviour might be common in lightcurve inversion. A similar occurrence is present in asteroid orbit computation, where the uncertainty of the orbital parameters decreases sharply as the observational time interval increases to a certain point (e.g., [START_REF] Muinonen | Asteroid orbits using phase-space volumes of variation[END_REF].

The number of lightcurves alone, however, is not a sufficient measure in the cases where the total time span of observations is short or the number of apparitions is low and we do not consider the case here. However, the uncertainty does not increase to infinity for limited data sets, but settles to a few percent of the rotation period for a single lightcurve covering one revolution. This is based on the period uncertainties obtained for targets with only a few lightcurves in [START_REF] Muinonen | Spins, shapes, and orbits for potentially hazardous near-Earth objects by NEON[END_REF]. It is left for future studies to find out if the sharp increase in period uncertainty is typical for lightcurve inversion, or was it just due to the relatively small number (eight) of sample objects. Figure 3 For comparison to our results, there is a summary of previous solutions in the Database of asteroid spin vector determinations (http://vesta.astro.amu.edu.pl/Science/Asteroids/).

The database was originally created by Per Magnusson, at that time working at the Uppsala Observatory (Sweden), and is currently maintained by Agnieszka Kryszczynska at the Poznan Observatory. Recent results of some targets included in this paper are found in [START_REF] Durech | Physical models of ten asteroids from an observers' collaboration network[END_REF]. Differences to the results in other publications are due to the use of different modelling techniques. Naturally, it enhances the reliability of the results if they are in

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
agreement with the ones obtained using other methods. On the other hand, if the results are different from others, further observations are needed to find the correct solution.

For the asteroid sizes, there are varying estimates in the literature. Our data for sizes and spectral classes have been collected from UAPC [START_REF] Lagerkvist | Asteroid Photometric Catalogue, 5 th update[END_REF].

We assessed the need for weighting lightcurves according to the sampling density of lightcurve points (Eq. 2). However, adding weights had no effect to the solutions and, to be consistent with earlier solutions, we omitted them. This analysis showed again the stability of the convex inversion method. However, there may be cases, where the densely observed lightcurves dominate so strongly that they prevent finding a good fit for sparsely observed curves of high quality. In such cases, applying weights may stabilize the problem.

(31) Euphrosyne, C, 247 km

We used ten lightcurves, one of which was new (Table 1). The other lightcurves were published in [START_REF] Kryszczynska | CCD Photometry of Seven Asteroids: New Spin Axis and Shape Determinations[END_REF], [START_REF] Barucci | Rotational properties of ten main belt asteroids -Analysis of the results obtained by photoelectric photometry[END_REF], [START_REF] Mccheyne | Visible and near-infrared lightcurves of eight asteroids[END_REF], and [START_REF] Schober | The remaining large minor planets with unknown rotational properties -31 Euphrosyne and 65 Cybele[END_REF]. The data spanned a total of 30 years. This is an example of a target with data sufficient to provide only a clearly restricted number of solutions instead of an unambiguous solution, although the amount of data is extensive (Fig. 4). The (P,rms)-plot

shows that the best fits to the data are achieved with the period range of P = 5.5300-5.5340

h, but periods up to 5.5375 h are possible. The (λ, β) plot shows that two distinct spin axis ranges, one around λ = 200 • , β = 0 • and another one around λ = 330 • , β = -30 • . Plotting λ and β against rms did not help to choose the correct pole and further observations should be planned according to the (λ, β)-plot to eliminate spurious solutions. In the Poznan database, there are also a wide range of solutions for this target none of which is, however, in agreement with our solution. Figure 4 (54) Alexandra, C, 160 × 135 km

We used a data set of 20 lightcurves, two of which were not published before (Table 1).

The other data were from [START_REF] Belskaya | Physical studies of asteroids XXVII: Photoelectric photometry of asteroids 14 Irene, 54 Alexandra and 56 Melete[END_REF], [START_REF] Tancredi | A comparison of two pole determination methods for asteroids[END_REF], Di Martino et al. (1987), [START_REF] Haupt | Mitt.der Univ. Sternwarte[END_REF], and [START_REF] Van Houten-Groeneveld | Photoelectric photometry of seven asteroids[END_REF]. The data spanned a total of 42 years. We obtained three possible solutions for the spin state, namely 1) λ = 307 • , β = 20 • , and P = 7.02263 h, 2) λ = 122 • , β = -36 • , and P = 7.02334 h, and 3) λ = 325 • , β = -37 • , and P = 7.02334 h, the latter two constituting a double solution.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
The shape solutions corresponding to the first and the double solution were different, and are not displayed here. There is a previous spin solution in the Poznan database, namely λ = 160/290 • , β = 50 • , and P = 7.026384 h, which is different from ours.

(79) Eurynome, S, 68 km

We used 33 lightcurves, one of which was not published before (Table 1). The other data were from [START_REF] Michalowski | Photoelectric photometry, parameters of rotation and shapes of asteroids 22 Kalliope and 79 Eurynome[END_REF], Di [START_REF] Martino | Photoelectric photometry of 14 asteroids[END_REF], [START_REF] Angelis | Lightcurves and pole determinations for the asteroids 69 Hesperia, 79 Eurynome and 852 Wladilena[END_REF], [START_REF] Scaltriti | A photometric study of the minor planets 192 Nausikaa and 79 Eurynome[END_REF], [START_REF] Schober | The period of rotation and the photoelectric lightcurve of the minor planet 79 Eurynome[END_REF], and [START_REF] Shevchenko | Magnitude-phase dependencies for three asteroids[END_REF]. The total time span of the observations was 32 years. An unambiguous spin solution was obtained at λ = 64 • , β = 15 • , and P = 5.97857 h. The shape solution is shown in Fig. 5. Our pole solution lies within ten degrees of one of the solutions in the Poznan database. Figure 5 (93) Minerva, C, 150 km 14 lightcurves were used, five of which had not been published before (Table 1). The rest of the data are published in [START_REF] Millis | The occultation of AG+29 • 398 by 93[END_REF], Debehogne et al. (1982), and Denchev We used 58 lightcurves of Elektra, four of which were new (Table 1), and the rest are published in [START_REF] Danforth | CCD Photometry of 130 Elektra[END_REF] Figure 7 (376) Geometria, S, 36 km

We used 35 lightcurves of Geometria, one of which was new. The old data were published in [START_REF] Kryszczynska | CCD Photometry of Seven Asteroids: New Spin Axis and Shape Determinations[END_REF], [START_REF] Barucci | Rotational rates of very small asteroids -124 Brunhild, 376 Geometria, 437 Rhodia, and 1224 Fantasia[END_REF] 8

. The pole longitude is similar to earlier solutions, but the latitude is different. Figure 8 (471) Papagena, S, 140 km

We used 13 lightcurves, one of which was new (Table 1). The older data has been published in [START_REF] Lustig | The period of rotation and the photoelectric lightcurve of the minor planet 471 Papagena[END_REF], Di [START_REF] Martino | Photoelectric photometry of 14 asteroids[END_REF], Scaltriti and Zappalà (1978),

and [START_REF] Surdej | Rotation period and photoelectric lightcurves of Asteroid 471 Papagena[END_REF]. The total time span of the observations was 31 years. We

Characterization of shape solutions

We used a set of 87 convex inversion shape solutions, five of which are published in this paper, and the rest in K01, Kaasalainen et al. (2002a,b;2003[START_REF] Tonon | Explicit exact formulas for the 3-D tetrahedron inertia tensor in terms of its vertex coordinates[END_REF], [START_REF] Slivan | Spin vectors in the Koronis family: Comprehensive results from two independent analyzers of 213 rotation lightcurves[END_REF], T03, [START_REF] Michalowski | Photometry and models of selected main-belt asteroids (I): 52 Europa, 115 Thyra, and 382 Dodona[END_REF][START_REF] Michalowski | Photometry and models of selected main-belt asteroids (II): 173 Ino, 376 Geometria, and 451 Patientia[END_REF][START_REF] A C C E P T E D M A N U S C R I P T Accepted Manuscript Michalowski | Photometry and models of selected main-belt asteroids (III): 283 Emma, 665 Sabine, and 690 Wratislavia[END_REF], [START_REF] Marchis | Shape, size and multiplicity of main-belt asteroids[END_REF], and [START_REF] Durech | Physical models of ten asteroids from an observers' collaboration network[END_REF]. Also the Database of Asteroid Models from Inversion Techniques (http://astro.troja.mff.cuni.cz/projects/asteroid was used in collecting the data. In the following, we discuss the ranges of axis ratios as well as correlations between shape, size, rotation period, and taxonomic type. It has to be kept in mind that, especially for non-elongated objects, the vertical dimension of the model shapes is somewhat ill-defined when using only relative photometric observations. Otherwise, the features and dimensions of the shapes have been shown to be stable, e.g., against using various scattering laws in lightcurve inversion (K01). Taxonomic types are defined according to the SMASS classification [START_REF] Bus | Phase II of the small main-belt asteroid spectroscopic survey. A feature based taxonomy[END_REF]. Most of the asteroids included in this study belong to the three main categories C, S, and X, and two are determined as D-type.

Properties of all the sample asteroids are collected in Table 4.

Axis ratios

If we define the overall dimensions and axis lengths of the characteristic ellipsoids as c < b < a, the corresponding principal moments are

C c > C b > C a .
To make comparison of these quantities easier, we discuss the axis ratios that lie in the range [0,1]; in the former case this to vary, we define the c-axis as the one closest to the original rotation axis.

In many other contexts, shapes are characterized with elongation (1b/a) and flatness

(1c/b) but, in asteroid research, axis ratios have been in common use. Thus, an elongated asteroid has a low b/a-ratio, while a flat asteroid has a low c/b-ratio. Table 3 In Table 3, we show the mean, standard deviation, and maximum and minimum values for the axis ratios obtained with the four different methods. In general, the methods give similar mean values, except for the C b /C c -ratio of PMI, whose mean is lower than the other corresponding ratios. The value of b/a is, by definition, constrained to the range [0,1] for all the methods, but c/b is allowed to exceed unity for the OD and PMI techniques. In Figs.

10a and b, we see that the first three approaches (OD, E1, and E2) give similar results for individual targets, whereas the results of the PMI technique are clearly different (Fig. 10c).

Thus, the axis ratios from the first three methods can be compared with one another, while the principal moments should be considered separately. Figure 11 There is an interesting feature of S-and D-type objects achieving significantly smaller b/a ratios than C-and X-types. From [START_REF] Bus | Phase II of the small main-belt asteroid spectroscopic survey. A feature based taxonomy[END_REF], it is evident, that C-and X-type asteroids are spectrally close to one another, while the spectra of S-and D-types deviate from each other and from those of the C-and X-types. The S-types also have higher albedos than the C-and X-types. The C-types are known to be of lower density than the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

other spectral types, which could explain the more spherical shapes, since the energy of a collision, strong enough to fracture, e.g., an S-type asteroid, is absorbed into the interior by breaking the structure of the porous material inside a C-type asteroid, thus producing a global crater on the surface of the target (e.g., asteroid (253) Mathilde), but leaving it otherwise unfractured. Further studies with more observations of different spectral types are called for to shed light on this possible correlation between the spectral types and shape.

Shape, rotation period, and size

Having at hand methods to characterize the shapes of convex bodies, we can study the correlations between shape and other physical properties. In the following, we show distributions with dimensions obtained using only the E2 technique, since other characterization techniques give similar results.

There is dependence between shape and rotational period as seen in Fig. 12, since all of the most elongated objects have periods between 4.5 -10 h, while the ones with axial ratios b/a > 0.7 are slow or fast rotators. It is also interesting, that only C-and S-type asteroids have periods over 11 h. There is no evident explanation for these results, and more data are required to make any further conclusions. Especially, long-period targets (P > 15 h) lack proper lightcurve data. Also, data of short-period asteroids (P < 2 h) would be of interest since, for these, the centrifugal forces begin to override the internal forces, and loose rubble pile asteroids may begin to separate into parts. This should be seen in the elongation of these objects. Figure 12 The shape vs size plot (Fig. 13) also shows non-random features: as would be expected, the small objects reach more elongated shapes than the large ones while spherical shapes are represented by the entire size range. Observations of the largest asteroids are called for to provide information about the limit above which the asteroids are spherical. This limit should also depend on the taxonomic type. In our sample, the largest asteroid, (2)

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Pallas (D = 523 km), is spherical, but the second largest, (10) Hygiea (D = 429 km), is clearly elongated with b/a = 0.75. Both are C-type asteroids. Taxonomic and shape data of asteroids smaller than 30 km in diameter is now almost completely lacking, and thus more observations of this population are also called for. Figure 13

Principal moments

The principal moments of the shape models are not suitable for the determination of the axis ratios, since they give systematically lower values than the techniques giving the actual dimensions (Fig 10c). They can, however, be used to estimate the goodness of the convex shape model.

The axis of the largest inertia of the convex shape model does not usually coincide with the rotation axis of the model. One plausible explanation for this is that the lightcurve features caused by non-convex shape features and possible albedo markings in the target asteroids are reproduced by the convex model shapes with homogenous albedo distribution, distorting the shape somewhat, and causing deviation from principal-axis rotation. However, even slight changes in shape of nearly spherical models may change the directions of the principal axis significantly. Principal axes are also unstable for models with b/a-ratio close to unity since the vertical dimension for such objects is ill-defined. Thus, knowing the angle θ between the axis of largest inertia and axis of rotation of the model shape may provide us information about the scale of non-convexities or albedo features of target shapes with b/a-ratio significantly different from unity.

As seen in Fig. 14, θ reaches values up to 45 • . For elongated objects it seems to be smaller, increasing with increasing b/a-ratio. A large range of θ values would be expected for objects with b/a > 0.9 but, curiously, it decreases sharply with b/a increasing from 0.9 to 1. Along with increasing c/b-ratio, the deviation increases steadily towards spherical shapes.

As an example, elongated shapes with b/a-ratio below 0.7 usually have θ < 20

• . A few such A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
objects, however, show larger θ values, suggesting strongly non-convex shape, or significant variegation of albedo on the surface. Figure 14 In Fig. 11d, we see that two of the shape models have principal moment ratios C b /C c above unity, which implies an unstable rotation state. Such solutions exist, since the shapes have been inspected only by eye, and it is impossible to decide if an irregular object has C b /C c = 1 or a bit larger. Slight excesses can be tolerated for two reasons. First, the dimension along the spin axis of the shape solution is sensitive to the scattering law. For determining the shape solutions used in this paper, a combination of Lommel-Seeliger and Lambert scattering laws has been used, and the weight of the Lambertian part has been about 10%. However, this scattering law is inaccurate for expressing the scattering behaviour of asteroid surfaces, and we do not know what the optimal weight for each asteroid would exactly be. Without altering the shape otherwise, the vertical dimension of the model can be moderately altered by changing the weight of the Lambertian part.

If the vertical stretching in principal-moments sense was only due to the ill-defined scattering law, there would be basis for changing the scattering parameters to produce a more spherical shape. However, another reason for too large C b /C c -ratios is the same as that for the tilt angle of axis of largest inertia to the rotation axis of the model -the convexity of the shape model. For these reasons we should be careful when considering refinement of the vertical axis dimension of the shape solution. This dimension is anyway the least constrained quantity of the shape model, when only relative photometry is available. If the vertical stretching is small, not visible to an eye, the solution can be accepted.

Irregularity of shapes

One way to demonstrate the irregularity of a shape is to compute the deviation from the best-fit ellipsoid. Here we use the root-mean-square value (rms) of the ellipsoid fit from the E2 technique

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT rms = χ 2 e n v , (9)
where χ 2 e is obtained from Eq. 7, and n v is the number of vertices. As seen in Fig. 15, rms lies within range [0.02,0.06]. Curiously, all the most elongated objects with b/a < 0.65 have rms < 0.04. There are only a few such objects, however, and more observations are needed to see if this is a real phenomenon. Towards flatter shapes, i.e., along with decreasing c/b-ratio, the shapes seem to become more irregular in general. Otherwise, there does not seem to be a correlation between rms and elongation or flatness. Between size and rms, there is a slight negative correlation with Pearson product-moment correlation coefficient r = -0.32. This seems natural, since large objects can be expected to be more regular in shape than small ones. Between rms and taxonomic type, there is no correlation. Figure 15 The targets with small rms have lightcurves characteristic for ellipsoidal shapes, i.e., two equal minima and maxima per lightcurve period, and they might as well be modelled as ellipsoids. Most asteroids, however, have lightcurves with more than two minima and maxima or unequal maxima or minima. If such objects are modelled with an ellipsoid, a lot of information is lost.

Another way do describe the irregularity of a polyhedron is to calculate how large a portion of the model surface is covered with large planar areas. When this portion increases, the irregularity of the objects increases as well. Since the distribution of the facet areas cannot be guaranteed to be independent on the shape model, using the distribution of facet areas as a classification method requires that all the models are obtained with the same modelling technique and the same number of shape parameters. In Table 4, we show, as an example, quantities f 1 , the fraction of the facets with area more than 20% of the largest facet area, and f 2 , the fraction of the largest facet area of the total surface area (we give these quantities for only 52 shapes, since the exact modelling technique for the rest was

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
unclear. f 1 and f 2 are somewhat correlated with Pearson's r = -0.57, i.e., when the number of large facets decreases, the fraction of the total surface area covered by the largest facet increases. In Fig. 16, f 1 is plotted against the E1 ellipsoid axis ratios. We see that, Figure 16 towards more spherical shapes, the number of large facets increases, which means a more even distribution of facet sizes, and implies a more regular shape. Results are in agreement with the observation that, on elongated shapes, only a few facets often comprise a large fraction of the total surface area. According to simulations and test analysis, these large facets can be interpreted to show potential locations of global-scale concavities (KT01, K01).

Within our sample, only C-type asteroids have f 1 > 0.12, while the most dense cluster of Xtypes is located in f 1 < 0.05. This implies that, among the C-type asteroids, there are more regular bodies than among other spectral types. This is in agreement with the assumption that the C-type asteroids, consisting of more porous material than the other types, are also less subjected to fractionation than other types. In Fig. 17 17c shows that a large part of the surface area is covered by only a few facets, while the model of (2) Pallas (Fig. 17a) has more evenly distributed facet areas.

This implies that Berbericia is irregular, while Pallas has more soft features. The case of (6) Hebe lies in between these. The linear-exponential shape of the distribution is typical for models for which facets were combined for Minkowski minimization after the convex inversion procedure. The sharp decrease of size in the small-end of the facets distribution is caused by the combination procedure, which leaves some facets unchanged. In Table 4, we give the slope k of the linear part for the 34 models, and in Fig. 18 we plot k against size.

Decreasing value of k means more evenly distributed facet areas, and thus, a more regular Figure 18 shape. According to the plot, the C-type asteroids get more irregular with decreasing size, while other types show more shattered k values.

In Fig. 19 we see that, as the dominance of the largest facets increases, the deviation θ

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
of the axis of the largest moment of inertia from the rotation axis increases as well. Objects with both large θ and small f 1 are likely to contain global scale nonconvexities, locations of which might be incicated by the locations of the largest facets. Figure 19 We also calculated the shapes of the largest facets as the axis ratios of the best-fit ellipse.

There was no evident correlation with facet shape and model shape or taxonomic type. Thus, facet shapes cannot be used to define the existence of global scale circular craters, expected to exist mainly on C-type asteroids.

Table 4 5 Conclusions and future work

Three approaches are presented for calculating the global dimensions of convex polyhedron shape models of asteroids, and the results from all the techniques are in agreement. The relatively good ellipsoid fits to model shapes suggest that the surface irregularity of asteroids might be described by multiplying the best-fit ellipsoid by a spherical harmonics series, as was done for potato tubers by Torppa et al. (2006). This would be a continuation to the work of [START_REF] Muinonen | Inversion of shape statistics for small solar system bodies[END_REF], who modelled asteroids statistically with Gaussian spheres. And, after all, asteroid and potato tuber shapes are not that far from each other.

The principal moments of inertia may be used to describe the dynamical state of the model, providing information about how well it describes the shape of the target asteroid.

Notable deviation of the principal axis of the model shapes from their rotation axis, especially for elongated objects and along with the existence of large planar areas, may reflect presence of large-scale concavities or strong albedo variegation on the surface. To evaluate the reliability of this interpretation, however, we have to carry out a large number of simulations for various shapes and scattering properties.

In addition to the OD, E1, E2, and PMI techniques, we calculated the principal components, i.e., the Karhunen-Loeve transform, for each model shape by sampling evenly a grid

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
of points on their surfaces. Axis ratios, however, differed significantly from those of the other methods, and it is left for future studies to find a physical interpretation for these quantities, and possibly a more suitable way to sample the points within the shape.

The statistical study of various physical properties of asteroids shows interesting features:

S-type asteroids achieve smaller b/a ratios than C-and X-types, large asteroids, as well as those with long rotation periods have less elongated shapes, and C-type objects dominate the population of large asteroids. Original facet information from the convex inversion procedure should be preserved since, it seems, that the distribution of facet areas is in somewhat binded to other physical properties and thus, could be used to classify asteroids. Combining facets for Minkowski minimization seems to change the original facet area distribution, and produce a logarithmic distribution for the intermediate size range, while preserving the spike for the largest facets (Fig. 17), showing those shapes in which large planar areas exist. Additional shape data are required, however, to make further conclusions about classifying asteroids according to the slope of the linear part of the logarithmic distribution. Additional shape, size, and taxonomic data of the extreme populations (those with diameters D < 30 km and D > 300 km, and with periods P < 2 h and P > 15 h) are, however, called for in order to perform a more thorough analysis of the asteroid population. Also, different orbital populations should be separated in statistical studies, but also this requires a lot more observational data. Especially, radar-derived shapes would increase the number of small, kilometer-sized near-Earth asteroids available for characterization. The shape solutions derived from the radar observations or spacecraft flybys, however, cannot be included in the analysis at the moment, since none of our techniques allows processing of non-star-like shapes. Thus, we plan to develop the existing techniques and, in the future, include also non-star-like shapes in the analyses.

(

  2000, in Lagerkvist et al. 2001). The total time span of the observations was 26 years. Two solutions were accepted, namely 1) λ = 49 • , β = -40 • , and P = 5.98313 h and 2) λ = 216 • , β = 21 • , and P = 5.98088 h. The shape solutions for both pole directions were similar, and the model is shown in Fig. 6. The previous spin solution in the Poznan database is λ = 196 • , β = 13 • , and P = 5.9784 h.Figure6(130) Elektra, G, 215 × 155km

  et al. (1987),[START_REF] Shevchenko | Magnitude-phase dependencies for three asteroids[END_REF][START_REF] Durech | Physical models of ten asteroids from an observers' collaboration network[END_REF].The total time range of the observations was 26 years. The solution that fitted the data the best was λ = 160 • , β = -85 • , and P = 5.22466 h. The corresponding shape model is shown in Fig.7. The solution is the same, within error limits, as the one by[START_REF] Durech | Physical models of ten asteroids from an observers' collaboration network[END_REF] (λ = 65 • , β = -88 • , and P = 5.224664 h.

  ,[START_REF] Hainaut-Rouelle | Lightcurves of selected minor planets[END_REF]), Zeigler (1987), and[START_REF] Michalowski | Photometry and models of selected main-belt asteroids (II): 173 Ino, 376 Geometria, and 451 Patientia[END_REF]. Observations ranged a total of 21 years. The best fit solution was λ = 68 • , β = 2 • , and P = 7.71002 h. The shape model is shown in Fig.

  obtained three acceptable spin solutions: 1) λ = 29 • , β = 41 • , and P = 7.7.11364 h, 2) λ = 222 • , β = 40 • , and P = 7.11246 h, and 3) λ = 236 • , β = 56 • , and P = 7.11512 h. The corresponding shape solutions are different, and are not displayed here. The first of the pole solutions is close to one of the two earlier solutions in the Poznan database.(776) Berbericia, C, 390 km 33 lightcurves were used. Two of the curves were unpublished (see Table1), and the rest published in[START_REF] Debehogne | 804, 814[END_REF], Di[START_REF] Martino | Photometric photometry of 17 asteroids[END_REF],[START_REF] Schober | 87 Aquitania and 776 Berbericia -Two slow spinning asteroids with rotation periods of nearly one day[END_REF][START_REF] Hainaut-Rouelle | Lightcurves of selected minor planets[END_REF][START_REF] Durech | Physical models of ten asteroids from an observers' collaboration network[END_REF]. The observations spanned a total of 28 years. Two solutions that fitted the data equally well were 1) λ = 170 • , β = +59 • , and P = 7.66701 h and 2) λ = 347 • , β = +11 • , and P = 7.66701 h. The shape solutions were but mirror images through horizontal reflection. The shape corresponding to the solution λ = 170 • , β = +59 • is shown in Fig.9.[START_REF] Durech | Physical models of ten asteroids from an observers' collaboration network[END_REF] obtained a single solution λ = 347 • , β = +12 • , and P = 7.66701 h, in agreement with one of our solutions.Figure9

  corresponds to c/a, c/b and b/a, while in the latter case to C a /C c , C b /C c , and C a /C b . For the OD and E1 techniques, the definition of the c-axis is straightforward, since it is fixed to be aligned with the rotation axis. For the other two methods, where all the axes are allowed

  Figure 10 We use the first three methods to define the actual dimensions of the asteroids. In Fig. 11a we show the distributions of dimensional ratios (c/b, b/a) from the E2 technique (the OD and E1 techniques produce similar distributions). We see, that c/b and b/a most commonly lie within the range [0.7,1], b/a getting slightly lower values. There are a few objects that are significantly elongated with b/a below 0.6. These can be binary asteroids, or Eros and Itokawa -like elongated objects. The c/a-ratio usually lies within the range [0.5,1]. The results are in accordance with laboratory and numerical simulations (see Holsapple et al. 2002), which suggest that, for fractured shapes, b/a = [0.6, 0.7] and c/a = [0.45, 0.5].

  we show how the facet sizes are Figure distributed in decreasing order on the surfaces of three different model shapes. The model of (776) Berbericia in Fig.

Figure 1 :

 1 Figure captionsFigure 1: Finding the overall dimensions of the shape. a) extent (2c) in spin axis (z) direction. b) largest extent (2a) in xy-plane and largest extent (2b) in xy-plane, perpendicular to a. The images are from equatorial and polar viewing directions, respectively.

Figure 2 :

 2 Figure 2: Previously unpublished lightcurve data of a) (31) Euphrosyne, b) (54) Alexandra, c) (79) Eurynome, d) (93) Minerva, e) (130) Elektra, f) (376) Geometria, g) (471) Papagena, and h) (776) Berbericia. The horizontal axes refer to the time (UT/hours) of observation; the observing date is given on each plot.

Figure 3 :

 3 Figure 3: Period uncertainty of the eight convex inversion solutions published in this paper, versus the number of observed lightcurves.

Figure 4 :

 4 Figure 4: Distribution of possible spin states of (31) Euphrosyne. a) rms of the model fit versus the rotation period (P ) and b) heliocentric ecliptic pole longitude (λ) versus pole latitude (β).

Figure 5 :

 5 Figure 5: Shape model of (79) Eurynome from two perpendicular equatorial viewing directions.

Figure 6 :

 6 Figure 6: Shape model of (93) Minerva from two perpendicular equatorial viewing directions.

Figure 7 :

 7 Figure 7: Shape model of (130) Elektra from two perpendicular equatorial viewing directions.

Figure 8 :

 8 Figure 8: Shape model of (376) Geometria from two perpendicular equatorial viewing directions.

Figure 9 :

 9 Figure 9: Shape model of (776) Berbericia from two perpendicular equatorial viewing directions.

Figure 10 :

 10 Figure 10: Comparison of different methods for characterizing the shapes of the objects. a) OD and E1, b) E1 and E2, and c) E2 and PMI.

Figure 11 :

 11 Figure 11: Distributions of axis ratios of 87 convex inversion shape solutions from the a) E2 and b) PMI techniques. Different taxonomic types are plotted with different symbols as =C, =D,+ =S, =X, and × =undefined.

Figure 12 :

 12 Figure 12: Comparing b/a to the rotation period of the asteroid. Dimensions are from the E2 technique. Different taxonomic types are plotted with different symbols as =C, =D,+ =S, =X, and × =undefined.

Figure 13 :

 13 Figure 13: Comparing b/a to the size of the asteroid. Dimensions are from the E2 technique. Different taxonomic types are plotted with different symbols as =C, =D,+ =S, =X, and × =undefined.

Figure 14 :

 14 Figure 14: Angle θ between the rotation axis of the model and the axis of the largest moment of inertia. θ is plotted versus b/a and c/b (from the E2 method). Different taxonomic types are plotted with different symbols as =C, =D,+ =S, =X, and × =undefined.

Figure 15 :

 15 Figure 15: rms of the E2 fit as a function of a) b/a and c/b, and b) size. Different taxonomic types are plotted with different symbols as =C, =D,+ =S, =X, and × =undefined.

Figure 16 :

 16 Figure16: Fraction of the facets with area more than 20% of the largest facet area (f 1 ) compared to the dimensions of the asteroid. Different taxonomic types are plotted with different symbols as =C, + =S, =X, and × =undefined.

Figure 17 :

 17 Figure 17: Distribution of facet areas on model of (2) Pallas, (6) Hebe, and (776) Berbericia. The logarithms of the fractions of the largest facets from the total surface area are plotted in decreasing order.

Figure 18 :

 18 Figure 18: Slope k of the log(A/A tot ) vs. size distribution. Different taxonomic types are plotted with different symbols as =C, + =S, and =X.

Figure 19 :

 19 Figure19: Fraction of the largest facets (f 1 ) versus the deviation of the rotation axis from the axis of the largest moment of inertia.

  Torppa, J., Valkonen, J. P. T.,[START_REF] Muinonen | Inversion of small-particle silhouettes for Gaussian-sphere parameters[END_REF]. Three-dimensional stochastic shape modelling for potato tubers. Potato Res., 49, 109.Weidenschilling S.J., Chapman C.R. Davis D.R., Greenberg R.,Levy D.H., and Vail S., 1987. Photometric geodesy of main-belt asteroids. I -Lightcurves of 26 large, rapid

	rotators, Icarus 70, 191.
	Zeigler K.W., 1987. Photoelectric photometry of asteroids 213 Lilaea and 376 Geometria.
	MPBu, 14, 11.

Table 2 .

 2 Summary of the spin solutions. Columns are as follows: 1) asteroid name, 2) pole direction as ecliptic longitude and latitude, 3) rotation period, 4) rms-residual of the fit, 5) number of lightcurves, 6) number of apparitions, and 7) total time span of the observations.

	Asteroid	λ, β( • )	P (h)	rms	N cur N app T(yr)
	(31) Euphrosyne distribution 5.534 ± 0.004 0.012-10	5	30
				0.0125
	(54) Alexandra	307,20	7.02263	0.016 20	5	42
		122,-36	7.02334	0.017
		325,-37	7.02334	0.017
	(79) Eurynome	64,15	5.97857	0.014 33	5	32
	(93) Minerva	49,-40	5.98313	0.015 14	5	26
		216,21	5.98088	0.014
	(130) Elektra	160,-85	5.22466	0.019 58	12	26
	(376) Geometria 68,2	7.71002	0.021 35	10	21
	(471) Papagena 29,41	7.11364	0.021 13	3	31
		222,40	7.11246	0.021
		235,56	7.11512	0.021
	(776) Berbericia 170,59	7.66701	0.014 33	7	28
		347,11	7.66701	0.013

Table 3 .

 3 Comparison of the results of shape characterization using the four different methods. Colums contain: 1) method, 2) axial ratio, 3) mean value of the axial ratio, 4) standard deviation of the axial ratios, 5-6) maximum and minimium values of the axial ratios.

	method	ratio	σ	mean max	min
	Overall dimensions	b/a	0.1119 0.8463 0.9989 0.4903
		c/b	0.1258 0.8709 1.2102 0.5203
	Ellipsoid fit (facets) b/a	0.1061 0.8478 1.0000 0.5556
		c/b	0.1131 0.8790 1.0000 0.4762
	Ellipsoid fit (radius) b/a	0.1232 0.8348 0.9945 0.4343
		c/b	0.0942 0.8670 1.0000 0.6550
	Principal moments	C a /C b 0.1191 0.8588 0.9977 0.4978
		C b /C c 0.1320 0.8124 1.3494 0.4788

Acknowledgments

The study has been funded by Vilho, Yrjö and Kalle Väisälä foundation, the Academy of Finland, the Finnish Graduate School in Astronomy and Space Physics, and the University of Helsinki. JT thanks all the Finnish amateur astronomers who took part in carrying out the observations.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT tables

Table 1. Information about the photometric lightcurve observations published in this study. Columns contain the following information: 1) asteroid name, 2) date of observation, 3) phase angle, and 4) observatory code (TH-Taurus Hill, SR-Slope Rock, JK-Jakokoski).

Asteroid Obs date α(


Keywords: Mars, Surface, Water

In this paper we attempt to answer the question, how formation of a small-scale trench in the Martian regolith affects local distribution of the subsurface ice. We are especially interested in the consequences of digging a trench to search for buried ice, as has been done during the Phoenix Mars Lander mission. However, the results may be also applicable for natural troughs, or cracks. We present results of simulations of diurnal exchange of water between the regolith and the atmosphere. Our model includes the heat and vapor migration in the regolith surrounding the trench, as well as formation of diurnal frost. We take into account scattering of light in the atmosphere and on the trench facets, as well as changes of atmospheric humidity on diurnal and seasonal time scales. Our calculations show, that the measurements of ice content in a sample obtained within one, or two days from the beginning of digging should not be affected. However, on somewhat longer time scale at the south facing site of the trench the regolith can be significantly depleted from ice. This effect should be taken into account if the excavation and taking samples from different depths will be performed in stages separated in time by a month, or more.

A C C E P T E D M A N U S C R I P T 1 Introduction

Our work is intended to simulate evolution of the surface frost in a small trench, like these dug by the Phoenix Mars Lander. We calculated also evolution of the ground ice beneath the trench. We performed our simulations before landing and the successful digging in Martian regolith. Thus, we consider properties of the regolith deduced by several authors using remote observations. Below, we present short review of the estimates regarding composition of the regolith.

The opinions about the amount of water present on Mars significantly evolved during past decades. Theoretical simulations performed by many authors [START_REF] Leighton | Behavior of carbon dioxide and other volatiles on Mars[END_REF][START_REF] Fanale | Global distribution and migration of subsurface ice on Mars[END_REF][START_REF] Mellon | The distribution and behavior of Martian ground ice during past and present epochs[END_REF][START_REF] Schorghofer | Stability and exchange of surface ice on Mars[END_REF][START_REF] Aharonson | Subsurface ice on Mars with rough topography[END_REF] led finally to the conclusion, that the subsurface ice may be present in wide range of latitudes. The neutron and gamma ray spectrometers on board Mars Odyssey spacecraft confirmed, that some subsurface ice is present almost everywhere on Mars except low latitudes [START_REF] Boynton | Distribution of hydrogen in the near surface of Mars: Evidence for subsurface ice deposits[END_REF][START_REF] Mitrofanov | Soil water content on Mars as estimated from neutron measurements by the HEND instrument on-board the 2001 Mars Odyssey Spacecraft[END_REF]. This is in agreement with other observations [START_REF] Kossacki | Surface temperature of Martian regolith with polygonal features: influence of the subsurface water ice[END_REF][START_REF] Titus | Exposed water ice discovered near the south pole of Mars[END_REF]. Direct search for the subsurface ice was planned for the failed Mars Polar Lander mission. Now, similar experiment involving excavation of a trench is performed by the Phoenix lander. The regions initially considered for the landing of Phoenix were at the latitude about 70 o N. The final landing site is at 68.16 o N, 233.35 o E.

At the latitude as high as 70 o regolith was commonly expected to be rich in ice. However, high resolution observations of the surface temperature indicated that the subsurface concentration of ice on depths up to ten centimeters can be nonuniform and locally very low. This made various theoretical studies important for the preparation of the digging experiment. [START_REF] Sizemore | Effects of soil heterogeneity on Martian ground-ice stability and orbital estimates of ice table depth[END_REF] simulated possible non-uniformity of the ground-ice distribution in heterogeneous Martian subsurface. They found, that presence of buried decimeter scale rocks should create some depressions in the ice table over a horizontal range of 1 -2 rock radii. They also found that, dust lenses of the same scale may locally reduce depth to the ground ice by few centimeters. Small-scale heterogeneity of the regolith, as described above, and the related non-uniformity of the ground ice distribution should affect thermal inertia. This problem was recently investigated by Putzig and Mellon (2007b). The authors calculated thermal inertia based on day time and nigh-ttime observations performed using the Mars Global Surveyor Thermal Emission Spectrometer (three Martian years). At low latitudes the derived apparent thermal inertia is only moderately dependent on local time, but in the polar regions diurnal differences locally exceed 600 J m -2 K -1 s -1/2 . Based on the comparison between the results of modeling and the seasonal maps of the thermal inertia Putzig and Mellon (2007b) concluded that in the polar regions ground ice is typically covered by ice-free sand. Similar result was obtained by [START_REF] Bandfield | High-resolution subsurface water-ice distributions on Mars[END_REF], who compared observations performed using the Thermal Emission Imaging System (THEMIS) on the Mars Odyssey spacecraft with results of thermal simulations.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
Another problem is local inclination of the surface. [START_REF] Fergason | High-resolution thermal inertia derived from the Thermal Emission Imaging System (THEMIS): Thermal model and applications[END_REF] compared the night time temperatures of slopes at different orientations. For low and moderate inclinations up to 30 o the temperatures of the opposite facing slopes appeared to be the same to within few degrees. However, [START_REF] Fergason | High-resolution thermal inertia derived from the Thermal Emission Imaging System (THEMIS): Thermal model and applications[END_REF] did not analyze areas at high latitudes, where the zenith angle of the Sun is much larger. Thus, their result is not applicable for Phoenix mission. More recently Putzig and Mellon (2007a) compared apparent thermal inertia calculated for surfaces either horizontal, or partially sloped. The authors state, that their result are applicable to mixtures where the material properties vary at scales larger than a thermal skin depth. Analysis of high resolution HiRISE stereo images (down to 2 m length scale) indicates that in the area selected for the Phoenix landing site rocks abundances are low and slopes are typically < 3 o even at fine scales [START_REF] Arvidson | Overview of Mars Exploration Program 2007 Phoenix Mission Landing Site Selection[END_REF]. Thus, in our work dealing with the trough up to meter scale we assume, that the surface is horizontal and initially flat.

Digging a trench disturbs temperature field in the regolith. It should also result in direct contact of the uncovered layers with the atmosphere. Both effects may lead to significant changes of ice content in the regolith surrounding the trench and hence of ice content of the sample analysed by the Phoenix instruments. In this work we attempt to answer two questions: (i) how fast are changes of the subsurface ice content in the regolith induced by digging, and (ii) how much frost can form in the trench due to condensation of the atmospheric vapor.

Model

We use the LWMRM (Lindau-Warsaw-Mars-Regolith-Model) model which was initially developed to simulate thermal conditions for the formation of frost in the trench to be excavated during Mars Polar Lander mission [START_REF] Kossacki | Effect of Surface Roughness on Ice Distribution in the South Subpolar Region of Mars[END_REF]. In subsequent applications of LWMRM we improved its accuracy, allowed for more complex geometries and implemented more physical processes. The most recent versions were used to investigate the possibility of seasonal ice melting in the Martian gullies [START_REF] Kossacki | Seasonal melting of surface water ice condensing in Martian gullies[END_REF], and to simulate the escape rate of subsurface ice in the southern part of the Elysium region [START_REF] Kossacki | Possible remnants of a frozen mud lake in southern Elysium, Mars[END_REF]. In the current version we modified the way of calculating the condensation of atmospheric vapor in the regolith.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Below we briefly describe the main features of the LWMRM model. Farther details can be found in the references listed above.

The LWMRM model describes diurnal and seasonal condensation and sublimation cycle of H 2 O and CO 2 ices within a Martian trench, or crater, and in the surrounding regolith. In the current work we are interested in a trench. The regolith is composed of equal size dust grains mantled by ice if this is present. The thickness of ice covering a dust grain is uniform and is determined by the local volume of ice in the regolith. In real situation ice can migrate along the surface of a grain due to local differences of curvature of the surface. However, including this effect significantly reduces speed of calculations. Thus, we consider only instantaneous thickness of ice mantling the grains. The model includes the time dependent solar illumination, infrared emission from the atmosphere toward surface, infrared emission from the surface, heat and vapor transport within the regolith and the heat exchange due to sublimation and condensation. The flux of absorbed light depends on the local orientation of the facets (wall or bottom) and on the current position of the Sun. The latter evolves due to orbital motion and rotation of Mars. The model includes absorption of direct solar light, as well as light scattered once within the considered trench. Infrared emission from the facets is considered as well, including their relative orientations. For the absorption of the infrared radiation from the atmosphere we take into account fraction of the sky visible in a given location in the trench. The total IR radiation emitted by the atmosphere toward surface can be either parametrized, or taken from results of a general circulation model such as GCM LMD/Oxford (Forget et al., 1999;[START_REF] Lewis | A climate database for Mars[END_REF], available at http://www-mars.lmd.jussieu.fr/.

In the current work we use the latter option. The seasonal and diurnal variations of atmospheric pressure and humidity are also taken from the results of LMD/Oxford GCM. Our model allows including both seasonal and diurnal cycles of atmospheric IR emission. However, during polar night when the atmospheric IR radiation has highest relative influence the GCM results have low accuracy. Thus, we include only seasonal changes of the atmospheric IR emission. Formation of ice due to condensation of water vapor from the atmosphere was in the previous versions of our model restricted to the surface. Now, as long as the surface ice is thinner than the size of grains composing regolith, ice condenses simultaneously on the surface (including trench facets) and in the near surface layer of the regolith.

Our model, as all models describing water exchange between regolith and atmosphere, has several parameters. For present purpose most important is sensitivity of the model to the distribution of the subsurface ice. In the section 2.4 we show sensitivity of the calculated surface temperature and the frost thickness to the local content of the subsurface ice and to the atmospheric humidity.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

Equations

The main equations are those describing: local illumination, the rates of condensation and sublimation of ices on the surface, the diffusion of heat and vapor in the regolith, and the instantaneous local thickness of ice on the regolith grains.

The flux F s of the ice subliming from the facets of the trough, or from the external surface is determined by the pressure p sat at the local surface temperature. For the flux of condensing water vapor we use the modified equation proposed by [START_REF] Ivanov | The Role of Sublimation for the Formation of the Northern Ice Cap: Results from the Mars Orbiter Laser Altimeter[END_REF]. Here,

F s,H 2 O = u β μ H 2 O R T surf (P s,H 2 O -P H 2 O ) , ( 1
)
where u and P H2O are the horizontal wind speed and the partial pressure of vapor at the height H above the surface. If it is not stated otherwise, H = 1 m. P s,H2O is the pressure of phase equilibrium between solid and gas H 2 O at the current temperature of ice surface T surf . For the drag factor β we applied the value 0.002 suggested by [START_REF] Paterson | The Physics of Glaciers[END_REF] for Earth glaciers. The speed of wind can not be calculated without high resolution 3D simulation of the atmosphere, including topographic effects and the presence of the lander. Thus, we take u as a parameter.

The diffusion of heat and vapor in the regolith, as well as the surface temperature strongly depend on the structure and the thermal conductivity of the regolith. Thermal conductivity itself depends on the regolith structure and the temperature. The regolith can be dry, or filled with some amount of ice. Thus, its thermal conductivity and the gas permeability can vary in wide ranges.

In our work the regolith is a matrix of sand grains, with the void space partially filled by ice mantling mineral grains. The model accounts for sizes of the contact areas between the grains. These areas are determined by the local content of ice and the porosity. There are various formulas in the literature for the thermal conductivity of granular media. Here,

λ = λ m + h(v i )λ i v i , ( 2
)
where λ m is the conductivity of the ice-free matrix (dust), λ i is the bulk thermal conductivity of ice, v i is the volume fraction of ice, and h is the Hertz factor. The latter is defined as the ratio between the grain-to-grain contact area and the effective cross-section of a grain including the ice coating of the grain if it exists. The functional form of h(v i ) depends on the microscopic model of the regolith. We assume that the ice mantle covering a mineral grain has

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
uniform thickness, as it was described in the section 2. This allows calculation of the instantaneous value of h, determined by the time dependent v i .

Complete set of equations included in the previous versions of the model is described in [START_REF] Kossacki | Seasonal melting of surface water ice condensing in Martian gullies[END_REF], and [START_REF] Kossacki | Possible remnants of a frozen mud lake in southern Elysium, Mars[END_REF].

Initial and boundary conditions

In the current work we consider trench long when compared to its cross-section.

The heat and vapor diffusion equations are solved in the plane perpendicular to the trench axis, hence in two dimensions. However, illumination of the trench is calculated in three dimensions, this makes the simulations 2.5 D. Boundary conditions on the facets of the trench and on the regolith surface external to the trench are variable. The temperature satisfies heat balance through the energy equation. The vapor flux is determined by the temperature, the atmospheric humidity and the wind speed. The latter is a free parameter. The numerical grid is large enough that the fluxes through the walls and the bottom of the grid are negligible. Thus, we assume zero values for the vertical heat and vapor fluxes at the bottom of the grid and zero horizontal fluxes at the sides of the grid.

Before simulating excavation we run the program for ten full Martian years. Thus, at the time of digging the temperature versus depth is not affected by the initial conditions. Therefore, we can assume uniform initial temperature of the regolith. The initial ice content v i 0 increases linearly with depth from zero at the surface to v i,max at some depth d. At larger depths v i 0 = v i,max .

Parameters

We consider trench excavated in an area similar to the landing site of Phoenix. We performed simulations when the exact latitude of the landing site was still unknown. We decided for the latitude 70 o N. The finally selected landing site is at the latitude 68.16 o .

We look at simple, two stage scenario of excavating the trench. In the first stage we dig straight a way to the depth D = 20 cm, and width at the bottom L = 40 cm, with the side walls inclined by 45 o (like the trench dug in Beacon Valley Antarctica [START_REF] Mellon | The presence and stability of ground ice in the southern hemisphere of Mars[END_REF]). In the second stage of digging, 3 sols later, we dig more to reach the depth of 30 cm. The final trough is 20 cm wide at the bottom. The simulated trough is oriented either N-S, or E-W. The former orientation is approximately the same as planned for the trench to be excavated, while the latter is added to investigate conditions at the ends of

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
the N-S trench. It should be noted, that there is some uncertainty of landing orientation, and the robotic arm has ability to pivot around the vertical axis. Thus, the exact orientation of the trench is a priori not known. The ends of the real trench will have complex shape, but will face approximately toward North and South, respectively. The simulated excavation starts at L s = 77.4 o . The simulations are performed with spatial resolution of 5 cm. The time step is 5 minutes.

In the regions considered as possible landing sites the regolith is two layered. The upper layer is 3 -6 cm thick and has the thermal inertia slightly above 200 J m -2 K -1 s -1/2 , while the lower layer has 20 -50 % of ice, and has thermal inertia 400 -1100 J m -2 K -1 s -1/2 [START_REF] Titus | Thermal Inertia Characterization of the Proposed Phoenix Landing Sites[END_REF]. Locally, around rocks thickness of the upper dry layer of the regolith may be significantly larger than 10 cm [START_REF] Sizemore | Effects of soil heterogeneity on Martian ground-ice stability and orbital estimates of ice table depth[END_REF]. In our model the regolith is composed of single size spherical grains whose radii r g are 5 μm, or 0.5 mm. The regolith density is 1300 kg m -3 , and the specific heat c = 820 J kg -1 K -1 . The porosity of the regolith, when it is free of ice, is 0.47. The initial volume fraction of ice v i 0 increases with depth from zero at the surface to v i,max at the depth d. At any depth below d, ice volume v i 0 = v i,max . In our canonical case v i,max = 0.38. This value corresponds to porosity of about 0.1, while the percolation threshold is typically 0.2 -0.3. Thus, at v i = 0.38 most pores are closed and we can consider the regolith as effectively saturated with ice. We consider two values of d: 0.2 m and 0.5 m. If it is not stated otherwise d = 0.5 m. Thermal conductivity k(v i ) of the regolith is calculated accounting for the contact areas of the grains, which in turn is related to the porosity. Condensation and sublimation of ice within this matrix modifies the conductivity according to how much ice is present at any particular time and place. The thermal conductivity of the regolith free of ice k(v i = 0) is 0.01 W m -1 K -1 . Thus, in our model thermal inertia linearly increases versus depth from 100 W m -1 K -1 at the surface to about 1000 W m -1 K -1 at the depth d. In addition, we performed simulations for the regolith of low ice content (v i,max = 0.09).

The albedo of ice covered surface strongly depends on the purity of ice, and the ice granulation. In our model albedo has one of two values: A for the unfrosted surface, or the surface covered by ice thinner than 0.1 mm, and A i for the surface covered by thicker ice. We consider A = 0.25 and A i = 0.35.

The seasonal and diurnal variations of the atmospheric pressure and humidity are taken from the results of LMD/Oxford GCM. In both cases we have 4 data points per Martian day. Thus, we perform interpolation to get values of atmospheric pressure and humidity for each step of integration. The atmospheric pressure in the considered region is about 1000 Pa. According to [START_REF] Tamppari | Atmospheric Characteristics Expected at the Phoenix Landing Season and Location[END_REF] it should be within the range 700 -1100 Pa. GCM LMD/Oxford used in our work predicts similar range 730 -940 Pa. According to this model

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

the atmospheric humidity expressed by the water vapor mixing ratio is ranging from 3 10 -9 to 2 10 -3 . In Fig. 1 we show water vapor mixing ratio versus time at the altitude 5 m above the local surface. The corresponding condensation temperature is plotted as well. The profile of the mixing ratio is drawn for the whole Martian year. Thus, for the diurnal oscillations only highest and lowest values can be distinguished. The plotted profile of the condensation temperature shows, that at the time of digging the trench water frost should condense in locations colder than 195 -215 K. Diurnal minimum of the condensation temperature is predicted for morning hours. During the day time, when the surface temperature rises, the condensation temperature rises as well.

The IR radiation emitted by the atmosphere toward surface is also taken from the results of GCM. In this case we use GCM data averaged over diurnal and seasonal time scales to get 16 points per Martian year. The wind speed u at elevation 5 meters above the surface is either 5 m s -1 , or 20 m s -1 . The latter is the highest wind speed considered as acceptable for the Phoenix landing site [START_REF] Arvidson | Overview of Mars Exploration Program 2007 Phoenix Mission Landing Site Selection[END_REF].

Sensitivity of the model

In this subsection we present briefly sensitivity of our model to most important parameters. The current paper is intended to present results of simulations for the particular location on Mars, where they can be compared with the in-situ observations. Thus, it is not necessary to discuss influence of latitude and elevation. They are determined by the location of the Phoenix landing site.

The most important is sensitivity of the results to the volume fraction of ice beneath the trench.

In Fig 2 we demonstrate sensitivity of the model to the local content of the subsurface ice, and to the atmospheric humidity. Plotted are; the temperature versus time, and the frost thickness versus time. The profiles are drawn for the middle of the bottom of the 30 cm deep trench, starting from the second stage of digging (increase of the depth from 20 cm to 30 cm). The volume fraction of ice v i just beneath the trench bottom is either 0.04, or 0.17. The atmospheric humidity is taken from the results of GCM LMD/Oxford (thin black lines), or modified (thick grey lines). The grey curves correspond to the situation when the volume fraction of ice beneath the trench is 0.17 and the humidity is: 2x higher, and 2x lower than predicted by the GCM.

The calculated profiles of the temperature are significantly different. Amplitude of the diurnal oscillations is: 60 K and 30 K respectively. The highest diurnal temperatures are, 200 K and 220 K. The condensation temperature oscillates in diurnal cycle from 200 K to 215 K. Close to noon, when the

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
bottom of the trench is warmest, condensation temperature is about 210 K. Hence, periodical disappearing of the surface frost may be determined by the concentration of the subsurface ice beneath the trench. Comparison of the frost thickness calculated using different profiles of the atmospheric humidity versus time shows high importance of the local near-surface humidity. The average condensation rate of the surface frost changes proportionally to the average humidity. When the ground ice content beneath the trench is high and the humidity is reduced, the frost thickness increases with the average rate as predicted for low subsurface ice content and the GCM humidity.

Results

In Fig. 3 we show diurnal cycle of the surface temperature within the N-S oriented trench, after the first phase of digging (the trench is 20 cm deep). In addition we show temperature profile for the surface external to the trench, 20 cm away of the trench edge. Wind speed is taken to be 5 m/sec. Ice volume in the regolith is high (v i = 0.38, upper panel), and low (v i = 0.09, lower panel). Wall 1 is the east facing west wall and wall 2 is the west facing wall. According to our simulations, just after excavation ice volume v i beneath the bottom of the trench is: 0.11 and 0.03, respectively. When the ice content is high, the temperature in the trench never falls below 165 K, even in the most shadowed locations. The situation changes when the regolith beneath the trench is almost free of ice and has low thermal inertia. In this case the lowest temperatures are slightly above 150 K. During the first evening after digging the trench, temperature of the east facing wall allows even formation of the CO 2 frost.

In Fig. 4 we show profiles of the surface temperature and of the frost thickness in the E-W oriented trench. The north facing wall is coldest, as it should be expected. It is illuminated only at night, when the zenith angle of the Sun is large. The angle between the direction to the Sun and the direction normal to the wall is also large. At midnight it is 40 o , but in the morning and in the evening it exceeds 60 o . This may indicate, that the north facing wall should be warmest at midnight. However, when the Sun is very close to the horizon light is strongly attenuated. Thus, in the evening and in the morning the north facing wall is about 10 K warmer than at midnight. Nevertheless, it is always cold. In the selected point 10 cm beneath the upper edge, amplitude of the diurnal oscillations is 25 K and the minimum of the temperature is about 165 K. For the south facing wall minimum of the temperature is about 75 K.

The temperatures indicate, that the north facing wall should start collecting frost immediately after digging the trench. This conclusion is confirmed by the profiles of the frost thickness discussed further.

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT

In Fig. 5 we show the temperature cycles in three places on the bottom of the trench after the second stage of digging to 30 cm depth. Ice volume in the regolith is high (v i = 0.38). The profiles are for the middle of the bottom of the trough and for two locations 5 cm away of the middle. All three curves almost overlap. The highest temperatures are about 200 K. This value is 10 degrees lower than predicted for the shallower trench. At night time the temperature falls down to about 170 K. The differences between temperatures in the trenches of different depths are due to different illumination, as well as due to different ice content beneath the trenches. The latter affects the thermal inertia and hence amplitude of the temperature oscillations.

In Fig. 6 we present profiles of the frost thickness versus time, after the second phase of digging. The trench is oriented N-S, or E-W. The regolith has high ice content (v i,max = 0.38), or low ice content (v i,max = 0.09). The profiles are for two points on the walls, 10 cm beneath the upper edges, and in the middle of the bottom of the trench. When the trench is oriented N-S, in the selected locations on the walls frost disappears when the temperature locally reaches diurnal maximum. However, at the bottom of the trench frost accumulates from day to day. Both results are in agreement with the profile of the condensation temperature shown in Fig. 1. This profile shows, that at the time of digging the trench water frost should appear in locations colder than 195 -210 K. Hence, in the middle of the bottom of the trench the temperature is the whole time below the frost point. Content of ice beneath the trench has moderate influence. In the considered cases v i beneath the trench is: 0.17 when v i,max = 0.38, but only 0.04 when v i,max = 0.09. In the former case frost cover on the bottom of the trench increases in thickness approximately 2 times faster than in the latter one. This is because higher volume of ground ice results in higher thermal inertia and hence in lower day time temperatures. When the trough is oriented E-W, profiles of the frost thickness for the bottom and the north facing wall of the trough almost overlap. It should be noted, that the thickness of the frost is in all cases very small reaching a value of 11 micron after 10 sols. Increasing wind speed to 20 m/sec results in the maximum frost thickness reaching about 45 microns, still a small value. This should not affect measurement of the ice content in the regolith samples taken from the trench short time after digging. Fig. 7 is analogous to the upper panel of Fig. 6. Different is the depth d, now it is only 0.2 m instead of 0.5 m. It can be seen, that after 8 sols frost starts accumulate on the walls, where previously it was present only periodically. However, accumulation is slow and may cease due to some decrease of the atmospheric humidity. and on the bottom (E-W trench). For each of the considered walls the curves are for a point 10 cm beneath the upper edge (point A) and a point 7 cm beneath the wall, in the direction perpendicular to the wall (point B). For the bottom of the trough we show content of ice in a point in the middle of the trench width. For this point are drawn two profiles: for the regolith rich in ice, and for the regolith of low ice content. In Fig. 9 we show changes of the ice content at different depths in the regolith beneath the south facing wall. The profiles are drawn versus distance perpendicular to the wall from the point A at the wall. The curves are for L s = 77.4 o (just after digging the trough) and for L s = 81 o . It can be seen, that the dessication of the regolith is visible within a few centimeters beneath the wall already after this short time.

Conclusions

In this work we investigated influence of digging the decimeter scale trench, as now performed by Phoenix lander, on two processes: (i) formation of the surface frost within the trench, and (ii) changes of the ice content in the regolith surrounding the trench. Fast formation of frost in the shadowed locations, as well as quick changes of the subsurface ice distribution may have significant effect on the planned measurements. Small-and medium scale troughs on Mars were described in several papers dealing with polygonal features (e.g. [START_REF] Mellon | Small-scale polygonal features on Mars: Seasonal thermal contraction cracks in permafrost[END_REF]; [START_REF] Kuzmin | Mars: Detaching of the Free Water Signature (FWS) Presence Regions on the Base of HEND/ODYSSEY Data and Their Correlation with Some Permafrost features from MOC Data[END_REF]; [START_REF] Kossacki | Surface temperature of Martian regolith with polygonal features: influence of the subsurface water ice[END_REF]; [START_REF] Mangold | High latitude patterned grounds on Mars: Classification, distribution and climatic control[END_REF]; Page (2007); [START_REF] Haltigin | Comparative Morphometric Analysis of Polygonal Terrain at Potential Mars Phoenix Landing Sites[END_REF]; [START_REF] Mellon | Hirise Team[END_REF]). However, authors focused attention on the troughs large enough to be identifiable in satellite images. Most commonly investigated topics seem: formation of the troughs [START_REF] Mellon | Small-scale polygonal features on Mars: Seasonal thermal contraction cracks in permafrost[END_REF], correlation of the polygonal troughs with other features indicating presence of the subsurface ice [START_REF] Kuzmin | Mars: Detaching of the Free Water Signature (FWS) Presence Regions on the Base of HEND/ODYSSEY Data and Their Correlation with Some Permafrost features from MOC Data[END_REF][START_REF] Kossacki | Surface temperature of Martian regolith with polygonal features: influence of the subsurface water ice[END_REF], and determining past climate on Mars [START_REF] Mangold | High latitude patterned grounds on Mars: Classification, distribution and climatic control[END_REF]. To our knowledge short time changes of ice content in and around small (decimeter scale) troughs is

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT
for the first time simulated in details in the present work.

Our simulations show, that the surface frost should accumulate on the bottom of the considered trough independently of the conditions (depth of the trough, wind speed and the ice content beneath the trench. When the trough is oriented E-W, or when the depleted in ice upper layer of the regolith is thin, frost should accumulate also locally on the walls of the trough. However, the amount of frost is not large. When the wind speed is 5 m s -1 frost thickness should not increase faster than about 1 micron per sol. Increasing the wind speed to 20 m s -1 may result in the formation of about 0.1 mm of frost during 20 sols. This should be already taken into account. South facing wall can be significantly depleted of water within a few sols. Thus, taking into account low amount of frost accumulated during first sols, the best locations to take samples of the regolith short time after digging the trench is its bottom.

If the regolith surrounding the trench is almost free of ice, surface frost should be taken into account. The depth of the trench is 30 cm. The volume fraction of ice just beneath the trench bottom is either 0.04 (dashed line), or 0.17 (solid line). The atmospheric humidity is taken from the results of GCM LMD/Oxford (thin black lines), or modified (thick grey lines). The grey curves correspond to the situation when the volume fractions of ice beneath the trench is 0.17 and the humidity is: 2x higher (upper grey line), and 2x lower than predicted by the GCM (lower grey line).

Fig. 8

 8 Fig.8presents changes of the ice content in the regolith at and beneath facets in both N-S and E-W oriented trenches. The profiles are drawn for selected points at: the east wall (N-S trench), the north and south walls (E-W trench),

Fig. 1 .Fig. 2 .

 12 Fig. 1. Water vapor mixing ratio versus time at the altitude 5 m (upper panel), and the corresponding condensation temperature (lower panel).

Fig. 4 .

 4 Fig.4. As Fig,3but for trench oriented E-W. Wall 1 is the north facing south wall and wall 2 is the north wall. The profiles are drawn for the selected locations: on the walls, 10 cm beneath the upper edge; and the middle of the trench bottom.

Fig. 5 .

 5 Fig. 5. As Fig. 3 but after the second stage of digging to 30 cm depth and for three locations on the trench bottom.

Fig. 7 .

 7 Fig. 7. As the upper panel of Fig.6, but for small depth d to the lower (filled with ice) layer of the regolith: 0.2 m instead of 0.5 m.

  The wind speed has largest of the considered values, 20 m s -1 . This makes the calculated desiccation rate of the regolith an upper estimate. Decreasing the wind speed proportionally decreases the desiccation rates. The plotted volume fractions of ice are normalized to the local values of ice volume v i,ref at L s = 75 o , few sols before digging the trench. Drawn are profiles of v i /v i,ref .In the selected locations at the walls (points A) ice concentration in the regolith increases, or decreases depending of conditions. In both cases changes of v i /v i,ref are significant. At the south facing wall of the E-W oriented trench (within the first numerical cell) ice can disappear from the regolith within about 20 sols. At the bottom of the trough ice content in the regolith increases in all studied cases. The ice volume changes by about 20% just beneath the bottom, but already 5 centimeters beneath the trough changes are small.

Acknowledgments

This work was partially supported by the Polish Ministry of Education and Science (grant 0576/H03/2007/32).

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT Fig. 3. Temperature within the N-S oriented trench. The depth is 20 cm. The ice content in the regolith is either high (upper panel), or low (lower panel). In the former case v i,max = 0.38, while in the latter v i,max = 0.09. The profiles are for two points on the walls, 10 cm beneath the upper edge (points A), and in the middle of the trench bottom. Wall 1 is the east facing west wall and wall 2 is the west facing wall. In addition, profile for the surface external to the trench, 20 cm away of the trench edge is also shown.

A C C E P T E D M

A N U S C R I P T ACCEPTED MANUSCRIPT Fig. 6. Frost thickness as function of time in trenches with both orientations, when the ice content in the regolith is high and low. The profiles are for two points on the walls, 10 cm beneath the upper edge, and in the middle of the trench bottom. Walls are defined as in the figures above. Wind speed is 5 m/sec. The value of v i,max is 0.38 (upper and lower panels) and 0.09 (middle panel).

A C C E P T E D M A N U S C R I P T ACCEPTED MANUSCRIPT