Cenozoic Dynamic Topography of Madagascar

It has been proposed that Oligo‐Miocene regional uplift of Madagascar was generated and is maintained by mantle dynamical processes. Expressions of regional uplift include flat‐lying Upper Cretaceous‐Paleogene marine limestones that crop out at elevations of hundreds of meters along the western seaboard and emergent Quaternary coral‐rich terraces that rim the coastline. Here, we explore the history of subcrustal topographic support through a combined analysis of four sets of observational constraints. First, we exploit published receiver function estimates of crustal thickness and spectral admittance between gravity and topography. An admittance value of ∼+40 ± 10 mGal km−1 at wavelengths >500 km implies that ∼1 km of topography is supported by subcrustal processes. Secondly, new apatite fission‐track and helium measurements from 18 basement samples are inverted, constraining temperature and denudation histories. Results suggest that 0.5–1.6 km of regional uplift occurred after ∼30 Ma. Thirdly, we calculate a history of regional uplift by minimizing the misfit between observed and calculated longitudinal river profiles. Results suggest that topography was generated during Neogene times. Finally, inverse modeling of rare earth element concentrations in Neogene mafic rocks indicates that melting of the asthenospheric source occurred at depths of ≤65 km with potential temperatures of 1300–1370 °C. Melting occurred at higher temperatures beneath Réunion Island and northern Madagascar and at lower temperatures beneath the Comores and southern Madagascar. These inferences are consistent with shear wave velocities obtained from tomographic models. We conclude that Madagascar is underlain by thinned lithospheric mantle and that a thermal anomaly lies within an asthenospheric layer beneath northern Madagascar.

