Stéphanie C Bodin
email: s.bodin90@gmail.com

Christophe Vaschalde
email: christophe_vaschalde@yahoo.fr

David Ollivier
email: dollivier@mmsh.univ-aix.fr

First insights into the wood management for the production of lime given by the anthracological study of a 19 th century lime kiln from Martinique, Lesser Antilles

Keywords:

The colonial economy of Martinique, French West Indies, has long been based on sugar cane exploitation, which started in the 17 th century and decreased during the 19 th century. Sugar economy led to woodland exploitation for the production of fuel. In particular, this fuelwood was needed for lime production, which was used for the purification of cane juice. As lime kilns are numerous in Martinique, their functioning implied an intense wood exploitation. To date, no archaeobotanical studies were conducted to assess how lime kiln workers used and managed wood for the production of fuel. In this study, we analyzed 207 wood charcoals fragments from a 19 th century lime kiln at the Sucrerie plantation, Les Anses-d'Arlet, Martinique. This led to the identification of 50 taxa, showing that wood collection was mainly diffuse and unfocused. But the occurrence of taxa known as being used for fuelwood or charcoal production suggests that lime kiln workers were looking for wood with good fuel properties. Our results also show that the collection of wood was made locally, as revealed by the presence of a xero-to mesophilic vegetation in the anthracological record, which is in accordance with botanical and ethnological data. The presence of pioneer and heliophilic taxa also suggest that the vegetation was already impacted by wood exploitation i.e., local forest associations were largely disturbed by this time and secondary regrowth was conspicuous. Nevertheless, the high proportion of small-gauge wood charcoals fragments (36%) and the presence of edible taxa in the anthracological record suggest that saplings and pruning waste may have been preferred to preserve useful and mature trees. This study allows for the first time to assess which fuelwood was used by Antillean lime kiln workers during the colonial Period which makes it possible to better understand one of the most widespread handicrafts during the sugar exploitation in the Caribbean.

Introduction

For several decades, the development of archaeology for the colonial period has led to an interest in the plantation dwellings of the West Indies [START_REF] Moussette | Archéologie de l'Amérique coloniale française[END_REF][START_REF] Kelly | Bitation, archéologie des habitations/plantations des Petites Antilles[END_REF]. In Martinique, the inventory of sugar refineries and distilleries [START_REF] Mousnier | Atlas historique du patrimoine sucrier de la Martinique (XVIIè-XXè s.). L'Harmattan[END_REF]) allowed several research themes to emerge: the habitat of slaves [START_REF] Kelly | Historical Archaeology in the French Caribbean: An Introduction to a Special Volume of the Journal of Caribbean Archaeology[END_REF]2013) the production of ceramics [START_REF] England | Acculturation in the Creole context : a case study of La Poterie Martinique[END_REF]1997;[START_REF] Amouric | Projet collectif de Recherche Interrégional : Poteries des îles françaises de l'Amérique, productions locales et importées, XVIIe-XXe siècles. SRA MArtinique ; LA3M APG[END_REF][START_REF] Amouric | Projet Collectif de Recherche Interrégional « Poteries des îles françaises d'Amérique : XVIIe-XXe siècles. SRA Martinique, service des patrimoines de la Guadeloupe[END_REF][START_REF] Amouric | Projet Collectif de Recherche Interrégional « Poteries des îles françaises d'Amérique : XVIIe-XXe siècles. SRA Martinique, service des patrimoines de la Guadeloupe[END_REF][START_REF] Amouric | Projet Collectif de Recherche Interrégional « Poteries des îles françaises d'Amérique : XVIIe-XXe siècles. SRA Martinique, service des patrimoines de la Guadeloupe[END_REF]2013;[START_REF] Amouric | Projet Collectif de Recherche Interrégional « Poteries des îles françaises d'Amérique : XVIIe-XXe siècles[END_REF][START_REF] Thiriot | Fouiller les encyclopédistes : transfert de modèles aux Antilles françaises[END_REF]Cloarec-Quillon et al. 2015) and that of lime [START_REF] Vidal | Tout autour du four. La fabrication artisanale de la chaux en Martinique, structures et savoir-faire[END_REF]Verrand and Vidal 2004;[START_REF] Vaschalde | Fours à chaux et chaufourniers en Martinique[END_REF][START_REF] Louis-Alexandre | A new subfamily classification of the Leguminosae based on a taxonomically comprehensive phylogeny[END_REF]2019b). Currently, archaeobotany is developing thanks to Cultural Resources Management archaeology archaeology, but no detailed study of wood charcoal from the colonial period has ever been carried out in the Lesser Antilles. Yet, the economy of large plantations has impacted the Caribbean forested landscape. As early as 1635, new settlers had to clear the land in order to benefit from advantageous settlement conditions granted by the royal power (Blérald 1986, p25). An export economy developed rapidly, first based on tobacco cultivation, then on sugar cane [START_REF] Doumenge | Experimental anthracology: Evaluating the role of combustion processes in the representivity of archaeological charcoal records in tropical forests, a case study from the Maya Lowlands[END_REF]. In 1670, large land estates appeared and enabled the French West Indies to export sugar to Europe on a massive scale. The sugar plantation was the basis of the colonial economy. It concentrates both the cultivation of cane, which implies the acquisition of land, and the means to transform it into sugar (mill, lime kiln...). For nearly two centuries, no innovation allowed production techniques to progress. The growth of the sugar economy was based solely on the increase in the slave population, the multiplication of plantations and the extension of cultivated lands (Blérald 1986, p27).

In the 18 th century, sugarcane cultivation often occupied less than half of the surface area of sugar estates (Blérald 1986, p27;Bégot 2008, p51;Bonnet 2008, p167), and a large part of the estates were covered with woods and savannahs, which were essential for livestock and firewood (Bégot 2008, p50). However, the lack of management of forest resources led to a shortage of timber, which was then imported (Bégot 2008, p50). In the 19 th century, profound changes affected the sugar economy: the modernization of the means of production (adoption of steam mills, shareholding), the abolition of slavery in 1848, and the competition from European beet sugar. However, the areas cultivated with cane continued to increase thanks to the development of the rum industry. Between 1855 and 1895, these areas increased by 28%, and Martinique became the first world exporter of rum (Blérald 1986, p60 and 74-79). During World War I, the number of distilleries doubled, as the production of sugar and alcohol was stopped in France because of the conflict. From the middle of the 20 th century onwards, many factories closed down, and the cultivated areas decreased.

Thus, between the end of the 17 th and the 20 th century, the slave sugar economy, then the capitalist one, led to the clearing of large areas of land. Moreover, the supply of fuel was indispensable to the West Indian economy. Before the arrival of industrial techniques and of coal at the end of the 19 th century, wood and bagasse were the only fuels available. Fire was used to heat sugar boilers [START_REF] Vidal | Le fonctionnement et l'évolution du système de chauffe des habitations-sucreries traditionnelles, à la Martinique, entre le XVIIe et le début du XIXe siècle[END_REF], potters' kilns that produced sugarloaves and molasses pots [START_REF] England | Acculturation in the Creole context : a case study of La Poterie Martinique[END_REF]1997;[START_REF] Thiriot | Fouiller les encyclopédistes : transfert de modèles aux Antilles françaises[END_REF], as well as lime kilns. Many kilns provided the lime necessary for the purification of cane juice. At the end of the 18 th century, Martinique had nearly 80 lime kilns (Verrand and Vidal 2004, p41). An archaeological inventory has allowed to identify more than 130 of them on the island (Vaschalde et al. 2019b). The impact of this activity on the vegetation has thus been significant, and may have required the development of techniques for woodland management.

The colonial archives do not provide any clues on this woodland management, and only anthracology is able to enlighten us. In 2014, an archaeological excavation was carried out on a 19 th lime kiln at the Sucrerie plantation, Les Anses-d'Arlet, and allowed the extraction of wood charcoal. Anthracology in the tropics is challenging because of the scarcity of the studies compared with those from temperate or Mediterranean environments; this makes comparisons difficult. In the Neotropics, anthracological studies are most common in South America (e.g., [START_REF] Tardy | Feux, sols et écosystèmes forestiers tropicaux[END_REF]Scheel-Ybert 2001;[START_REF] Beauclair | Fire and ritual: bark hearths in South-American Tupiguarani mortuary rites[END_REF][START_REF] Bianchini | Plants for life and death: evidence of use of plant resources in funerary activities of shellmound builders through the anthracological analysis of Jabuticabeira-II site[END_REF][START_REF] Fernandes Caromano | Revealing Fires and Rich Diets: Macro-and Micro-archaeobotanical Analysis at the Hatahara Site, Central Amazonia[END_REF][START_REF] Scheel-Ybert | The forest people: landscape and firewood use in the Araruama region, southeastern Brazil, during the late Holocene[END_REF][START_REF] Bachelet | Landscape and firewood selection in the Santa Elina rock shelter (Mato Grosso, Brazil) during the Holocene[END_REF][START_REF] Bodin | Unraveling pre-Columbian occupation patterns in the tropical forests of French Guiana using an anthracological approach[END_REF] but are still rare in the Antilles [START_REF] Newsom | ): the sites of Norman Estate, Anse des Pères and Hope Estate with a contribution to the "La Hueca problem[END_REF]. To our knowledge, the only anthracological study conducted in Martinique is that by Christophe Tardy on the pre-Columbian site of Vivé [START_REF] Tardy | Rapport anthracologique du site archéologique de Vivé (Martinique)[END_REF]. In this paper, we assess how lime kiln burners managed woodland to produce charcoal for the functioning of a 19 th century lime kiln, in the context of the sugar cane exploitation in Martinique. Specifically, we ask these questions to better understand the fuel economy:

