Thermal Structure of Eastern Australia's Upper Mantle and Its Relationship to Cenozoic Volcanic Activity and Dynamic Topography

Spatio‐temporal changes of upper mantle structure play a significant role in generating and maintaining surface topography. Although geophysical models of upper mantle structure have become increasingly refined, there is a paucity of geologic constraints with respect to its present‐day state and temporal evolution. Cenozoic intraplate volcanic rocks that crop out across eastern Australia provide a significant opportunity to quantify mantle conditions at the time of emplacement and to independently validate geophysical estimates. This volcanic activity is divided into two categories: age‐progressive provinces that are generated by the sub‐plate passage of mantle plumes and age‐independent provinces that could be generated by convective upwelling at lithospheric steps. In this study, we acquired and analyzed 78 samples from both types of provinces across Queensland. These samples were incorporated into a comprehensive database of Australian Cenozoic volcanism assembled from legacy analyses. We use geochemical modeling techniques to estimate mantle temperature and lithospheric thickness beneath each province. Our results suggest that melting occurred at depths ≤80 km across eastern Australia. Prior to, or coincident with, onset of volcanism, lithospheric thinning as well as dynamic support from shallow convective processes could have triggered uplift of the Eastern Highlands. Mantle temperatures are inferred to be ∼50–100°C hotter beneath age‐progressive provinces that demarcate passage of the Cosgrove mantle plume than beneath age‐independent provinces. Even though this plume initiated as one of the hottest recorded during Cenozoic times, it appears to have thermally waned with time. These results are consistent with xenolith thermobarometric and geophysical studies.

