Rémy Indebetouw

T Wong
email: wongt@illinois.edu

C.-H Rosie Chen

Amanda Kepley
email: akepley@nrao.edu

Vianney Lebouteiller
email: ney.lebouteiller@cea.fr

Suzanne Madden
email: suzanne.madden@cea.fr

Joana M Oliveira
email: j.oliveira@keele.ac.uk

STRUCTURAL AND DYNAMICAL ANALYSIS OF 0.1PC CORES AND FILAMENTS IN THE 30 DORADUS-10 GIANT MOLECULAR CLOUD

High-resolution (<0.1 pc) ALMA observations of the 30Dor-10 molecular cloud 15 pc north of R136 are presented. The 12 CO 2-1 emission morphology contains clumps near the locations of known midinfrared massive protostars, as well as a series of parsec-long filaments oriented almost directly towards R136. There is elevated kinetic energy (linewidths at a given size scale) in 30Dor-10 compared to other LMC and Galactic star formation regions, consistent with large scale energy injection to the region. Analysis of the cloud substructures is performed by segmenting emission into disjoint approximately round "cores" using clumpfind, by considering the hierarchical structures defined by isointensity contours using dendrograms, and by segmenting into disjoint long thin "filaments" using Filfinder. Identified filaments have widths ∼0.1 pc. The inferred balance between gravity and kinematic motions depends on the segmentation method: Entire objects identified with clumpfind are consistent with free-fall collapse or virial equilibrium with moderate external pressure, whereas many dendrogramidentified parts of hierarchical structures have higher mass surface densities Σ LT E than if gravitational and kinetic energies were in balance. Filaments have line masses that vary widely compared to the critical line mass calculated assuming thermal and nonthermal support. Velocity gradients in the region do not show any strong evidence for accretion of mass along filaments. The upper end of the "core" mass distribution is consistent with a power-law with the same slope as the stellar initial mass function.

INTRODUCTION

One major open question in star formation is how gas is accreted from the clump and cloud scale (∼10 pc) onto nascent stars. For many years there has been debate about whether isolated cores (∼0.1 pc) form which then collapse into individual stars and multiples, or whether a significant fraction of the eventual stellar mass is accreted from further away in the cluster gravitational potential well (often referred to as competitive or collaborate accretion). A related question is whether the stellar initial mass function is already set by cloud fragmentation into the core mass function, before most of the matter has accreted onto protostars (e.g. Motte et al. 2018, and refs therein). Over the last decade, it has become clear that in solar neighborhood clouds, filamentary structures (∼0.1 pc×1 pc) are ubiquitous, contain most of the dense cores in such clouds, and likely play a role in transferring matter from the cloud to core scales (e.g. André 2015;Fernández-López et al. 2014).

To develop any universal understanding of star formation it is important to determine whether the same core and filament structures exist in a wider range of molecular clouds than exist in the solar neighborhood. Distant regions in the Milky Way disk suffer from both significant distance (and thus source luminosity) uncertainty, and confusion along the line of sight. The Large Magellanic Cloud (LMC) is not much farther, but being close to face-on suffers from little distance uncertainty or line of sight confusion. The 30Dor-10 molecular cloud (Johansson et al. 1998) is located 15pc to the NE of the rich star cluster R136 in the LMC. The ultraviolet (UV) radiation field affecting 30Dor-10 is 3500× that of the solar neighborhood (Werner et al. 1978), and the reduced dust abundance at the region's ∼1/2 solar metallicity (Russell & Dopita 1992;Peimbert 2003) makes the interstellar medium more permeable to that radiation (Poglitsch et al. 1995). This makes 30Dor-10 a good laboratory in which to study star formation in the presence of strong external radiative feedback. 30 Doradus is also remarkable compared to typical star formation regions in the Magellanic clouds or Milky Way, and perhaps signatures of what caused such an intense star formation event are still evident in the structure of the residual molecular gas. Low and intermediate-mass star formation is actively ongoing in 30Dor-10 (Rubio et al. 1998;Walborn et al. 2013;Sabbi et al. 2016) and now with the Atacama Large (sub)Millimeter Array (ALMA) we can obtain a detailed picture of the molecular gas down to 0.1pc scales.

In ALMA Cycle 0 we mapped the 30Dor-10 cloud in 12 CO 2-1, 13 CO 2-1, C 18 O 2-1, and 1.3 mm continuum at a resolution of 2.3 ×1.5 (2011.0.00471.S, Indebetouw et al. 2013, hereafter Paper I), equaling = 0.56 pc×0.36 pc at the adopted distance of 50 kpc (e.g. de Grijs et al. 2014). We found that resolved parsec-scale structures have >3× larger velocity linewidths at those scales than other star formation regions in the Milky Way and LMC (Indebetouw et al. 2013;Nayak et al. 2016). This can be explained by an external pressure of P e /k>10 6 cm -3 K, either from the bubble and ionized regions around R136, or merely from the weight of the molecular cloud envelope in which the observed clumps are embedded. The slope of the size-linewidth relation agreed with other regions within uncertainties, and there were no strong trends of clump property with distance from R136, so the question of whether star-forming clumps in 30Dor-10 are affected by feedback, or that the entire region is simply very turbulent, was unresolved. This study extends the previous work to much smaller angular resolution, now resolving the 0.1 pc-sized cores that are expected to be actively participating in individual or multiple star formation. The greater spatial dynamic range permits more robust determination of size-linewidth-mass relations, and a core mass function more directly relevant to star formation can be measured, as described in the following sections.

OBSERVATIONS AND DATA REDUCTION