1. Did craftspeople favor some kind of vegetation? 2. Did they collect wood locally?

To answer these questions, we identified wood charcoals found in the remains of the lime kiln. This study, combining anthracology with an ethnoarchaeological approach, is the first one for a Caribbean colonial site.

Material and methods

Regional setting

The study was undertaken in Les Anses-d'Arlet (14°29′26″, 61°04′50″), in Martinique, a French overseas department located in the Lesser Antilles (Fig. 1). The climate is tropical humid (Peel et al. 2007) with a mean annual rainfall of 1300 mm and a mean annual temperature of 27°C recorded in Les Anses-d'Arlet between 1951and 1970[START_REF] Portécop | Phytogéographie, cartographie écologique et aménagement dans une île tropicale[END_REF]. The dry season extends from January to April and the rainy season occurs between June and November. In the South, the geological substratum is older and the topography lower than in the North-Western part of the island where the Montagne Pelée stands. The climate is also drier in the South [START_REF] Portécop | Phytogéographie, cartographie écologique et aménagement dans une île tropicale[END_REF]. Martinique belongs to the Caribbean hotspot of biodiversity, which hosts more than 11,000 species of native seed plants with a species endemism rate of about 72% [START_REF] Myers | Biodiversity hotspots for conservation priorities[END_REF]Conservation International 2020). The most recent update of the floristic checklist counts about 3200 angiosperm species in Guadeloupe-Martinique [START_REF] Fournet | Progrès récents dans la connaissance de la flore de Guadeloupe et de Martinique[END_REF]. In 1994, Fiard reports 396 tree species of which 20% are endemic to the Lesser Antilles. According to [START_REF] Portécop | Phytogéographie, cartographie écologique et aménagement dans une île tropicale[END_REF], the flora of Martinique has a strong affinity with that of Puerto Rico and, in a lesser extent, with that of the North of South America. In Les Anses-d'Arlet area, the vegetation is xerophytic, with a Tabebuia-Bursera forest, a Ceiba-Ficus sandy-clay facies and a Lonchocarpus-Inga riparian forest [START_REF] Portécop | Phytogéographie, cartographie écologique et aménagement dans une île tropicale[END_REF]).

Historical context

During the second half of the 18 th century, the Sucrerie plantation was one of the most important sugar plantation in southwest Martinique. Around 1770, the map of Moreau du Temple indicates that the domain was called "Lebrun" [START_REF] Bégot | La Martinique de Moreau du Temple, 1770. La carte des ingénieurs géographes[END_REF]. In 1783, the Delpech family built a sugar factory, and in 1788, the plantation was bought by Jean-François Hayot (1733-1788) (Bruneau-Latouche et al. 2002, p569;Elisabeth 2005, p1). During the 19 th century, the Hayot family established itself throughout southern Martinique and diversified its economic activity (sugar, rum, pottery, building materials, etc.). It increased the extent of the Sucrerie plantation, which reached 116 ha in 1788, 177 ha in 1820, and 284 ha in 1828 (Bruneau-Latouche et al. 2002, pp569-573;[START_REF] Ollivier | Habitation La Sucrerie : Fouille du four à chaux. La Martinique, les Anses-d'Arlet[END_REF]. Starting in the 1840s in Martinique, the creation of modern factories equipped with steam boilers disrupted the operation of the Sucrerie plantation. Around 1845, Michel Hayot (1786-1852) put an end to the pressing of sugar on the domain, and sent the cane to the new sugar factory at Pointe-Simon (Fort-de-France). This factory was founded by Joseph Quennesson (1819-1890), which bought the plantation in 1878 [START_REF] Ollivier | Habitation La Sucrerie : Fouille du four à chaux. La Martinique, les Anses-d'Arlet[END_REF]. After the closure of the factory, Simon Hayot bought the Sucrerie plantation in 1899. Taking advantage of the development of the rum industry, he set up a distillery in 1918 [START_REF] Louis-Alexandre | A new subfamily classification of the Leguminosae based on a taxonomically comprehensive phylogeny[END_REF]. In 1926, the map of the Compagnie française aérienne indicates that sugar cane plantations occupy large areas in the plain around the house, and that the hills are forested. The distillery closed definitively in 1950 [START_REF] Mousnier | Atlas historique du patrimoine sucrier de la Martinique (XVIIè-XXè s.). L'Harmattan[END_REF]Flohic et al. 1998, p29).

Remains associated to the Sucrerie plantation

The plantation is organized around a hill overlooking the Ravine des OEillets and the village of Les Anses-d'Arlet. The Sucrerie plantation is composed of a master's mansion probably dating from 1760-1770 AD [START_REF] Charlery | Maisons de maître et habitations coloniales dans les anciens territoires français de l'Amérique tropicale[END_REF]) installed on the hilltop, and of various buildings related to its past economic activity. The only vestiges that can be linked to the sugar activity of the plantation in the 18th-19th centuries are those of the sugar mill that has been identified south of the ruins. The 20th century distillery, built on the sugar factory, stretches from east to west on the northern piedmont of the hill. It consists of a large building, restored in 2011 and containing several vats or vat supports, a steam boiler and the mechanical installation still in place, also restored. Other ruins are located to the north and west. In 2014, the archaeological excavation uncovered walls of basalt rubble bound with lime mortar and partially plastered (Fig. 2; MR 1002 and MR 1025). Dating from the second half of the 18 th century, their function remains unknown.

Lime kiln structure and functioning

The lime kiln is built to the North-East of the distillery. It is a circular building measuring approximately 5 m in diameter for a preserved height of 2.20 m (Fig. 2). The external facing is built in irregular apparatus of large basalt rubble bound with a lime mortar, the whole being covered with a smooth lime plaster, like the walls previously described.

To the West is the mouth of the kiln, facing a ruined building that may have served as a "lime hut". The mouth is composed of a basalt sill and refractory stone pedestals. The lintel topping the whole has disappeared. The mouth is preceded by a large splay, which sheltered the working area, partially filled by alternating layers of lime, containing fragments of charred coral and wood mixed with ash, which testify to the activity of the kiln (US 1015). Inside the kiln, the firing chamber has a diameter of 1.50 m at the base. It has the shape of an inverted truncated cone. The floor consists of a paving laid on a bed of lime mortar and the walls are built with cut refractory stones.

To the north of the kiln, a masonry ramp allowing access to the top in order to feed it with corals has been identified despite the poor state of conservation.

The architecture of the structure indicates that it is a short-flame kiln (Vaschalde 2018, pp91-97), a process that is very widespread in Martinique [START_REF] Vidal | Tout autour du four. La fabrication artisanale de la chaux en Martinique, structures et savoir-faire[END_REF]Verrand and Vidal 2004;[START_REF] Vaschalde | Fours à chaux et chaufourniers en Martinique[END_REF][START_REF] Louis-Alexandre | A new subfamily classification of the Leguminosae based on a taxonomically comprehensive phylogeny[END_REF]2019b). The raw material, composed of coral and shell fragments, is charged by forming beds alternating with fuel (wood charcoals). At the end of each firing, the lime is discharged through the mouth at the bottom of the kiln, and part of the charcoal is spread on the floor of the work area.

Dating of the lime kiln

The dating of the lime kiln is based on stratigraphy, textual sources and artifacts found during the excavation. After the walls MR 1002 and MR 1025 were levelled, the ground was backfilled and the kiln's foundations were then dug. An inventory drawn up in 1820 does not mention the lime kiln, indicating that it was built after that date. In 1878, the lime kiln appears in the deed of sale of the estate to Joseph Quennesson [START_REF] Ollivier | Habitation La Sucrerie : Fouille du four à chaux. La Martinique, les Anses-d'Arlet[END_REF]. It can be assumed that it was in good condition at this time. Thereafter, it is probable that Joseph Quennesson did not have lime produced at the plantation anymore, since the canes were brought to the factory of Pointe-Simon. Finally, the excavation of the work area uncovered two fragments of glass bottle bottoms, dated to the mid-19th century. The lime kiln was therefore in operation in the mid-19th century, and seems to have been abandoned after a few decades of use.

Charcoal sampling and identification