for provincial names). Red/white/blue contours = positive/zero/ negative free-air gravity anomalies at 10 mGal intervals from DIR-R5 database, low-pass filtered for wavelengths >800 km (Bruinsma et al., 2014). Black polygons = loci of Cenozoic volcanism (Raymond et al., 2012). (b) Shearwave velocity averaged between 100 and 200 km depth from global tomographic model together with principal age-independent volcanic provinces (see Table 1 for province names; Schaeffer & Lebedev, 2013). Shear wave velocity shown relative to a reference model (Schaeffer & Lebedev, 2013). Red/white/blue contours = negative/zero/positive ΔV s at 0.05 km s −1 intervals. (c) Location of dated intraplate volcanic samples, colored by age. Triangles = age-progressive volcanic samples; circles = age-independent samples (Data Set S1); dashed lines = flow lines constructed by tracking motion of Australian plate relative to Pacific hotspot reference frame . (d) Radiometric dates as function of latitude for Cosgrove hotspot track. Large/small colored symbols = 40 Ar-39 Ar/ 40 K-40 Ar and 14 C dates; gray symbols = Cenozoic radiometric dates for eastern Australia. (e) Same as panel D for Tasmantid hotspot track. son, 1989;Wellman & McDougall, 1974). Volcanic provinces are placed into these categories based upon a variety of spatial, temporal, and petrologic factors. The most subjective of these factors is whether or not volcanic provinces fit within predetermined North-South age-progressive chains of volcanoes or whether their eruption is independent of these pathways. Age-progressive activity forms chains of volcanoes that appear to track the passage of mantle plumes, such as the Cosgrove and Tasmantid chains, beneath the horizontally translating plate (Figures 1c-1e; Davies et al., 2015;Wellman & McDougall, 1974;Sutherland et al., 2012). Offshore seamount chains, central volcanoes and leucitites lie along age-progressive trends (Davies et al., 2015;Johnson, 1989).
Lava fields are defined as age-independent volcanic provinces, which consist exclusively of mafic material (Johnson, 1989;Wellman & Mc-Dougall, 1974). This age-independent activity does not obviously coincide with the passage of mantle plumes beneath the plate. Instead, it is attributed to a variety of processes that include convective upwelling at locations with steep gradients of lithospheric thickness, influx of fertile slab-derived material from subduction of the Pacific plate along the Tongan-Kermadec Trench, offshore lithospheric stretching, transtensional mantle decompression, and heat transfer from the adjacent Pacific plate mantle (Figure 1c; Demidjuk et al., 2007;Mather et al., 2020;Rawlinson et al., 2017). However, numerous lava fields are active at the same time as putative mantle plumes pass adjacent to them and so it is unclear how useful this categorization is ( Figure 1c). Here, we refer to offshore seamounts, central volcanoes and leucitites as age-progressive provinces. We refer to lava fields as age-independent provinces.
Mantle plumes are vertical upwellings of hot mantle material which originate at thermal boundary layers, such as the core-mantle boundary (Morgan, 1972). When a plume reaches the base of a plate, melting can occur. Variations in melt fraction are reflected in the volume and chemical composition of mafic lavas erupted at the surface (Klein & Langmuir, 1987;McKenzie & O'Nions, 1991). Here, we characterize changes in the thermal structure of the upper mantle during Cenozoic times by exploiting the geochemical composition of these rocks. Temperature variations that we infer from analysis of volcanic rocks can be compared with independent observations. Our principal goal is to determine asthenospheric temperature and lithospheric thickness variations along the eastern seaboard of Australia and to show how this structure relates to the timing and distribution of age-progressive and age-independent volcanic provinces. Regional topography is moderated by changes in upper mantle thermal structure (Hoggard et al., 2016). In the final section, we discuss how observed mantle temperature and lithospheric thickness variations contributed to the formation of the Eastern Highlands during Late Cretaceous and Cenozoic times.

Upper Mantle Structure
During a field campaign to Queensland in 2017, 78 igneous rock samples were collected ( Figure 2). This campaign was designed to target a combination of age-independent and age-progressive volcanic provinces where limited published geochemical information exists (e.g., Chudleigh, Sturgeon, Nebo). Samples were analyzed for major, trace and rare earth elements using XRF and ICP-MS techniques at Geoscience Australia, at the University of Cambridge, and at the University of Edinburgh (Supporting information; Data Set S2). These measurements helped to augment an extensive database of legacy geochemical measurements for eastern Australia that also includes isotopic analyses (Data Set S3; Table 1). The combined database is used to estimate potential temperature, T p , the projected temperature that a packet of mantle has if adiabatically decompressed to surface pressures, and lithospheric thickness, a, beneath each volcanic province. By estimating T p and a at a range of locations, spatio-temporal variations in mantle structure can be established. BALL ET AL.   (Cohen et al., 2017;Coventry et al., 1985;Gibson, 2007;Griffin & McDougall, 1975;Johnson, 1989;Jones et al., 2020;Stephenson & Griffith, 1976;Stephenson et al., 1998;Whitehead et al., 2007;Wyatt & Webb, 1970). Black polygons = surface expression of volcanism; circles/triangles = loci of age-independent/progressive whole-rock geochemical analyses; red/white fill = newly acquired/legacy data (Jones et al., 2020;O'Reilly & Zhang, 1995;Stephenson et al., 2000;Zhang et al., 2001). Inset map shows location along eastern seaboard of Australia.

Trace Element Modeling
Concentrations of rare earth elements (REEs) within a basaltic rock are sensitive to the conditions present when it's primary melt is generated (Gast, 1968;Kay & Gast, 1973). REEs are incompatible which means that their concentration within the melt phase reduces as melt fraction increases. Heavy REEs, such as Yb, are more compatible in garnet than light REEs, such as La. Since garnet is only present at depths that are approximately ≥63 km, REEs can be used to determine depth and temperature of melting (Jennings & Holland, 2015;Kay & Gast, 1973;McKenzie & O'Nions, 1991). If mantle composition is fixed and if melting ceases at the base of the lithosphere, depths to base and top of the melt column are predominantly controlled by T p and a, respectively (McKenzie & Bickle, 1988). Final melt fraction is determined by the difference between these two depths. Here, we exploit the INVMEL-v12 forward model to generate a suite of REE profiles by varying T p and a in 1°C and 1 km increments between 1,250 and 1,500°C and between 40 and 80 km, respectively (Ball et al., 2021;McKenzie & O'Nions, 1991;White & O' Nions 1992). Averaged REE profiles from each volcanic province are compared with these modeled profiles (e.g., Figures 3a, 3c and 3e). A grid search is carried out to identify models with the smallest root mean squared (rms) misfits between calculated and average observed concentrations of La, Ce, Sm, Dy, and Yb normalized by their respective standard deviations (e.g., Figures 3b, 3d and 3f). These elements are selected because they span an appropriate range of REE partition coefficients and they are commonly measured. Models are deemed acceptable if they yield rms misfit values that are <1.5 times the misfit value of the best-fitting model.
Prior to modeling, Data Set S3 is screened in order to identify the most primitive samples from each province and thus mitigate the effects of fractional crystallisation. If possible, we exploit samples with MgO > 8.5 wt%. In provinces where <5 samples have MgO > 8.5 wt%, this MgO threshold was reduced by increments of 0.5 wt% until there are ≥5 samples with a lower limit imposed at 7 wt% of MgO. Back-calculation of olivine loss is carried out for acceptable samples until they are deemed to be in equilibrium with forsterite-90 olivine (Lee et al., 2009). Observed REE concentrations are then corrected for the amount of olivine addition, assuming that olivine contains no REEs. Note that in several instances, particularly when the MgO threshold is ≤ 8 wt%, clinopyroxene fractionation could have occurred (Jones et al., 2020). This unaccounted-for fractionation of clinopyroxene leads to an overestimate in REE concentrations for primary melts. Hence T p values obtained by modeling samples that have experienced significant clinopyroxene fractional crystallization are necessarily underestimates.
BALL ET AL.
Partition coefficients for olivine, orthopyroxene and spinel required by the INVMEL-v12 algorithm are listed by McKenzie and O'Nions (1995). Coefficients for clinopyroxene and garnet are calculated using parameterizations described by Wood and Blundy (1997) and by Westrenen et al. (2001), respectively. Mineral compositions together with mantle mineralogy as a function of pressure and melt fraction are provided by McKenzie and O'Nions (1991) and by McKenzie and O'Nions (1995). We assume that the spinel-garnet transition occurs between 63 and 72 km (Klöcking et al., 2018). Source composition is calculated by matching the average ɛNd value for each volcanic province. This calculation is carried out using a linear mixture BALL ET AL.
10.1029/2021GC009717 6 of 22 of primitive and depleted peridotite mantle sources, which are assumed to have ɛNd values of 0 and +10, respectively (McKenzie & O'Nions, 1995). In volcanic provinces where no ɛNd measurements are available, the average ɛNd value for Cenozoic volcanic provinces of the eastern seaboard of Australia is assumed (i.e., +4.18; Table 2). The hydrous melting model of Katz et al. (2003) is used to calculate melt fraction as a function of depth for each T p -a pair. Note that values of various parameters exploited by this melting model have been revised to honor more recent experimental constraints (Shorttle et al., 2014). Source water content is set relative to the Ce content of the source (H 2 O/Ce = 200; Michael, 1995). Assuming that mid-oceanic ridges are underlain by depleted mantle (i.e., ɛNd = +10), this specific melting parameterization predicates a T p value of 1,312°C in order to generate average oceanic crustal thickness (6.9 km; Richards et al., 2018). We assume that this value represents the ambient temperature of sub-plate mantle.

Results
Figure 3 summarizes results for three volcanic provinces: McBride, an age-independent volcanic province of Northern Queensland that has been active since ∼8 Ma; Peak Range, an age-progressive volcanic province linked to the Cosgrove plume track that records volcanism between ∼26 and 47 Ma; and the Newer Volcanic Province, an age-independent province in Victoria that has been active since ∼6 Ma ( Figure 1). In each case, an excellent fit between observed and calculated REE concentrations is obtained. Light REEs (e.g., La) typically exhibit greater standard deviations than heavy REEs (e.g., Yb). Consequently, the ratio of melting within the spinel and garnet stability fields, which principally controls heavy REE concentrations, has a greater effect upon misfit than melt fraction, which controls light REE concentrations. Therefore, a negative trade-off between T p and a exists whereby a hotter asthenosphere but thinner plate, or a cooler asthenosphere but thicker plate, can match the observed concentrations within acceptable limits ( Figure 3f). Results for individual provinces are presented in Supporting Information Figures S1-S4 and summarized in Table 2. The highest potential temperatures are recorded beneath the age-progressive provinces of central Queensland (CQ: 1348°C; Figure 4a). These values are ∼50-100°C hotter than those values obtained for the age-independent provinces of New South Wales and Victoria (NSW and V: 1286°C). Temperature estimates for the age-independent volcanic provinces of Northern Queensland lie between these two extrema (NQ: 1312°C). Beneath the eastern seaboard of Australia, the upper boundary, at which melting ceases, is equivalent to a lithospheric thickness of 45-70 km.
T p estimates for NSW and Victoria are up to 50°C cooler than the ambient mantle value (e.g., 1,312°C). However, the exact values of T p and a determined using the INVMEL-v12 algorithm can change with respect to this ambient mantle value if model parameters are varied (Ball et al., 2019). For example, compositional heterogeneity of the upper mantle can influence both temperature and plate thickness. Increasing the mantle water content by 0.01 wt% shifts the value of T p by up to +10°C and the value of a by up to −2 km. Changing mantle composition from primitive to depleted mantle decreases T p by −20°C and increases a by +3 km. Deepening the spinel-garnet transition zone by 5 km increases T p and a by +15°C and +5 km, respectively. The INVMEL-v12 algorithm does not take into account the presence of pyroxenitic or harzburgitic lithologies within the upper mantle. If pyroxenite or harzburgite occur within the source region, recovered values of T p will be over-or under-estimates, respectively (Matthews et al., 2021). In summary, accurately constrained source compositions are required to improve the reliability of absolute T p and a estimates. We acknowledge that the assumed mantle source composition is unlikely to represent the actual complexity of eastern Australian mantle. However, isotopic analyses do imply that similar mantle sources exist beneath the Cenozoic volcanic provinces of this region (Jones et al., 2020). Even though the recovered values of T p and a have significant uncertainties, we contend that inter-provincial differences in T p and a are robust.
Mantle melts can become enriched in incompatible elements, such as REEs, by lithospheric interactions during ascent toward the surface (Foley, 1992). Instances of contamination caused by melting of metasomatized lithospheric mantle have been reported for lavas from various Australian Cenozoic volcanic provinces (e.g., Jones et al., 2020;O'Reilly & Griffin, 1988;Shea & Foley, 2019;Zhang et al., 2001). This contamination increases REE concentrations within the melt, which causes calculated values of T p to be underestimated (Ball et al., 2019;Klöcking et al., 2020). In locations where lithospheric contamination is reported, values of T p should be regarded as lower estimates. In extremal cases, lithospheric contamination can enrich melts to such an extent that the INVMEL-v12 algorithm no longer successfully matches observed concentrations, Note. n s = Number of samples; XTB = lithospheric thickness estimates from xenolith/xenocryst thermobarometry ; NVP = Newer Volcanic Province; OVP = Older Volcanic Province  Figure 1a). We infer that lavas from these particular provinces could have either experienced significant lithospheric contamination, or have been sourced directly by melting of metasomatized lithosphere, which would yield extremely enriched incompatible element concentrations.

Comparison With Thermobarometric and Tomographic Observations
It is instructive to compare our results with estimates of T p and a estimated by thermobarometric studies of mantle xenoliths and by analyzing calibrated seismic tomographic models. Mantle xenoliths occur in many of the Cenozoic volcanic provinces (Kay & Kay, 1983;O'Reilly & Griffin, 1985). Compositions of minerals within these xenoliths can be used to estimate the final temperatures and pressures at which xenoliths equilibrated (Nickel & Green, 1985;Nimis & Taylor, 2000). Hoggard et al. (2020) fit geothermal profiles through xenolith and xenocryst thermobarometric estimates at different locations across Australia (Figure 4c). In each case, the lithosphere-asthenosphere boundary is estimated to occur at the depth at which temperature exceeds 1,175°C. In oceanic regions, this isothermal surface coincides with a peak change in azimuthal anisotropy which demarcates the transition between conductive and convective mantle (Burgos et al., 2014;Richards et al., 2018). Values of a estimated in this way can be compared with those calculated by REE modeling ( Figure 5). In each case, these estimates agree to within an average of ±12 km.
Shear-wave velocities, V s , within the mantle are sensitive to variations in temperature, T. For a given shearwave tomographic model, temperature structure of the mantle can be determined by calibrating the relationship between V s and T as a function of depth. A range of global and regional tomographic models together with calculated estimates of T p and a, exist for the Australian continent (Fishwick & Rawlinson, 2012;Hoggard et al., 2020;Kennett et al., 2013;Schaeffer & Lebedev, 2013). Here, we compare values of T p that are calculated by REE modeling with a global mantle temperature model (Figure 4b; Hoggard et al., 2020). This particular temperature model was constructed by exploiting an empirical anelastic parameterization which describes V s as a function of T Schaeffer & Lebedev, 2013;Yamauchi & Takei, 2016). The parameterization includes unconstrained constants and it is calibrated by varying these constants to minimize the misfit between observations of mantle temperature structure and the calculated temperature model. Our chosen temperature model is underpinned by the SL2013sv tomographic model and it was principally calibrated by comparing V s (T) with an oceanic plate cooling model between depths of 0 and 125 km and to an average asthenospheric temperature profile calculated between depths of 225 and 400 km (Richards et al., 2018Schaeffer & Lebedev, 2013;Shorttle et al., 2014).
Along the eastern seaboard of Australia, Cenozoic volcanism is focused in regions where slow shear-wave velocity anomalies and ambient or elevated temperatures lie beneath the plate. In Figure 4b, we compare tomographic estimates of T p that are averaged between depths of 100 and 200 km with REE-modeled estimates. In order to compare these different estimates of T p , results from each modeling approach are normalized with respect to the ambient mantle value, which yields temperature anomalies, ΔT p . Both REE-modeled and tomographic estimates of T p are verified by their ability to generate the average oceanic crustal thickness using an identical melting parameterization (Shorttle et al., 2014). However, the ambient value of T p for the tomographic model is slightly higher because, in the case of the tomographic model, the mantle source is assumed to be anhydrous (1,333°C compared to 1,312°C; Ball et al., 2021;Hoggard et al., 2020;Richards et al., 2018Richards et al., , 2020. Since the Australian plate is rapidly translating northward at ∼7 cm/yr, REE-modeled and tomographic ΔT p estimates are only strictly comparable for the youngest volcanic provinces (<10 Ma; Northern Queensland and Newer Volcanic provinces; DeMets et al., 2010). With the exception of Silver Plains, Piebald and McLean (SPM), both ΔT p estimates agree to within ± 20°C (Figure 4d). Note that when these techniques are applied to a global database of intraplate volcanic rocks, a systematic offset of 30-50°C is observed between REE-modeled and tomographic estimates of ΔT p (Figure 4d; Ball et al., 2021). This offset could explain why REE-modeling estimates in NSW and Victoria are lower than ambient mantle values by ∼30°C.
Changes of lithospheric thickness can also be estimated from a global temperature structure obtained for the SL2013sv model by tracking the 1,175°C isothermal surface . Lithospheric thicknesses calculated using this V s -to-T calibration, which is based upon an oceanic plate cooling model, have BALL ET AL.

10.1029/2021GC009717
greater uncertainties in continental regions than in oceanic regions (Priestley & McKenzie, 2013). Since Australia is dominated by thick and depleted lithosphere, it may be more appropriate to use a lithospheric thickness model calibrated to continental temperature observations. Here, we adopt a lithospheric thickness model generated by calibrating the regional FR12 tomographic model to geothermal profiles calculated from xenolith thermobarometry (Figure 4c; Fishwick & Rawlinson, 2012;Hoggard et al., 2020). This regional lithospheric thickness model predicts thinner lithosphere beneath cratonic Australia compared with the global model. In areas of thin lithosphere relevant to our study, these models broadly agree. Both regional and global tomographic estimates of a are consistent with geochemical estimates within an average BALL ET AL.

10.1029/2021GC009717
10 of 22  Schaeffer & Lebedev, 2013). Red/white/blue contours = positive/zero/negative ΔT p at 25°C intervals. Circles = age-independent provinces <10 Ma colored by ΔT p ; dashed red line between C and C′ = locus of transects shown in Figure 6, which tracks relative motion of Cosgrove plume (Davies et al., 2015). (c) Lithospheric thickness calculated by tracking 1,175°C isothermal surface from temperature-calibrated FR12 model (Fishwick & Rawlinson, 2012;Hoggard et al., 2020). Diamonds = lithospheric thickness estimates constrained by xenolith/xenocryst thermobarometry  of ±12 km along the length of the eastern seaboard ( Figure 4e). Therefore, seismic tomographic models provide a useful means of extrapolating lithospheric thickness between point-wise xenolith thermobarometric constraints. Tomographic estimates of a are <80 km for all volcanic provinces, with the exception of the east Australian leucitites (e.g., El Capitan (E), Griffith (G), Cosgrove (Co); Figure 1a). Elevated REE concentrations indicative of lithospheric contamination reported for Cosgrove and El Capitan, coupled with a lack of BALL ET AL.
xenolith thermobarometric observations, mean that these tomographic estimates of lithospheric thickness cannot be geochemically corroborated.

Mechanisms for Generating Cenozoic Volcanism
All three techniques used to estimate upper mantle structure along the eastern seaboard of Australia are subject to significant uncertainties. Nevertheless, the self-consistency of our results is notable, which supports the inference that chemical composition of Cenozoic volcanic rocks is a function of upper mantle conditions at the time of melt generation. Here, we discuss potential mechanisms for generating both age-progressive and age-independent volcanism in this region.
The REE modeling results presented here include the first attempt to estimate the temperature of the putative Cosgrove mantle plume. To reproduce REE concentrations of the Central Queensland age-progressive volcanic rocks, potential temperatures, that are ∼50-100°C hotter than those needed to produce age-independent volcanism of New South Wales, are required. Assuming that age-independent volcanism of NSW was generated at ambient mantle temperatures, we estimate that the Cosgrove plume had an excess temperature of 50-100°C when Central Queensland translated over it. REE concentrations and absolute temperatures estimated for Central Queensland are similar to those calculated using the same techniques applied to the present-day Réunion mantle plume (Ball et al., 2021). These temperatures are amongst the hottest recorded on Earth during Neogene-Quaternary times using this modeling scheme (Ball et al., 2021).
The Older and Newer Volcanic provinces are geographically proximal to one another within the state of Victoria (Figure 1a). The Older Volcanic Province (OVP) was active between ∼90 and 14 Ma, whilst the Newer Volcanic Province (NVP) has been active since ∼6 Ma (Gray & McDougall, 2009;Price et al., 2014). Volcanic activity of the NVP is coincident with passage of the Cosgrove mantle plume beneath Victoria (Davies et al., 2015). Since, at the time of eruption, the Cosgrove plume was present beneath the NVP and absent beneath the OVP, we would expect to see a difference in T p estimates between these provinces. However, we obtain similar T p results in both cases, which indicates that the plume may have thermally waned by the time it passed beneath Victoria ( 1308°C, respectively; Figure 6a). Significantly, a discrete low velocity anomaly at the predicted location of the present-day plume head is generally absent from tomographic models (C' in Figure 4b). We note, however, that recent seismic imaging using observations from the WOMBAT and BASS arrays, reveals radial anisotropy that could be generated by radial flow away from a plume conduit (Bello et al., 2019).
The Eastern Australian leucitites are also thought to have been generated by passage of the Cosgrove mantle plume (Davies et al., 2015). These K-rich, leucite-bearing lava flows crop out between Central Queensland and the Newer Volcanic Province (e.g., El Capitan, Griffith, Cosgrove volcanic provinces; Figure 1a). Using either global or regional tomographic models, this area is predicted to have lithosphere >100-200 km thick (Figure 6b). Thick lithosphere largely inhibits melting and so, based upon these models, no magmatism would be expected to occur beneath regions where leucitites crop out. However, a regional P-wave tomographic model based upon observations from the WOMBAT array, which is positioned south of ∼29°S, suggests that these low-volume leucititic melts could form over localized regions of thin lithosphere that cannot be resolved by large-scale tomographic models (Davies et al., 2015). Unfortunately, regional tomographic models are difficult to calibrate with a view to extracting temperature and lithospheric structure. Nonetheless, we can estimate values of lithospheric thickness from the WOM-BAT-based model by assuming that minimum and maximum P-wave velocities along the Cosgrove track coincide with the thinnest and thickest lithosphere, respectively (Davies et al., 2015). By linearly scaling P-wave velocities with respect to values of lithospheric thickness between these extreme values, regions of lithosphere <80 km thick are predicted to occur beneath all leucitite outcrops (Figure 6b). In these regions, the lithosphere is sufficiently thin to permit small degrees of partial melting (Figure 6c). Note that, if the leucitites are formed by melting of metasomatized lithosphere, which is a reasonable conclusion given their K-rich composition, elevated asthenospheric temperatures are required to initiate melting and to form age-progressive volcanism.
Although our results suggest thermal waning of the Cosgrove plume, we cannot constrain its timing. Mapping of volcanic provinces across Central Queensland reveals a southward decrease in melt volume, BALL ET AL.

10.1029/2021GC009717
which is attributed to a waning plume (Jones et al., 2017). In this region, it is likely that lithospheric thickness increases southward since volcanic activity occurs closer to the thick cratonic interior (Figure 4c). Therefore, a reduction of melt volume within Central Queensland could, instead, be due to the plume beginning to track beneath thicker lithosphere (Davies et al., 2015). However, exploiting current tomography-based lithospheric thickness models, a pronounced increase in lithospheric thickness is not observed until south of the Central Queensland volcanic provinces, and lithospheric thickness variations cannot account for T p differences between Central Queensland and Victoria (Figure 6; Hoggard et al., 2020). Gradual, rapid, or even delayed, thermal waning of the Cosgrove plume between Central Queensland and Victoria can each account for formation of all volcanic provinces along the plume track ( Figure S5). Unfortunately, a paucity of REE compositions for age-progressive volcanism in Southern Queensland, for NSW, and for the offshore seamounts limits the spatial coverage of our study. Further analysis of mafic rocks from age-progressive volcanic provinces would be helpful to constrain in more detail the number, location and temperature excess of mantle plumes that have passed beneath Australia during Cenozoic times.
Results of REE modeling suggest that generation of age-independent volcanic provinces across NSW requires a combination of either ambient or modestly elevated asthenospheric temperatures and thinned lithosphere. This combination could be consistent with the hypothesis that age-independent volcanism is generated by either edge-driven convection and/or shear-driven upwelling of mantle material. In this hypothesis, the juxtaposition of thick cratonic lithosphere with thinner coastal lithosphere facilitates generation of local edge-driven convective cells that sustain intermittent volcanic activity over protracted periods BALL ET AL.   Figure 4b which tracks the Cosgrove plume (Davies et al., 2015). (a) Temperature profile of Cosgrove plume estimated using results of geochemical modeling. Circles/triangles = ageindependent/-progressive estimates of T p from rare earth element (REE) modeling; red line = estimated temperature profile for the Cosgrove plume. (b) Red solid/dashed lines = estimates for lithospheric thickness from FR12/FR12 modified by linearly scaling lithospheric thickness using average P-wave velocities over depth range of 50-200 km from WOMBAT database (Davies & Rawlinson, 2014;Hoggard et al., 2020); circles/triangle/diamonds = age-independent/progressive/xenolith thermobarometric estimates of lithospheric thickness; dotted line = depth of solidus, assuming that mantle H 2 O concentration = 210 ppm (Shorttle et al., 2014); NE = Nebo; PR = Peak Range; SP = Springsure; BU = Buckland; BY = Byrock; EC = El Capitan; BH = Begargo Hill; GR = Griffith; NV = New Volcanic Province. (c) Red solid/dashed lines = melt fraction estimates obtained by assuming adiabatic upwelling to base of lithosphere from FR12 and modified FR12 models (Shorttle et al., 2014); circles/triangles = age-independent/-progressive melt fraction estimates.
of time (Davies & Rawlinson, 2014;Demidjuk et al., 2007;Rawlinson et al., 2017). Rapid northward translation of Australia could modify the flow field in the mantle beneath, initiating shear-driven upwelling and altering the distribution and extent of volcanic activity at the surface (Davies & Rawlinson, 2014;Rawlinson et al., 2017). Since eruptions have not been continuous throughout Cenozoic times, it is possible that minor variations of mantle temperature, lithospheric structure, mantle flow and/or asthenospheric composition act to modulate the melting process (Demidjuk et al., 2007;Mather et al., 2020;Rawlinson et al., 2017).
Our REE modeling results suggest that ambient to modestly elevated asthenospheric potential temperatures exist at the present day beneath Northern Queensland. These results are corroborated by seismic tomographic estimates of T p (Figures 4b and 4d). As such, it is unlikely that age-independent volcanism of Northern Queensland and NSW are generated by identical processes. Northern Queensland lies at the leading edge of the translating Australian plate, which is a dynamically unfavorable locus for melting to have occurred either by edge-driven convection or by shear-driven upwelling (Davies & Rawlinson, 2014;Rawlinson et al., 2017). Instead, intraplate volcanism of Northern Queensland could have been generated by the recent arrival of a mantle plume (Kennett & Davies, 2020). Beneath the Coral Sea, a deep-rooted plume is inferred to exist within the lower mantle that has associated seismic attenuation within the upper mantle (French & Romanowicz, 2015;Kennett & Davies, 2020). However, no bathymetric swell is observed at the predicted location of this plume. Instead, a significant negative residual depth anomaly is recorded throughout the Coral Sea where stratigraphic evidence suggests that this subsidence occured during Neogene times (Czarnota et al., 2013;Hoggard et al., 2017;Müller, Gaina, & Clark, 2000). Absence of a clearly defined plume within the upper mantle, combined with the geographic distribution of volcanism, could result from interaction between a rising plume and subducted oceanic lithosphere within the mid-mantle (Kennett & Davies, 2020). Further investigation is required to determine the mechanism that generates the volcanism of Northern Queensland, although the modestly elevated temperatures estimated by REE modeling are indeed compatible with the presence of a mantle plume.

Evolution of the Eastern Highlands
We have used a database of newly acquired and legacy geochemical measurements to constrain values of T p and lithospheric thickness beneath the eastern seaboard of Australia. Our results agree to within an average of ±20°C and ±12 km of tomographic and xenolith thermobarometric estimates. They show that lithosphere beneath the eastern seaboard is < 80 km thick throughout Cenozoic times. This thin lithosphere is underlain by asthenospheric mantle whose potential temperature varies as a function of time and space. The bulk of Cenozoic volcanism occurs within the Eastern Highlands, which sit at present-day elevations of 500-1,500 m. Here, we explore the geodynamical implications of our results. We suggest that a combination of lithospheric thinning and asthenospheric flow could be responsible for the formation and longevity of the Eastern Highlands.
To assess the upper mantle's contribution to formation of the Eastern Highlands, we must first constrain the timing and extent of uplift. In Northern Queensland, Late Cretaceous marine sedimentary rocks capped by lateritic deposits are distributed along the western edge of these highlands (Coventry et al., 1985;Wellman, 1987). Laterite deposits build up during prolonged spells of sub-aerial weathering in hot tropical environments where topographic gradients are minimal. These combined stratigraphic and geomorphic observations suggest that sub-aerial exposure of Northern Queensland occurred between 100 and 45 Ma (Wellman, 1979(Wellman, , 1987. We infer that Cenozoic volcanism occurred onshore and commenced when topographic relief was low. To the south, post-Triassic marine sedimentary rocks are absent which means that the commencement of sub-aerial exposure is poorly constrained. Within the Monaro volcanic province, basaltic lavas that erupted at ∼45 Ma flowed along deeply incised paleovalleys, which implies that a component of regional topography existed during Mesozoic times (Taylor et al., 1985;Vandenberg, 2010). Forward and inverse modeling of longitudinal river profiles from the Eastern Highlands reveals a two phase uplift history (Czarnota et al., 2014;Salles et al., 2017). Within the southeastern highlands, where erosion rates are best-constrained, these phases occurred in Late Cretaceous (i.e., 120-80 Ma) and Cenozoic (i.e., 60-30 Ma) times, respectively (Czarnota et al., 2014;Young & McDougall, 1993). Apatite (U−Th)/He thermochronologic analyses of sedimentary BALL ET AL.

10.1029/2021GC009717
rocks from the Drummond basin of southern Queensland suggest that cooling and exhumation of this basin began at ∼90 Ma (Zhang et al., 2019). A subsequent period of Cenozoic uplift and denudation is consistent with other geomorphologic, fission track and stratigraphic studies (Faiz et al., 2007;Holdgate et al., 2008;O'Sullivan et al., 1995O'Sullivan et al., , 2000Wellman, 1987). This second phase of regional uplift coincides with onset of intraplate volcanism. Sub-aerial lava flows with radiometric dates of ∼30 Ma crop out at elevations of ∼30 m along the coastline adjacent to the southeastern and central highlands (Young & McDougall, 1982;Young & Wray, 2000). These flows demonstrate that most of the regional uplift occurred prior to Oligocene times. Although this two phase uplift model for the Eastern Highlands is widely accepted, the relative contributions of each uplift phase to present-day topography remains debated (Czarnota et al., 2014;Hill, 1999;Holdgate et al., 2008;Salles et al., 2017;Vandenberg, 2010;Webb, 2017).
Some previous studies have attributed regional uplift of the Eastern Highlands to crustal, rather than mantle, processes. These processes include tectonic shortening, rifting during opening of the Tasman Sea, magmatic underplating, and erosion of pre-existing topography (Holdgate et al., 2008;Kohn et al., 2002;Lister & Etheridge, 1989;Stephenson & Lambeck, 1985;Van der Beek et al., 1999). These alternative proposals generally assume that present-day topography is supported by crustal thickness and density variations. For example, in the southernmost Eastern Highlands, topography does indeed correlate with crustal thickness and with present-day seismic activity. Therefore, on local scales, this assumption may hold (Braun et al., 2009;Kennett et al., 2011;Sandiford et al., 2004). Models of the SE Australian stress regime indicate that an E-W compressive stress field developed at 5-10 Ma (Sandiford et al., 2004). Crustal shortening attributable to these stresses could account for a small proportion (i.e., < 200 m) of uplift during or after the second phase of Eastern Highlands formation (Braun et al., 2009). However, the Eastern Highlands coincide with positive long wavelength free-air gravity anomalies (Figure 1a). Thus it is less likely that the whole of the Eastern Highlands is supported by crustal isostasy alone. Instead, admittance analysis of eastern Australia implies that a combination of mantle density variations and normal stresses, probably generated by asthenospheric flow, support the Eastern Highlands (Czarnota et al., 2014).
It is unlikely that the Eastern Highlands have been consistently supported by thermally buoyant asthenosphere since our geochemical analyses suggest that the generation of age-independent volcanism does not require excess mantle temperatures. Rapid northward translation of the Australian plate during Cenozoic times also negates this possibility. Numerous studies have advocated lower mantle density variations as a key control of eastern Australian topography. Regional uplift of the Eastern Highlands could conceivably be produced by migration of the Australian continent away from a sinking slab within the mantle and toward the Pacific Large Low Velocity Province (Gurnis et al., 1998;Müller et al., 2016;Salles et al., 2017). The predicted surface rebound of ∼1,000 m generated by lateral motion over these deep mantle features encompasses both the Eastern Highlands together with the Coral and Tasman seas, which lie to the east. In contrast to the Eastern Highlands, the Tasman and Coral seas are associated with negative long wavelength free-air gravity anomalies, negative residual topography, and hundreds of meters of Neogene-Quaternary subsidence (Figure 1a; Czarnota et al., 2013;Müller, Lim, & Isern, 2000). Accordingly, whilst this hypothesis can match the magnitude and timing of uplift observed for the Eastern Highlands, it is difficult to reconcile it with proximal offshore observations. Finally, modest predictions of dynamic topographic support at the present day (i.e., −100 m to +200 m) are at odds with the existence of a positive free-air gravity anomaly and with the requirement for sub-crustal support across the Eastern Highlands (Czarnota et al., 2014). These observations imply that additional processes are necessary to generate and maintain both the Eastern Highlands and the dynamic topographic low of the Coral and Tasman seas.
An alternative (or additional) mechanism for generating regional uplift is thinning of lithospheric mantle (Bird, 1979). Our results show that the lithosphere beneath the eastern seaboard of Australia was thin prior to, or coincident with, onset of volcanism. Thin lithosphere persists to the present day. Using a simple isostatic calculation, thinning of lithosphere from 120 to 70 km by removal of lithospheric mantle material will generate either 0.72 km or 1.26 km of regional uplift, depending upon whether the remnant lithosphere is thermally dis-equilibrated or re-equilibrated (Appendix A). Therefore, lithospheric erosion can provide sufficient buoyancy to support the Eastern Highlands, which have a present-day elevation of 0.5-1.5 km. Any hypothesis that accounts for support of the Eastern Highlands must reproduce its two phase uplift history (Czarnota et al., 2014;Müller et al., 2016;Salles et al., 2017). These highlands could initially have been uplifted by erosion of lithospheric mantle prior to the opening of the Southern ocean and of the Tasman sea at ∼90 Ma (Gaina et al., 1998). The second phase of uplift coincides with onset of volcanism and with an increase in plate speed as the Australian plate begins to separate from Antarctica and translate toward the Melanesian subduction zone at 60-30 Ma (DeMets et al., 2010;Müller, Gaina, & Clark, 2000). Along the eastern seaboard of Australia, thin coastal lithosphere abuts the thick cratonic interior at the present day ( Figure 4c). As plate speed increases, local convection cells that develop at marked changes in lithospheric thickness are conceivably enhanced (Conrad et al., 2010;Davies & Rawlinson, 2014;Rawlinson et al., 2017). These putative cells could destabilize overlying lithospheric mantle and could also contribute hundreds of meters of dynamic support at wavelengths of order 10 2 -10 3 km, as a consequence of normal stresses imparted by shallow mantle flow (Colli et al., 2016;Kim & So, 2020). Return convective flow could also generate observed subsidence of the Coral and Tasman seas. The surprising longevity of thinned lithosphere along the eastern seaboard and, therefore, of the Eastern Highlands, could result from some combination of shear-driven upwelling, edge-driven convection and/or periodic small-scale melting, which act to inhibit lithospheric rethickening.

Conclusions
A revised and augmented geochemical database of Cenozoic volcanism along the eastern seaboard of Australia has been compiled and analyzed. This database includes newly acquired samples and analyses from a series of under-represented locations across Queensland. Geochemical modeling of this database suggests that continental lithosphere is < 80 km thick and that ambient, or modestly elevated, mantle temperatures existed beneath age-independent volcanic provinces during Cenozoic times. Both thinning of the lithospheric mantle and development of small-scale convection cells within the uppermost asthenosphere prior to, or coincident with, the onset of volcanism could have contributed to regional uplift of the Eastern Highlands. Excess mantle temperatures of 50-100°C are required to produce the age-progressive volcanic provinces of Central Queensland. These excess temperatures mark the passage of the Cosgrove mantle plume beneath the horizontally translating plate. We observe a thermal waning of this plume as it passes beneath the plate. Results of modeling of REE concentrations agree with independent constraints obtained from xenolith thermobarometric and seismic tomographic analyses, which validates our proposed lithospheric architecture. In the future, quantitative geochemical modeling of older volcanic provinces should yield useful insights into the way in which mantle structure of the Australian plate has evolved through geologic time.

Appendix A: Isostatic Calculations
The amount of regional uplift, U, generated by removal of dense lithospheric mantle is given by where z and ρ refer to thicknesses and densities, respectively. Subscript labels are as follows: c indicates crust; l indicates lithosphere; lm indicates lithospheric mantle; and a indicates asthenosphere. Subscripts o and n refer to original and modified lithospheric columns, respectively. Lithospheric and asthenospheric densities are functions of pressure and temperature. Here, we use identical equations and parameter values to those described by Ball et al. (2019). We assume T p = 1,330°C, an asthenospheric temperature gradient of 0.6°C km −1 , a pressure gradient of 0.033 GPa km −1 , a reference mantle density of 3.33 Mg m −3 , a thermal expansion coefficient of 4 × 10 −5 °C −1 , and a bulk modulus of 115.2 GPa. In the thermally equilibrated case, the base the of lithospheric column is assumed to have a T p of 1,330°C with temperature that linearly decreases to zero at Earth's surface. In the disequilibrated case, temperature profile remains fixed so that the base of the lithosphere has T p < 1,330°C and temperature linearly decreases to zero at the Earth's surface. For these calculations, we assume that z c = 35 km and that the lithosphere is not depleted. If lithosphere with a thickness of 120 km is instantaneously thinned to 70 km, U = 0.72 km. If this perturbed lithospheric column is thermally re-equilibrated, it becomes less dense and U increases to 1.25 km. The effects of mantle flow and regional denudation can act to amplify U.

Data Availability Statement
Newly acquired analyses are provided in Supporting information, through Geoscience Australia and in the EarthChem data repository (https://doi.org/10.26022/IEDA/111889). INVMEL-v12 software was developed and is owned by D. McKenzie to whom requests for access should be directed. For access to other software please contact authors.