. This Late Cretaceous flood basalt province is often attributed to the Marion hotspot (e.g., Storey et al., 1995). During Cenozoic times, Madagascar was bordered by passive margins and marine deposition occurred along the western and northern coastlines. Small amounts of magmatic activity occurred between 51.4 and 23.5 Ma in the northern portion of the island, where basaltic tuffs within dolomitic limestones and scattered igneous intrusive rocks occur (Cucciniello et al., , 2017Estrade et al., 2014). Lithospheric extension commenced on the eastern seaboard during Neogene times (Bertil & Regnoult, 1998;de Wit, 2003). Global positioning system (GPS) campaigns combined with earthquake focal mechanisms and mapping of active fault scarps demonstrate that east-west crustal extension is ongoing at rates of several mm per year (Figure 2c; Bertil & Regnoult, 1998;Kusky et al., 2010;Rindraharisaona et al., 2013;Saria et al., 2014;Stamps et al., 2018Stamps et al., , 2021. This phase of extension is manifest by northeast-southwest oriented troughs filled with Neogene sedimentary rocks, which are now regarded as part of the East African Rift system (Kusky et al., 2010;Saria et al., 2014;Stamps et al., 2021). In post-Oligocene times (i.e., ∼23 Ma-Recent), mafic magmatism occurred in central, southern, and northern parts of Madagascar ( Figure 2d; Bardintzeff et al., 2010;Cucciniello et al., 2017Cucciniello et al., , 2018Emerick & Duncan, 1982;Estrade et al., 2014;Melluso et al., 2011Melluso et al., , 2016Rasoazanamparany et al., 2021). Monogenetic lava fields, stratovolcanoes, and dyke swarms were first described by Lacroix (1915) and have since been attributed variously to rift magmatism resulting from lithospheric thinning to emplacement of a mantle plume (Cucciniello et al., 2017;de Wit, 2003).
There is evidence that Madagascar sits within a region whose bathymetry is supported by mantle processes. Oceanic residual depth measurements from the Mozambique channel to the west, and from the Indian STEPHENSON ET AL. 10.1029/2020GC009624 2 of 38 Figure 1. Regional setting of Madagascar. (a) Long-wavelength (i.e., >800 km) free-air gravity field taken from GGM03C model (Tapley et al., 2005). (b) Average shear-wave velocity anomalies between 100 and 200 km depth that were extracted from tomographic model of Schaeffer and Lebedev (2013). Anomalies are relative to their version of AK135 model. (c) Topography and bathymetry of region encompassing Madagascar extracted from SRTM and ETOPO1 databases, respectively (Becker et al., 2009) Ocean to the east, have values of +0.25 to +1.5 km (Figure 2a; Hoggard et al., 2016Hoggard et al., , 2017. We have augmented these measurements along the Mozambique margin by exploiting the seismic wide-angle survey of Mueller and Jokat (2017), following the methodology described by Hoggard et al. (2016Hoggard et al. ( , 2017. Significant positive anomalies encompass Madagascar and demonstrate that the entirety of the region is dynamically supported. Although the relationship between gravity and topography fields can be complex, this conclusion is corroborated by the existence of a +40 mGal long-wavelength free-air gravity anomaly that protrudes northward from the southern ocean superswell (Figure 2b).
Several lines of evidence indicate that the southern Horombé, central Antananarivo, and northern Maromokotro Plateaux underwent regional uplift during Neogene times. First, Eocene marine sedimentary rocks fringe the coastline at elevations that range between sea level and ∼850 m which exceeds STEPHENSON ET AL.   (Besairie, 1964); colored circles = water-loaded residual depth measurements corrected for both sedimentary and crustal loading ; circles with bold black rims = residual depth measurements calculated using seismic survey of Mueller and Jokat (2017); colored upward-/downward-pointing triangles = upper/lower bounds for residual depth measurements where only sedimentary loading is corrected; colored linear tracks = selected shiptrack bathymetric measurements corrected for sedimentary thickness using global grid . (b) Map of long-wavelength free-air gravity anomalies filtered between 730 and 13,000 km taken from GGM03C model (Tapley et al., 2005). Red/black/blue lines = positive/zero/negative contours at 10 mGal intervals. (c) Beach balls = focal mechanisms of large (M w > 5) earthquakes taken from Centroid Moment Tensor Catalog (Dziewonski et al., 1981;Ekström et al., 2012); large black circles = events where M w > 5 from Mozambique Channel or events with no focal mechanisms elsewhere; small black circles = events where 3 < M w < 5 (Bondár & Storchak, 2011;Rakotondraibe et al., 2020); arrows = vectors showing horizontal displacement with respect to GPS station at Antananarivo determined from GPS measurements (Stamps et al., 2018); open ellipses = uncertainties for GPS measurements; red circle = Antananarivo. (d) Red polygons = Cenozoic igneous rocks; pink polygons = Cretaceous (i.e., Turonian-Santonian) igneous rocks. (e) Colored circles = loci of Upper Cretaceous and Cenozoic marine fauna and ichnofossils (see color bar for age); up/down vertical arrows = emergent/submerged Quaternary shoreline deposits (e.g., beachrock and coral reef) with respect to mean sea level (Battistini, 1959(Battistini, , 1965(Battistini, , 1977Robequain, 1953;Stephenson et al., 2019). glacio-eustatic amplitudes (Besairie, 1964;Delaunay, 2018;Stephenson et al., 2019). For example, Nummulites, a marine protozoan, is found in Eocene limestones of the Antsiranana basin at the northern tip of the island. Nummulites are also found within the Morondava basin in the southwestern corner together with scleractinian corals at elevations of 450 and 850 m, respectively (Table 1; Figure 2e; Alloiteau, 1958;Besairie, 1964;Stephenson et al., 2019). A variety of Miocene shipworm, Kuphus, and Mammalia sirenia occur on Nosy Makamby in association with fossil turtles and crocodiles at an elevation of ∼20 m (Collignon & Cottreau, 1927;Samonds et al., 2009;Samonds & Fordyce, 2019). Quaternary emergent coral reef terraces rim the northern and southern coastlines (Battistini, 1959(Battistini, , 1965(Battistini, , 1977Stephenson et al., 2019). Secondly, a switch from carbonate to clastic sedimentation and a marked increase in sedimentary flux onto the western margin have been identified that commence at 30 Ma. This rapid increase of eroded material coincides with stratigraphic evidence for several hundred meters of uplift of the coastline (Delaunay, 2018). Thirdly, the Glorieuse Seamount, which is located ∼150 km northwest of Madagascar, underwent ∼400 m of rapid uplift according to the chronology of submerged terraces. This putative uplift event started at 15 Ma (Leroux et al., 2020). Finally, qualitative analysis of lateritized erosional surfaces has been used to suggest that an originally subdued landscape was uplifted in a series of phases during Cretaceous and Neogene times (Delaunay, 2018;Dixey, 1960).
Exactly what constitutes dynamic topographic support is often debated in the literature (e.g., Davies et al., 2019;Molnar et al., 2015). Here, we use the generally accepted definition first described by Hager and Richards (1989), although we acknowledge that alternative definitions are sometimes employed. Deflections of the Earth's surface are driven by viscous stresses, which are generated by density contrasts within the convecting mantle. These deflections constitute dynamic topography and are given by where δa lm is dynamic topography of the Earth's surface and A l (r) is the normalized surface response kernels as a function of radius, r, for a putative density anomaly located at different depths within the mantle. a is the radius of the Earth, Δρ a is the density contrast at the surface boundary, and δρ lm (r) represents mantle density anomalies. Superscripts l and m refer to spherical harmonic degree and order, respectively. A l (r) depends upon the radial viscosity structure of the mantle. It is important to appreciate that A → 1 within the STEPHENSON ET AL. uppermost mantle for all values of l regardless of viscosity structure, which means that δa lm → δρ/Δρ a (e.g., Colli et al., 2016;. Consequently, the sum of the flow and isostatic components of dynamic topography within, say, an asthenospheric channel immediately beneath the lithospheric plate can be confidently approximated using a thermal isostatic relationship . This standard equation is predicated upon the fact that density variations and convective flow within the mantle produce radial stresses which are balanced by vertical deflections of the Earth's surface. These stresses are generated in three interlinked ways. First, changes in the thickness of the upper boundary layer (i.e., the lithospheric mantle) are isostatically compensated by surface deflections. This process gives rise to the well-known oceanic plate cooling trend but its continental expression is more complex (e.g., Guerri et al., 2016). For clarity, the definition of dynamic topography does not usually extend to including oceanic plate cooling (pace . Secondly, asthenospheric temperature anomalies that reside immediately beneath the plate must be isostatically compensated at the surface, even in the absence of flow (e.g., Richards et al., 2020). Finally, viscous flow itself creates pressure gradients that impose vertical tractions on the overlying plate (e.g., Lithgow-Bertelloni & Silver, 1998;Spasojevic & Gurnis, 2012). It is important to emphasize that all three of these components are automatically embedded within the definition of dynamic topography since they are inextricably interlinked from a fluid dynamical perspective. Hence, in this contribution, we use the term "dynamic topography" to mean topography arising from flow within the mantle, from the isostatic consequences of static thermal anomalies within the subplate mantle, and from changes in the depth to the base of the thermal boundary layer (i.e., the base of the lithospheric plate). A primary distinction is between vertical stresses acting within the mantle and horizontal stresses generated by plate tectonic processes.
In this contribution, we present and analyze four independent sets of observations with a view to investigating the extent to which Cenozoic regional uplift, denudation, and magmatism of Madagascar are caused by mantle processes. First, we exploit a database of published receiver function estimates of velocity structure and crustal thickness to establish crustal constraints on lithospheric structure and residual topography. This approach is combined with analysis of the spectral admittance between the free-air gravity and topographic fields to gauge the limit of flexural support for Malagasy plateaux. Secondly, we present and model a suite of newly acquired apatite fission-track and apatite helium measurements to place constraints upon the spatio-temporal history of exhumation. Thirdly, a calibrated inverse model is used to calculate a regional uplift history from a revised and augmented inventory of ∼2000 longitudinal river profiles. Finally, geochemical modeling is carried out upon a suite of published and newly acquired Neogene basalt samples to constrain the depth and extent of isentropic melting beneath the Malagasy lithosphere. In this way, we aim to elucidate geodynamic processes that link uplift, erosion, and magmatism of Madagascar.

Receiver Function Analyses
Positive residual bathymetry of oceanic lithosphere surrounding Madagascar implies that at least some fraction of its topography might not be supported by crustal isostasy. To estimate this fraction, it is instructive to compare the relationship between observed crustal thickness, density, and elevation with that relationship predicted by isostatic calculation. Andriampenomanana et al. (2017) jointly inverted receiver functions and surface waves to constrain the one-dimensional shear wave velocity structure at different locations across the island (Figure 3a). Crustal thickness decreases from 40-45 km beneath the Antananarivo Plateau to 30-35 km beneath the Maromokotro Highlands and Horombé Plateau and to <30 km in the Morondava and Antsiranana basins. Beneath the northern Antsiranana basin, the crust is as thin as 18 km. We note that these estimates broadly agree with those reported in other studies (e.g., Paul & Eakin, 2017;Rindraharisaona et al., 2017Rindraharisaona et al., , 2013).
An unknown component of topography is supported by density differences within the crust and within the sedimentary cover which obscures the signal caused by differences in crustal thickness variations and subcrustal processes. Here, we estimate crustal density by converting shear-wave velocity as a function of depth, v s (z), into P-wave velocity as a function of depth, v p (z), exploiting the empirical relationship presented by Brocher (2005). v p (z) is then converted into density as a function of depth, ρ c (z), using the Nafe-Drake relationship (Brocher, 2005;Ludwig et al., 1970 , where ρ cr = 2.8 Mg m −3 is a reference crustal density and t cc is crustal thickness. We assume an asthenospheric density of ρ a = 3.18 Mg m −3 . The isostatic correction accounts for differences in bulk density between regions where, for example, low-density sedimentary cover is present. If the uncertainty in conversion from seismic velocity to density is ± 0.05 Mg m −3 , we obtain an uncertainty for the corrected elevation of ±0.3 km if the crust is 18-km thick and ±0.5 km if the crust is 30-km thick. For transparency, the uncorrected measurements are shown in Figure S1. There is a demonstrable relationship between corrected elevation and crustal thickness (Figures 3c  and 3d). These corrected observations can be compared with predicted elevations obtained by balancing an idealized column of continental lithosphere against the density structure of continental lithosphere whose surface lies at sea level.
The elevation, e, of a column of continental crust with thickness t cc is given by where t cr = 35 km and ρ cr = 2.8 Mg m −3 are the thickness and density of reference continental crust at sea level. The difference between observed and calculated elevation is defined as residual topography, which can be generated either by changes in lithospheric mantle thickness or by subplate convection (Figure 3b). Our positive residual topographic estimates suggest that the lithospheric mantle is anomalously thin, that some component of topography is supported by subplate mantle processes, or a combination of the two STEPHENSON ET AL.
10.1029/2020GC009624 6 of 38  Hoggard et al. (2020). Colored circles = selected crustal thickness estimates, t cc , determined by joint inversion of receiver functions and surface wave dispersion measurements (Andriampenomanana et al., 2017). (b) Residual elevation, e r for low-pass filtered (i.e., wavelengths > 30 km) topography after accounting for crustal isostatic support given by reference column with crustal thickness, t cr = 35 km and density ρ cr = 2.8 Mg m −3 . (c) Surface elevation plotted as function of crustal thickness for different lithospheric thicknesses. Colored circles = elevation corrected for density differences as function of crustal thickness, assuming ρ cc = 2.8 Mg m −3 for low-pass filtered (i.e., wavelengths > 30 km) topography, where color refers to lithospheric thickness. Dotted lines = predicted elevation as function of crustal thickness for equilibrated lithosphere where a 1 = 50, 100, and 150 km, and ρ cc = 2.8 Mg m −3 , t cr = 35 km, and a r = 120 km; horizontal bar of cross at right-hand side = average uncertainty for crustal thicknesses determined from receiver function analysis; vertical bar of cross = average uncertainty of corrected elevation if uncertainty in density conversion, Δρ cc = 0.05 Mg m −3 for t cc = 30 km. (d) Same for disequilibrated lithosphere after instantaneous removal of 50 km from base of plate without thermal re-equilibration.
( Figure 3b). To place bounds upon the thickness of the lithospheric mantle beneath Madagascar, we calculate elevation as a function of crustal and lithospheric thickness for ρ cr = 2.8 Mg m −3 (Figures 3c and 3d). In the absence of subplate topographic support, this elevation is given by where a r = 120 km and a 1 are the thickness of the lithospheric plate for the reference and predicted columns, respectively. Values of lithospheric mantle densities, ρ mr and ρ m , are calculated by assuming a linear geothermal gradient within the lithospheric plate and an ambient asthenospheric temperature of 1370 °C at the base of the plate (i.e., the temperature of asthenospheric mantle at 100 km given potential temperature, T p = 1330 °C). The density of mantle rock at standard temperature and pressure is assumed to be ρ • = 3.33 Mg m −3 and density as a function of temperature is given by ρ • (1 − αT) where α = 3.28 × 10 −5 °C −1 (Section S1). Two end-member states are considered for the lithospheric mantle beneath Madagascar. First, we calculate the elevation of a thermally equilibrated column of continental lithosphere ( Figure 3c). Secondly, we assume that the lithosphere has been instantaneously thinned (i.e., delaminated) but that the plate has not been reheated ( Figure 3d). In the latter case, the lithosphere is therefore both cooler and denser, which gives rise to a lower than expected elevation compared with thermally equilibrated lithosphere of a given thickness (M c Nab et al., 2018). Lithospheric thickness estimates, a 1 , beneath Madagascar are subject to considerable uncertainty. For example, Rocco et al. (2017) infer a local maximum lithospheric thickness of 60 km based upon the absence of garnet-bearing peridotitic xenoliths from magmatic rocks in northern Madagascar. This value is ∼30 km thinner than that calculated by Priestley and M c Kenzie (2013), but in closer agreement with the results of Rajaonarison et al. (2020) and Hoggard et al. (2020), who calculate thicknesses of <80 km beneath the northern and central highlands. Note that these three studies exploited different surface wave tomographic models to obtain lithospheric thickness estimates. In Figures 3c and 3d, we show the expected elevation for a lithospheric column that is between 50 and 150 km thick for both end-member thermal states. If the lithosphere is in thermal equilibrium, a plate that is 80-120 km thick is required to bound the observed elevations. This range of thicknesses lies at the upper end of independent lithospheric thickness estimates Priestley & M c Kenzie, 2013;Rajaonarison et al., 2020;Rocco et al., 2017, this study). If Malagasy lithosphere is thermally disequilibrated following, say, rapid loss of 50 km of lithospheric mantle from the base of the plate then, for the same range of lithospheric thicknesses, topography would be elevated by several hundred meters with respect to equilibrated lithosphere. In this case, a thinner plate of about 60-90 km is required to bound the observed elevations (Figures 3d and S1). We point out that these results are highly dependent upon the density structure of the reference continental column, the density of lithospheric mantle, and uncertainties in the crustal density correction. The effects of depletion can reduce density by up to 60 kg m −3 , which in turn yields greater predicted elevation (e.g., Jordan, 1978;Lee, 2003). Thus increasing depletion progressively increases the lithospheric mantle thickness required to fit the observations (e.g., M c Nab et al., 2018). The presence of localized residual melt can further reduce lithospheric density. Nevertheless, our results suggest that some component of topographic support of the Malagasy Plateaux exists beneath the crust. We note that a significant consequence of thinning the lithospheric plate by removal of mantle material from its base is to generate regional rock uplift since the cold lithospheric mantle root is replaced by warm asthenospheric mantle.

Admittance Analysis
Madagascar is characterized by a dramatic long-wavelength (i.e., 800-1300 km) free-air gravity anomaly of +40 mGal, which implies that its overall topography is not exclusively supported by changes in crustal thickness and density. Notwithstanding edge effects, it is useful to investigate the spectral relationship between the free-air gravity and topography fields to place constraints on the relative significance of flexural and mantle dynamical support ( Figure 4). The admittance, Z(k), is the ratio between coherent components of the gravity and topography fields as a function of wavenumber, k = 2π/λ, where λ is wavelength. We calculate Z(k) for the SRTM30_PLUS digital elevation and EIGEN-6C free-air gravity models using a two-dimensional multitaper approach for a 600 × 1500 km 2 box (Becker et al., 2009;Förste et al., 2012; (Figures 4d and 4f). This value is broadly consistent with independent studies of free-air admittance and Bouguer coherence, which yield values of T e <20 km (e.g., Audet, 2014;Rakotondraompiana et al., 1999).   (Becker et al., 2009). Box = region of admittance analysis. (b) Free-air gravity anomalies taken from EIGEN-6C database (Förste et al., 2012). Box as in panel (a). (c) Observed and calculated admittance as function of wavenumber, k. Open circles with error bars = observed values of admittance ±1σ, calculated from spectral analysis of panels (a) and (b), used to estimate elastic thickness; solid circles = observed values at longest wavelengths not used to estimate elastic thickness; black line = best-fitting flexural model with The length scale over which topographic loads are supported by flexural isostatic processes is controlled by the value of T e . The half-wavelength (i.e., the distance from the center of a load to the crest of the flexural bulge), and E = 70 GPa is Young's modulus, g = 9.81 m s −2 is acceleration due to gravity, σ = 0.25 is Poisson's ratio, ρ a = 3.18 Mg m −3 is the density of the asthenospheric mantle, and ρ s = 2.0 Mg m −3 is density of the infilling material. If topography is air-loaded, we obtain x km for the calculated range of T e . Since this value is smaller than the length scale of the Antananarivo and Maromokotro highlands, it is reasonable to infer that Malagasy topography is not maintained by plate flexure. At long wavelengths, Z → 40 ± 10 mGal km −1 which cannot easily be accounted for by flexure. Instead, this value probably manifests deeper-seated density variations such as those generated by mantle convective processes and/or by removal of lithospheric mantle (e.g., Ball et al., 2019;Colli et al., 2016;Ebinger et al., 1989;M c Kenzie, 2010). Given the long-wavelength free-air gravity anomaly over Madagascar is about 40 mGal, then this value of Z implies that ∼1.0 ± 0.25 km of topography is supported by subplate processes. We emphasize that the relationship between free-air gravity and subplate support is necessarily ambiguous since the depth, size, and contrast of a putative anomaly cannot be uniquely determined (e.g., Colli et al., 2016;Molnar et al., 2015).

Apatite Fission-Track and (U-Th-Sm)/He Measurements
Constraining the history of Malagasy uplift is challenging within the central and northern highlands where stratigraphic constraints are largely absent. However, regional uplift commonly triggers denudation, the result of which is indirectly manifest in the geologic record. Thermochronologic techniques provide a useful way for investigating the spatial and temporal evolution of denudation. Apatite fission-track analysis is a well-known radiometric dating approach that exploits spontaneous nuclear fission of 238 U (e.g., Wagner, 1968). Charged fission products ejected from the nucleus create a damage trail within the crystal lattice which becomes visible when etched by acid. These tracks anneal (i.e., heal and shorten) as a function of temperature and time such that they eventually disappear (e.g., Green et al., 1985Green et al., , 1986Ketcham et al., 1999;Laslett et al., 1987). The temperature range over which the annealing rate significantly varies is called the partial annealing zone and occurs between ∼110 and 60 °C (e.g., Carlson, 1990;Gleadow et al., 1986;Ketcham et al., 2007Ketcham et al., , 1999Laslett et al., 1987). The complementary apatite helium system is a radiometric dating technique predicated upon alpha decay of U, Th, and Sm. He is lost geologically instantaneously from the crystal by volume diffusion at temperatures greater than ∼70 °C. It is retained over long periods of time at temperatures lower than 40-50 °C. Over this temperature range, which is known as the He partial retention zone, He is lost at a rate that increases with temperature in accordance with the Arrhenius relationship (Farley, 2000;Wolf et al., 1998). The sensitivity of both apatite fission-track and He systems to temperatures <110 °C means that these thermochronometers are a useful means for probing the thermal history of the upper crust.
We collected and analyzed samples along two transects that cross the edge of the Antananarivo Plateau and the Maromokotro Highlands, where erosional histories across escarpments are likely to vary as a function of elevation ( Figure 5). Fission-track ages range from 111 ± 20 to 280 ± 20 Ma with mean track lengths STEPHENSON ET AL.
10.1029/2020GC009624 9 of 38  (Emmel et al., 2006;Seward et al., 2004); black circles = locations of other published samples; two black boxes = loci of central and northern transects, respectively (see Figures 6 and 7). of 11.02 ± 0.45 to 13.78 ± 0.21 μm (Table 2). A smaller selection of apatite He analyses were also carried out at the London Geochronology Center. The resultant single grain He ages range from 10.5 to 586.6 Ma and there is significant variation both within and between samples (Supporting Information Database). A complete description of the analytical protocol followed for fission-track and He analyses is provided in Section S2.

Inverse Modeling
The existence of a crustal geothermal gradient means that the thermal evolution of a sample can, in the absence of other processes, be related to burial and denudation. Here, we are primarily interested in recovering the most parsimonious thermal history that best fits the thermochronologic observations obtained for each sample. We have tackled this problem using a threefold strategy that exploits the QTQt and HeFTy modeling packages of Gallagher (2012) and Ketcham (2005), respectively. First, we invert fission-track observations to determine a thermal history for each sample. Secondly, we jointly invert fission-track and apatite He observations using only He ages that are consistent with fission-track data within the same sample. Thirdly, we carry out a series of guided forward modeling tests with a view to investigating alternative thermal histories. Results obtained using the QTQt package are presented here and those obtained using STEPHENSON ET AL.  The QTQt package consists of a transdimensional Bayesian inverse algorithm that calculates thermal histories from thermochronologic observations (Gallagher, 2012). Fission-track length is calculated using the multikinetic annealing model described by Ketcham et al. (2007). A Bayesian framework seeks a range of models within a posterior distribution, given both prior assumptions about model parameters and the likelihood (i.e., a measure of the fit to observational constraints). We exploit the log-likelihood measure of fit described by Gallagher (1995). A uniform prior distribution across temperature-time space between 400 and 0 Ma and between 0 °C and 140 °C is assumed for all models. A Markov chain Monte Carlo (i.e., MCMC) approach is taken, whereby a new model is proposed at each iteration based upon perturbing the current model. This new model is then either accepted or rejected based upon the difference in posterior probability between it and the current model (Gallagher, 2012;Gallagher et al., 2009;Sambridge et al., 2006). The transdimensional aspect of the modeling approach permits the number of parameters (i.e., the number of time-temperature nodes) to be constrained by the observations. This scheme ensures that simpler models are favored, unless adding an extra dimension markedly improves the model fit between observation and prediction (e.g., Gallagher, 2012;Sambridge et al., 2006).
For each sample, 5 × 10 5 iterations of the MCMC inverse algorithm were carried out, following an initial "burn-in" of 10 5 iterations that are discarded from the posterior distribution. This "burn-in" ensures that the search algorithm adequately explores the T-t model space and converges on that part of the misfit space which best fits the observational constraints. The algorithm yields a suite of thermal histories, which each have an associated posterior probability. According to Bayes' theorem, this probability is proportional to the product of the prior and likelihood functions. From this suite of thermal histories, we extract a mean model that is weighted by the posterior probability of each model. Furthermore, credible intervals are calculated which bound that portion of the model space within which 95% of acceptable models reside. These intervals effectively represent the uncertainty in the calculated mean thermal history: where the thermal history is closely constrained by observations, credible intervals tend to coalesce and the mean model closely reflects the suite of accepted models in the posterior distribution; where the thermal history is less well-defined, credible intervals tend to diverge and the mean model smooths the complexity of individually acceptable models (Gallagher, 2012). In this way, the MCMC inverse algorithm identifies the smoothest thermal histories that best fit observed fission-track observations.

Thermal Histories
The results of inverting fission-track measurements from the central and northern transects are presented in Figures 6 and 7. Four results from the central transect demonstrate that both apparent ages and distributions of fission-track lengths have been adequately matched within uncertainty. They suggest that two significant episodes of cooling occur. The first episode occurs between ∼275 and 200 Ma and is apparent, although poorly resolved, in samples 2-4 (Figures 6c-6e). The second episode occurs between ∼30 Ma and the present day and is clearly manifest by sample 2 from the central transect, which is located at the foot of the erosional escarpment ( Figure 6c). Significantly, both cooling episodes are apparent where the HeFTy algorithm is used, although the older episode is better resolved (Figures S2-S4). The more densely sampled northern transect yields similar results ( Figure 7). Here, the second episode is clearly manifest by samples 7b and 11 (Figures 7c and 7d).
Consideration of all samples from both transects suggest that, despite some discrepancies, there are broadly two categories of results ( Figures S2-S4). Thermal histories from higher elevations show broadly monotonic cooling since Triassic times (e.g., samples 4, 5, 13, and 16). In contrast, histories from locations at or beneath the erosional escarpments show a rapid increase in cooling rate after ∼30 Ma (e.g., samples 1, 2, 7b, 10, and 11). A smaller number of samples have hybrid thermal histories whereby either monotonic or twophase cooling is permitted (e.g., samples 3, 8, and 15). These inferences are corroborated by independently calculating thermal histories with the HeFTy algorithm, although this approach tends to favor the two-step cooling history ( Figures S2-S4).
Notwithstanding the important role that inverse modeling plays in identifying smooth thermal histories that best fit fission-track observations, it is instructive to explicitly test which features of the recovered thermal histories are specifically required to match fission-track observations and to investigate the extent to which the observational constraints are blind to thermal complexity (Fox & Carter, 2020;Green & Duddy, 2020). We have used two complementary approaches. In the initial approach, we have run a series of guided forward models with a view to investigating two aspects of recovered thermal histories. First, we test the requirement for an episode of Neogene cooling (Figure 8). Our results demonstrate that omission of a discrete Neogene cooling event significantly degrades the fit between observed and calculated ages and track length distributions (Figures 8c and 8f). Continuous (i.e., monotonic) cooling from 300 Ma until the present day also degrades this fit (Figures 8b and 8e). Secondly, we test whether or not the inclusion of a STEPHENSON ET AL.
10.1029/2020GC009624 12 of 38    we invert these synthetic data to identify which features of a particular thermal history are resolvable. Significantly, we find that cooling episodes which occur within the last 50 Ma are likely to be well resolved. Older (i.e., Mesozoic) episodes are poorly resolved and intermediate (i.e., Cretaceous) episodes can only be resolved if their amplitude exceeds ∼20 °C ( Figure S7).
Apatite He measurements were carried out on a subset of samples, which means that joint inverse modeling of fission-track and He measurements can be undertaken ( Figure 9; Supporting Information). In principal, lower closure temperature of the complementary apatite He system should enable the younger portions of thermal histories to be better constrained even though He ages can, in practice, be dispersed and difficult to fit by joint inverse modeling (e.g., Brown et al., 2013). We acknowledge that, even after screening for broken grains, the apatite He data are clearly dispersed, which limits the interpretability of some of our modeling results. Nonetheless, our results show that, for a subset of samples, fission-track and He ages are mostly in agreement, bolstering the requirement for a discrete and rapid episode of Neogene to present-day cooling along both transects (Figure 9). There are, however, three results that must be treated with considerable caution. Figures 9a-9f show that joint inverse modeling of samples 2 and 3 yields a significant residual misfit between observed and calculated track length distributions. These failures combined with evidence for significant dispersion of He ages, especially for sample 10, indicate that the calculated thermal history shown in Figure 9s is spurious (Figure 9t).
In summary, the recovered thermal histories suggest that two episodes of cooling affected central and north- more recent exhumation occurs around northwestern and western coastlines. Note that many of our thermal histories are poorly constrained during Mesozoic times and that they could be consistent with reheating until Cenozoic times, followed by cooling to surface temperatures during the Neogene period (Figures 6c,  6d, 7c, 7d, S6a, and S6b). Such a history is consistent with progressive burial beneath sedimentary and volcanic rocks that fill the Mahajanga basin. This basin started to form in Late Permian to Early Triassic times. It continued to subside throughout the Mesozoic Era (e.g., Razafindrazaka et al., 1999). Thermal histories differ in the eastern central highlands, where more rapid cooling appears to have occurred (Jöns et al., 2009). Finally, Emmel et al. (2012) report that a suite of AFT and He analyses from the central highlands did not undergo Neogene cooling, which is consistent with the results presented here. This inference is consistent with evidence for low highland erosion rates obtained from catchment-averaged 10 Be exposure dating (Cox et al., 2009).

Uplift and Denudation
A straightforward way to cool a parcel of rock sitting within the upper crust is to remove its overburden and exhume it along the geothermal gradient. If so, calculated thermal histories at lower elevations suggest that a period of significant denudation started at ∼30 Ma. This denudation may be a consequence of erosional retreat of escarpments that occurs when drainage networks are externally forced by an episode of regional uplift, which could have been triggered in a number of different ways, including by tectonic or mantle dynamical processes. This forcing causes a kinematic wave of fluvial erosion to sweep through the landscape as the steepened portions of river channels retreat. The magnitude of denudation, D, required to generate the observed cooling histories can be gauged by assuming a geothermal gradient of 25 ± 5 °C km −1 . For example, the thermal history of sample 7b is consistent with a Neogene cooling episode of up to ∼50 °C, which yields (Figure 7). We note that a cooling episode with a similar amplitude is supported by inverse modeling of a number of samples ( Figures S3-S5). Assuming that sample 16, which is at a higher elevation, experienced a background cooling rate, the differential cooling between these two samples is 40 °C, which yields a value of differential denudation, This value is about equal to the height of the escarpment. If there has been no isostatic compensation in response to this denudation, then we can take this value of ΔD to be equal to the amount of regional uplift, U. Alternatively, we can assume that denudation results as a consequence of both regional uplift and isostatic amplification. In this case, regional uplift generated by tectonic or mantle processes is less than both the amount of denudation and the total amount of rock uplift (England & Molnar, 1990). If the total cooling is 50 °C, the amount of regional uplift required to generate denudation is now given by U = (1 − ρ o /ρ a )D where ρ a = 3.18 Mg m −3 is the density of the asthenosphere and ρ o = 2.4 Mg m −3 is the average density of the overburden (Brodie & White, 1995). This relationship yields regional uplift, U ≈ 0.5 ± 0.1 km. Given that some degree of isostatic amplification of denudation is likely to have occurred over the last 30 Ma, this lower estimate of the magnitude of regional uplift may be closer to the true value. Nonetheless, the end-member range of U ≈ 0.5-1.6 km is consistent with offshore estimates of positive residual depth anomalies and with the elevation of Eocene marine limestones ( Figure 2a; Table 1). Significantly, there is independent evidence for an erosional event that is manifest by a switch from carbonate to clastic deposition along the western margin at ∼30 Ma (Delaunay, 2018). This event coincides with our inferred increase in denudation. An increase in clastic flux between Neogene times and the present day is also recorded on borehole-calibrated seismic reflection profiles on the same margin (Delaunay, 2018

Malagasy River Networks
Tangible evidence for regional uplift of the Malagasy landscape is provided by the scattered distribution of elevated Upper Cretaceous and Cenozoic marine strata. Nevertheless, the actual spatial and temporal history of uplift is more difficult to constrain. Since it is likely that fluvial landscapes are conditioned by tectonic forcing, patterns of regional rock uplift should be recorded by the geometry of longitudinal river profiles (i.e., elevation as a function of distance along a river channel). When landscapes reach equilibrium, such that the rate of rock uplift is matched by the rate of vertical incision, river profiles are expected to achieve smooth concave-upward geometries (e.g., Whipple & Tucker, 1999). When the rate of uplift exceeds the rate of vertical incision, river profiles are characterized by long-wavelength convexities (i.e., knickzones ;Perron & Royden, 2013;Whittaker et al., 2007). Fluvial disequilibrium has been successfully exploited to extract histories of regional rock uplift (Conway-Jones et al., 2019;Fox et al., 2014;Gallen, 2018;Glotzbach, 2015;Goren et al., 2014;M c Nab et al., 2018;Pritchard et al., 2009;Richards et al., 2016;Roberts & White, 2010;Roberts, Paul, et al., 2012;Rodriguez Tribaldos et al., 2017;Rudge et al., 2015;Stephenson et al., 2014). Here, we apply this inverse modeling strategy to the drainage network of Madagascar.
First, a network of 1833 rivers is extracted from the 1-arc-second SRTM digital elevation model (Becker et al., 2009). Anomalous spikes and holes are excised and filled to avoid extracting spurious drainage planforms. Secondly, a suite of flow-routing algorithms from the ArcGIS software package is used to extract complete river profiles. A minimum drainage area threshold for overland flow of 64 km 2 helps to automatically locate river sources and drainage divides. The fidelity of the recovered drainage network was checked using Landsat satellite imagery and any extracted rivers that deviated significantly from observed channels were excised. For example, a minor catchment from northern Madagascar was erroneously routed into a larger catchment and so it was removed, and upstream drainage area of tributaries from the larger catchment were corrected. Figure 10 shows recovered river profiles from four typical catchments. These profiles are evidently disequilibrated since they display steep knickzones with length scales of 200-300 km and heights of 200-800 m. They also have numerous knickpoints (i.e., short wavelength steps and discontinuities) that are probably caused by a combination of lithologic variation, minor faulting, and man-made dam construction. Unlike knickzones, these knickpoints are unlikely to have been generated by long-wavelength rock uplift. Roberts, Paul, et al. (2012) showed that changes in bedrock lithology of Madagascar generally do not coincide with significant knickzones at wavelengths greater than ∼2 km. Kirby et al. (2003) came to a similar conclusion by analyzing Tibetan drainage patterns.

Inverse Modeling Strategy
It is generally recognized that the geometry of a river channel is controlled by the competing influence of the rate of rock uplift, U, and the rate of erosion, E, along the channel such that where x is the distance along the river from the mouth and t is time before present day. At large scales, erosion can be parametrized using a simplified version of the empirical stream-power law where A(x) is the measured upstream drainage area. v, m, and n are erosional parameters whose values can be determined by a combination of optimization and independent calibration (e.g., Perron & Royden, 2013;Roberts & White, 2010;Whipple & Tucker, 1999). v is equivalent to the erodibility coefficient of geomorphic literature, K. The value of v is poorly known as a function of space and time but it is probably controlled by a combination of factors that include climate, precipitation and bedrock erodability (e.g., Lague, 2014;Whipple, 2004). Coefficients m and n probably depend upon drainage basin geometry and hydrology (Howard, 1994;Whipple & Tucker, 1999). Inverse modeling demonstrates that significant inventories of river profiles can be optimally fitted if n = 1 (Paul et al., 2014;Roberts, Paul, et al., 2012;Rudge et al., 2015;Stephenson et al., 2014). We acknowledge that this result does not necessarily mean that n = 1 on timescales and length scales that are significantly shorter than those which we consider (e.g. Harel et al., 2016). Rudge et al. (2015) exploit the method of characteristics to pose and solve the inverse problem for Equation 6. Here, we apply this algorithm to calculate the smoothest history of regional rock uplift that best fits an inventory of observed river profiles. The solution can be written in the form z(x, t) = z(x(t), t) since which can be expressed as two ordinary differential equations: where the boundary conditions are The boundary conditions of Equation 10 refer to the present-day river profile, which has a height of z* at a distance x* along the channel measured from its mouth (Rudge et al., 2015). The landscape response time or Gilbert Time, τ G , is the time taken for a kinematic wave to propagate upstream from the mouth to any point x on the profile (e.g., Weissel & Seidl, 1998). The boundary condition of Equation 11 represents a time in the past when the characteristic curve intersects the mouth of the river (i.e., x = 0, z = 0). τ G is given by and represents the maximum amount of time over which uplift can be recorded along a given river. Solutions of Equations 8-12 take the form These equations govern the evolution of a single river profile. Here, the goal is to simultaneously invert hundreds of river profiles to determine the temporal and spatial pattern of regional uplift (e.g., Fox et al., 2014;Roberts & White, 2010;Rudge et al., 2015). Uplift is specified at a set of spatial and temporal nodes within a finite mesh such that uplift is free to vary in both spatial dimensions as a function of time. The value of τ G is determined by integration using the known values of A and x*, which means that it is straightforward to calculate characteristic curves from Equation 13. Integration of uplift rate at each point, x(t), along the river using Equation 14 yields the height, z, of the point, x, as a function of time, t. Following discretization, this approach can be written down in matrix form such that where z is the elevation at each point along any one river profile, M is the model matrix, and U is the uplift history at each spatial and temporal node (Rudge et al., 2015). M contains the information required to translate an uplift history into river profiles at each temporal node. Required information includes the x and y coordinates of river profiles, upstream drainage area, A, and erosional parameters v and m. U is calculated by inverting Equation 15. A damped nonnegative least squares scheme is used to minimize The nonnegativity constraint ensures that oscillations between positive and negative values of U are suppressed (Roberts & White, 2010;Rudge et al., 2015). λ S and λ T are spatial and temporal smoothing parameters that control regularization of the problem by penalizing high values of gradient for spatial and temporal components of the uplift history. Spatial and temporal roughness are penalized using and In this way, many hundreds of river profiles are simultaneously inverted to yield a smooth history of regional uplift (

Regional Uplift History
This optimization approach is used to calculate a regional rock uplift history of Madagascar from the inventory of 1833 river profiles. This analysis augments that presented by Roberts, Paul, et al. (2012), who inverted 98 rivers using an older, nonlinear, numerical scheme. The fit between observed and calculated river profiles is shown for selected catchments in Figure 10. This fit is predicated upon the cumulative rock uplift history shown in Figure 11. To obtain these results, it is necessary to determine values of both the erosional and smoothing parameters: m and v; λ S and λ T . It is important to emphasize that our inverse strategy makes two assumptions. First, the topology of the drainage network remains fixed through time which means that we do not account for drainage capture or divide migration. Drainage simulation tests show that this assumption is probably a reasonable one since stream power depends upon a small fractional power of upstream drainage area (O'Malley et al., 2020;Paul et al., 2014;Stephenson et al., 2014;Wilson et al., 2014). However, we acknowledge that drainage migration and capture may have occurred in parts of the island that are experiencing extension (Kusky et al., 2010;Seward et al., 2004). Secondly, we initially assume that v does not change as a function of time and space.
Values of m, λ S , and λ T influence both the misfit between observed and calculated river profiles as well as model roughness (Figure 12). Here, this parameter space is explored to locate optimal values of these three parameters that yield the smoothest model which best fits observed river profiles in accordance with the approach of Parker (1977). In general, significant trade-offs exist between m, λ S , and λ T (M c Nab et al., 2018;Richards et al., 2016). A parameter sweep is carried out where m, λ S , and λ T are varied between 0.1 and 0.75, 10 −3 and 10 3 , and 10 −3 and 10 3 , respectively (Figures 12a and 12b). A weak global minimum misfit occurs at m = 0.35. This value is consistent with independent studies which have shown that 0.35 ≤ m ≤ 0.65 if n = 1 (e.g., Kirby & Whipple, 2012;Perron & Royden, 2013;Stock & Montgomery, 1999;Whipple & Tucker, 1999). Values of λ S and λ T are less well determined since multiple combinations yield good fits. Nevertheless, since we are primarily interested in the smoothest model that yields an optimal match between observed and calculated river profiles, λ S = 2 and λ T = 1 are appropriate ( Figure 12). Note that this form of calibration STEPHENSON ET AL. was not carried out by Roberts, Paul, et al. (2012) owing to the punitive computational expense of their optimization scheme at that time.
The value of v sets the pace of advective erosion across the drainage network and thus constrains the timescale for regional uplift. It is important to note that v does not affect the misfit between observed and calculated river profiles. The timescale is determined by comparing a range of calculated uplift histories to independent rock uplift constraints (Figure 12c). For example, scattered outcrops of elevated Eocene limestones provide convincing evidence for regional uplift (Besairie, 1964). We find that there is a reasonable correlation between observed and calculated uplift rates if v ≈ 4.0 ± 1 m 0.3 Ma −1 . A significant consequence of this value is that τ G ≈ 90 Ma, which implies that the drainage network is principally sensitive to tectonic forcing on Cenozoic timescales.
The calculated cumulative uplift history of Madagascar suggests that the modern landscape began to develop at ∼30 Ma, which broadly agrees with previous results (Roberts, Paul, et al., 2012; Figure 11). Regional uplift initiates in the southern half of the island. By ∼20 Ma, significant uplift has commenced at the northernmost tip of the island. By 10 Ma, the entire landscape has become emergent with the possible exception of the north western seaboard. The average rate of uplift between ∼30 and 20 Ma is ∼0.01 mm yr −1 . It is noticeable that during the last 10 Ma, the rate of uplift increases to almost 0.1 mm yr −1 within the central Antananarivo Plateau and 0.05 mm yr −1 within the Maromokotro Highlands. In the south, uplift does not see such an acceleration in Miocene-Recent times, remaining at ∼0.01-0.02 mm yr −1 . We infer that the interior of the Antananarivo Plateau and the Maromokotro Highlands have been uplifted by about 0.8-1.2 km since ∼30 Ma, while the southern Horombé Plateau experienced only ∼0.6 km of uplift during this time.
The value of τ G ≈ 90 Ma implies that Malagasy drainage can potentially record regional uplift over the Cenozoic Era. In practice, however, this value is likely to be an overestimate since it assumes that uplift is always inserted at the mouths of river channels which is probably incorrect. The spatial resolving power of the model with respect to the changes of uplift rate at any given node within the mesh as a function of time is more accurately gauged by summing the number of non-zero values within the model matrix, M (Rudge et al., 2015). Figure 11 indicates that spatial resolution is initially low and that only uplift signals inserted adjacent to river mouths are preserved within present-day longitudinal profiles. The highest resolution is obtained over the last 10-15 Ma as model coverage extends across most of the island.

Mafic Magmatism
Minor volumes of intraplate Neogene magmatism occur through Madagascar and on adjacent islands (Figure 2d;Cucciniello et al., 2011Cucciniello et al., , , 2017Cucciniello et al., , 2018Emerick & Duncan, 1982;Melluso & Morra, 2000;Melluso et al., 2016). There are two significant provinces: the Ankaratra and the neighboring Itasy volcanic fields located on the Antananarivo Plateau; and the northern volcanic zone, which comprises Nosy Bé, the Ampasindava Peninsula, Ankaizina, Amparihy, Massif d'Ambre, and the Babaomby region at the northernmost tip. Small amounts of volcanism are recorded elsewhere, notably throughout the northern part of the island, along the east coast and at Ankililoaka in the southwest corner (Cucciniello et al., 2018). Offshore, significant volcanism occurs on Réunion, on Mauritius and on the Comoro islands, which include Grand Comore, Mohéli, Anjouan and Mayotte.
A database of 1,049 samples has been assembled that includes 1,039 published samples and 10 samples which were collected in June 2012 and October 2015 from Nosy Be, Massif d'Ambre, and from the Babaomby peninsula (see Supporting Information Data). New samples were analyzed for major, trace, and rare earth elemental (REE) compositions using inductively coupled plasma-mass spectrometry (ICP-MS) and X-ray fluorescence (XRF) at the universities of Cambridge and Edinburgh, respectively. For modeling purposes, we have divided the database into six geographic provinces: the Comores, Massif d'Ambre, Nosy Bé, Ankaratra, Ankililoaka, and the Mascarene Islands ( Figure 13). The bulk of samples under consideration consist of basalts and basanites with minor numbers of low alkali foidites and picritic basalts. Samples from the Comores have the lowest SiO 2 content. To account for samples that have undergone either high levels of fractionation or lithospheric contamination, analyses were screened to ensure 8.5 < MgO < 14.5 wt% (see Supporting Information Data). High MgO values are consistent with low degrees of fractionation since olivine crystallization removes MgO from the melt. However, MgO > 14.5 wt% could indicate assimilation of peridotitic xenoliths.

Modeling Strategy
We wish to explore the quantitative relationship between basaltic magmatism, dynamic topography, and mantle convective processes that occur beneath the lithospheric plate. Our starting point is the acknowledgment that the combination of tectonic setting and observed isotopic signatures rule out a subduction-related explanation and instead favor an intraplate source (e.g., Cucciniello et al., 2011;Melluso et al., 2016). It is possible that decompression melting of a subplate source could be driven by anomalously high asthenospheric temperature, by passive upwelling caused by lithospheric thinning, by mantle flow, or by some combination of all three mechanisms. It is also possible that one of these mechanisms could mobilize metasomatic material within the lithospheric mantle. Here, we explore these different possible mechanisms by applying two independent geochemical approaches that enable the temperature and pressure of asthenospheric melting to be calculated from major and trace element observations. The first approach exploits an empirical thermobarometric scheme developed by Plank and Forsyth (2016) which determines the pressure and temperature of melt equilibration from major element measurements. The second approach calculates melt fraction as a function of depth from REE concentrations using an inverse modeling strategy developed by M c Kenzie and O'Nions (1991) that is predicated upon the way by which partition coefficients control the composition of polybaric melting beneath the plate.
Since we are primarily interested in gauging whether the subplate temperature beneath Madagascar differs from that of the ambient convecting mantle, it is useful to choose a reference value for mantle potential temperature, T p . The value of T p used in different studies varies (e.g., Jennings & Holland, 2015;Katz et al., 2003;M c Kenzie & Bickle, 1988). Here, we adopt the melting model of Katz et al. (2003)

Major Element Thermobarometry
Partitioning of major elements between melt and solid is pressure-and temperature-dependent (Lee et al., 2009). Building on the work of Lee et al. (2009), by compiling a comprehensive suite of melt equilibration observations from laboratory experiments, Plank and Forsyth (2016) parameterized major element compositions as a function of pressure and temperature. This approach is predicated upon carefully correcting the concentration of each major element for olivine fractionation to determine the primary melt composition. This procedure is carried out by adding in olivine, that is in equilibrium with the melt composition, until the melt reaches equilibrium with mantle rock of an assumed forsterite content. Here, the approach of M c Nab et al. (2018)  Following fractionation correction, which yields the composition of the primary melt, temperature, T, and pressure, P, are calculated using the following pair of empirical equations: The change in T that results from water in the melt is given by  (Lee et al., 2009;Plank & Forsyth, 2016).
Temperature estimates are dependent upon the value of Fe 3+ /ΣFe and upon sample water content. Following Plank and Forsyth (2016) and M c Nab et al. (2018), Fe 3+ /ΣFe is assumed to be 0.2, which is halfway between the values for mid-ocean ridges and island arcs. This assumption is in agreement with geochemical observations from western North America (Plank & Forsyth, 2016 (Dixon et al., 2002). Note that higher values of H 2 O/Ce will yield lower temperatures (M c Nab et al., 2018).
The strategy of M c Nab et al. (2018), who adapted the approach of Plank and Forsyth (2016), is used to determine values of T p from equilibration temperatures and pressures. Melting paths are fitted to thermobarometric estimates for all samples from a given province. First, we assume that the melting region is dry so that melting paths can be taken from the anhydrous parametrization of Katz et al. (2003). A best-fitting melt path is obtained by calculating a combined least-squares misfit value between these melt paths and the suite of equilibrated thermobarometric estimates. Observed and calculated values of P and T are first normalized by their uncertainties (i.e., ±39 °C and ±0.24 GPa, respectively). Next, a least-squares misfit is calculated for both P and T before both estimates are summed to yield a joint misfit value (Ball et al., 2019;M c Nab et al., 2018). The locus of intersection between a best-fitting melt path and the anhydrous solidus corresponds to an adiabatic decompression gradient, which also intersects the solidus in the same place. A parcel of melt that decompresses along this adiabatic gradient will evolve and equilibrate along the best-fitting melt path. The value of T p is calculated by extrapolating this adiabatic gradient to the surface to obtain the temperature. The uncertainty for T p is given by values that yield a root-mean-squared misfit which is double the value at the global minimum. Samples that plot beneath the solidus are not included in the misfit calculation. In practice, this approach assumes that each sample represents a batch melt that equilibrated at a given pressure and temperature before extraction without any further change. More than half of the thermobarometric equilibration estimates from Ankaratra and Ankililoaka lie beneath the solidus curve. In contrast, all of the estimates from Nosy Bé fall beneath this curve. This difference can partly be explained if water within the source region is incorporated into the melting parametrization. The presence of water increases the depth (i.e., pressure) of the solidus as a function of temperature, which means melting starts at greater depths (e.g., Katz et al., 2003). The effect of adding 0.01 wt% water to the source region for each province is shown in Figure S8. For a wetter source, the recalculated solidus curve shows that the majority of samples plot within the melting region. Pressure and temperature estimates for Nosy Bé also fall within the melting field, corresponding to a potential temperature of path if a dry source region is assumed. It is important to state that the presence of CO 2 in the source region would have a similar, but larger effect. We note that the potential temperature estimates calculated using this empirical thermobarometric approach probably represent upper bounds since values obtained for mid-oceanic ridge basalts are significantly hotter than the expected 1330 °C, by perhaps as much as 70 °C (Ball et al., 2019;Lee et al., 2009). Regardless of the source, the shallowest equilibration depths are between 50 and 60 km.

Inverse Modeling of Rare Earth Element Concentrations
Modeling of REE concentrations can be used to constrain the degree and depth of melting of mantle rocks. REEs are incompatible during peridotitic melting and they preferentially partition into the melt phase, becoming increasingly diluted by other elements as melting proceeds. The exact mineralogy of the source region plays a major role in controlling REE compatibility. Over a particular depth range within the mantle source, the aluminous phase changes from spinel to garnet. When garnet is present in the source, heavy REEs partition much less readily into the melt phase than in the presence of spinel since heavier REEs are preferentially incorporated into garnet. The spinel-garnet transition is primarily pressure-dependent and so relative depletion of heavy REEs within a melt is indicative of the depth of melting. This systematic behavior is exploited by the inverse modeling strategy developed by M c Kenzie and O'Nions (1991).
The INVMEL(v12.0) algorithm attempts to match the REE composition of primitive mafic igneous rocks by varying melt fraction as a function of depth. Melt fraction is discretized at regular depth intervals. At each interval, the partition of REEs into the melt phase is calculated and their total concentration is determined by integrating melt fraction as a function of depth. The least-squares misfit between observed and calculated REE concentrations is minimized by systematically varying melt fraction with depth using a standard conjugate direction search routine. Mantle potential temperature is separately calculated by comparing the recovered melt fraction as a function of depth with mantle melting models (e.g., Katz et al., 2003). The depth of the garnet-spinel transition zone is instrumental in determining the depth of onset of melting since any sample that is highly depleted in heavy REEs is likely to have undergone melting within the garnet field. Here, a spinel-garnet transition zone depth of 63-72 km is used (Jennings & Holland, 2015). Following Klöcking et al. (2018), the base of this transition zone has been slightly deepened to stabilize the inverse algorithm.
Two end-member source compositions are specified: primitive and depleted MORB mantle. For each province, an average source composition is estimated using the observed value of ϵNd. ϵNd = 0 and ϵNd = 10 correspond to primitive and depleted sources, respectively. Partition coefficients were set using the parametrization of M c Kenzie and O'Nions (1995) for olivine, orthopyroxene, and spinel. REE partition coefficients for clinopyroxene and garnet are calculated for mineral compositions listed in M c Kenzie and O'Nions (1991) using the parametrizations of Wood and Blundy (1997) and Van Westrenen et al. (2001), respectively. Melt fraction is discretized at fixed 3-km depth intervals. Melting is assumed to cease at a depth, d top , which represents the base of the lithospheric column. d top is systematically varied to identify the optimal value. Although the INVMEL algorithm minimizes the misfit between observed and calculated REE concentrations, it is useful to use the recovered melt fraction as a function of depth to predict, by forward modeling, other trace element concentrations, which acts as a useful independent constraint. Figure 15 shows the results of geochemical inverse and forward modeling. For the Mascarene Islands, a cumulative melt fraction of 6.8% is required to match REE concentrations and approximately one-third of this melting is required to occur within the garnet stability field. The best-fitting cumulative melt fraction with depth pathway almost exactly falls along the 1365 °C adiabatic curve and d top = 50 km. In general, Malagasy and Comorean samples have similar concentrations of heavy REEs but greater concentrations of light REEs, which is consistent with lower melt fractions with some degree of melting within the garnet field. Cumulative melt fractions are 3.2% for Massif d'Ambre, 1.5% for Nosy Bé, 1.8% for the Comores, 1.3% for Ankaratra, and 0.9% for Ankililoaka. Massif d'Ambre, Nosy Bé, Ankaratra, and Ankililoaka all require a small amount of melting within the garnet field, but at low temperature, which results in a non-adiabatic pathway for melt fraction with respect to depth. Nevertheless, at depths of <70 km, melt fraction more closely follows an adiabatic relationship, which corresponds to T p = 1348 °C for Massif d'Ambre, T p = 1325 °C for Nosy Bé, T p = 1330 °C the Comores, and T p = 1315 °C for Ankaratra and Ankililoaka. In other words, the recorded potential temperature is taken to be that at the top of the melt column. These results suggest that temperatures are elevated beneath the Mascarene Islands and Massif d'Ambre, but are ambient to cool beneath other provinces. Melting ceases at d top ≈ 60 km for all Malagasy provinces.
For all locations with the exception of Ankililoaka, REE concentrations are successfully matched with RMS misfits <1. Other trace element concentrations have also been broadly successfully matched. At Nosy Bé and Ankililoaka, light REE concentrations are somewhat underpredicted, possibly because small melt fractions have interacted with, or are partly sourced from within the lithospheric mantle. Indeed, authors have previously noted that Ankaratran and Ankililoakan melts are likely to have interacted with the lithosphere (e.g., Cucciniello et al., 2018;Melluso et al., 2016). Calculated concentrations of Cs, Rb, and K for Nosy Bé, the Comores, Ankaratra, and Ankililoaka are too great. These three elements have elevated solubilities, which means that their concentrations are susceptible to surficial weathering processes. The melt fraction tails required to fit observed REE concentrations at Ankaratra and Massif d'Ambre, Nosy Bé, and Ankililoaka indicate that small amounts of partial melting are required within the garnet stability field that cannot be explained by adiabatic melting. Such melting could result from the presence of minor amounts of volatiles within the source region that deepen the onset of melting but are rapidly exhausted (e.g., H 2 O or CO 2 ). Alternatively, fluids originating from within the garnet field may ascend to shallower depths and contaminate melts within the spinel stability field. Such contamination has been previously invoked to explain the composition for Malagasy Cenozoic volcanism (e.g., Cucciniello et al., 2011Cucciniello et al., , 2018Melluso et al., 2007Melluso et al., , 2016. Low melt fraction tails have been reported elsewhere and their existence is supported by major element thermobarometric calculations (Gibson & Geist, 2010;M c Nab et al., 2018).

Discussion
A range of geologic and geophysical observations have been used to investigate the timing and amplitude of regional epeirogenic uplift for a region that encompasses Madagascar. The most tangible evidence comprises emergent marine sedimentary strata that are patchily distributed along the western seaboard. These fossiliferous deposits demonstrate that ∼850 m of uplift occurred along the southern rim of the island and that ∼450 m of uplift occurred along the northern rim between Eocene times and the present day. Idealized isostatic calculations show that regional elevation is broadly consistent with crustal thickness measurements and crustal density estimates provided 50-75 km of lithospheric mantle has been removed. Admittance analysis of the topographic and gravity fields indicates that this regional elevation is unlikely to be supported by crustal flexure. At wavelengths of greater than ∼300 km, an observed admittance value of Z = +40 ± 10 mGal km −1 is indicative of subcrustal (i.e., mantle convective) support. It is significant that Neogene magmatism occurs along the length of Madagascar and in its surroundings. Inverse and forward modeling of the rare earth elemental composition of mafic igneous samples suggests that isentropic melting occurred at depths as shallow as 50 km and at temperatures which are neutral to 50 °C higher than that of ambient asthenospheric mantle.
Youthful regional uplift has evidently helped to sculpt the modern fluvial landscape which is characterized by escarpment retreat. Inverse modeling of thermochronologic measurements from a series of regional transects that cross prominent escarpments shows that a resolvable episode of cooling commenced at the foot of these escarpments in Neogene times. This cooling represents 1.6 km of denudation in response to ∼0.5-1.6 km of regional tectonic uplift. This episode of erosion is consistent with escarpment retreat caused by fluvial incision, with a switch from carbonate to clastic deposition at ∼30 Ma, and with an increase in sedimentary flux. Longitudinal river profiles that drain the Central and Northern Highlands contain prominent knickzones that indirectly record regional rock uplift. Calibrated inverse modeling of an island-wide inventory of these profiles suggests that regional uplift commenced during Oligocene times in the southern part of the island and during Miocene times in the Northern Highlands at rates of ∼0.01 mm yr −1 , before accelerating to ∼0.1 mm yr −1 during the last 10 Ma. We conclude that inferred rates of fluvial erosion and thermochronologic measurements are broadly self-consistent. Nevertheless, it has not escaped our notice that there is remarkably limited evidence for a regional cooling episode during Late Cretaceous times when significant basaltic magmatism and lithospheric rifting occurred. This conclusion is consistent with other thermochronologic observations that appear to preclude significant Late Cretaceous and Paleogene burial and/or denudation in central, northern, and western regions (e.g., sample 17; Emmel et al., 2012;Jöns et al., 2009;Seward et al., 2004). However, published fission-track analyses from the foot of the eastern escarpment do indeed have Late Cretaceous fission-track ages which perhaps hint at rapid Cretaceous cooling (Emmel et al., 2006;Seward et al., 2004). While Emmel et al. (2012) present a limited number of inverse models which require Cretaceous cooling from west of the central highlands, our results suggest that this significant period of cooling is not routinely recoverable within central, northern and western regions (Emmel et al., 2006(Emmel et al., , 2008Seward et al., 2004).
These different observations and inferences about Neogene regional uplift can be tested by using seismic tomographic models, which provide a present-day snapshot of subplate structure, to further investigate the causes of epeirogeny. Shear wave speed is sensitive to temperature, which means that seismic tomographic imaging provides an independent constraint on the thermal structure of subplate mantle. Here, we examine the results of four tomographic studies ( Figure 16).  carried out a regional study of the upper mantle beneath Africa that includes Madagascar. In their model, the average shear wave velocity between 100 and 200 km shows that subplate mantle beneath the northern portion of Madagascar is anomalously slow with respect to their reference model (Celli, Lebedev, Schaeffer, & Ravenna, 2020; Figure 16a). Faster velocities occur beneath the central part of the island and there is a band of slow shear wave velocities beneath the southwestern corner. These results are in general agreement with a detailed local study of (Pratt et al., 2017). This study suggests that there are three patches of slow shear velocity anomalies beneath northern, central, and southern Madagascar, which coincide with the distribution of Cenozoic volcanic rocks (Figure 16b). Finally, Fishwick (2010) and Mazzullo et al. (2017) each present regional tomographic models of the western Indian Ocean, which show that slow shear wave velocity anomalies principally occur beneath northern and central Madagascar and coincide with the Massif d'Ambre and Ankaratra volcanic fields (Figures 16c and 16d).
We determine the temperature structure of the subplate mantle by calibrating the AF2019 tomographic model of . To convert shear wave velocity into temperature, we exploit the thermodynamic calibration developed by Yamauchi and Takei (2016) and subsequently revised by Richards et al. (2020) and Klöcking et al. (2020), both of which build upon the empirical parametrization of Priestley and M c Kenzie (2013). It is worth noting that away from Africa and the South Atlantic Ocean, where coverage of great circle paths is denser, the AF2019 model closely resembles the global tomographic model of Schaeffer & Lebedev (2013), which Richards et al. (2020) used to parametrize a revised calibration. Figure 17 shows a north-south transect through the calibrated AF2019 model together with independent geologic observations. Beneath northern Madagascar, this transect shows that a significant thermal anomaly with a potential temperature of 1350°-1360 °C sits between 100 and 300 km depth (Figure 17e). The size of this anomaly is comparable to that determined by inverse modeling of REE concentrations (Figure 18). Beneath central and southern Madagascar, this anomaly appears to shallow with potential temperature decreasing to between 1280 °C and 1310 °C. These values are slightly cooler than those determined by inverting REE concentrations (Figure 18). Along the entire transect, the depth to the 1200 °C isothermal surface is generally less than 100 km, in broad agreement with melting depth estimates determined by major element thermobarometry and inverse modeling of REE concentrations. We note that mantle xenoliths erupted with Cenozoic basalts contain spinel but not garnet, which indicates that melts ascended through lithospheric mantle whose base is shallower than the garnet-spinel transition (Rocco et al., 2017).
In accordance with the definition of dynamic topography, it is straightforward to calculate the amount of uplift that is generated by anomalously hot subplate temperatures, by lithospheric thinning, or by some combination of both. At the northern end of the transect, excess temperatures of up to about ∼50 °C within a 150-km-thick asthenospheric channel are consistent with ∼300 m of air-loaded uplift, which accords with the height of Eocene limestone outcrops, with the amount of uplift inferred from thermochronologic measurements, and with offshore residual depth anomalies although there is a small systematic discrepancy (Figure 17a; Richards et al., 2020;Rudge et al., 2008). The temporal evolution of this prominent temperature anomaly is probably manifest by regional tilting of the shoreline since marine isotope stage 5e (i.e., the last interglacial period; Stephenson et al., 2019). At the southern end of the transect, the observed temperature anomalies are smaller and less able to account for offshore residual depth anomalies and for the presence of Cenozoic basalts.
There are a number of ways to account for this discrepancy. The simplest explanation is predicated upon geologic evidence that supports a significant change in lithospheric thickness. Stratigraphic observations demonstrate that much of Madagascar lay close to mean sea level during Late Cretaceous and Paleogene times. If the lithospheric mantle was thermally equilibrated during Paleogene times then, in the absence of subplate topographic support, isostatic considerations show that the lithosphere should have been ∼100-135 km thick for crust that is 35-40 km thick with a density of 2.85 Mg m −3 . A combination of geochemical, tomographic, and isostatic observations suggest that the present-day lithosphere is as thin as 60 km. Hence if Malagasy lithosphere was initially thermally equilibrated and, say ∼120-km thick, then instantaneous thinning of the lithosphere generates regional uplift where km depth taken from regional tomographic model of Africa . Anomalies are with respect to their version of AK135 model. Black polygons = Cenozoic volcanic rocks (Besairie, 1964). (b) Local tomographic model of Madagascar (Pratt et al., 2017). (c) Regional tomographic model of Western Indian Ocean (Mazzullo et al., 2017). (d) Regional tomographic model of Africa and Madagascar (Fishwick, 2010). Note. Separate scale bar for panel (d).
Figure 17. Regional epeirogeny of Madagascar. (a) Observed and calculated regional uplift along SSW-NNE oriented transect shown in inset of panel (c). Solid/dashed lines = air-loaded uplift and its uncertainty calculated using inferred subplate thermal anomaly shown in panel d for ambient asthenospheric potential temperature of 1330 ± 30 °C. Large black circles with error bars = air-loaded residual depth measurements and their uncertainties, where both sedimentary and crustal corrections are applied, taken from Hoggard et al. (2017) and projected onto transect as necessary; upward/downward pointing triangles = lower/upper estimates of residual depth for which only sedimentary corrections are applied (sign of crustal correction inferred from regional constraints); small black circles = air-loaded residual depth measurements determined from inventory of ship-track bathymetry that is binned every 1°; blue square = elevation of Eocene nummulitic limestone that crops out in northernmost Madagascar; green diamonds with vertical bars = regional uplift ±1σ determined by inverse modeling of apatite fission-track measurements. (b) Same as panel (a) but solid/dashed lines = air-loaded uplift and its uncertainty calculated from combined effects of subplate thermal anomaly and of instantaneously thinning lithosphere from 120 to 60 km. (c) Topography along transect x to x′ (see inset). (d) Vertical slice through earthquake tomographic model of  where red/blue colors indicate negative/positive shear wave velocity anomalies relative to their version of AK135 global reference model adjusted for variable crustal structure. Gray polygon = lithospheric plate where calculated temperature <1200 °C; dashed line = locus of 1200 °C isothermal surface; (e) Vertical slice showing thermal structure determined by converting absolute shear wave velocities into temperatures  a 0 and a 1 are the original and present-day lithospheric thicknesses, α = 3.28 × 10 −5 °C −1 is the thermal expansivity of mantle rock, and T 1 = 1370 °C is the temperature beneath the lithospheric plate (i.e., the temperature of asthenospheric mantle at 100 km given T p = 1330 °C). If the lithosphere is thinned from 120 to 60 km, U = 0.7 km. Figure 17b shows that the combined effects of anomalously hot asthenosphere and lithospheric thinning yield a more satisfactory agreement between observed and calculated regional uplift. Note that higher asthenospheric temperatures lead to slightly higher values of U in response to lithospheric thinning. This combination simultaneously accounts for stratigraphic and thermochronologic observations, for offshore residual depth measurements, and for the geochemical composition of Cenozoic basalts (Figure 19). We infer that this combination accounts for the elevation of the central Antananarivo and northern Maromokotro Highlands, where the timing and amplitude of regional uplift cannot at present be constrained (Dixey, 1960).

Conclusions
A combination of geophysical, geomorphic, and geochemical observations has been used to investigate the dynamic topography of Madagascar. Thermochronologic measurements imply that significant escarpment retreat occurred during Neogene times which is consistent with the cumulative uplift history calculated by calibrated inverse modeling of an STEPHENSON ET AL.
10.1029/2020GC009624 32 of 38 Figure 18. Comparing different temperature estimates. Relationship between seismologic temperature estimates determined from AF2019 model and geochemical estimates determined by inverse modeling of rare earth element concentrations. Blue/red/green/dark red/purple/yellow circles = Ankaratra/Mascarene Islands (i.e., Réunion and Mauritius)/ Massif d'Ambre/Nosy Bé/Comores/Ankililoaka. Gray squares = potential temperature estimates determined by major element thermobarometry. inventory of river profiles. We suggest that regional uplift commenced in southern Madagascar at 30-40 Ma and spread into the central and northern parts of the island. This history is consistent with stratigraphic observation and with the age and height of emergent marine terraces that rim the northern and southern coastlines. Geochemical modeling of major and REE concentrations for Cenozoic basalts distributed along the central spine of Madagascar and throughout adjacent oceanic basins indicates that isentropic melting of an asthenospheric mantle source occurred at ambient to elevated asthenospheric temperatures beneath a lithospheric plate as thin as 60 km. Four different seismic tomographic models of a region encompassing Madagascar confirm that slower than ambient shear wave velocity anomalies lie immediately beneath the thinned lithosphere. Calibrated anelastic parametrization of one of these models suggests that a significant thermal anomaly sits beneath northern Madagascar and that the island is characterized by thin lithosphere. Some combination of modest thermal anomalies beneath a thinned plate yields a satisfactory agreement between observed and calculated regional uplift. Finally, we acknowledge that we have not accounted for the possibility that the viscous flow of subasthenospheric mantle plays a role in generating and maintaining some proportion of dynamic topography. The significant difficulty is how the putative geometry of this flow can be observationally constrained. Notwithstanding this limitation, we argue that Neogene uplift, erosion, and magmatism of Madagascar are primarily explained by subcrustal processes.

Data Availability Statement
Crustal thickness and velocity profiles exploited in Section 2.1 are available from Andriampenomanana et al. (2017). Oceanic residual depth measurements are available from Hoggard et al. (2017). Additional residual depth measurements provided in the Supporting Information of this contribution were calculated by M. Holdt. GPS data are available from Stamps et al. (2018). Earthquake location data are available from Rakotondraibe et al. (2020) and the ISC catalog. Tomographic data are available from the authors and lithospheric thickness grids are in the Supporting Information of Hoggard et al. (2020). Thermochronologic data are available in the Supporting Information and are archived at Geochron (http://www.geochron.org/ dataset/html/geochron_dataset_2021_03_23_vYSnw). Thermochronologic modeling software is available from Richard Ketcham and Kerry Gallagher. Longitudinal river profile inverse modeling software is available from Gareth Roberts. INVMEL v12.0 is available from Dan M c Kenzie. The mafic igneous geochemical database was assembled with the aid of GEOROC and is available in the Supporting Information and is archived with EarthChem (https://doi.org/10.26022/IEDA/111941). Department of Earth Sciences contribution number esc.6028.