Sampling was done in a deposit of mixed charcoal, coral and stone fragments in front of the opening of the kiln from a single layer, the stratigraphic unit 1015 (Fig. 2), which probably corresponds to the last firing. The 2-liters deposit was collected in clods with a trowel to avoid charcoal breaking. It was then water-sieved in 4 and 2 mm meshes to retrieve charcoal fragments. Charcoal fragments were then dried away from the direct sunlight to prevent them from bursting.

Only fragments > 4 mm were considered for anthracological analysis [START_REF] Chabal | L'anthracologie[END_REF]Scheel-Ybert 2001). Below this size, the identification is harder because fewer anatomical features are visible. This is especially true for the tropical material where there is a strong anatomical homogeneity. Charcoals were observed under a reflected light microscope with dark and bright field, at x100, x200, x500 and x1000 magnifications. In the absence of a reference collection of charred wood from the Lesser Antilles at the Institut des Sciences de l'Evolution de Montpellier, identification was performed with the help of the charcoal reference collection from French Guiana (Scheel-Ybert 2016), as many taxa are common between the two regions. Specific literature [START_REF] Détienne | Evaluating methods and results for the application of anthracology to high diversity and high endemism environments: Case study in the Tiwaka Valley, north-eastern Grande Terre of New Caledonia[END_REF][START_REF] Cirad ; Quillon | Production locale de formes à sucre pour l'industrie sucrière de la Martinique, l'exemple de la poterie Dalençon au Marin[END_REF][START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]) and identification tools (InsideWood 2004-onwards; [START_REF] Bodin | CharKey: An electronic identification key for wood charcoals of French Guiana[END_REF] were also used to complement the comparisons with the reference collection. Many tropical woods are characterized by a homogeneous anatomy, making them difficult to distinguish, especially within the same family. Identification at the genus level is possible in the case of a monogeneric family or when genera present distinguishable anatomical features. Similarly, identification at the species level is only possible in the case of a monospecific genus or when species are anatomically heterogeneous. Poor preservation of charcoal masking diagnostic features or the presence of features common to several taxa led to the designation of morphological types (e.g., type Lonchocarpus or Leguminosae type Lonchocarpus). A quasi-secured identification is preceded by "cf." for confer (e.g., cf. Spondias or Spondias cf. mombin). When a taxon has an anatomical affinity with another taxon but differs from it in major aspects, it is preceded by the mention "aff." for affinis. Finally, when two species may match with the observed features, they are both mentioned (e.g., Spondias dulcis/mombin). The nomenclature used follows that of the Angiosperm Phylogeny Group (2016) and that of the Legume Phylogeny Working Group (2017) for the Legume family.

Results

The anthracological analysis resulted in the identification of 183 charcoal fragments out of 207, distributed in 50 taxa including 16 morphological types and 7 indeterminate (Fig. 3). Among the 207 charcoal pieces, 24 were indeterminable due to a very poor preservation (e.g., vitrification) or because of a distorted anatomical structure (knots). The anthracological assemblage contains an important proportion of Leguminosae (45%) distributed in 19 types of which 6 are assigned to genera or species which have a close anatomy. Among the 207 analyzed charcoal fragments, 29% come from twigs and 6% come from knots.

Discussion

Palaeological significance: a firewood collected locally

The high number of types in our anthracological spectrum reflects the need of a specific reference collection for study areas with peculiar vegetation composition: even if the reference collection of charred woods from French Guiana helped through the identification process, many taxa are endemic to the Lesser Antilles or simply do not occur in French Guiana. However, building a reference collection is a long process which is not compatible with short projects. This is an issue that was previously underlined during a pioneering research in New Caledonia, another tropical region where anthracological studies are still scarce (Dotte-Sarout 2011). Hence, due to the high number of morphological types, the following interpretation relies on a reduced list of taxa that may shed light on the collect of wood used as fuel for the lime kiln. This short list already makes it possible to formulate hypotheses.

First, some taxa found in the anthracological assemblage may indicate that wood collection was made in a xerophilic forest. These are Spondias mombin, cf. Zygia latifolia, an Apocynaceae resembling Plumeria sp. and the types Acacia macracantha and Senna bicapsularis. To this list can be added the types Psidium guajava and Citharexylum which have xero-mesophilic affinities and Coccoloba/Sapindaceae which may also correspond to a xerophilic facies [START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]. It is noteworthy that the porknut Acacia macracantha, which was the most abundant taxon in the charcoal assemblage, and the golden apple Spondias dulcis, an introduced species from Oceania, are used as firewood and for the production of charcoal [START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]. Some species of the genera Citharexylum, Coccoloba and Plumeria are also used as firewood and to produce charcoal in the French West Indies, such as the fiddlewood Citharexylum spinosum, the seagrape Coccoloba uvifera and the frangipani Plumeria alba [START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]ONF 2015).

There are also taxa corresponding to a mesophilic vegetation, which suggest that wood collection was not restricted to the xerophilic forest. These are the breadfruit Artocarpus altilis, the amarillon Lonchocarpus cf. heptaphyllus, which is used as firewood in the Lesser Antilles [START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]ONF 2015) and Melastomatoideae. To these taxa can be added the type Cassipourea guianensis and two Inga morphotypes.

Second, other taxa may indicate a pioneer and/or heliophilic vegetation such as Spondias dulcis/mombin, Melastomatoideae, Salicaceae, cf. Zygia latifolia, Coccoloba/Sapindaceae and the type Acacia macracantha. This suggests that the fuelwood was collected in a secondary forest or in disturbed areas. Three morphological types found in the anthracological assemblage might support the hypothesis of a collection from a secondary forest association. These are the types Citharexylum, Cordia and Psidium guajava. The species of the genus Citharexylum are heliophilic, the common guava P. guajava is a pioneer and some Cordia species are short-lived pioneers, such as C. alliodora [START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]. Interestingly, some Cordia species are also used as firewood in the French Antilles, such as the white manjack C. sulcata [START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]ONF 2015). The presence of pioneer and heliophilic taxa may be the result of an intense landscape exploitation. Third, the presence of possible riparian taxa, such as the amarillon Lonchocarpus cf. heptaphyllus, cf. Zygia latifolia, Spondias mombin and the type Senna bicapsularis, may also suggest that some woods were collected near water courses [START_REF] Guppy | Plants, seeds, and currents in the West Indies and Azores; the results of investigations carried out in those regions between 1906 and[END_REF][START_REF] Portécop | Phytogéographie, cartographie écologique et aménagement dans une île tropicale[END_REF][START_REF] Rollet | Arbres des Petites Antilles. Tome 2 : Description des espèces[END_REF]ONF 2015).

Obviously, the uncertainty of some identifications (e.g., morphological types) requires caution in interpreting their presence. Nevertheless, the presence of xero-to mesophilic, secondary forest and riparian taxa in the anthracological record is coherent with the present-day flora found near the Sucrerie plantation. This suggests that the craftspeople collected the firewood in the vicinity of the plantation, and, given the high number of taxa represented in the anthracological record (50 potential taxa) without any preferential collection of certain species for fuel use. However, it is possible that the craftspeople favored woods with good fuel properties, as suggested by the presence of several taxa known to be used as firewood or for charcoal production.

The fact that the wood was collected in the near environment of the lime kiln is consistent with ethnological data from another lime kiln site in Martinique. In Le Lamentin, near Fort-de-France, a retired craftsman told that he used the wood of different mangrove tree species (Vaschalde et al. 2019b). In the bay of Fort-de-France, the mangrove vegetation is abundant, especially near Le Lamentin, which reinforces the idea of local supply. According to the craftsmen of Le Lamentin, mangrove wood was appreciated by the lime kiln workers. In Guadeloupe, this use is attested since the 19 th century (Archives de Médecine Navale 1869). In Les Anses-d'Arlet, the mangrove is absent, and no mangrove charcoal was found in the lime kiln. As in Le Lamentin, artisans have harvested wood around the site, and have adapted their practices to local taxa. 4.2 On the most abundant taxa: a reflect of the palaeoenvironment or a differential preservation?