The 30Dor-10 cloud was observed with ALMA in 12 CO 2-1, 13 CO 2-1, C 18 O 2-1, 1.3 mm continuum, the H30α recombination line, H 2 CO 3 0,3 -2 0,2 , 3 2,2 -2 2,1 and 3 2,1 -2 2,0 (218.22219, 218.47563, 218.76007GHz; formaldehyde data will be presented separately). Most of the cloud was mapped as part of project 2011.0.00471.S (Indebetouw et al. 2013) at a resolution of 2.3 ×1.5 = 0.56 pc×0.36 pc. The brightest parts of the cloud were observed in project 2013.1.00346, targeting the same lines at higher angular (<0.1 pc) and spectral resolution. (C 18 O 2-1 and the H 2 CO lines were observed at 122 kHz=168 m/s resolution, 12 CO 2-1 and 13 CO 2-1 at 61 kHz 80 m/s resolution.) The observation was executed 7 times between 2015/06/27 and 2015/09/24. The phase calibrator was J0635-7516 (0.46-0.53 Jy at 230.5 GHz during the time range of observations). J0635-7516 was also used as bandpass calibrator in all but one execution which used J0538-4405 (1.2 Jy at 230.5 GHz). The amplitude calibrator was J0519-454 (0.66-1.0 Jy at 230.5 GHz). Data were calibrated using the ALMA Calibration pipeline version Cycle3R1 included in Common Astronomy Software Applications (CASA; http://casa.nrao.edu) v. 4.3.1 (McMullin et al. 2007). Calibrated visibilities were subsequently continuum-subtracted in the uv domain and deconvolved using clean, tclean in CASA 4.5, and a prototype version of the auto-multithresh automasking routine now included in CASA (since v. 5.1.0 Kepley et al. 2020). Visibilities from both projects were included in the deconvolution. The 12 CO 2-1 and 13 CO 2-1 data were imaged at 0.4 km s -1 resolution with Briggs weighting, robust=0.5, and multi-scale deconvolution at 0,5, and 9 times the 0.032 arcsec pixel, achieving a beam of 0.31 ×0.22 = 0.08 pc×0.05 pc, and rms noise in line-free channels of 2.3mJy bm -1 . The 12 CO 2-1 image was feathered with APEX single dish data to recover all of the ∼25% of the large-scale emission that was resolved out by the interferometer; the combined image used for subsequent analysis has an rms of 6mJy bm -1 . Comparison of the ALMA and APEX data for 13 CO 2-1 indicates that the interferometer recovered all emission to within uncertainties, so the ALMA-only image is used for 13 CO 2-1 analysis. The C 18 O 2-1 data were imaged at 0.4 km s -1 resolution with natural weighting, achieving a beam of 0.44 ×0.28 and rms = 1.6 mJy bm -1

CO MORPHOLOGY AND ASSOCIATED STAR FORMATION

Figure 1 shows the extent of 12 CO 2-1 emission relative to optical and near-infrared emission imaged with HST for the HTTP project (Sabbi et al. 2013), and the location of features that will be discussed in this section. Figure 2 shows the peak and integrated emission of 12 CO2-1 and 13 CO2-1, and Figure 3 shows three zoomed-in subregions in F160W∼H, F110W∼J, and 1mm continuum.

Nearest R136 in the southwest, the CO-traced molecular gas is characterized by photodissociated pillars (P1,2,3 in Figure 1) each of which has embedded 1mm sources visible in Figure 3, third panel. Interestingly, the pillar "P1" closest to R136 in projection happens to be the largest sized pillar, and contains two very embedded YSOs, only detected in the millimeter continuum. Just in "front" (SW, i.e. on the R136-side) of the pillar are two near-infrared-detected, more exposed or evolved, stars. This sequence of more embedded sources being further back from the photodissociating source has been studied in the Milky Way, and cited as evidence for triggered star formation by compression of the pillar heads, but causality is challenging to establish (Fukuda et al. 2002;Dale et al. 2015). A 1720GHz OH maser was detected at the pillar position (Figure 3, third panel; Brogan et al. 2004), but those data had insufficient angular resolution to associate the maser with the embedded protostars, or the outer part of the pillar compressed by the HII region. 1720GHz OH masers are often associated with supernova remnant-molecular cloud interactions, and a similar shock could be being driven into the pillar, but such masers are also associated with massive YSOs (e.g. Gray et al. 1992). Pillar 1 is also associated with one of three "massive YSOs" identified in this region with Spitzer (Whitney et al. 2008;Gruendl, & Chu 2009;Walborn et al. 2013). However, even with IRAC at ≤8µm, Spitzer's resolution is nearly a half-parsec, so its not possible to distinguish a small cluster from an individual massive protostar. With HST's resolution of ∼15000AU, resolved red near-infrared sources at the positions of the red Spitzer sources can be much more convincingly called single or multiple MYSOs, and now with ALMA the small-scale structure of their associated molecular gas is resolved. The Pillar 1 "massive YSO" is clearly a small group including the more embedded millimeter continuum sources and the more exposed NIR sources in the bright photoionized rim.

Slightly further away from R136 in projection (∼3pc) are two other pillars "P2,3" which each contain a YSO visible in the NIR. Both embedded sources are slightly extended in NIR and 1mm continuum.

The main part of the region, between the south-west pillars and the central cluster and clumps, is characterized by multiple narrow filaments pointing approximately at the central stellar cluster R136. Photoionization can erode clouds leaving shadowed structures and pillars pointing towards the radiation source (Gritschneder et al. 2010), but it is also possible that this structure is instead related to the formation of 30 Doradus: Rahner et al. (2018) model the region as an older burst of star formation, followed by recollapse and formation of R136, at high molecular cloud mass and density. Although their model is only 1-dimensional, repeated feedback-driven recollapse and non-spherical re-expansion could naturally lead to dense radial structures.

In the center of the mapped region is a small (∼1.5pc diameter, see Figure 1) intermediate-mass cluster associated with diffuse Hα emission, and at least one very massive star (O3-6; Walborn, & Blades 1997). Remnant filamentary CO-traced molecular gas is coincident with the eastern side of the cluster, likely in front of the cluster, since the CO filaments correspond to shadows in the diffuse Hα. These filaments may be under compression by the expanding ionized region, but there is no evidence in millimeter continuum or NIR of embedded protostars. Just east of the cluster is the brightest CO clump in our data, now resolved into a southern core with an embedded near-infrared source, and a northern "hub" with filaments leading into it (although there are no clear kinematic signatures of accretion along those spokes). The northern source has an extended 1mm source, not detected in the NIR, and is likely the young MYSO that dominates the MIR emission from the region (Spitzer cannot resolve these two MYSOs). That more embedded northern source is associated with a H 2 O maser (Imai et al. 2013), thought to be collisionally pumped, and most commonly associated with massive YSO outflows (Elitzur et al. 1989;Walsh et al. 2011). The "hub-and-spoke" morphology of the CO emission is also commonly seen in Galactic MYSOs forming in molecular clouds. Thus, although the source is suggestively located on the edge of the small HII region, it is difficult to unambiguously claim that its formation was triggered by that HII region. On the western side of the small cluster is another CO clump coincident with two detected H 2 O masers. One is at the position of a 1mm continuum source, the other might coincide with a NIR-detected star.

On the eastern side of our mapped area, furthest from R136, is another very massive clump, with at least one deeply embedded massive protostar detected in 1mm continuum. The kinematic structure is complex, and suggestive of multiple outflows. It is not clear whether the near-IR source is an illuminated outflow cavity from the massive protostar(s) or an unassociated source.

4. CO 2-1 STRUCTURE ANALYSIS 4.1. Methods Molecular clouds have complex hierarchical structure on a large range of spatial scales, but there is insight to be gained by segmenting the emission and analyzing it as discrete structures. The most common technique has been segmentation of position-position-velocity cubes into approximately round entities a few times the spatial resolutionthese are typically called clumps or cores, for ∼parsec-sized, and ∼0.1 pc-sized entities, respectively. The segmentation process begins by finding local maxima, with two critical parameters: δI the intensity difference between a local maximum and highest connected saddle point, and I min the minimum intensity to consider. Neighboring pixels are F110W∼J,and F658N (Hα). The CO beam is the small black dot to the left of the "12CO2-1" label. The dominant sources of ionizing radiation are the stars to the southwest including R136. Pillar structures P1,2,3, and a small cluster in the center of the cloud, are labeled in blue, and the three Spitzer-identified massive YSOs in this region are marked in cyan -see text for discussion. Figure 2 shows the CO structure in more detail. Figure 3 shows zooms of the three regions in magenta.

assigned to each local maxima to define a set of clumps. Popular implementations are clumpfind (Williams et al. 1995(Williams et al. , 2011)), which assigns all emission down to I min , and cprops (Rosolowsky & Leroy 2006, 2011), which only assigns emission down to the lowest isointensity surface that contains a single local maximum.

An alternate analysis considers isointensity surfaces as a set of hierarchical entities, instead of assigning emission to disjoint clumps. This is more naturally suited to the hierarchical structure in molecular clouds, but a given emitting pixel gets plotted and analyzed multiply, as part of multiple isointensity structures. The smallest structures associated with local intensity maxima are called "leaves", and the largest isolated regions of emission called "trunks" or "islands". The most commonly used implementation is dendrograms (Rosolowsky et al. 2008). For clump segmentation, we use the quickclump python implementation (https://github.com/vojtech-sidorin/quickclump/) to which we added a I minpk parameter, the minimum peak intensity required for a valid clump (https://github.com/indebetouw/ quickclump). This allows the fainter envelopes of bright clumps to be included without keeping faint noise peaks. ), magenta contours: 13 CO 2-1 beak brightness at 0.04,0.06,0.08,.1 Jy bm -1 (cube rms = 2.2mJy bm -1 ). H2O masers are marked with yellow circles of radius=1 , the quoted positional uncertainty (Imai et al. 2013). An OH maser is marked similarly with a cyan circle (Brogan et al. 2004).

Properties (moments, fitted sizes, etc) of clumps are calculated with a python translation of the moments calculations in cpropstoo (https://github.com/akleroy/cpropstoo). For dendrograms we use the python implementation at https://github.com/dendrograms/astrodendro.

Somewhat more recently it has become popular to identify elongated regions of emission and refer to them as filaments. We use filfinder (Koch & Rosolowsky 2015) and its python implementation (https://github.com/ indebetouw/FilFinder) which allowed for easy modification, testing, and incorporation into other analysis scripts. Filaments are initially identified in the (2D) peak intensity image -if there is not an overabundance of sightlines with multiple velocity components, this should work well, and indeed for this data, the method produces a visually very satisfactory result (Figure 4). Manual examination and searching of this data cube only finds one position at which there are two bright structures at different velocities, and fewer than 10 positions with a bright structure at one velocity, and structure wings or faint emission at another velocity. filfinder performs pruning of small and multiply-connected filament branches. We modified filfinder to do that pruning in 3D as well as the default 2D. That requires choosing a metric to calculate length in position-position-velocity space, so we used 1km s -1 ∼0.036pc. Our trials found that the final filament skeletons differ only in minor ways between the default 3D pruning and pruning based on 3D lengths, so here we present only the default 2D results using the publicly available code.

It is important to note the bias of the various algorithms: clump segmentation will divide emission into approximately round entities a few times the beam size. Filament finders will identify filament skeletons for any distribution of emission, whether visually filamentary or not. Dendrograms do not impose such geometric constraints, but multiplycount most emission. To leverage the objective nature of dendrograms in filament analysis, we match filament branches to elongated dendrogram entities. We identify and place greater confidence in those filaments for which there is an isointensity surface that surrounds most of the filament, which contains most of only one filament, and is elongated (in the end a cutoff aspect ratio > 2.5:1 was used, but higher cutoffs up to 5:1 yielded similar results). Figure 4 shows the filfinder filaments identified in the peak intensity image, and the elongated isointensity surfaces that match some of the filament branches. Filaments are found using a global threshold of 0.17 and 0.07 times the peak intensity for 12 CO and 13 CO, respectively. The image is flattened at 90% of peak, smoothed to 3 times the beamwidth, and size and adaptive thresholds of 24 times the beam area and 14 times the beamwidth used (see filfinder documentation for parameter descriptions). Branches were pruned with thresholds of 5 beams and 10 pixels. Only isointensity structures more elongated than 2.5:1 with areas between 0.2 and 0.5 square parsecs were matched to filament branches. Overall, there is good but not perfect agreement between identified filaments and elongate isointensity contours, so a filament analysis will reveal different information relative to analyzing all dendrogram structures (defined by isointensity surfaces). In regions where filaments are relatively isolated, the same filaments are identified using 12 CO and 13 CO, although sometimes the 13 CO intensity is low enough to cause gaps or breaks in the 12 CO-identified filament. In dense clumps such as those in the center and east of the region, containing the most massive Spitzer-identified protostars (see § 3), 12 CO becomes very optically thick, and the relatively flat spatial intensity profile over the entire clump causes the filament finder to break the clump up into loops or cells (see especially the central clump with cyan loops in the 12 CO image, first panel of Figure 4.

Mass calculation

Physical analysis of structures identified in clouds requires calculating the mass of each structure. We calculate the mass per PPV pixel from 12 CO and 13 CO brightness using what is sometimes referred to as the LTE or "standard" method (Bourke et al. 1997;Indebetouw et al. 2013). The 12 CO excitation temperature is derived from the 12 CO 2-1 brightness temperature. The 13 CO excitation temperature is assumed to be the same as 12 CO 2-1, to calculate the 13 CO optical depth and column density.

To assess the accuracy of this method, we ran a grid of non-LTE excitation models with Radex (van der Tak et al. 2007), spanning N( 12 CO)∈[10 16 ,10 21 ]cm -2 , T K ∈[2,100]K, n(H 2 )∈[10 2 ,10 7 ]cm -3 , and N( 12 CO)=75 N( 13 CO) (Heikkilä -1 e 10.6/2.7 -1 -1 N ( 13 CO) = 1.5 × 10 14 T ex e 5.3/Tex τ 13 v dv 1 -e -10.6/Tex

We considered only models with N( 13 CO)<10 14 n(H 2 ), corresponding to line-of-sight pathlength <1pc for an abundance H 2 / 13 CO =5×10 5 . The results of this section are insensitive to path lengths and H 2 / 12 CO/ 13 CO abundances differ by up to a factor of 3 in either direction. Figure 5 compares N( 13 CO) calculated with the LTE method to the actual value input to each model. The method works fairly well, with a tendency to overestimate the true column density in brighter regions. There is a modest effect that the calculated column density (and mass) of bright structures may be overestimated by up to a factor of 2 relative to faint structures.

Structure analysis result: filament velocity structure

Part of understanding how stars accrete mass, and how that relates to molecular cloud structure, is determining the extent to which mass is accreted along filaments. To begin to quantify filamentary accretion in this region, we calculated the intensity-weighted mean velocity at each point in the map (the first moment of the intensity cube), and then calculated the 2-dimensional gradient of that velocity field. Figure 6(L) shows the 12 CO2-1 moment 1 velocity map, with filaments overlaid. One diagnostic is whether the velocity gradients are predominantly aligned with filaments, as would be the case if accretion along relatively long-lived filaments were dominant, or across the filaments, as would be the case if the filaments were the result of a turbulent velocity field, and not very dynamically important in the cloud. Figure 6(R) shows the distribution of angles between the local gradient in the moment 1 velocity, and the filament skeleton direction. No strong trend is evident in the alignment, although there is a weak trend for filaments with a velocity rms larger than 1 km s -1 to have a velocity gradient more across than along the filament. This result suggests that accretion along filaments, although it may be present, does not dominate over the stochastic turbulent velocity motions in the cloud.

Structure analysis result: filament stability

The stability of cylindrical shapes has been calculated for numerous cases including infinite homogeneous (Chandrasekhar & Fermi 1953;Ostriker 1964), polytropic (Ostriker 1964), and magnetized (Stodólkiewicz 1963;[START_REF] Tilley | [END_REF]. Unmagnetized isothermal filaments with line mass i.e. mass per unit length M l > 2σ 2 /G are unstable to gravitational collapse, where σ is the velocity dispersion and G the gravitational constant. Other geometries and equations of state change the critical value by a factor of order unity. We calculate the critical line mass at each point of each filament, for only thermal support using the CO excitation temperature i.e. σ 2 = kT ex /µ where µ is the mean molecular mass 2.36m H , and for thermal and turbulent support σ 2 = kT ex /µ + σ 2 v where σ v is the second velocity moment of the filament, calculated over the part of the cube within 5 km s -1 of the filament peak velocity.

We calculate the actual line mass three ways: for each point on the filament, filfinder extracts the profile perpendicular to that point. We fit a Gaussian of width σ r to the profile of 12 CO 2-1 intensity, and calculate M l = σ r √ 2πN pk using the LTE column density in the center of the filament at that point, N pk . The fitted widths and peak column densities are calculated from the average of a 3-pixel neighborhood (smaller than the angular resolution) to increase signal-to-noise. In the second method, we fit the width and amplitude of the perpendicular profile in N , and use M l = σ r √ 2πN amp . In the third method, we simply integrate N across each point of the filament. When structures are close together in the direction perpendicular to the filament axis, the width-fitting methods do not have as many points to fit that are unambiguously associated with the given filament, and thus can underestimate σ r and M l , but on the other hand, the integrated N will be an over estimate because it includes unassociated emission. Figure 7 shows the three line masses and two critical masses along one example filament. The more visually fragmented portion of the filament ∼1pc from the bottom (the "gap") has M l << M l,crit and that material is likely not gravitationally bound. By contrast, the bright knot ∼1.7pc from the bottom of the inset (from the left in the profile) appears to be unsupported against collapse, with M l > M crit,nonthermal . Support for other portions of the filament are more ambiguous, with M crit,thermal < M l < M crit,nonthermal .

The gravitational stability of 0.1pc scale structures in 30Dor-10 analyzed as filaments is shown in Figure 8. All points for all filaments are now plotted together, after analyzing all filaments as shown in Figure 7. Points with M l > M crit have more mass per unit length than thermal or kinetic support can prevent from collapsing (the nonthermal critical mass is shown). The points cluster at M l = M crit , but there are many with M l > M crit which might be unstable, or supported by magnetic fields. There is a systematic uncertainty in M l depending on the choice of abundance ratio n(H 2 )/n( 13 CO). 5×10 6 was used here; a lower ratio, or higher 13 CO abundance, by a factor of 5 would move most points into the stable part of the plot. Alternately, a factor of 3 lower ratio combined with a systematic overestimate of N LT E by a factor of two (Figure 5, §4.2) would also bring most points into the stable regime. Given these systematic uncertainties it may be more illuminating to examine the trends -clearly, locations with brighter 12 CO 2-1, or those (L) points colored by 12 CO 2-1 brightness temperature from lower (blue) to higher (red). "x" symbols are those filaments that correspond well to elongated isointensity contours (Figure 4). (R) points colored by 8µm brightness from lower (blue) to higher (red).

with brighter 8µm emission (from Spitzer/SAGE Meixner et al. 2006), are statistically less gravitationally stable. Both trends are present when all filament branches are considered (dots) as well as when only considering filaments that are associated with elongated isointensity contours (crosses, see Figure 4 and §4.1. It is not unexpected that less stable parts of the cloud might have more significant associated star formation, and hence brighter 8µm emission, but the Spitzer resolution of >0.5pc makes it difficult to unambiguously associate emission with 0.1pc scale molecular structures or to draw strong conclusions.

In Figure 9 we present the width distribution of filaments. The spatial moment of 12 CO intensity perpendicular to the filament was calculated at each position, as well as the spatial moment of calculated column density. The weighted mean of those moments is used as the best estimator of the spatial sigma σ r , and the deconvolved width is (2.354σ r ) 2 -0.075 2 . There is no significant difference in the width distribution between points on filaments that have line mass greater than critical line mass including nonthermal (turbulent) support. The median width is 0.09pc, with standard deviation of 0.06pc, so the filaments in 30 Dor-10 have consistent widths to the 0.1pc width found in solar neighborhood clouds (Arzoumanian et al. 2019). Distribution of widths of positions along filaments, where "width" = 2.534 times the fitted spatial moment of intensity, with the beam size of 0.075pc subtracted in quadrature. All points with fitted width > 2 times the width uncertainty are plotted, as well as the subset with M l > M crit,nonth . The dotted line indicates the beamsize.

Structure analysis result: size-linewidth-flux relations

The relations between size, velocity dispersion, and surface density (or luminosity as a proxy for mass) have long been studied to yield insight into the structure and gravitational stability of molecular clouds (e.g. Heyer et al. 2009;Goodman et al. 2009;Ballesteros-Paredes et al. 2011;Bolatto et al. 2013). Figure 10 shows the relation for this cloud, for structures segmented with dendrograms and clumpfind, and compared to other molecular clouds. The "radius" plotted is the historical value of 1.91 times the second moment in the plane of the sky, and both that size and the velocity dispersion have had the instrumental resolution deconvolved. The differences between algorithms are clearthe dendrogram analysis classifies all emission, down to very small structures, and multiply-counts that emission, as part of very large structures. By contrast, clumpfind is highly biased towards selecting clumps in a narrow range of sizes. The relation is fit to a power-law with the Levenberg-Marquardt least-squares fitter in the Kapteyn package [START_REF] Terlouw | Kapteyn Package version 2.3[END_REF]. The slope of 0.75±0.05 is somewhat steeper than the theoretical slope of 0.5 for a medium dominated by turbulent shocks (Ballesteros-Paredes et al. 2011, and refs therein), and on the higher side but within 3σ of values measured in other clouds (see references and discussion of 30Dor-10 at lower angular resolution in Nayak et al. 2016).

N159, a pair of massive star forming clouds just to the south of 30 Doradus, has been observed with ALMA at similar angular resolution to the data presented here (Saigo et al. (2017) 2019): resolution of 0.28 ×0.25 ) dendrogram structures in N159 have a size-linewidth slope of 0.7±0.1, consistent with 30Dor-10, but the linewidth at a given scale is 0.33±0.05 dex lower. The molecular mass of N159 is significantly higher than 30Dor-10, so the enhanced linewidths in 30Dor cannot be explained by global equilibrium between gravitational and kinetic energy. The enhanced kinetic energy is either a result of feedback, or an equilibrium including external pressure on the molecular cloud. There is embedded star formation within 30Dor-10, but the total infrared luminosity and inferred star formation rate is again significantly less than in N159. Thus, internal sources of mechanical energy are unlikely the source of elevated turbulence, and instead this is either caused by thermal pressure from the ionized gas, or large-scale dynamics of the region. Interestingly, Lee et al. (2019) analyze high-J CO and far-infrared line emission at low spatial resolution in 30 Doradus, and find that they require significant energy input from low-velocity shocks to explain the CO line spectral energy distribution. They also conclude based on the distribution of main sequence and protostellar wind sources that local kinetic energy injection is unlikely to dominate, but instead they favor kpc-scale energy injection related to the overall dynamics of the region. This kpc-scale energy input is very likely related to kpc-scale colliding filaments which intersect at exactly the center of 30 Doradus, and may well be the reason that the super-star cluster could form there in the first place (Fukui et al. 2017).

Also shown in Figure 10 are the relation fit to Galactic molecular clouds δv=0.72 R 0.5 (Heyer et al. 2009;Solomon et al. 1987, SRBY), and structures in the Perseus A cloud in the solar neighborhood (Ridge et al. 2006). The fitted slope of 0.5±0.05 and intercept are consistent with the SRBY relation, with dispersion at a given scale 0.65±0.05 lower than in 30Dor-10.

Figure 11 shows the relation between LTE and virial surface density for structures in 30Dor-10. Σ=M/πR 2 where 1987;Heyer et al. 2009). Data from all regions were re-analyzed the same way, with clumps segmented using dendrograms.

R is the cloud "radius" 1.91 √ σ x σ y calculated from the weighted spatial moments σ x and σ y . One can also use the "exact area" or full spatial extent of the pixels assigned to each structure -this area is 1.2-1.9 times larger, decreasing both surface densities accordingly. Σ LT E is calculated from the LTE mass, and Σ vir from the virial mass 5σ 2 v R/G. This plot is sometimes referred to as a "boundedness" plot because it probes the degree to which kinetic energy is gravitationally bounded. Structures in virial equilibrium between gravitational and kinetic energy would lie along Σ vir =Σ LT E , and those in free-fall (hierarchical) collapse driven by gravity along Σ vir =2Σ LT E (Ballesteros-Paredes et al. 2011). Alternately clumps might be in virial equilibrium between kinetic and gravitational energy with external pressure (Field et al. 2011), if they are sufficiently long-lived to achieve that virialization (Bonnell et al. 2006). The resulting curves Σ vir = Σ LT E + 20 3π20.9 nT Σ LT E are marked on Figure 11 for external pressures of nT = (10 5 , 10 6 ) cm -3 K.

Structures segmented using clumpfind fall in a similar locus to much larger (1-10pc) molecular clouds in the Milky Way disk and Galactic central molecular zone (Heyer et al. 2009;Oka et al. 2001). These all fall somewhat above the unity relation expected from virial equilibrium, and more consistent with either external pressure confinement or free-fall collapse. Structures in 30Dor-10 lie in between the disk and central molecular zone (CMZ) clouds; if the relevant model is confinement by external pressure P e , then that pressure is higher in 30Dor-10 than typical Milky Way, but lower than the CMZ. For the large Milky Way clouds, P e should be interpreted as pressure from the neutral interstellar medium surrounding the molecular cloud, but for the sub-parsec structures plotted in 30Dor-10, P e is the pressure external to a core, which is a combination of any warm ISM pressure acting on the entire cloud, and the effective pressure of the diffuse cloud acting on the core. In Paper I we used the formula of Bertoldi & McKee (1992) to estimate that interclump pressure to be 3×10 6 cm -3 K, consistent with the location of core and clumps in Figure 11.

Structures identified by dendrogram have significantly higher Σ LT E than those identified with clumpfind. This conclusion is unchanged if either the "exact" cloud area is used to calculate Σ, or if observed quantities (without the velocity and angular resolution deconvolved) are plotted. Dendrogram structure properties are calculated in Figure 10 by assuming all emission above each isointensity surface is part of the structure (bijection), but one could also subtract the value of that lower bounding surface (clipping), as discussed in Rosolowsky & Leroy (2006). Clipping decreases Σ vir by only 0.02dex, whilst decreasing Σ LT E by 0.3dex -the mean difference of Σ LT E between dendrogram and clumpfind structures is 0.8dex without clipping, and decreases to 0.5dex when clipping is used. Rosolowsky & Leroy (2006) conclude that clipping underestimates clump masses for their solar neighborhood clouds, but for the data presented here, clipping would result in structures with Σ LT E closer to those identified by clumpfind.

There is some tendency for larger dendrogram structures to be evidently less gravitationally stable: The second panel of Figure 11 shows Σ vir and Σ LT E colored by structure size. The virial mass was calculated for a constant density sphere, but molecular cloud substructures are clearly not spheres. Generalizing from spheres, the gravitational energy of an ellipsoid can be written in closed form for various radial density profiles (Neutsch 1979), including constant density,

U = 3 5 GM 2 l w(η) w(η) = η sinh -1 η 2 -1 η 2 -1 ,
where M is the mass, l the semi-major axis, and η = l/R the aspect ratio (Lee et al. 2017). For the same mass, a more elongated cloud has greater gravitational energy by a factor of (2,3) for a cloud aspect ratio of ∼(3.5,10). Thus the virial mass required for gravitational energy to balance a given kinetic energy is lower by a factor of a few for an elongated cloud compared to the assumed spherical shape. If we applied this correction, we would conclude that the clouds are even less gravitationally stable, and located even lower down in Figure 11. We measured the elongation or aspect ratio of each structure in the Figure, but there is no clear trend with boundedness, and applying an elongation correction does not decrease the scatter of points in the plot. Another source of systematic uncertainty was mentioned in §4.4: A factor of 5 higher assumed 13 CO abundance would move most dendrogram points into the stable part of the plot. Alternately, a factor of 3 lower ratio combined with a systematic overestimate of N LT E by a factor of two (Figure 5) would also bring most points into the stable regime. However, either of those corrections would cause the clumpfind clumps to disagree significantly with previous measurements of molecular clouds in this and many other regions.

The most likely cause for the difference between the boundedness of clumpfind and dendrogram segmented structures is that the dendrogram structures by design do not include entire clumps. Consider a model core with power-law density distribution outside of a core radius r c . The density distribution ρ ∝ r -1 for example was used in Solomon et al. (1987) to derive the commonly used relation between the second spatial moment σ x and the "edge" or "effective radius" of the core R e = 1.91 σ x . If one assumes that the velocity dispersion σ v follows a power-law size-linewidth relation σ v ∝ R α , where R is the projected radius in the plane of the sky, then one can calculate the virial surface density Relation between virial surface density 5σ 2 v /(GπR) and LTE surface density MLT E /(πR 2 ) for a model spherical core. The core has constant density within rinner, and a power-law radial density distribution outside of that. Each part of the core within a projected radius R is assumed to follow the size-linewidth relation σv ∝ R α . Calculation of the two surface densities for a partial core will yield significantly higher Σ and moderately lower Σvir compared to the values calculated for the entire core.

Σ vir (R) = 5σ 2 v πGM (R)R measured
as a function of projected radius R, where M (R) = R 0 4rρ(r)dr. The measured radius R measured is what one could measure from the data within a radius R, and could be simply R itself (i.e. the radius of the assignment area area/π), or more commonly the "Solomon effective radius" calculated from the two spatial second moments of the emission, 1.91 √ σ x σ y . Figure 12 shows the virial surface density Σ vir and the actual measured surface density Σ = M/π/R 2 measured , as a function of R for a model core with two different size-linewidth relations α and two different power-law density profiles. Clearly, if one calculates a surface density for only the brightest part of a core, one will measure a significantly higher surface density Σ than for the entire core. At the same time, the calculated virial surface density Σ vir will be somewhat lower than the whole-core value. This underscores the need to use the same segmentation method when comparing different datasets. It also suggests that if parts of clouds segmented with dendrograms are analyzed, the relative boundedness between two different clouds is a robust comparison, but the absolute value of that boundedness relative to the theoretical lines of stability should be interpreted cautiously. As a final test, we calculated the sizelinewidth-mass relations for structures identified with clumpfind but using a higher and higher noise floor or cutoff, raising it gradually up to 20× the noise level in the cube. Raising the cutoff floor causes the assigned structures' properties to smoothly move over to the location of the dendrogram-assigned structures, as expected. 4.6. Structure analysis result: Core mass function A fundamental question in star formation is whether the stellar initial mass function (IMF) is predetermined by the molecular cloud structure. Similarity between the dense core mass function and the IMF would support that premise. We analyze compact structures identified using clumpfind (as implemented in cprops). The distinction between calling a molecular cloud structure a clump or core is inhomogeneous in the literature, with ∼1pc structures usually called clumps and ∼0.1pc structures usually called cores; These structures are 0.1-0.2pc in diameter, so we call them cores for brevity, without intending any implication about stability or concentration.

Figure 13 shows the core luminosity and mass distributions. Different values of the clumpfind segmentation parameters are shown -dT and T min , the difference between each local maximum and nearest saddle, and minimum level to analyze. The main effect of varying those parameters is to increase the number of faint cores when T min is decreased. Fitting power law distributions has well-studied uncertainties (e.g. Clauset, Shalizi, & Newman 2009;Maíz Apellániz 2009); we show both the Hill maximum likelihood estimator with its statistical uncertainty (dashed lines), as well as a simple linear fit to the log number n, weighted by n -0.3 (dotted lines; the exponent of the weight makes little difference to the result). It is evident that even with the mathematical uncertainties due to fit method, the uncertainties due to how the emission is segmented into structures are even larger.

12 CO luminosity is proportional to cloud mass on large (>pc) scales (e.g. Bolatto et al. 2013). Although one would expect the luminosity of an optically thick line to be a less reliable mass tracer on small scales when emission fills the beam, it is still interesting to consider the shape of the 12 CO luminosity function (Figure 13 left). The bright end is fit with a power law of slope α=1.7±0.2. If luminosity were proportional to mass, this would imply a core differential mass function N (> M ) ∝ M -α slope of somewhat steeper than the stellar initial mass function (IMF) slope -1.35. However, the 12 CO optical depth increases systematically with mass τ ∼M 0.3 LT E , so the observed 12 CO luminosity function is expected to be steeper than the mass function by a slope of about 0.3, and the inferred mass function from 12 CO alone would then be consistent with the IMF, within uncertainties. The second plot of the figure shows the differential LTE mass distribution (from 12 CO and 13 CO), for different core segmentation parameters. Slopes are fitted directly to the differential mass distribution as well as using the aforementioned Hill maximum likelihood estimator. The systematic biases of direct fitting are more evident than when fitting the cumulative distribution. The best estimate slope of -1.3±0.15 is consistent with the stellar IMF slope -1.35.

The upper end of the "core" mass distribution is consistent with a power-law with the same slope as the stellar initial mass function. A shallower mass function (more massive stars) has been found for the main-sequence massive stars in 30Doradus (Schneider et al. 2018). The next generation of star formation in 30Dor-10 will either not have that massive stellar excess, or the stellar mass distribution is not predetermined by the current core mass function. The discussion of stability above raises the natural question of whether the clumps and cores whose mass distribution is being analyzed should have any correspondence with the stellar mass function, especially if many of those cores are not gravitationally bound or collapsing. To test this, we analyzed the distribution of only cores with Σ virial /Σ less than a threshold (we tried thresholds between 2 and 4; see Figure 11 and discussion in section 4.5 about the normalization of Σ virial /Σ). Interestingly, the shape of the mass distribution, and the fitted slope of the upper end are the same for the more gravitationally bound subset(s) of cores. If the mass distribution of cores is independent of core stability, and is related to the stellar IMF, then this result supports the notion that the core and cloud structure predetermines the stellar IMF, before gravitational forces become dominant.

CONCLUSIONS

High-resolution (<0.1pc) observations of the 30Dor-10 molecular cloud 15 pc north of R136 are presented. The CO emission morphology contains clumps near the locations of known mid-infrared massive protostars, as well as a series of parsec-long filaments oriented almost directly towards R136, most of which don't show signs of embedded star formation. The aligned filaments could possibly be "pillars" left behind by photoionization, or could be a relic of the initial collapse and formation of 30Doradus and the R136 star cluster.

Analysis of the cloud substructures is performed by segmenting emission into disjoint approximately round "cores" using clumpfind, by considering the hierarchical structures defined by isointensity contours using dendrograms, and by segmenting into disjoint long thin "filaments" using Filfinder. The balance between gravitational and kinetic energy of "cores" and dendrogram branches are analyzed with formulae appropriate to spheres and ellipsoids, and the balance is analyzed for filaments using formulae appropriate to infinite cylinders. We find that the filaments have widths of ∼0.1pc, similar to those in solar neighborhood clouds.

There is elevated kinetic energy (linewidths at a given size scale) in 30Dor-10 compared to other LMC and Galactic star formation regions. The slope of the size-linewidth relation is also a bit steeper than on those other regions, although not dramatically so. A steeper slope (more energy at larger scales), and a lack of correlation with local sources of kinetic energy (stellar winds and protostellar outflows), suggests that energy is injected on large (100s of pc to kpc) scales. This agrees with the analysis of high-J CO and far-infrared line emission at low spatial resolution in 30 Doradus Lee et al. (2019), and the existence of kpc-scale colliding filaments in the region (Fukui et al. 2017).

Clumps and cores when analyzed as entire objects with clumpfind lie in a similar part of size-linewidth-mass parameter space as Milky Way and other molecular clouds, and are consistent with free-fall collapse or virial equilibrium with moderate external pressure. A significant fraction of dendrogram structures and filaments have mass surface densities Σ LT E or line masses M l in excess of the corresponding quantities if gravitational and kinetic energies were in balance. The discrepancy can be resolved if the 13 CO/H 2 abundance is a factor of 5 or more higher than the value that was assumed here. Alternately, one physical explanation could be that the small scale structures and filaments in this cloud have significant magnetic support against gravity.

The upper end of the "core" mass distribution is consistent with a power-law with the same slope as the stellar initial mass function. Slopes are fitted directly to the differential mass distribution as well as using the aforementioned Hill maximum likelihood estimator. The systematic biases of direct fitting are more evident than when fitting the cumulative distribution. The best estimate slope of -1.3±0.15 is consistent with the stellar IMF slope -1.35. The shape of the mass distribution including the fitted slope does not change if only a subset of cores that have large Σ/Σ vir are considered; the cores that are more likely to collapse and form stars according to our stability measurement do not have a statistically different mass distribution than those that are less likely to collapse. This fact, and the fact that the most reliably measured part of the mass distribution has the same slope as the stellar IMF, support the notion that the stellar IMF is predetermined by the characteristics of turbulent fragmentation in the pre-collapse molecular cloud. A shallower mass function (more massive stars) has been found for the main-sequence massive stars in 30Doradus (Schneider et al. 2018). The next generation of star formation in 30Dor-10 will either not have that massive stellar excess, or the stellar mass distribution is not predetermined by the current core mass function.

Figure 1 .

 1 Figure1. The black contour shows the 3σ 12 CO 2-1 peak intensity 30Dor-10, overlaid on a RGB composite of HST F160W∼H, F110W∼J, and F658N (Hα). The CO beam is the small black dot to the left of the "12CO2-1" label. The dominant sources of ionizing radiation are the stars to the southwest including R136. Pillar structures P1,2,3, and a small cluster in the center of the cloud, are labeled in blue, and the three Spitzer-identified massive YSOs in this region are marked in cyan -see text for discussion. Figure2shows the CO structure in more detail. Figure3shows zooms of the three regions in magenta.

Figure 2 .

 2 Figure 2. Peak and integrated 13 CO 2-1 and 12 CO 2-1 at 0.1-pc resolution in 30Dor-10 (beams are small white dots to the left of each label).

Figure 3 .

 3 Figure3. Zooms on three sections of 30Dor-10, with the same image stretch and contour levels. Colors are red: 1mm ALMA continuum in square root scale, green: F160W HTTP HST in log scale, blue: F658N (Hα) HTTP HST in log scale, white contours: 12 CO 2-1 peak brightness at 0.125,0.2,0.275 Jy bm -1 (cube rms = 3.6mJy bm -1 ), magenta contours: 13 CO 2-1 beak brightness at 0.04,0.06,0.08,.1 Jy bm -1 (cube rms = 2.2mJy bm -1 ). H2O masers are marked with yellow circles of radius=1 , the quoted positional uncertainty(Imai et al. 2013). An OH maser is marked similarly with a cyan circle(Brogan et al. 2004).

Figure 4 .

 4 Figure 4. Filaments identified in peak intensity maps of 12 CO 2-1 (L) and 13 CO 2-1 (R). Filaments that match with elongated dendrogram structures are shown in blue, with their corresponding dendrogram structures in red. Filaments that don't match elongated dendrogram structures are marked in green. Elongated dendrogram structures that don't match filaments are shown in yellow.

Figure 5 .

 5 Figure5. Assessment of validity of the commonly used "LTE method": N( 13 CO) calculated from modeled 12 CO and 13 CO brightness temperatures, compared to the true values input to each model. The top panel is the average NLT E /Ntrue for all models with given 12 CO 2-1 and 13 CO 2-1 brightness temperatures. The contours are a histogram of the observed values in our cubes. The lower panel shows the value of the error for which 90% of the models are closer to the true value. The method works fairly well, overestimating the true column density by up to a factor of 2 at 90% confidence.

Figure 6 .

 6 Figure 6. (L) The intensity-weighted first moment map of 12 CO 2-1, with filaments overlaid. At each point, the local gradient in the 2-dimensional velocity field was measured. (R) A histogram of the angles between each filament branch's skeleton and the local velocity gradient -an angle of zero would result of the velocity gradient were completely aligned with the filament branch. Histograms are shown for all filament branches, and for those that are matched to elongated dendrogram structures. No strong trend is evident in the alignment.

Figure 7 .

 7 Figure 7.Top: A representative filament, showing 12 CO and 13 CO peak intensity as blue and green. Bottom: Line mass calculated three different ways along the filament (green, blue, orange dots). Critical mass for only thermal support (red), and thermal+nonthermal support (magenta), are also plotted. The more visually fragmented portion of the filament ∼1pc along the filament (the "gap") has M l << Mcrit. The bright knot ∼1.7pc along appears to be unsupported against collapse, with M l > M crit,nonthermal .

Figure 8 .

 8 Figure 8. Filament stability: points with M l > Mcrit have more mass per unit length than thermal or kinetic support can prevent from collapsing.(L) points colored by 12 CO 2-1 brightness temperature from lower (blue) to higher (red). "x" symbols are those filaments that correspond well to elongated isointensity contours (Figure4). (R) points colored by 8µm brightness from lower (blue) to higher (red).

  Figure 9.

  : resolution of 1.21 ×0.84 , Tokuda et al. (2019): resolution of 0.29 ×0.25 , Fukui et al. (2015):resolution of 1.3 ×0.8 , Fukui et al. (

Figure 10 .

 10 Figure 10. (L) Size-linewidth relation in 30Dor-10 for dendrogram (green) and clumpfind (magenta) structures. Overlaid is the slope fitted to dendrogram structures log(σv) = a log(r) + b.(R) Comparison of size-linewidth relations between 30Dor-10, the nearby less evolved LMC massive star formation region N159, Galactic molecular cloud PerseusA, and the "historic" Milky Way disk relation (dashed line; the data from which that was derived are all off the right-hand side of the plotSolomon et al. 1987;Heyer et al. 2009). Data from all regions were re-analyzed the same way, with clumps segmented using dendrograms.

Figure 11 .

 11 Figure 11. Relation between virial surface density 5σ 2 v /(GπR) and LTE surface density MLT E /(πR 2 ). Virial equilibrium between gravitational and kinetic energy Σvir=ΣLT E and gravity-driven collapse Σvir=2ΣLT E are marked as dashed lines. Virial equilibrium with external pressures nT =(10 5 ,10 6 )cm -3 K are solid green curves. (L) the different segmentation methods in this region (green=cprops, magenta=dendrogram) are compared to previous measurements of the galactic center (CG; Oka et al. 2001) and Milky Way disk (MW; Heyer et al. 2009). (R) Only dendrogram structures in this region are shown, now colored by size from smaller=blue to larger=red.

Figure

  Figure 12. Relation between virial surface density 5σ 2v /(GπR) and LTE surface density MLT E /(πR 2 ) for a model spherical core. The core has constant density within rinner, and a power-law radial density distribution outside of that. Each part of the core within a projected radius R is assumed to follow the size-linewidth relation σv ∝ R α . Calculation of the two surface densities for a partial core will yield significantly higher Σ and moderately lower Σvir compared to the values calculated for the entire core.

Figure 13 .

 13 Figure 13. (L) cumulative luminosity distribution of 12 CO 2-1 clumpfind-identified cores with two different sets of segmentation parameters. Power-law fits are shown using a maximum likelihood estimator (Clauset, Shalizi, & Newman 2009, eq 3.1; dashed) and a simple linear fit to the log number (dotted). (R) differential LTE mass function of cores with three different segmentation parameters.

Software: CASA (v4.3.1; v4.5; v5.1.0 McMullin et al. 2007), Kapteyn (Terlouw & Vogelaar 2015), filfinder (Koch & Rosolowsky 2015), quickclump (https://github.com/vojtech-sidorin/quickclump/), clumpfind (Williams et al. 1995(Williams et al. , 2011)), cprops (Rosolowsky & Leroy 2006, 2011), dendrograms (Rosolowsky et al. 2008) This paper makes use of the following ALMA data: ADS/JAO.ALMA#2011.0.00471.S, ADS/JAO.ALMA#2013.1.00346, ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. This research has made use of NASAs Astrophysics Data System Bibliographic Services. RI was partially supported during this work by NSF AST-1312902.


J Hinderer

B Hector

U Riccardi

Séverine Rosat

J-P Boy

M Calvo

F Littel

J-D Bernard

Hinderer

F & Bernard

A study of the monsoonal hydrology contribution using a 8 year record (2010-2018) from superconducting gravimeter OSG-060 at Djougou (Benin, West Africa)

Keywords: time variable gravity, Africa, hydrology, loading, monsoon

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A study of the monsoonal hydrology contribution using a 8-year record (2010-2018) from superconducting gravimeter OSG-060 at Djougou (Benin, West Africa) Introduction The superconducting gravimeter OSG-060 from GWR Instruments was installed in Nalohou, close to Djougou in northern Benin (West Africa) in July 2010. This installation was done in the frame of the GHYRAF (Gravity and Hydrology in Africa) program that lasted from 2008 to 2012 [START_REF] Hinderer | Land water storage changes from ground and space geodesy: first results from the GHYRAF (Gravity and Hydrology in Africa) experiment[END_REF]). This gravimeter is part of the IGETS (International Geodynamics and Earth Tides Service, (http://igets.u-strasbg.fr) (see also [START_REF] Boy | Superconducting Gravimeter Data from Djougou -Level 1[END_REF], under the umbrella of IAG (International Association of Geodesy). Fig. 1 shows the location of the Djougou (DJ) station (the geographical coordinates of the station are 9.7424 °N and 1.6056 °E) among the other stations of the international network of superconducting gravimeters (SG). Fig. 1. A map of the superconducting gravimeter stations belonging to IGETS (International Geodynamics and Earth Tides Service, http://igets.u-strasbg.fr ). The operational stations in 2018 (about 20) are in red and former stations now stopped are in blue. The location of Djougou in Benin makes this station unique for investigating the hydrology contribution of the West African monsoon. The choice of this location for installing the SG was driven by the fact that Nalohou is part of the Ara catchment (see Fig. 2) , one of the key sites in the Sudanian region monitored by AMMA-CATCH hydro-meteorological observatory (http://www.amma-catch.org/, AMMA-CATCH, 1990) with a high density of data available over a long time span of several decades [START_REF] Galle | AMMA-CATCH a Critical Zone Observatory in West Africa Monitoring a Region in Transition[END_REF]. The right hand side of Figure 2 (Fig. 2d) shows the location of the various hydro-meteorological instrumentation deployed (piezometers, neutron probe boreholes, rain gauge, Parshall, trenches) over the catchment as well as the gravimeter (FG5, SG-060, CG5) stations and MEP (Multi-electrode Electrical Profiling) profiles. In this study we use an absolute gravimeter FG5#206 from Micro-g Solutions Inc., and a superconducting gravimeter SG-060 from GWR Instruments Inc. .The huts for episodic absolute gravity measurements (FG5#206) and continuous superconducting gravimeter monitoring (GWR OSG-060) are very close to each other (a few meters) and located on the summital part of the catchment. Since 2011 there is also a network of relative gravity stations that has been regularly repeated with a Scintrex, model CG-5 (serial number 9379) gravimeter.  To observe non-invasively a signature of the internal redistribution of water within the critical zone for hydrological processes identification (proxy)  To provide additional constraints for hydrological modeling  To inform ultimately on groundwater recharge and annual storage/discharge that is critical for resource managers. The question of the hydrological signature in gravity has been extensively studied by many groups possessing superconducting gravimeters. It would be long to cite them all and we refer the reader to the review done in [START_REF] Hinderer | Superconducting gravimetry, in treatise on geophysics[END_REF], as well as in the more recent one by [START_REF] Van Camp | Geophysics from terrestrial time-variable gravity measurements[END_REF]. These studies mostly concern Europe (e.g. Creutzfeld et al. 2010;[START_REF] Weise | The Improved Hydrological Gravity Model for Moxa Observatory, Germany[END_REF][START_REF] Camp | Hydrogeological investigations at the Membach station, Belgium and application to correct long periodic gravity variations[END_REF][START_REF] Boy | Diurnal atmospheric tide and induced gravity variations[END_REF]) and hence the Djougou station is interesting by its unique location in Africa and because the investigation of monsoonal effects in West Africa has shown a clear intensification of rainfall and hydrological cycles that might be linked to climate changes (see [START_REF] Taylor | Frequency of extreme Sahelian storms tripled since 1982 in satellite observations[END_REF][START_REF] Galle | AMMA-CATCH a Critical Zone Observatory in West Africa Monitoring a Region in Transition[END_REF][START_REF] Nkrumah | Recent Trends in the Daily Rainfall Regime in Southern West Africa[END_REF]. Gravity data processing The gravity record has been processed jointly with air-pressure in a classical way starting with a preprocessing to remove the major disturbances (spikes, offsets, gaps, offsets) of instrumental and/or manmade origin. We refer the reader to [START_REF] Hinderer | Superconducting gravimetry, in treatise on geophysics[END_REF] for more details on the processing steps. The gaps due to earthquakes are filled with a synthetic local tide coming from a precise tidal analysis of the 8-year data set [START_REF] Hinderer | A study of the solid Earth tides, ocean and atmospheric loadings using an 8-year record[END_REF]. The cleaned pressure and gravity data sets are shown on Fig. 4 for the investigated period (2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018). The reduction for lunisolar tides (solid Earth + ocean tidal loading), as well as for atmospheric effects and Earth's rotation leads to the gravity residuals (in black on Fig. 5). We refer the reader to [START_REF] Hinderer | A study of the solid Earth tides, ocean and atmospheric loadings using an 8-year record[END_REF] for more details on the solid Earth tidal model as well as on the ocean tidal model and atmospheric corrections for the OSG-060 gravity record. A further correction for the instrumental drift of the gravimeter (in red on Fig. 5) leads then to the corrected residual signal (in blue on Fig. 5) that is reportedly of hydrological origin and characterized by an annual periodicity (see [START_REF] Hector | Hydro-gravimetry in West-Africa: First results from the Djougou (Benin) superconducting gravimeter[END_REF]). The drift assessment through least-square fitting shows that the instrumental drift of the GWR OSG-060 is composed of an initial exponential term followed by a linear term that remains stable over the years. This is very classical for SGs and confirms the earlier results based on a 1.5 year record of SG-060 (Hinderer et al. 2014a). Fig. 4. Gravity (upper panel, in nm s -2 ) and barometric pressure (lower panel, in hPa) observed at the Djougou station from 2010 to 2018. The fit to an exponential function g = g 0 + A 1 *exp(-(x-x 0 )/T 1 ), where g 0 and x 0 are the starting values of gravity and time, leads to an amplitude A 1 of -147.6 ± 1.2 nm s -2 and a time decay T 1 = 13.8 ± 0.3 days. This exponential drift estimate was done on the residual gravity first corrected for MERRA global hydrology model [START_REF] Reichle | Assessment of MERRA-2 land surface hydrology estimates[END_REF]) (see below for this choice) because this contribution mixes up with the instrumental part and must be removed. MERRA model is here the version 2 of The Modern-Era Retrospective analysis for Research and Applications. This hydrology model is available at the EOST loading service (http://loading.u-strasbg.fr). The linear drift is + 213.6 ± 0.1 nm s -2 /year which is rather large for SGs [START_REF] Hinderer | Superconducting gravimetry, in treatise on geophysics[END_REF], but linearity is very strong and apparently stable in time from our tests.

Fig. 5. Residual gravity signal (in black) after correction for solid Earth and ocean tidal loading, atmospheric effects and Earth's rotation; corrected gravity residuals (in blue) after subtracting a linear instrumental drift (in red). Solid Earth tides, ocean tidal loading and atmospheric loading

We analyzed the whole dataset of gravity observations (nearly 8 years) using the [START_REF] Hartmann | The HW95 tidal potential catalogue[END_REF] tidal potential catalogue. Amplitude and phase tidal parameters for groups combining inseparable constituents as well as the air pressure admittance factors were adjusted in the long period, diurnal and sub-diurnal frequency bands by means of ET34-ANA v7.1 software [START_REF] Schüller | Theoretical basis for Earth tide analysis and prediction[END_REF]Ducarme & Schüller 2018). A detailed discussion of the tidal (solid Earth + oceans) effects, as well as the atmospheric contribution to gravity, is done in [START_REF] Hinderer | A study of the solid Earth tides, ocean and atmospheric loadings using an 8-year record[END_REF]. The equatorial location of Djougou station leads to strong pressure waves of thermal origin modulated in amplitude (e.g. [START_REF] Gegout | Influence of atmospheric pressure on the Free Core Nutation, precession and some forced nutational motions of the Earth[END_REF][START_REF] Schindelegger | Surface Pressure Tide Climatologies Deduced from a Quality-Controlled Network of Barometric Observations[END_REF]. The air pressure also exhibits long period features among which the annual component Sa (one cycle per year) is the largest and the other terms are at Ssa (2 cycle per year) and Sta (3 cycle per year). The gravity loading due to atmospheric pressure at Djougou takes into account the gravitational attraction of the air masses as well as the elastic deformation of the Earth and hence cannot be reduced to a simple function of the local pressure like in the frame of the classical concept of barometric admittance [START_REF] Crossley | Effective barometric admittance and gravity residuals[END_REF]Hinderer et al. 2014b) but is more complicated because the non-local part involves the convolution of the pressure field worldwide with the atmospheric loading Green's function (see e.g. [START_REF] Boy | Reduction of surface gravity data from global atmospheric pressure loading[END_REF][START_REF] Boy | Diurnal atmospheric tide and induced gravity variations[END_REF]. Different atmospheric gravity loadings were computed using ECMWF (European Centre for Medium Range Weather Forecasts) operational and reanalysis (ERA interim) pressure data, assuming either an inverted barometer response of the oceans (without any induced pressure effect at the ocean bottom see e.g. [START_REF] Wunsch | Atmospheric loading and the oceanic "inverted barometer" effect[END_REF] or a dynamic response using TUGO-m [START_REF] Carrère | Modeling the barotropic response of the global ocean to atmospheric wind and pressure forcing -comparisons with observations[END_REF]. The best pressure reduction was found using an hybrid approach computed by retaining the MERRA modelled non-local component and replacing the local one with the observed air pressure multiplied by the nominal admittance coefficient coherently with model resolution. In the MERRA model the resolution is 0.625 ° and we used a local cell with a radius of 0.1 ° together with a nominal admittance of -2.21 nm.s -2 /hPa. The hybrid method enabled us to account for both very local atmospheric effects, as probed by the surface pressure observations, and larger scale contributions to gravity originating from the planetary thermal pressure waves. Hydrological loading There are several global hydrology models available at http://loading.u-strasbg.fr. In this study we used GLDAS/Noah [START_REF] Rodell | The global land data assimilation system[END_REF], ECMWF operational and reanalysis (ERA interim) and MERRA2 [START_REF] Reichle | The MERRA-Land Data Product[END_REF][START_REF] Reichle | Assessment of MERRA-2 land surface hydrology estimates[END_REF]. Global hydrology models are based on a land surface model (for example, Noah for GLDAS) forced by atmospheric parameters, such as precipitation, temperature, long and short wavelength radiations, etc. They can be run fully coupled with an atmospheric model in case of the ERA interim and MERRA2 reanalysis, or independently as for GLDAS using variables from various origins (atmospheric model, satellite data, etc.). The other main discrepancy is the difference between operational and reanalysis model. In the first case, the model and/or the forcing can change in time, whereas the entire system is stable for reanalysis. All the hydrology models used in this study are providing estimates of soil moisture content in the first meters, and also canopy water in the case of GLDAS/Noah. We compare all these different models as they can show some discrepancies (see Figure 6), due to differences in the land surface models, but also in the forcings. Fig. 6. Comparison of gravity changes due to hydrology according to global models (ERA, GLDAS, MERRA2) with superconducting gravity residuals at Djougou. For sake of comparison gravity loading from these three hydrology models, computed for the period 2008-2018, are superimposed onto the SG gravity residuals (Fig. 6), which are only available after July 2010 (installation of gravimeter). The beginning of the SG residuals in summer 2010 is clearly different from the hydrology predictions because of the initial instrumental drift of the SG after installation. By computing the discrepancy between the observed gravity and the hydrology models we are able to find out the preferred model that leads to the best fit. This is shown by the histograms in Fig. 7 depicting the differences between gravity observations and model predictions. The RMS is 42.9, 40.4, 36.8 and 26.1 nm s -2 for models GLDAS/Noah, ERA Interim, ECMWF and MERRA2, respectively. Hence, the best model results to be clearly MERRA2. Fig. 7. Histograms of the differences between observed gravity and 4 global hydrology models. As already shown for shorter duration in previous papers [START_REF] Hector | Hydro-gravimetry in West-Africa: First results from the Djougou (Benin) superconducting gravimeter[END_REF][START_REF] Hector | Water storage changes as a marker for baseflow generation process in a tropical humid basement catchment: insights from hybrid gravimetry[END_REF], there is a fair agreement between the continuous superconducting gravimeter observations, episodic absolute gravity measurements and space gravimetry GRACE data that all superimpose onto global hydrology models. The comparison of MERRA hydrology model with OSG-060, FG5#206 and GRACE data is shown in Fig. 8 where the GRACE solutions come from iterated global "mascons" [START_REF] Luthcke | Antarctica, Greenland and Gulf of Alaska land ice evolution from an iterated GRACE global mascon solution[END_REF]. We use the latest version (v02.4) of the 1-degree equal area global iterated "mascon" solution provided by NASA Goddard Space Flight Center from January 2003 to July 2016. A mascon (mass concentration) represents an excess or a lack of surface mass compared to the a priori mean gravity field in a predefined region, directly inverted from the GRACE K-band range rate measurements. The excess or lack of mass is represented as a uniform water layer, expressed in centimeters of equivalent water over this area. These "mascon" solutions differ from the classical constrained or unconstrained spherical harmonic solutions, as they allow to introduce spatial constraints in the inversion of the Kband range rate residuals [START_REF] Rowlands | Global mass flux solutions from GRACE: A comparison of parameter estimation strategies -Mass concentrations versus Stokes coefficients[END_REF], and to increase the spatial resolution of the solution, as post-processing filtering is no longer required. Fig. 8

. Comparison of MERRA2 hydrology model with superconducting gravimeter (OSG-060), absolute gravimeter (FG5#206) and GRACE data.

One has to keep in mind that the spatial resolution of these data sets is quite different. Surface gravity is a point measurement, MERRA has a 0.625 ° grid resolution and GRACE data merely are representative of much larger scales (typically 300 x 300 km). The fact that there is a fair agreement between surface gravity and GRACE in Fig. 8 mostly comes from the fact that the hydrology loading of monsoonal origin is coherent over a wide region in West Africa, and probably also because of the specific uphill location of the gravimeter. Fig. 9

. Comparison of gravity observed with the superconducting gravimeter (OSG-060) (black line) and predicted by MERRA global hydrology model (green line), soil moisture (red line) , water table depth (blue line) and rain data in mm/day.

Our last absolute gravity data unfortunately was measured in 2013 and GRACE data stopped in 2016. In Fig. 9 we show for the period 2010-2015 the comparison of MERRA hydrology model with gravity (OSG-060) (like in Fig. 8) but superimposed onto local hydrological data such as soil moisture (available between 0 and 1.2m using 6 probes at depths of 5, 10, 20, 40 and 60 cm, AMMA-CATCH, 2005), water table depth [START_REF] Rain_Od | Soil dataset (soil moisture, temperature, and succion profiles),within the Donga watershed (600 km2), Benin[END_REF] and rain data(AMMA-CATCH 1999) in the Ara catchment available from the data base of the AMMA-CATCH observation system [START_REF] Galle | AMMA-CATCH a Critical Zone Observatory in West Africa Monitoring a Region in Transition[END_REF], http://www.ammacatch.org/).

Fig. 10. Correlation coefficient as a function of time between observed gravity (SG) and soil moisture (SM) and water table (WT).

In this ground context, mainly characterized by direct water infiltration, soil moisture reacts more quickly to the rain than the deeper water table. Fig. 10 shows the correlation coefficient computed using a 2-month sliding window (to filter out the higher frequency contributions and keep only the major processes during the hydrological cycle) between SG and both water table and soil moisture. At the onset of the rain season, the level of the water table continues to decrease, the discharge being driven by the hydraulic gradient at the catchment scale [START_REF] Hector | Water storage changes as a marker for baseflow generation process in a tropical humid basement catchment: insights from hybrid gravimetry[END_REF][START_REF] Hector | Hydrological functioning of West-African inland valleys explored with a critical zone model[END_REF]. This results in a phase lag inducing a decorrelation of the water table and the gravity signal, due to the time offset between the yearly minima of the two series. At the same time, the rain water is being stored in the shallow soil and the correlation between soil storage and gravity takes over until the core of the rainy season, where the soil storage compartment is full (soil moisture time variability shows a plateau), while still maintaining a downward infiltration flux, recharging the water table and increasing the gravity. Thereafter the gravity correlates well to the water table, until the end of the rainy season, where gravity starts to decrease due to the system discharge by evapotranspiration, while the water table still rises in response to the slowly percolating water. During the core of the dry season, all the system discharges and both water table heights and soil moisture values are well correlated to gravity. Fig. 11. Modeling the WSC using a simulation of evapotranspiration (ET), streamflow and induced gravity effect based on Parflow-CLM model [START_REF] Hector | Hydrological functioning of West-African inland valleys explored with a critical zone model[END_REF]) and its validation through the observed gravity. This means that the transfer function between gravity and water storage changes is variable in time because several processes occur in the water table and soil moisture during the hydrological cycle driven by the monsoon; moreover the WSC are not uniform over the entire catchment as already shown by the analysis of the repeated microgravimetric network observations that could identify zones with preferential infiltration [START_REF] Hector | Water storage changes as a marker for baseflow generation process in a tropical humid basement catchment: insights from hybrid gravimetry[END_REF]. Finally in Fig. 11, the observed gravity is superimposed onto the WSC originating from a simulation [START_REF] Hector | Hydrological functioning of West-African inland valleys explored with a critical zone model[END_REF]) of a small headwater catchment (720 m x 300 m) encompassing the Nalohou catchment (Fig. 2) with the help of the physically-based critical zone model ParFlow-CLM [START_REF] Maxwell | Development of a coupled land surface and groundwater model[END_REF][START_REF] Kollet | Capturing the influence of groundwater dynamics on land surface processes using an integrated, distributed watershed model[END_REF] run at a 20 m lateral resolution. ParFlow solves the water transfer in porous media on a gridded domain using the 3D Richards equation for the saturated and nonsaturated zone, and simultaneously solves the overland flow using the kinematic wave and Manning equation on the surface cells which exhibit a positive pressure [START_REF] Kollet | Integrated surface-groundwater flow modeling: A free-surface overland flow boundary condition in a parallel groundwater flow model[END_REF]. ParFlow is tightly coupled to CLM, the Common Land Model [START_REF] Dai | The Common Land Model[END_REF][START_REF] Maxwell | Development of a coupled land surface and groundwater model[END_REF], a land surface model which solves the water and energy budget through the user-defined N first cells. The ParFlow model is forced at a 0.5hr time step for 11 years with atmospheric and vegetation variables from the AMMA CATCH observational service; further details on the physical parametrization of the model domain as well as on the model set up can be found in [START_REF] Hector | Hydrological functioning of West-African inland valleys explored with a critical zone model[END_REF]. Three additional years are used as compared to the paper by [START_REF] Hector | Hydrological functioning of West-African inland valleys explored with a critical zone model[END_REF]. Simulation outputs have been extensively tested against multiple field data (streamflow, evapotranspiration, soil humidity, water storage). The SG is seated on the South-East edge of the catchment, on a saprolite unit. Therefore, only the saprolite simulation from [START_REF] Hector | Hydrological functioning of West-African inland valleys explored with a critical zone model[END_REF] is shown on Fig. 11. As the SG is located on the edge of the simulation domain, no attempt to calculate the gravity effect of the simulated water storage change (WSC) is done. This would have required a finer resolution close to the gravimeter, and the explicit simulation of the shelter effect as described in [START_REF] Deville | On the impact of topography and building mask on time varying gravity due to local hydrology[END_REF] or in the detailed study by Reich et al. (2018). Instead, the simulation represents the typical behavior of a small headwater catchment composed by saprolite rock cover, an elementary hydrological unit in the area. To convert the water storage changes (in mm of water equivalent) to gravity (nm s -2 ), the best fit (leading to the smallest differences in time between observations and predictions) yields a factor of 0.37 nm s -2 /mm that is lower than the classical Bouguer plate coefficient of 0.42 nm.s -2 /mm of water for a flat water layer of density 1000 kg/m 3 . This difference cannot be explained by the topography effect, as computed by Hector et al., 2015 (Their Fig. B1), which would imply a slightly higher coefficient of 0.44 nm.s -2 /mm. However, this is in accordance with the range of coefficients found when computing the effect of the shelter (Hector et al., 2014, Fig. 5) for an average depth of water storage changes between 2 and 3 meters, in agreement with the actual depth range comprised between the surface and the water table. Furthermore, this value is still higher than the 'rainfall admittance', the mean observed ratio of gravity increase over precipitation amount computed from a representative diversity of individual rain events, of 0.22 nm.s - 2 /mm. This value is the lower bound for fast water storage changes occurring close to the surface and for which the gravity measurement is very sensitive to the shelter extent and the sensor height above ground [START_REF] Hector | Hydro-gravimetry in West-Africa: First results from the Djougou (Benin) superconducting gravimeter[END_REF]. This observation clearly calls for a higher resolution simulation of the surroundings of the gravimeter, explicitly simulating the shelter effect to be able to calculate directly the gravity effect. Such a model would also require the simulation of the whole catchment to constrain the lateral groundwater flux, responsible for the groundwater discharge. The simulated water storage changes match to a reasonable extent the observed gravity changes behaviour (R² = 0.77), and while the annual cycle and the inter-annual variations are overall satisfying, there are clear discrepancies around specific events. This may encourage the use of the current model as a correction model for data gap-filling, for long term geodetic applications, but not for high frequency signal analysis, for which a higher resolution model, taking into account the local substratum variability in the vicinity of the gravimeter, together with the shelter effect, is needed (see [START_REF] Reich | Reducing gravity data for the influence of water storage variations beneath observatory buildings[END_REF]).

Conclusions

We report on an 8-year long continuous gravity monitoring (2010-2018) by the superconducting gravimeter OSG-060 at Djougou (Benin, West Africa). First, a tidal analysis using ET34-ANA v7.1 software enables us to remove the solid Earth tides, ocean tidal loading and Polar motion effects in the gravity record. Atmospheric pressure effects are reduced by means of a hybrid model using the global atmospheric model MERRA everywhere except in in the local zone, where the model pressure was replaced by the observed pressure accounting for the specific admittance factor. The observed residual gravity signal of hydrological origin was then investigated. We first compared this signal to several global hydrology models (ERA, GLDAS, MERRA) and found that the best model (i.e. the one leading to the smallest discrepancy between observations and model) is MERRA2 with a standard deviation of 26.1 nm s -2 . The MERRA hydrology model was also compared to absolute gravity data available from 2008 to 2013 and to GRACE mascons solutions between 2008 and 2016 when GRACE data stopped. In a further step, the gravity signal is compared to local hydrological data like soil moisture in the very superficial layer (0-1.2 m), water table depth and rainfall. The temporal evolution of the correlation coefficient between the gravity observation and both the soil moisture and the water table is well explained by the direct infiltration process of rain water together with the lateral transfer discharging the water table. An attempt to model water storage changes was finally done using a simulation based on Parflow-CLM model of the catchment. This model that is forced by rain is in agreement with evapotranspiration and stream flow data and leads to simulated water storage changes. These changes nicely superimpose onto gravity signal using a 0.37 nm s -2 /mm conversion factor. The next step will be to use a finer resolution model in the vicinity of the gravimeter, together with the explicit simulation of the shelter effect, in order to accurately simulate the observed gravity changes and use it as a gap-filling support tool. The long term monitoring of WSC will help identify possible effects of land use intensification currently occurring in the area (lower duration crop rotation cycles, increased use of pesticides resulting in longer periods of soil denudation, village growth...) on groundwater recharge through increase of surface runoff. Furthermore, a major science question still unresolved in the area is whether the observed groundwater discharge during the dry season is driven by evapotranspiration only, or also by deeper drainage through permeable zones like faults. A comparative analysis of the gravity and water table decreases together with local measurements of evapotranspiration may help to further constrain this issue [START_REF] Descloitres | The contribution of MRS and resistivity methods to the interpretation of actual evapo-transpiration measurements: a case study in metamorphic context in north Bénin[END_REF]. Acknowledgments We acknowledge the financial support of INSU-CNRS for operating SG 060 in Djougou as part of the French Observatory Service (SNO Gravimétrie). The AMMA-CATCH regional observing system (which is part of the OZCAR research infrastructure in critical zones) was set up thanks to an incentive funding of the French Ministry of Research that allowed pooling together of various pre-existing small-scale observing set-ups. The continuity and long-term perennity of the measurements are made possible by continuous IRD funding since 1990 and by continuous CNRS-INSU funding since 2005. All (or most of) the computations presented in this paper were performed using the Froggy platform of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (GRANT CPER07_13 CIRA), the OSUG@2020 labex (reference ANR10 LABX56) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d'Avenir supervised by the Agence Nationale pour la Recherche. Special thanks to Théo Ouani (IRD) and Saré for their continuous help on site.

Fig. 2 .

 2 Fig. 2. Location of the Ara catchment inside the Upper Ouémé catchment in North Benin and available hydro-meteorological and geophysical instrumentation (Hector et al. 2015); 2a) shows the upper Ouémé catchment in northern Benin, 2b shows the Ara catchment inside the upper Ouémé catchment, and 2c shows the location of the Nalohou catchment inside the Ara catchment. 2d) indicates all available instrumentation in hydrology and geophysics in the Nalohou catchment (WGS 84 / UTM zone 31N projected coordinate system). NP means Neutron probe, MEP Multi-electrode Electrical Profiling, Parshall a Parshall flume for measuring the flow of water and Bas-Fond is the French and local name for an inland valley wetland.

Geophysical Journal International