The most abundant taxa of the anthracological record are the type Acacia macracantha, Lonchocarpus cf. heptaphyllus and the type Psidium guajava which count respectively for 14%, 10% and 8.2% of the total number of charcoal fragments, representing together 32.2% of the anthracological record. As there is a large amount of other taxa represented, reaching a total of 50, it is unlikely that these three taxa were preferred as fuelwood and specifically selected by the craftspeople. The dominance of these taxa in the anthracological record may simply reflect their abundance in the environment rather than a wood preference. Similar observations have been made in Europe for potters and lime kiln workers from the historical period [START_REF] Chabal | Forêts et sociétés en Languedoc (Néolithique final, Antiquité tardive). L'anthracologie, méthode et paléoécologie[END_REF]Vaschalde 2018). Another explanation for their dominance could be a differential preservation leading to an over-fragmentation of charcoals. Several studies have shown that post-depositional processes affect charcoal preservation and can induce a bias in the proportions between taxa [START_REF] Labat | Temperature, compression and fragmentation: an experimental analysis to assess the impact of taphonomic processes on charcoal preservation[END_REF][START_REF] Théry-Parisot | Anthracology and taphonomy, from wood gathering to charcoal analysis. A review of the taphonomic processes modifying charcoal assemblages, in archaeological contexts[END_REF][START_REF] Chrzazvez | Impact of post-depositional processes on charcoal fragmentation and archaeobotanical implications: experimental approach combining charcoal analysis and biomechanics[END_REF][START_REF] Doumenge | Experimental anthracology: Evaluating the role of combustion processes in the representivity of archaeological charcoal records in tropical forests, a case study from the Maya Lowlands[END_REF][START_REF] Henry | Old taphonomy issues, new charcoal data for Mesolithic contexts: Impact of fragment size and sampling context on the assemblages of Escabasses cave (SW France)[END_REF]. To assess the likelihood of the hypothesis of a differential preservation, we would need this kind of data for each of our taxa, which would require time-consuming experiments. Furthermore, this would be applicable only for those which have been securely identified, not for the morphological types. Testing this hypothesis is therefore very difficult.

Collection of cultivated and/or edible species as fuelwood

As pointed above, the large amount of potential taxa (i.e. including the morphological types) found in charcoals suggest that there was no selection of a specific taxon for fuel use. Even edible species were collected: the breadfruit Artocarpus altilis, a species introduced in the late 18 th century [START_REF] Leakey | Breadfruit reconnaisance study in the Caribbean region[END_REF], and the hog plum or the golden apple (Spondias mombin or S. dulcis), also both introducedprobably much earlier concerning the hog plum, by natural seed dispersion [START_REF] Guppy | Plants, seeds, and currents in the West Indies and Azores; the results of investigations carried out in those regions between 1906 and[END_REF]. To these taxa might be added the common guava Psidium guajava. The use of edible species as fuel for lime kilns in Martinique has already been reported by Father Labat in the late 17 thearly 18 th century with the orange tree Citrus sp., introduced in the 16 th century [START_REF] Labat | Temperature, compression and fragmentation: an experimental analysis to assess the impact of taphonomic processes on charcoal preservation[END_REF][START_REF] Bellec | A la découverte des fruits des Antilles[END_REF], and more recently by [START_REF] Vidal | Tout autour du four. La fabrication artisanale de la chaux en Martinique, structures et savoir-faire[END_REF] with the mango tree Mangifera indica, introduced in the 18 th century (Le [START_REF] Bellec | A la découverte des fruits des Antilles[END_REF]). Therefore, it is not surprising to find edible species in the wood used for the functioning of lime kilns. The charcoal from the lime kilns can thus testify to the taxa exploited and/or cultivated within an exploitation specialized in sugar production. In Europe, the presence of exploited/cultivated taxa in anthracological assemblages results from the pruning of fruit trees, as it is the case with vine for example (Ros and Vaschalde, in press). At the Sucrerie plantation, pruning waste from Spondias dulcis/mombin, Artocarpus altilis and Psidium guajava could thus have been harvested for fuel use.

4.4 On the presence of charcoals of twigs in the fuel About 30% of the charcoal fragments came from twigs, as revealed by the non-parallel arrangement of the rays in the transverse section. If we consider that the knots also come from twigs or branches, this percentage reaches 36%. It is therefore likely that small-gauge wood was favored by the craftspeople. Several reasons could explain this. First, the abundance of small-gauge wood in the environment, such as in a semi-open forest with small trees and shrubs as suggested by the presence of pioneer and heliophilic taxa (see above). This type of vegetation would be in agreement with the intensive land exploitation by the sugar economy of the 19 th century. This has been observed with Mediterranean lime kilns fueled by wood from the garrigue or the maquis [START_REF] Vaschalde | Charcoal analysis of lime kiln remains in Southern France: an original process of mediaeval and modern traditional lime burning[END_REF]2019a;Vaschalde 2018). Second, this could be due to wood management practices. The managers of the plantation could have wanted to preserve the useful parts of the trees (trunks, main branches) by forbidding them to be cut down. Artisans would have used only small branches, dead wood, or brushwood. This type of practice is attested in the south of France during the modern period (Vaschalde 2018), and would be in agreement with the recovery of supposed pruning waste from fruit trees found in the lime kiln of the Sucrerie plantation.

Conclusion

The study of wood charcoals from the lime kiln of the Sucrerie plantation is one of the firsts from the colonial period in the West Indies. It shows that the lime burners used wood harvested from local vegetation composed of a great diversity of species. This vegetation was open to semi-open, and appears to have been highly impacted by the exploitation of resources. However, tree management practices may have been applied, explaining the important use of small-gauge wood. This management would notably concern cultivated and/or exploited trees and food crops in the plantation. Anthracology makes it thus possible to address issues that cannot be dealt with on the basis of a study of the archives. Thanks to this archaeobotanical approach, we get a better understanding on how lime burners managed woody vegetation during the colonial period. Further research is however needed on this kind of activity in the Antilles, and it is necessary to study a greater number and diversity of archaeological sites from the colonial period to able comparisons between practices used in other islands of the archipelago. Finally, it is also necessary to carry out floristic surveys of the current vegetation around the sites in order to highlight the impact of the colonial sugar economy on the presentday environment.

Acknowledgements

We acknowledge an Investissement d'Avenir grant from the Agence Nationale de la Recherche (CEBA: ANR-10-LABX-25-01).

The archaeological excavation and the project was granted by the DAC Martinique and the Fondation Clément. We are grateful to Arkaeos association for supporting the project. We also thank L. Paradis and A. Affouard for their help in the realization of Figure 1.

Declarations of interest: none. SB, CV and DO wrote the manuscript. CV did fieldwork (charcoal extraction and sieving), DO led the archaeological excavation and SB did charcoal identification.


Laurent Beauregard

Annafederica Urbano

Stéphanie Lizy-Destrez

Joseph Morlier

Multidisciplinary design and architecture optimization of a reusable lunar lander

With renewed interest in lunar exploration and the upcoming deployment of the lunar space station, the Lunar Orbital Platform-Gateway (LOP-G), the scientific community is focusing on the design of a lander to bring people back to the lunar surface. This work focuses on optimizing two aspects of the lunar lander concurrently: the mission architecture and the vehicle design, often treated independently in the literature. A methodology is introduced to enumerate and preliminarily

bring people back onto the lunar surface. The international community has agreed to deploy an orbital station around the Moon, named the Lunar Orbital Platform-Gateway (LOP-G) in the 2025 timeframe. This station is planned to be used as a hub. The crew would first travel to the LOP-G, there they would transfer to a lander to reach the lunar surface. Their return would pass by the LOP-G to transfer back into the Orion capsule, allowing a return to Earth. While it has been decided that the Space Launch System (SLS) and Orion will bring people to and from the LOP-G, the vehicle that will bring the crew unto the lunar surface has not been finalized yet. Several architectures have been proposed in the past, including by Lockheed Martin [START_REF] Cichan | Concept for a Crewed Lunar Lander Operating from the Lunar Orbiting Platform-Gateway[END_REF] which proposed a reusable one stage lander with liquid hydrogen as fuel. Another proposal was by NASA with the Altair [START_REF]America's Spacecraft for a New Generation of Explorers The Altair Lunar Lander[END_REF] lander: a two-stage lander similar to Apollo with the landing stage using hydrogen and the upper stage using storable propellants. Due to the high cost of bringing materials and people around the Moon, it is vital that the vehicles are optimized with a metric that correlates strongly with the total cost of the mission. To obtain an optimal vehicle architecture according to this metric, Multidisciplinary Design Optimization (MDO) will be employed. MDO techniques have gained significant traction in the last decades and been applied in the field of aerospace engineering, including airplane design [START_REF] Gray | Aero-propulsive Design Optimization of a Turboelectric Boundary Layer Ingestion Propulsion System[END_REF] with aeropropulsive design optimization, expendable launcher design [START_REF] Castellini | Global and Local Multidisciplinary Design Optimization of Expendable Launch Vehicles[END_REF] with global and local optimization approach, reusable launchers design [START_REF] Brevault | Multi-Objective Multidisciplinary Design Optimization Approach for Partially Reusable Launch Vehicle Design[END_REF] with multilevel, multimission, decoupled MDO formulation and satellite design [START_REF] Gray | Automatic Evaluation of Multidisciplinary Derivatives Using a Graph-Based Problem Formulation in OpenMDAO[END_REF] with millions of design variables. In this work a MDO tool developed by NASA Glenn Research Center called OpenMDAO will be used.

OpenMDAO has seen several successful applications [START_REF] Gray | OpenMDAO: An open-source framework for multidisciplinary design, analysis, and optimization[END_REF]. A vast amount of work has been done on lunar landers/orbiters, analyzing in detail their components such as the propulsion system, structure, life support system, reusability, maintenance, GNC and so on [START_REF] Ntrs | Lunar Lander Conceptual Design[END_REF][9] [START_REF] Cichan | Concept for a Crewed Lunar Lander Operating from the Lunar Orbiting Platform-Gateway[END_REF]. Recent work has focused on a MDO approach to reusable lunar cargo lander design [START_REF] Robertson | A Conceptual Design Study for an Unmanned , Reusable Cargo Lunar Lander[END_REF].

In the context of this article, a mission architecture is defined as: Number of vehicles, sequence of transfers and number of reuses. The choice of these variables impacts massively the lunar lander design process. Work on simultaneous trajectory and vehicle optimization for a lunar supply chain have been done before [START_REF] Taylor | Concurrent Trajectory and Vehicle Optimization: A Case Study of Earth-Moon Supply Chain Logistics[END_REF]. However, for lunar lander missions, there is a lack of literature that treats the optimization of the mission architecture and the vehicle design concurrently. Most literature on lunar landers consider either fixed mission architecture, or only a few alternative choices [START_REF] Ntrs | Lunar Lander Conceptual Design[END_REF]. The present work aims at pushing an understudied aspect of the lunar mission design: the concurrent optimization of the mission architecture and vehicle design. This is done in a two-step optimization process. In the first step, a systematic approach to generate all possible mission architectures is introduced. Those mission architectures are then classified using a basic metric. This first step in known in the literature as the exhaustive enumeration (EE) method or combinatorial optimization [START_REF] Venkataraman | Applied Optimization with Matlab Programming[END_REF]. For the second step, the best mission architectures are combined with a multidisciplinary design optimization process to obtain a series of optimal vehicle designs. The workflow for this article is shown in Fig. 1. The result of this study is a ranking of the mission profiles and the number of reuses as a function of the average mass to perform the mission. The relevant results will be compared to NASA's proposed "Notional Human Landing System Reference Architecture" [START_REF]HUMAN LANDING SYSTEM BAA, INDUSTRY FORUM[END_REF]. The complete results are displayed in tables 10, 11 and 12. The goal of this work is to assist a system engineer in choosing the best mission architecture for a lunar mission.

This paper is composed of 7 chapters. The underlying theory and methods used throughout this work will be introduced in section II. Section III explores the first phase of the work: mission architecture enumeration. Section IV presents the engineering models for each subsystem considered for the lunar lander. The list of tools used in this article is shown in figure 1. Section V will consider the MDO process, the objective function and the optimization. Section VI will present the results obtained and, lastly, section VII will conclude this work. As seen in figure 1, this research is broken down into two parts. First, the mission architectures are exhaustively enumerated in a discrete way. A MDO framework is then built to analyze the vehicles. Finally the MDO framework is coupled with the mission architectures to produce a family of preliminary designs of optimized vehicles.

II. Theory and methods

To describe and optimize the lunar lander, several tools are used. First, the field of multidisciplinary design optimization (MDO) is described in II.A. Optimal control theory is introduced in II.B. Lastly, section II.C deals with the models of orbital mechanics used in this work. The tools of mission architecture enumeration and engineering models of the lunar lander subsystems will be described in their own sections in III and IV respectively.

A. Multidisciplinary design optimization

A multidisciplinary design optimization problem can be brought into the following canonical form: [START_REF] Martins | A Short Course on Multidisciplinary Design Optimization[END_REF]:

min x,y,z F(x, y, z) (1a)
Subject to

y = C(x, y, z) (1b) 0 = R(x, y, z) (1c) 0 ≤ g(x, y, z) (1d)
where x ∈ R N are the design variables, y ∈ R D are the coupling variables and z ∈ R S are the state variables.

N is the number of design variables of the system, D is the number of coupling variables of the system, S is the number of state variables necessary to describe the system. I is the number of inequality constraints.

F : U F ⊆ R N × R D × R S →
R is the objective function to be minimized.

C : U C ⊆ R N × R D × R S → R D represents the disciplines of the system. R : U R ⊆ R N × R D × R S → R S
are the residual of the state variables describing the system.

g : U I ⊆ R N × R D × R S → R I
represents the inequality constraints of the system.

To numerically solve these optimization problems, a nonlinear program software must be used. Examples of commonly used software include IPOPT [START_REF] Wächter | Line Search Filter Methods for Nonlinear Programming: Motivation and Global Convergence[END_REF], SNOPT [START_REF] Gill | SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization[END_REF] or SLSQP [START_REF] Kraft | A software package for sequential quadratic programming[END_REF]. In this work, the MDO framework used is NASA's python library OpenMDAO [START_REF] Gray | OpenMDAO: An open-source framework for multidisciplinary design, analysis, and optimization[END_REF].

B. Optimal control

Several instances of optimization problem must performed over a set of functions. A special case of those are called optimal control problems. Most optimal control problem can be brought into the form [START_REF] Bryson | Applied Optimal Control: Optimization, Estimation, and Control[END_REF]:

min u S (2a)
Where

S = ∫ t f t 0 L(t, x, u) dt + φ(t 0 , x 0 , t f , x f ) (2b) Subject to x = f (t, x, u) (2c) 0 = Ψ(t 0 , x 0 , t f , x f ) (2d)
Where t 0 is the initial time, t f is the final time, x 0 is the initial state and x f is the final state, S is the (scalar) cost function to be minimized over the trajectory, u(t) is the control of the system which must be found, x(t) is the state of the system, f is the dynamics of the system, L is the running cost, φ is the boundary cost and Ψ is the boundary constraint.

In this work, the software used to deal with optimal control problems is DYMOS, a python subpackage of OpenMDAO [START_REF] Falck | Optimal Control within the Context of Multidisciplinary Design, Analysis, and Optimization[END_REF]. It is a direct method that uses a Gauss-Lobatto [START_REF] Herman | Direct Optimization Using Collocation Based on High-Order Gauss-Lobatto Quadrature Rules[END_REF] or Radau Pseudospectral [START_REF] Garg | A Unified Framework for the Numerical Solution of Optimal Control Problems Using Pseudospectral Methods[END_REF] direct transcription method to convert the optimal control problem into a non-linear program, for which the optimization software described at the end of section II.A can be used. In this particular work, the Gauss-Lobatto method was used due to lower computational time.

C. Orbital mechanics

Two orbital models are used in this work depending on the location of the orbiter. At lunar launch and in low lunar orbit (LLO), the gravity of the Moon dominates the dynamics and so two-body dynamics is considered. For farther orbits, e.g. near the LOP-G, the dynamic of trajectories are modeled by an Earth-Moon circular restricted three-body problem (CR3BP) [START_REF] Ross | Dynamical Systems, the Three-Body Problem, and Space Mission Design[END_REF].

Low Lunar Orbits

For the lunar surface to low lunar orbit (LLO) trajectory, the gravity of the Moon dominates. During an ascent and descent it can be assumed that the motion is restricted to a plane [START_REF] Wilhite | Lunar Module Descent Mission Design[END_REF]. In an inertial frame around the Moon, the dynamics can be described by the following system of equations:

r = v r (3a) θ = v t r (3b) v r = - µ r 2 + v 2 t r + T m sin(α) (3c) v t = - v r v t r + T m cos(α) (3d) m = - T g 0 I sp (3e)

Fig. 2 Geometry of low lunar orbit dynamics

Where r is the radial position of the spacecraft, θ is the angular position of the spacecraft, v r is the radial component of the velocity, v t is the tangential component of the velocity, m is the mass of the spacecraft, T is the thrust of the spacecraft, g 0 = 9.80665 m/s 2 is the standard gravity, I sp is the specific impulse of the rocket engine and α is the angle of attack of the spacecraft.

Circular restricted three body problem

For objects that orbit farther than low lunar orbit, the gravity of the Earth must be accounted for. The model chosen in this work is the CR3BP which assumes that i) only the gravity of the Earth and the Moon are considered and ii) the Moon and the Earth have circular orbits. Usually the dynamic is expressed in the rotating frame of the Earth-Moon system whose origin is the center of mass of the system. In that case, the equation of motions are:

ì r = - µ E |ì r -ì r E | 3 (ì r -ì r E ) - µ M |ì r -ì r M | 3 (ì r -ì r M ) -ì ω × ( ì ω × ì r) -2 ì ω × ì r (4)
Where ì r is the vectorial position of the spacecraft, ì r E and ì r M are the vectorial positions of the Earth and Moon respectively (at rest in the rotating frame), µ E and µ M are the standard gravitational parameters of the Earth and Moon, ì ω is the angular rotation vector of the Earth-Moon system.

The rich dynamics of this system has been studied extensively and several periodic orbits of the CR3BP exist [START_REF] Gomez | Dynamics and mission design near libration points, advanced methods for collinear points[END_REF]. A class of well-known orbits is called halo orbits. Earth-Sun and Earth-Moon halo orbits have been successfully used for spaceflight missions such as LISA Pathfinder and Chang'e 4. The LOP-G will be placed on a subcategory of halo orbits called Near Rectilinear Halo Orbit (NRHO).

III. Mission architecture enumeration

A mission architecture is a specific sequence of transfers and maneuvers between different orbits and spacecrafts in order to achieve a particular objective. A large set of mission architectures exist which could be used to bring a crew from the LOP-G to the surface and back. In this chapter, a systematic way of considering all mission architectures is described.

While more sophistical optimization techniques could be used, such as branch and bound [START_REF] Morrison | Branch-and-bound algorithms: A survey of recent advances in searching, branching, and pruning[END_REF], combinatorial optimization was found to be sufficient for this application. Graphically, the architecture can be represented by graphs similar to Fig. 3.

In particular Fig. 3 represents a one stage architecture; the vehicle starts at the NRHO, performs a transfer to LLO, then to the lunar surface, then launch from the surface to reach LLO again and, lastly, returns to the NRHO. This diagram can also be shortened by the one on the right with the label "1" indicating the vehicle 1 is performing the transfers. To ensure an exhaustive search of the mission architectures, a general description of the framework of orbital transfers has to be constructed. The following 3 basic parameters will be used to parametrize the architectures: the number of valid locations, N S , the number of vehicles, N V and the number of allowed transfers, N F . In this work N S = 3, that is, only 3 locations are considered: the orbit of the LOP-G (NRHO), the low lunar orbits (LLO) and the lunar surface. While an arbitrary amount of vehicles and transfers could be considered, the search complexity grows exponentially with those variables: the number of possible architectures to check is bounded above by N N V (N F +1)

S

. This research will focus on the most promising cases; the number of vehicles N V will be varied from 1 to 3, (3 being the number Apollo used: service module, ascent and descent stages [START_REF] Orloff | Apollo by the Numbers[END_REF]). The number of allowed transfers N F will be varied from 4, the minimum to achieve the mission, represented in Fig. 3, to 6 transfers. Beyond 6, the search space is too large to compute with the available tools and techniques. Moreover, as orbital transfers are very expensive in terms of fuel, the optimal mission architectures will necessarily have a low number of transfers. Therefore, there is no drawback to searching only for a low number of transfers.

When N vehicles travel on the same trajectory there are exactly:

N n=1 N n n N -n (5)
number of different ways of arranging the transfers. The first few values are shown in Table 1. The explicit construction of the 10 cases for N = 3 is shown in the appendix 14. Because of the sheer number of architectures to consider, a simple and fast mass model to order the architectures had to be considered. The mass model, at this stage, is given by several linear relations between the variables [START_REF] Fortescue | Spacecraft Systems Engineering[END_REF]:

M T ot al = M P + M S + M T + M F (6a) M F = M T ot al 1 -e -∆v/(g 0 I s p ) (6b) M T = αM F (6c) M S = βM T ot al (6d)
Where M P is the payload mass of the spacecraft, M F is the fuel mass of the spacecraft, M T is the tank mass of the spacecraft, M S is the structure mass of the spacecraft, M T ot al is the total mass of the spacecraft, α is the tank to fuel mass ratio, β is the structural efficiency (structure/total mass) and ∆v is the velocity increment of the stage. The values of α = 0.05 and β = 0.15 were chosen according to the work of [START_REF] Fortescue | Spacecraft Systems Engineering[END_REF].

Due to the simplicity of the system, it is possible to write the solution of these equations in closed form as:

M T ot al = 1 (1 + α)e -∆v/(g 0 I s p ) -α -β M P (7)
Once the mission architectures are generated, a metric must be introduced to order them. At this stage, the metric chosen is not critical since this analysis is only to remove extremely inefficient architectures. Since the vehicles are meant to be reusable, the simple metric used is the amount of fuel per mission. This metric allows one to obtain a preliminary ranking of the various mission architectures. Examples of mission architectures are graphically represented in Fig. 4a to 4d. The number labeling refers to which spacecraft is performing the maneuver. Double and triple lines indicates vehicles moving together, with the first number being the vehicle that provides the thrust. For this analysis, the sole goal is to bring spacecraft "1" to the lunar surface and back at the NRHO. After taking into account combinatorial duplicates, there remain several mission architectures that, while theoretically valid, require strictly more fuel than another mission architecture. Example 4c shows a loop of vehicle "2" from LLO to NRHO, this loop is a strict waste considering the sole purpose of bringing spacecraft "1" on the surface and back. Example 4d has vehicle "3" coming back from the surface to LLO, this maneuver is strictly unnecessary in the context of bringing spacecraft "1" to the surface and back.

These architectures could have niche applications but, for the purpose of this study, any architecture which is strictly inferior in terms of mass than another one is discarded. Table 12 and 13 of the appendix show all mission architectures considered in this work. Table 13 and 14 enumerate the characteristics of each of the mission architectures.

IV. MDO model description

In this section, the physical models of the spacecraft will be introduced. Discipline models are needed to build the MDO tool. In particular, three disciplines are considered: propulsion, trajectory and structure as shown in an extended design structure matrix in Fig. 5. The MDO process also needs an objective function to optimize the vehicles; this will be introduced in the next section.

A. Propulsion

Since the vehicle is intended to carry humans and be refuelable in orbit, only liquid-fueled propulsion systems will be considered in this work. For these types of applications, 4 different types of propellant are usually considered [START_REF] Mills | Propellant Selection for the Lunar Lander Ascent Stage[END_REF] • LH2/LOX: Liquid hydrogen + Liquid Oxygen The propulsion module is modeled in three steps, first, the combustion analysis is performed using Chemical Equilibrium with Applications (CEA), the NASA developed state of the art combustion analysis software [START_REF] Gordon | Computer program for calculation of complex chemical equilibrium compositions and applications[END_REF]. Then the fluid flow is modeled by an isentropic one-dimensional flow through a convergent-divergent nozzle [START_REF] Sutton | Rocket Propulsion Elements[END_REF]. Lastly the mass of the engine is obtained via a multitude of Mass Estimating Relations (MER) [START_REF] Humble | Space propulsion analysis and design[END_REF]. While the software has been developed for all 4 propellants above, only the results for LH2/LOX will be shown. This allows a more direct comparison with the recent LH2/LOX Lockheed Martin lander. Figures 6 shows the Isp results of the propulsion module for a combustion pressure of 30 bar. The optimality line of best mixture ratio for a given expansion ratio is also shown. Humble [START_REF] Humble | Space propulsion analysis and design[END_REF], the tank mass can be based on a simple hoop analysis of thin pressure vessels with corrective factors.

B. Structural, thermal and mass analysis

Those assumptions allow the mass of the tank M t to be proportional to the mass of fuel (or oxidizer) inside the tank M F :

M t = N real 2 ρ t ρ F P σ M F (8)
with ρ t the density of tank material, ρ F the density of the fuel (or oxidizer), the operating pressure of the tank P and the allowable stress of the material σ. The allowable stress of the pressure vessel is usually half (safety factor of 2) of the ultimate strength of the material σ = 1 2 σ max . The actual mass of the tank is typically 2 to 2.5 times larger than the theoretical value, this is encoded in N r eal . To be conservative, the value of N real = 2.5 was chosen in this work. The choice of tank material is in itself a large subject. For simplicity, aluminum 2219 was chosen because of its history in the space sector [START_REF] Humble | Space propulsion analysis and design[END_REF]. Its relevant properties can be found in Table 5.

Table 5 Physical properties of aluminum 2219

Density (kg/m 3 ) Ultimate strength (GPa) 2800 0.413

As for the mass of the structure, landers and orbiters have different stress requirements [32] [9]. Order of magnitude estimate of the structural mass for orbiters and landers are: 15% of the dry weight for orbiters and 35% of the dry weight for landers.

For vehicles involving cryogenic propellant, like hydrogen, oxygen and methane, heat management is primordial to prevent excessive boiloff. Thermal control essentially falls into two categories, passive and active cooling. In passive cooling, the tanks are heavily insulated and the liquids are routinely vented to prevent pressure build-up inside the vessel.

While the boil-off rate (BOR) depends on the nature of the system, the quality of its insulation, geometry and exposure to the sun etc, some order of magnitude are available in the literature [START_REF] Rapp | Human Missions to Mars: Enabling Technologies for Exploring the Red Planet[END_REF], hydrogen BOR: 0.127% /day = 3.81% /month, oxygen BOR: 0.016% /day = 0.48% /month and methane BOR: 0.016% /day = 0.48% /month. When passive cooling is not sufficient for the mission, active cooling may be employed to keep the fuel and oxidizer from boiling away. The situation where all boil-off is prevented is called the zero boiloff (ZBO) condition. To achieve this state, a cryocooler must be used to remove heat leaks into the tank. At temperatures of liquid hydrogen, the efficiency of a cryocooler is in the order of 0.1%, requiring 1 kW of electricity to remove 1 W of heat. Oxygen and methane however are stored at a much higher temperature and so cryocoolers in that range have a 2% efficiency [START_REF] Radebaugh | Cryocoolers: The state of the art and recent developments[END_REF]. For a cylindrical vessel of constant radius, the heat leak scales proportionally to the mass of the fuel. Since the mass of the power system scales approximately linearly to the power required, this ultimately implies an approximately linear relationship between the mass of the cooling system and the mass of the propellant. The following relation is obtained:

M cooling = 2k I r ρ F M f uel = α cooling M f uel (9)
where k = 25 kg/kW is the amount of mass necessary to generate a given power with solar arrays [START_REF] Jpl | Solar Cell and Array Technology for Future Space Science Missions[END_REF], I = 0.4 W/m 2 is the heat leakage rate of an insulated tank [START_REF] Davis | Future orbital transfer vehicle technology study[END_REF], r ∼ 1 m is the radius of the cylindrical tank, ρ F is the density of the fuel and is the efficiency of the cryocooler. The proportionality constant can be estimated, α cooling = 0.3 for hydrogen and α cooling = 0.01 for methane/oxygen. In this work, only active cooling was considered. Further detailed comparison of thermal management systems should reconsider passive cooling as well. [START_REF] Whitley | Options for staging orbits in cislunar space[END_REF]. To consider the gravity losses due to design variation in the TWR and Isp, the trajectory analysis must be part of the optimization process.

C. Trajectory analysis

NRHO to LLO

The first phase of the transfer will move the spacecraft from the NRHO to the LLO. As a zeroth-order approximation, the transfer ∆v was obtained in this work via optimization of the Lambert problem [START_REF] Izzo | Revisiting Lambert's Problem[END_REF]. For the CR3BP, this value is approximately 750 m/s and agrees well with the literature [START_REF] Whitley | Options for staging orbits in cislunar space[END_REF]. However, the exact ∆v requirement for a particular mission depends on various components of the spacecraft, the most important factors being the thrust to weight ratio and the specific impulse. The gravity losses (or gravity drag) of a maneuver is the difference between the actually ∆v spent and the idealized ∆v inst with an instantaneous burn.

As a novel contribution, an analytical expansion of the gravity losses as a function of the maneuver time t f was found.

The spacecraft starts on a circular orbit and applies a prograde burn with the objective to change the apoapsis or periapsis.

This expansion is done in MATLAB symbolic toolbox and is truncated to order t 4 f :

∆v -∆v inst ∼ 1 48 µ r 3 A 0 t 3 f - 1 32 µ r 5 A 2 0 t 4 f + 1 96 µ r 3 A 2 0 t 4 f g 0 I sp ( 10
)
where ∆v is the actual cost to perform the maneuver, ∆v inst is the idealized cost assuming instantaneous maneuvers, A 0 is the initial acceleration of craft and r is the distance from the body at which the maneuver is performed. This formula is similar to Robbin's gravity loss formula [START_REF] Robbins | An analytical study of the impulsive approximation[END_REF]:

∆v -∆v inst ∼ 1 24 µ r 3 A 0 t 3 f (11)
Equation 10 used in this article contains higher order terms than Robbin's formula, hence it is more precise. Moreover, the boundary conditions assumed in deriving Robbin's formula are that the starting and ending orbits are fixed. Those conditions are stricter than the ones in this article. Thus, in this case, Robbin's formula overestimates the losses by a factor of 2.

LLO to Surface

The second phase of the transfer involves landing or taking off from the lunar surface. Considering first the takeoff, it is clear that the initial thrust must overcome the weight F > m µ r 2 . There are several ways to report the "cost to orbit" of a particular spacecraft, two of which are the fuel fraction ζ and the ∆v:

ζ = M i -M f M i = F m 0 t f g 0 I S P (12)
∆v = g 0 I S P ln m 0 m 0 -F g 0 I S P t f (13)
When ascending from the lunar surface, the ground must be cleared fast enough to prevent collision with the ground. In this work, the following simple geometric constraint was employed:

r(θ) ≥ R + θ θ h min + 1 k R ( 14
)
where R is the radius of the moon, h min is the altitude clearance and k is the initial slope of the constraint. Those constraints are shown geometrically in Fig. 7.

Fig. 7 Geometry of the ascent trajectory (Not to scale)

The optimal control of the trajectory described by equations 3a to 3e with the control variable α is solved using DYMOS with a Gauss-Lobatto transcription of order 3. Fixing the parameters k = 10 and h min = 6 km, the resulting fuel fraction can be calculated. A sample plot of the state and control variables is shown in 8. The solutions are calculated for a 20x20 uniform grid with 1 < TW R < 4 and 250 s< I sp < 500 s and is shown in Fig. 9. To interpolate the grid, the Python library Surrogate Modelling Toolbox (SMT) is used [START_REF] Bouhlel | A Python surrogate modeling framework with derivatives[END_REF].

V. MDO optimization specification A. Objective function

Several optimization metrics for aerospace vehicles have been used throughout the literature, including maximization of the payload at fixed total mass [START_REF] Balesdent | Multidisciplinary Design Optimization of Launch Vehicles[END_REF], minimization of cost [START_REF] Koelle | The transcost-model for launch vehicle cost estimation and its application to future systems analysis[END_REF], minimization of the dry weight [START_REF] Balesdent | Multidisciplinary Design Optimization of Launch Vehicles[END_REF] and minimization of the total mass [START_REF] Balesdent | Multidisciplinary Design Optimization of Launch Vehicles[END_REF]. While cost models are sometimes used as the metric, it is common, in the preliminary phase of the design, to use the mass of the vehicle as the objective function as a substitute for cost assessment. In this work, the objective function to minimize is the average mass M av necessary to accomplish the mission. The choice of this objective function assumes that one of the most important parameter to minimize for all lunar missions is the mass that must be sent to cis-lunar space, regardless of the nature of the mass itself. The objective function is thus given by:

F Objective = M av = M total + N reuse M expendable 1 + N reuse (15)
Where N reuse is the number of times the system can be reused, M total is the total mass of the system and M expendable is the mass necessary to refurbish the vehicles every mission, such as the fuel mass or expendable stages.

B. Optimization

For the optimization of the lunar vehicles, the following assumptions are made: all vehicles use hydrogen and oxygen as propellants and the number of reuses vary from 0 to 6. Each vehicle has its own thrust, expansion ratio and mixture ratio, which will be optimized, for a total of 3 optimization variables per vehicle. For 2 and 3 vehicles, all 20 mission scenarios are considered. The bounds on the design variables are shown in Table 8. In this model there is only one constraint: for vehicles that land on the lunar surface, the thrust has to be greater than the weight at touchdown and takeoff:

Thrust ≥ Weight Touchdown/Takeoff [START_REF] Gill | SNOPT: An SQP Algorithm for Large-Scale Constrained Optimization[END_REF] This constraint can be said to be "unrelaxable" [START_REF] Digabel | A Taxonomy of Constraints in Simulation-Based Optimization[END_REF], that is, the output of the model does not exist if this constraint is not respected. In this particular case, the ∆v would be infinite or undefined.

VI. Results and discussion

The best performance obtained is for a 2 vehicle architecture with 6 reuses, with an average mass of 46.0 t. In contrast, a one-stage expendable lander has a total mass of 80.2 t, an increase of 74% compared to the best one. The mass breakdown for selected mission architecture is shown in Table 9. The results of a one-stage lander can be found in Appendix Table 10 for a number of reuses from 0 to 6. For 2 and 3 vehicles and for the number of reuse varying from 0 to 6, the vehicles were optimized for each of the 20 mission architectures. Table 11 and 12 in the appendix shows all the possible cases with the ordering of the best 20 mission architectures with the percentages increase in average mass as compared to the best mission architecture (46 t). As a point of comparison, the Lockheed Martin Lunar Lander would be reused from 4 to 10 times [START_REF] Cichan | Concept for a Crewed Lunar Lander Operating from the Lunar Orbiting Platform-Gateway[END_REF] and has a wet mass of 62 t and a dry mass of 22 t. This can be compared to the results obtained in this article for a one-stage 6 times reusable lunar lander with a wet mass of 80.3 t and a dry mass of 25.3 t.

More generally, for one stage landers with multiple reuses, the average mass per mission is shown in Fig. 10a. As the number of reuse increases, the average mass necessary to complete the mission decreases. It is important to note the diminishing return behavior as the number of reuse increases. Reusing the system only once decreases the average cost by about 18%, however, the next reusability decreases the system mass by only 7%. Figures 10b-10c-10d show the values of the optimized design variables: thrust, expansion ratio and mixture ratio as a function of the reusability number. As the number of reusability increases, the optimization tends to favor more efficient vehicles even at the cost of higher dry weight. An increase in efficiency is primarily driven by two factors: lower gravitational losses and a gain in Isp. The three graphs are then explained as follows

• An increase of thrust diminishes the gravity losses at the cost of higher engine mass

• A higher expansion ratio increases the Isp at the cost of a higher engine mass

• A mixture ratio closer to the absolute optimal increases the Isp at the cost of an increase in tank mass Figure 11 shows architectures #4 #5 #7 and #9 for two vehicles and how the ranking of the architectures changes as the number of reuses increases. At zero reuse, the ranking is 4-5-7-9. However for three reuses, the ranking is 7-9-4-5.

The swapping around of the ordering of mission architectures is partially due to the fact that some architectures have expendable vehicles and some do not. In this particular case, architecture 4 and 5 both involve an expendable vehicle preferred whereas for higher reusability architectures 7-9 tends to wins in terms of efficiency. It is worth keeping in mind that the objective function only considered the average mass of the vehicles and did not include any penalty originating from the maintenance of the system or the refurbishment of the vehicles. For example, the mass/performances of the engines were modeled independently of the number of reuse. In reality, reusable space vehicles tend to be more complex and so would incur higher development costs than expendable ones. In addition, as systems are made to be reusable more and more times, the system requirements get stronger and harder to comply with, increasing the development and operating cost. Proper implementation of these considerations with a penalty term in the objective function would most likely result in an optimal design that occurs at a finite number of reuses. Moreover, no penalty was considered for the complexity of the mission architectures. Tables 13 and14 in the appendix show the main characteristics of the mission architectures for 2 and 3 vehicles. In future models, the objective function could have penalty terms related to the number of burns and the number of docking/undocking. Safety considerations should also be taken into account.

For example, surface rendezvous and docking may be extremely challenging and carry high risks, this should factor into the optimization. This work also did not include any cost associated with keeping material (like fuel) in a lunar orbit.

VII. Conclusion

In order to design a reusable vehicle connecting the LOP-G with the lunar surface, an innovative methodology has been developed and presented in this paper. The methodology comprehends a dual approach. Firstly, the The two methods are then merged to obtain the best preliminary design for a family of architectures. The MDO process involved the 3 disciplines of propulsion, trajectory and structural/thermal analysis. The optimization variables were thrust, mixture ratio and expansion ratio of the vehicles. The objective function was chosen to be the average mass expended per mission with a number of reuses of the vehicles varying from 0 to 6. The effects of the reusability number on the optimization variables were observed and their behavior interpreted. With the 20 most promising architectures for 2 and 3 vehicles, a large database of results was established, ranking the different architecture by their average mass usage as a function of the reusability number. Interestingly, the ordering of the best architectures depends on the number of reuses. As part of the trajectory analysis software, an analytical expression of the gravity losses was obtained. With these results, a system engineer designing a lunar lander can use this paper and his own mission requirements to choose the most appropriate mission scenario for his/her situation. Future work will improve engineering models and consider Ordering of the best architectures #1-#20, increase in relative mass %

5 0

 5 are then coupled with a multidisciplinary design optimization process by modeling the various components of the spacecraft and optimizing over a set of design parameters. The need for fast computational models, particularly in trajectory optimization, resulted in an analytical approximation of gravitational losses. This work resulted in a hierarchy of mission architectures that are ranked according to the average mass necessary to perform the mission. This work is intended to help a system engineer designing a lunar lander in choosing the best number of vehicles, the number of reuses and the mission profile for his/her mission requirements. years ago humanity stepped for the first time on the lunar surface. Nowadays a new space race is underway to

Fig. 1

 1 Fig. 1 Schematic of the workflow and tools used in this work.

Fig. 3

 3 Fig. 3 Schematic of the mission with a single vehicle.

Fig. 4

 4 Fig. 4 Mission architecture examples

Fig. 5

 5 Fig. 5 Extended design structure matrix for the lunar lander with the three disciplines of trajectory, propulsion and structure.

Fig. 6

 6 Fig. 6 LH2/LOX specific impulse (s) as a function of mixture and expansion ratio. Calculated using CEA [30].

Fig. 8

 8 Fig. 8 Sample optimal control for an ascent trajectory with path constraints. Initial TWR = 1.95, Isp = 309 s. Calculated using Dymos [19].

Fig. 9

 9 Fig. 9 Fuel mass fraction ζ for lunar take-off as a function of Isp and TWR. Calculated using Dymos [19], interpolated using SMT [41].

Fig. 10

 10 Fig. 10 Optimized objective function and design variables for a one stage vehicle

Fig. 11

 11 Fig. 11 Increase average mass for two vehicles architecture 4, 5, 7, 9. See Table 11

Fig. 13

 13 Fig. 12 Mission architectures for 2 vehicles

Table 1 Number of rearrangements

 1

	Number of vehicles traveling on the same path # rearrangements
	1	1
	2	3
	3	10
	4	41
	5	196

Table 2 Inputs/Outputs for propulsion module [28]

 2

	Inputs	Dimension	Outputs	Dimension
	Type of fuel	-	Specific impulse, I sp	s
	Mixture ratio	-	Combustion temperature	K
	Combustion pressure	Pa	Molecular weight	g/mol
	Expansion ratio	-	Adiabatic index	-
	Maximum thrust	N	Engine mass	t
	• UDMH/NTO: Hydrazine based compounds + Nitrogen tetroxide

Table 3 Inputs/Outputs for structural and thermal module

 3

	Inputs	Dimension	Outputs	Dimension
	Density of fuel	kg/m 3	Fuel tank mass	t
	Density of oxidizer	kg/m 3	Oxidizer tank mass	t
	Mixture ratio	-	Cooling system mass	t
	Tank material	See Table 5	Structural mass	t
	Total propellant mass	t

Table 4 Fixed mass for the two types of vehicles

 4

	Payload mass	Comments

Table 6 Inputs/Outputs for trajectory module

 6 One of the fundamental disciplines of any spacecraft is the analysis of its trajectory. In this work, the spacecraft must be able to operate between a NRHO, a low lunar orbit (LLO) and the surface of the Moon. Approximate values for the ∆v requirements are given in Table7

	Inputs	Dimension	Outputs	Dimension
	Thrust	N	Propellant mass	t
	Total vehicle mass	t	Maneuver ∆v	m/s
	Specific impulse, I sp	s	-	-

Table 7

 7 Approximate ∆v requirement

	Transfer	∆v (m/s)
	NRHO <-> LLO	750
	LLO -> Lunar surface	2100
	Lunar surface -> LLO	1900

Table 8 Bounds on the design variables

 8

	Variable	Min Max
	Thrust (kN)	10	400
	Expansion ratio	10	1000
	Mixture ratio	2	8

  On a Intel® Core™ i7-4800MQ CPU @ 2.70GHz × 8, it takes on average 8 seconds to optimize a one stage lander, 35 seconds for a two stage lander and 70 seconds for a three stage lander. It takes on average 2 hours to optimize all 41 mission architecture with 7 reuses (7 × 8 s + 7 × 20 × 35 s + 7 × 20 × 70 s ∼ 2 hours). To avoid local minima, the standard technique of multi-start is employed with 30 initial states randomly chosen in the design space. The total runtime was thus about 60 hours. The local optimization is performed by a SLSQP algorithm with a relative tolerance of 1e-6 as a termination criterion and maximum number of iterations of 100. If an initial guess does not satisfy the unrelaxable constraint 16, the optimization is started over with a newly randomized initial point.

Table 9 Selected architecture mission

 9

		One-stage	Most efficient architecture	NASA proposed architecture
		Fully expendable		6 reuses		3 reuses
	Average mass per mission (t)	80.2		46.0		58.6
	Vehicle number	1	1	2	1	2	3
	Total Mass (t)	80.2	50.8	16.1	36.0	26.1	13.3
	Mix Ratio	7.29	6.82	6.86	6.67	6.59	6.73
	Engine Thrust (kN)	216	189	44.0	139	246	52.5
	Expansion Ratio	128	115	904	112	68.3	523
	Dry Mass w/o payload (t)	13.0	10.3	1.71	8.64	3.43	1.41
	Dry Mass w/ payload (t)	25.0	22.3	2.11	20.64	3.83	1.81
	H2 Mass (t)	6.66	3.64	1.77	2.00	2.94	1.48
	O2 Mass (t)	48.6	24.9	12.2	13.4	19.4	9.99
	Fuel Tank + Coolant Mass (t)	2.95	1.62	0.786	0.890	1.30	0.658
	Oxygen Tank + Coolant Mass (t)	0.432	0.221	0.108	0.119	0.172	0.0889
	Engine Mass (t)	0.855	0.651	0.504	0.406	0.614	0.393
	Engine Isp (s)	463	465	491	466	458	486
	Structure Mass (t)	8.74	7.80	0.317	7.22	1.34	0.272

whereas architecture 7 and 9 have both vehicles reusable. Therefore for a low number of reuses, architectures 4-5 are

Table 12 Results for 3 stages architecture.

 12

	Reuse	0	1	2	3	4	5	6

Acknowledgments

The help of the whole SACLab team at ISAE-SUPAERO is greatly appreciated. Special thanks goes to Alberto Fossa and Giuliana Miceli for their contribution of the optimal ascent/descent trajectory.

Funding Sources

This work was possible due to the funding of Airbus DS, ArianeGroup and the SACLab at ISAE-SUPAERO

Appendix 25