Spatial and temporal patterns of Australian dynamic topography from River Profile Modeling

Despite its importance, the temporal and spatial evolution of continental dynamic topography is poorly known. Australia's isolation from active plate boundaries and its rapid northward motion within a hot spot reference frame make it a useful place to investigate the interplay between mantle convection, topography, and drainage. Offshore, dynamic topography is relatively well constrained and can be accounted for by Australia's translation over the mantle's convective circulation. To build a database of onshore constraints, we have analyzed an inventory of longitudinal river profiles, which is sensitive to uplift rate history. Using independently constrained erosional parameters, we determine uplift rates by minimizing the misfit between observed and calculated river profiles. Resultant fits are excellent and calculated uplift histories match independent geologic constraints. We infer that western and central Australia underwent regional uplift during the last 50 Myr and that the Eastern Highlands have been uplifted in two stages. The first stage from 120 to 80 Ma, coincided with rifting along the eastern margin and its existence is supported by thermochronological measurements. A second stage occurred at 80–10 Ma, formed the Great Escarpment, and coincided with Cenozoic volcanism. The relationship between topography, gravity anomalies, and shear wave tomographic models suggest that regional elevation is supported by temperature anomalies within the lithosphere's thermal boundary layer. Morphology and stratigraphy of the Eastern Highlands imply that these anomalies have been coupled to the base of the plate during Australia's northward motion over the last 70 Myr.


Introduction
Australia is an important natural laboratory for studying dynamic topography, generated by an interplay between plate motion and subplate convective circulation [Sandiford, 2007;DiCaprio et al., 2009;Heine et al., 2010]. Its importance stems from the continent's distance from active plate boundaries after Cretaceous times together with its rapid northward motion since Eocene times (70 km Myr −1 ). Offshore, the spatial and temporal pattern of dynamic topography is consistent with horizontal translation of Australia over a changing pattern of convective circulation (Figure 1a) [Czarnota et al., 2013]. Here our principal aim is to investigate how methods for extracting uplift rate histories from drainage networks can be used to constrain the onshore evolution of dynamic topography. We start by reviewing independent constraints on Australian dynamic topography.
The Australian landscape is often thought of as being old, flat, and stable [e.g., Pain et al., 2012]. Neo-Proterozoic landforms and Carboniferous weathering profiles indicate that in some places, modern topography has been inherited [Ollier et al., 1988;Stewart et al., 1986;O'Sullivan et al., 2000;Pillans, 2005]. However, landscape preservation can be attributed to successive periods of burial and exhumation rather than to prolonged exposure [O'Sullivan et al., 2000;Belton et al., 2004;Weber et al., 2005;Pillans, 2007]. , dynamic topography) of oceanic lithosphere relative to oceanic age-depth curve from Czarnota et al. [2013]. Colored circles and upward/downward triangles = accurate estimates and lower/upper bounds of residual bathymetry from seismic profiles in 2 • × 2 • bins. Colored squares = relative epeirogenic movement (color = magnitude of movement, number = onset age in Ma). Filigree of lines = residual bathymetry from ship tracks corrected for sediment loading. Red/black/blue contours = positive/zero/negative contours of dynamic topography calculated from long-wavelength gravity anomalies assuming a water-loaded admittance of 23 mGal/km for the oceans; contour interval = 0.5 km [Czarnota et al., 2013]. Colored polygons = youngest onshore marine and coastal strata [Langford et al., 1995]. Black polygons = laterites [Raymond et al., 2012]. (b) Shuttle Radar Topography Mission (SRTM) topography, seismicity, and Neogene faults. Black dots and white circles = epicenters of earthquakes M w 4-6 and >6, respectively (Incorporated Research Institutions for Seismology catalog, 1973Seismology catalog, -2012. Earthquake focal plane solutions M w or M L ≥ 4.5 from centroid moment tensor catalog and P wave polarity local network studies with ≥25 stations and quality rating ≥C [Zoback, 1992;Leonard et al., 2002]. Fault plane solutions of largest earthquakes in a cluster are shown. Red lines =

Geophysical and Geologic Constraints
There is a fragmented record of Cenozoic vertical motion preserve along the edges of the Western Plateau, Pilbara, Kimberley, and Arnhem Land, all of which currently sit on thick lithosphere ( Figure 2d). Rapid regional uplift is recorded by paleoshorelines in the Eucla Basin, on the southern edge of the Western Plateau [Sandiford, 2007]. These paleoshorelines record ∼ 300 m of uplift since ∼ 40 Ma [Hou et al., 2008]. The present-day tilt of these paleoshorelines matches that of the Western Plateau. Uplifted marine sedimentary rocks of a similar age and elevation occur at isolated localities along the western side of this plateau ( Figure 1a) [Collins and Baxter, 1984;Haig and Mory, 2003]. In contrast, subsidence analyses from the Northwest Shelf show that northwestern Australia has been drawn down by up to 700 m during the last ∼9 Ma [Czarnota et al., 2013]. Combining these observations with residual depth estimates from oceanic crust reveals a diachronous pattern of long-wavelength vertical motion. Southern Australia has steadily emerged from a broad region of drawdown which now coincides with the Australian-Antarctic Discordance. More recently, northern Australia has started to submerge as it drifts toward Southeast Asia [Lithgow-Bertelloni and Gurnis, 1997;Sandiford, 2007;Czarnota et al., 2013].
The uplift history of the Eastern Highlands is controversial [e.g., Bishop, 1982;Jones and Veevers, 1982;Lambeck and Stephenson, 1985;Wellman, 1987;Ollier and Pain, 1994;Persano et al., 2002]. Maximum age of uplift is constrained by elevated marine Cretaceous sedimentary rocks in the north and by uplifted Triassic deltaic sequences of the Sydney Basin in the south (Figure 2a) [Wellman, 1987]. The minimum age is constrained by the vertical extent of paleochannels infilled by Cenozoic basalts [Wellman and McDougall, 1974b;Wellman, 1987]. In southeastern Australia, subaerially erupted basalts, which are ∼30 Ma, occur at an elevation of ∼30 m along the coast. Their existence implies that the erosional Great Escarpment was established by this time and that negligible uplift has subsequently occurred [Young and McDougall, 1982]. Similar observation from the Central Highlands suggests that the escarpment developed before 30 Ma with a maximum of 150 m of subsequent uplift [Young and Wray, 2000]. The uplift history between Triassic/Cretaceous and mid-Oligocene times is largely unconstrained although thermochronological data show that a significant phase of rock cooling occurred between 80 and 100 Ma [Gleadow et al., 2002;Persano et al., 2002].

Admittance Analysis
At short (< 10 2 km) wavelengths, topography is supported by lithospheric flexure. At longer wavelengths, topography is usually supported isostatically (i.e., by lithospheric density variation) and/or dynamically (i.e., by mantle convective processes). The admittance, Z, is defined as the ratio of free-air gravity and topography anomalies in the spectral domain. Its variation as a function of wave number allows the mechanism of topographic support to be constrained, provided that topography and gravity anomalies are coherent [Watts, 2001;McKenzie, 2010]. Here we have calculated Z for a box, which encompasses Eastern Australia, CZARNOTA ET AL.
©2014. American Geophysical Union. All Rights Reserved. Admittance between gravity and topography in eastern Australia. (a) Onshore EIGEN6c gravity shaded by topography [Förste et al., 2011]. (b) Decimated SRTM topography. (c) Predicted dynamic topography shaded by topography. Filtered gravity data to exclude wavelengths <500-800 km divided by an admittance of 52 mGal km −1 . (d) Admittance, Z, plotted as function of wave number and wavelength, for data shown in Figure 3a. Solid circles with error bars = Z with 1 error within wave number bands. Blue solid line = predicted values of Z, calculated using an elastic model, which best fit observed values (elastic thickness, t e = 15 km for internal load fraction, F2 = 0.34) [McKenzie, 2003]. Crustal thickness = 38 km (mean for eastern Australia) [Salmon et al., 2013]; upper crustal thickness and density = 10 km and 2.45 Mg/m 3 , lower crustal thickness and density = 28 km and 2.85 Mg/m 3 . Note observed admittance of ∼52 mGal km −1 at longest wavelengths, indicative of dynamic support. Red solid line = predicted dynamic support values of Z, calculated using an analytical solution to Stokes flow equation, which best fits observed values (t e = 15 km, mechanical lithospheric thickness, a = 60 km) [McKenzie, 2010] (e) Coherence between free-air gravity and topography plotted as a function of wave number.
where Δg and t are the multitaper Fourier transforms of gravity and topography, ⟨⟩ represents annular averages over wave number bands, and ⋆ denotes the complex conjugate [McKenzie and Fairhead, 1997]. The coherence is given by Our results are summarized in Figure 3d. At wavelengths of <150 km, Z ∼ 100 mGal km −1 , which is a measure of upper crustal density (i.e., uc = 2.45 kg m −3 ). This value is consistent with the density of unconsolidated Cenozoic sedimentary rocks and weathered basement rocks ( Figure 2). Z decreases rapidly between wavelengths of 150 and 250 km which constrains the elastic thickness of the plate, T e . By minimizing the difference between observed and predicted admittance, we obtained T e = 15 +15 −5 km. This estimate accounts for the presence of internal crustal loads using the method of McKenzie [2003, Table 1]. It agrees with values of 15-35 km obtained by Zuber et al. [1989] who ignored the effects of internal loading and used a coherency method even though 2 fa = 0.4. In the absence of dynamic support, Z should decrease to zero with increasing wavelength. In Eastern Australia, Z ≈ 50 mGal km −1 at wavelengths > 500 km, which is indicative of dynamic support. Similar values have been estimated in the oceanic realm and in Africa, where long-wavelength topography is maintained by sublithospheric convective circulation [Crosby and McKenzie, 2009;Jones et al., 2012]. Numerical convective experiments predict Z ∼ 50 mGal km −1 [Anderson et al., 1973;McKenzie et al., 1980;Parsons and Daly, 1983]. We assume that a thermal anomaly, T(x, z), located below the base of the mechanical lithosphere can be described using where T • determines the amplitude of a thermal anomaly. The admittance at wavelengths >250 km is described by an analytical solution of the Stokes Equation (Table 1) [McKenzie, 2010]. This solution is sensitive to the value of elastic thickness and to the location of the base of the mechanical boundary layer. For T e = 15 +15 −5 km, the best-fitting distribution of Z is obtained for a mechanical boundary layer thickness of 60-80 km (Figure 3d). This range is consistent with the depth of the deepest mantle nodules from Cenozoic volcanic fields across the Eastern Highlands [Ferguson et al., 1979;O'Reilly and Griffin, 1985]. Thus, it is possible that a thermal anomaly located within or below the thermal boundary layer of the lithosphere maintains Eastern Highland topography [cf. Parsons and Daly, 1983;Demidjuk et al., 2007]. The amplitude of this support can be estimated by dividing long-wavelength free-air gravity by an admittance of ∼50 mGal km −1 . Figure 3c shows that the estimated dynamic topography closely matches the observed topography of the Eastern Highlands. Immediately to the west, a region of downwelling could account for the loci of the Murray, Darling, and Lake Eyre basins.
The amount of topography, e, supported by variations in crustal thickness and density can be predicted by balancing a column of continental lithosphere against a mid-oceanic ridge so that CZARNOTA ET AL.
©2014. American Geophysical Union. All Rights Reserved. Mean density of continental crust e Gridded in Figure 4 kg m −3 a Density of the asthenosphere 3.2 × 10 3 kg m −3 a For air-loaded t e calculations w is set to zero. b t u + t l = average crustal thickness of Eastern Australia from Australian Seismological Reference Model (AusREM) [Salmon et al., 2013]. c Density of topography which corresponds to the admittance at short wavelengths. d Average crustal density for Eastern Australia from AusREM model [Salmon et al., 2013]. e Mean density of crustal column along transect shown in Figures 4b and 4c.
where values of parameters are given in Table 1. We assume that the thickness and density of the lithospheric mantle are constant, and use the AusREM model of Salmon et al. [2013] to define crustal thicknesses and densities. In the Eastern Highlands, we find that there is considerable residual topography which matches long-wavelength gravity anomalies and which cannot be accounted for by crustal isostasy alone. If, for example, the Eastern Highlands are compensated at the base of the crust predicted and observed topography ought to be similar. Figures 4e and 4f show that predicted topography is fairly uniform. These observations can be accounted for by differences in the density structure of the mantle. It is noteworthy that lithospheric thickness variations determined from surface wave tomography broadly correlate with topography and with the pattern of volcanism (Figures 4a-4d).

Magmatism
Eastern Australia is peppered by Cenozoic volcanism, which tracks the drainage divide and supports the notion arising from the geophysical analysis that a significant thermal anomaly occurs beneath the Eastern Highlands (Figures 2a-2b) [O'Reilly and Zhang, 1995]. Volcanism can be divided into at least two categories: mafic lava fields and bimodal central volcanoes. The volcanoes progressively young southward and record Australia's ∼70 km Myr −1 northward translation within a hot spot reference frame [Wellman and McDougall, 1974a;Knesel et al., 2008]. In contrast, the lava fields have no systematic age progression. Instead, scattered eruptions occur along the length of the Eastern Highlands between ∼70 Ma and the termination of central volcano volcanism (Figure 5a). These lavas have ocean island basalt source affinities contaminated by subcontinental lithospheric mantle [O'Reilly and Zhang, 1995]. Mafic granulitic and rare eclogitic xenoliths suggest high geothermal gradients at the time of volcanism, and widespread intrusion of mafic material at the base of the crust (i.e., underplating) [Lovering and White, 1969;O'Reilly and Griffin, 1985;Griffin et al., 1987;O'Reilly et al., 1988;Roach, 2004]. These high geothermal gradients are consistent with > 90 mW m −2 heat flow measurements in areas of recent (< 10 Ma) volcanism [Purss and Cull, 2001].

Journal of Geophysical Research: Solid Earth
10.1002/2013JB010436 have been 400 km Gaussian filtered to remove high-frequency variations. Notice elevation of Eastern Highland cannot be accounted for by crustal isostasy, implying that it must be controlled by changes in mantle lithospheric thickness and density and/or subplate thermal anomalies.  [Wellman and McDougall, 1974a]. (b) Cumulative volume of volcanism as a function of age [Wellman and McDougall, 1974a].
©2014. American Geophysical Union. All Rights Reserved.  [Farr et al., 2007] (see Table A1). Solid lines = real rivers; dashed lines = interconnected lakes and/or paleodrainage; dotted lines = artifacts of the extraction process, localized to region of internal drainage around Lake Eyre, and karst drainage in Eucla Basin.
The age of magmatic underplating is unknown. However, equilibration temperatures and pressures from phenocrysts of Cenozoic Central Highland basalts are consistent with magma fractionation at Moho depths [Ewart et al., 1980]. A combination of seismic wide-angle and receiver function profiles from the southeastern highlands reveal thick (35-55 km) crust with Vp velocities of 7.35-8.08 km s −1 and Vs velocities of > 4 km s −1 [Finlayson, 1979;Tkalčić et al., 2012]. These velocities are consistent with magmatic underplating which may have caused some of the regional uplift [Christensen, 1982;Johnson, 1989;Lister and Etheridge, 1989]. However, the pronounced crustal root beneath the Snowy Mountains probably makes a negligible contribution to uplift because of its high density ( Figure 4c).

Drainage Patterns
Longitudinal river profiles are sensitive to spatial and temporal variations in uplift [e.g., Whipple and Tucker, 1999;Schoenbohm et al., 2004;Roberts and White, 2010;Roberts et al., 2012a]. Given that it is mostly arid, Australia would seem to be an unlikely place to exploit this sensitivity. However, extensive temperate and wet rain forests prevailed until late Miocene times [Kershaw et al., 1994;Macphail et al., 1994]. In early Pleistocene (0.78-2.59 Ma) times, a transition to arid conditions coincided with a change from freshwater to evaporitic and gypsiferous deposition within lakes of the Western Plateau, of the Central Ranges and of the Murray Basin [Zheng et al., 1998;Chen and Barton, 1991;English et al., 2001;Dodson and Lu, 2005;McLaren and Wallace, 2010;Fujioka and Chappell, 2010]. In late Pleistocene times, particularly during the last glacial maximum, dune fields were established [Herczeg and Chapman, 1991;Nanson et al., 1992;Hesse et al., 2004]. The youthfulness of this aridity makes it possible to identify and map out extensive drainage systems.
CZARNOTA ET AL.  We extracted 254 river profiles of Strahler order ≥ 4 from the 90 m SRTM data set ( Figure 6) [Strahler, 1952;Farr et al., 2007]. The fidelity of recovered drainage networks is ensured by a three stage process. First, local sinks and spikes were removed. Second, a drainage network was reconstructed by calculating flow directions ( Figure 7). Finally, sets of longitudinal river profiles were extracted from this data set which best match drainage patterns imaged by satellites (Figures 7a-7d and Table A1). However, planform discrepancies of > 1 km occur in anabranching river channels on very low gradient depositional slopes (e.g., the Murray-Darling Basin). Elevations of longitudinal river profiles located within 1 km of spot heights measured by the Global Navigation Satellite System are consistently 20-30 m lower [Brown et al., 2011]. Comparison with published river profiles derived from topographic maps are also in good agreement although our recovered profiles are ∼20% longer in headwater regions [Wellman, 1979b;Young and McDougall, 1993;van der Beek and Bishop, 2003]. We attribute this difference to the superior spatial resolution of SRTM data.
The drainage network is divided into three groups consisting of present-day river channels, paleoriver channels, and artifacts ( Figure 6). Paleoriver channels comprise interconnected lake systems partially infilled by alluvial sediments (e.g., Figures 7e-7f ). Southwest of the Central Ranges and between the Pilbara and the Kimberley, paleoriver channels often consist of depressions within dune fields which constitute a drainage pattern. There are two important sources of artifacts, both of which stem from the requirement that drainage from any point must reach the coast. Thus, internal drainage basins are not preserved and all rivers drain to the coast via interconnected lakes. An extreme example is the Lake Eyre Basin where rivers drain to an elevation of ∼ 80 m above sea-level, controlled by the divide, separating Lake Eyre and Lake CZARNOTA ET AL.

10.1002/2013JB010436
Torrens, which formed during Late Miocene times . A second form of artifact occurs when rivers are forced to reach the coast despite evidence for karstic drainage (e.g., the Eucla Basin). We have included some of these profiles since their upstream tracts contain useful information.
Many Australian river profiles have irregular, convex-upward shapes ( Figure 8). Knickzones of 20-50 km wide can separate segments with lower relief and these knickzones can contain multiple knickpoints < 5 km wide. Ollier [1982] suggested that prominent knickzones along eastward-draining rivers indicate erosion of the Great Escarpment. Knickzones with similar amplitudes occur at common elevations along westward-draining highland rivers. Elevations of both east-and west-facing knickzones match the amplitude of dynamic support (Figures 8a-8j). This spatial correlation suggests that the pattern of knickzones is an expression of the growth of dynamic topography along the east coast. On the southwest side of the continent, river profiles have consistent convex-up profiles with knickzone amplitudes of 200-300 m.
It is unlikely that lithology alone controls knickzone morphology ( Figure 8). Figures 8j and 8o show that changes in lithology and patterns of Neogene faulting can generate short-wavelength knickpoints with amplitudes of ∼100 m. However, these short-wavelength features do not control the overall form of a river profile. The correlation between changes in tensile strength of channel substrate, Δ t , a proxy for rock erodibility, and both slope and curvature along a profile is shown for three rivers in Figure 9. Our results suggest that there might be negligible lithologic control along the Tumut, Snowy, and Shoalhaven rivers. In their catchment-wide study of the Lachlan River, Bishop and Goldrick [1999] concluded that lithology plays a minor role in controlling profile morphology [cf. Bishop et al., 1985]. This conclusion agrees with the results of mapping subbasaltic river profiles throughout the Eastern Highlands whose morphology is dominated by the headward propagation of knickzones [Young and McDougall, 1993;Nott et al., 1996].

Uplift Histories From River Profiles
Recent work has shown that inverse algorithms can successfully determine the spatial and temporal evolution of uplift rate history from longitudinal river profiles [Pritchard et al., 2009;Roberts and White, 2010;Roberts et al., 2012aRoberts et al., , 2012b. The rate of change of elevation along a river profile, z∕ t, is given by where z and x are elevation and distance along a given river profile, U is the uplift rate which can vary through space and time, t, and E is the rate of erosion [e.g., Rosenbloom and Anderson, 1994;Whipple and Tucker, 1999]. River profiles can be modeled by bedrock channel incision where E is given by The left-hand, advective term describes how a kinematic wave travels upstream in accordance with the well-established detachment-limited stream power model [e.g., Howard et al., 1994;Whipple and Tucker, 1999]. m and n help to determine the transient shape of a river profile. A(x) is upstream drainage area and A(x) m is related to fluvial discharge [Hack, 1957]. If n = 1 and m = 0, then v is, by definition, the knickzone retreat velocity. If n ≠ 1 and m > 0, the advective velocity becomes a nonlinear function of local slope and upstream drainage area. The second term is an approximation of the transport-limited component of erosion, where is "erosional diffusivity" [e.g., Willgoose et al., 1991;Whipple and Tucker, 1999]. This downwearing term allows synchronous erosion along the length of a river [Rosenbloom and Anderson, 1994].
The challenge is to solve the inverse problem. Given a river profile, its upstream drainage area, and values of v, m, n, and , can U(x, t) be determined? Pritchard et al. [2009] and Roberts and White [2010] used synthetic and observational data sets to demonstrate that U is recoverable if upstream drainage area is invariant. Along the Eastern Highlands, the distribution of Cenozoic lavas suggests that the crest of the Great Dividing Range and the locations of individual rivers have remained largely fixed since Early Cenozoic times [e.g., Young, 1983;Bishop et al., 1985;Taylor et al., 1985;McDougall, 1985, 1993;Nott, 1992]. Paleocurrent data from North Queensland imply that this drainage divide may have persisted since Mesozoic times [Nott and Horton, 2000]. However, large-scale drainage reversal and reorganization does occur in regions with low relief (e.g., along northwestern edge of the Western Plateau) . In these areas, we have excluded anomalous river profiles from inverse modeling.

Landscape Response Time
Rivers act as tape recorders of uplift history [Pritchard et al., 2009]. The position of a knickzone along a river can be used to infer the time of an uplift event. However, once a knickzone has propagated to the head of a river, its position can only be used to infer the minimum age of uplift. If = 0 and n = 1, the time taken for a knickzone to propagate upstream, which we refer to as Gilbert Time, G , is given by CZARNOTA ET AL.
©2014. American Geophysical Union. All Rights Reserved.  Maps of G provide a useful estimate of maximum landscape response times ( Figure 10). Thus, Australia's drainage network potentially resolves Cretaceous uplift events. Pritchard et al. [2009] showed that if = 0 and if A can be approximated by x 2 , uplift rate history as a function of time can be determined analytically ( Figure 11). This rule-of-thumb analysis suggests that the Eastern Highlands have been affected by several uplift events. The youngest event occurs at 15-20 Ma and is visible at the northern end of the Eastern Highlands. Several distinct uplift events occur between 30 and 50 Ma. The oldest uplift event has a minimum age of 70 Ma, which is at the limit of resolution.

Erosional Parameters
On the western slopes of the southeastern highlands, paleoriver channels are preserved beneath ∼21 Ma basalt flows [Young and McDougall, 1993]. The shape and location of these paleoriver channels can be used to calibrate the four erosional constants. Despite a gap of > 20 Ma, Miocene and present-day river profiles have remarkable similar gradients (Figures 12a and 12e). This similarity of form suggests that highland landscape evolution is dominated by headward propagation of knickzones. Two large rivers mapped by Young and McDougall [1993], the Tumut River and Tumbarumba Creek, are suitable for modeling since they have large erosive signals.
Advective erosional parameters were calibrated by forward-modeling river profiles for the last 21 Ma. Paleoriver channels were used as starting solutions. A systematic sweep through n, m, and v was carried out by allowing the parameters to vary between 0-2, 0-1, and 10 −5 -10 3 m 1−2m Myr −1 , respectively. The misfit between modern observed and modeled river profiles was calculated. Finite-difference schemes that implement the Crank-Nicolson scheme to solve the diffusive term and the upwind scheme to solve the advective term were used to model river profile evolution. The schemes were combined using operator splitting techniques, and the time step was set by the Courant criterion [Press et al., 1992;Roberts and White, 2010]. A(x) was measured directly from digital elevation models. We set U = 0 and = 0.01. can vary by 7 orders of magnitude without affecting our results (see section 4.3.1). In Figure 12 we show orthogonal slices through the misfit function at n = 1 together with corresponding minimum misfit values of v and m. As noted by Roberts and White [2010], there are important trade-offs between all three parameters. The 2 misfit in the direction of trade-off is shown above each panel in Figure 12. Optimal fits are obtained when n is close to 1 and when m is 0.3-0.5. Trade-off functions for each river are slightly different. Given that the Tumut River has a larger drainage area and more extensive vertical and lateral paleoriver channel constraints, we have greater confidence in the robustness of this result.
©2014. American Geophysical Union. All Rights Reserved.  There is reasonable agreement between our results and those of both Stock and Montgomery [1999] and van der Beek and Bishop [2003] (Figures 12a and 12e). There is a consensus that optimal fits between modern and paleoriver channel profiles are obtained for n = 1. Their estimates of m and v lie along our minimum misfit trade-off, except for the Stock and Montgomery [1999] estimate from Tumbarumba Creek. This discrepancy probably arises from their choice of boundary conditions. They applied a constant gradient boundary at the downstream limit of the paleoriver channel, effectively simulating a constant rate of uplift for which there is little direct evidence. van der Beek and Bishop [2003] applied five different erosional models to smaller relief paleoriver channel profiles in the headwaters of the Lachlan River. They found that best fits are obtained for a detachment-limited stream power law and for an undercapacity model which includes a term for channel width. Values of v for the Tumut River are about 4 times smaller for m < 4.5 compared with those obtained for the Colorado River by Roberts et al. [2012a]. The implication is that Australian incision rates are slower [e.g., Bishop and Goldrick, 1999].

Uplift Rate as Function of Time
An important objective is to extract uplift rate histories from analysis of longitudinal river profiles. Given the parameter trade-off relationship for the Tumut River, how well can observed river profiles be fitted and how sensitive are extracted uplift histories to parameter variation? The simplest way to analyze this problem is to examine river profiles where it is reasonable to assume that uplift rate varies solely as a function of time (i.e., block uplift with knickzone propagation from the river mouth). Given the variation of dynamic topography across the Eastern Highlands, this assumption is reasonable for eastward-draining rivers where the drainage divide lies within 200 km of the coastline (Figure 6).
Uplift rate histories are obtained using the inverse approach developed by Roberts and White [2010]. If z(x) is known, equations (5) and (6) Stock and Montgomery [1999]; triangle = best fit parameters for the Lachlan River and Wheeo Creek from van der Beek and Bishop [2003]; white circles = parameters and range used by van der Beek and Braun [1998] and van der Beek and ; dashed white line = minimum misfit trade-off for Colorado River [Roberts et al., 2012a]. (e-h) As above, for Tumbarumba Creek. River locations are shown in Figure 6.
The first term calculates 2 misfit between observed, z o i , and calculated, z c i , river elevation at all positions, N, along the river, where the uncertainty in observed elevation is i . The second and third terms control the smoothness of U(t) by penalizing the slope and curvature of U(t), respectively. U(t) is subdivided into M number of discrete time steps, Δt. U ′′ k denotes the second derivative of U(t). The final term on the right-hand side is a positivity constraint, which is applied when U < 0 where f 1 = cosh(U k ) − 1. W 1 , W 2 , and W 3 are weighting coefficients and controls the model smoothness by covarying W 1 and W 2 . A conjugate gradient technique is used to minimize H [Press et al., 1992]. In our starting model, U(t) = 0. Other parameter values are listed in Table 2.

Sensitivity Analysis
The sensitivity of calculated uplift rate histories can be assessed by systematically covarying erosional parameters. We confine our analysis to the Snowy River but the inferences we make apply elsewhere. The headwaters of the Snowy River are located on the other side of the drainage divide from the Tumut River ( Figure 6). During optimization, the mean residual misfit, 2 , between the observed and calculated river profile decreased from 31.5 to 1.4 (Figure 13a). Model smoothness and data misfit trade-off against each other and it is important to determine the amount of model smoothness which yields a good fit [Parker, 1994]. Figures 13b-13h summarizes the effect that varying the regularization parameter, , has upon data misfit, 2 , and model roughness, ||m||, where   Variations in erosional parameters v, m, and n give rise to large differences in both calculated uplift histories and residual misfit. In Figures 14a-14g, we show the result of covarying m and v for n = 1 according to the Tumut River minimum misfit trade-off ( Figure 12b). Since this trade-off relationship was established in a region with a small upstream drainage area, a decrease in m produces older uplift ages and vice versa. The sensitivity of the calculated uplift history is shown in Figures 14d and 14e. Thus, varying m from 0.3 to 0.2 while covarying v appropriately shifts uplift events back in time.
Covarying all three parameters v, m, and n along the minimum misfit for the Tumut River shows that the best fit occurs for n ≈ 1 (Figures 14h-14n). Since river slopes are less than 45 • , a decrease in n produces faster erosion rates and vice versa. can vary by many orders of magnitude without affecting our results (Figures 14o-14u). More importantly, upstream drainage area can vary by ±30% and only shift peak uplift rate by small amounts (Figures 14v-14ab). This insensitivity to substantial variations in upstream drainage area suggests that putative river capture events have a modest effect upon our principal results. Furthermore, uncertainties in estimating upstream drainage area from SRTM data should not affect the robustness of our results [Hancock et al., 2006]. In Figures 14ac-14ai, we explore the consequences of changing river length. Along the east coastline, the shelf break generally lies within 50 km of the coast and estuaries extend less than 20 km inland. Changing river lengths by these amounts can shift peak uplift rates by up to 10 Myr.

Independent Constraints
Changes in v, m, and n can yield differing uplift rate histories. It is important to identify the combinations of erosional parameters that yield uplift histories that honor independent geologic constraints. We selected four rivers from the Eastern Highlands that have independently constrained uplift histories. The head of each of these rivers lies within 200 km of the coast (Figure 15). We also identified rivers, which drain southern and western Australia, where uplift histories can be independently constrained (Figures 16 and 17).
Uplifted marine strata provide an important constraint for maximum uplift age (e.g., Northern Queensland; Figures 15q and 15w). Bands of marine strata along the coastal plain may arise from the interplay between knickzone migration and eustasy (e.g., west coast of Australia; Figures 1a, 17h, and 17l). If uplift varies only as a function of time, then knickzone initiation should be the same age as, or younger than, the onset of rifting at a margin. Fault architecture within the Gippsland Basin constrains the age of rifting, which culminates in Tasman Sea floor spreading at 80-100 Ma [Rahmanian et al., 1990;Lowry and Longley, 1991;Norvick and Smith, 2001]. Similar constraints along the southern and western margin imply that rifting occurred between 165 and 140 Ma and 160 and 133 Ma, respectively [Totterdell et al., 2003;Longley et al., 2002]. The vertical extent of basalt flows constrains the minimum cumulative uplift at the time of basalt eruption (Figures 15e, 15k, 15q, and 15w). Incised basalt flows also provide spot estimates of paleoriver channel elevation (Figures 15b and 15h).
By varying n, m, and v along a minimum misfit relationship for the Tumut River, independent constraints are most closely honored when 1 ≤ n ≤ 1.05, 0.25 ≤ m ≤ 0.35, and 2.11 ≤ v ≤ 16.78 (Figures 15-16). The optimal values are n = 1, m = 0.3, and v = 5.96 m 0.4 Myr −1 . These values result in a residual misfit that is close to the global minimum for the Tumut River (Figures 12b and 14g). A formal goodness of fit is one of several considerations when selecting geologically valid erosional parameters.
In this way, we have attempted to honor all maximum age constraints and all minimum cumulative uplift constraints. For the Shoalhaven River, however, we underpredict cumulative uplift by ∼ 200 m at 30 Ma while satisfying the paleoelevation of 30 Ma basaltic flows preserved along the rim of the Shoalhaven River gorge (Figures 15h and 15k) [Nott et al., 1996]. This discrepancy probably arises because of the 60 km separation between the river and the cumulative uplift constraints and/or because of uncertainties in K-Ar dating of the basalts [Young and McDougall, 1982]. Our inability to match all constraints from the Eucla Basin stems from the spatial variation of uplift and from ∼ 70 m glacio-eustatic sea level variations, which modulate the uplift signal ( Figure 16) [Hou et al., 2008]. On the west coast, our cumulative uplift predictions are within 50 m of the independent constraints, which is within the amplitude of glacio-eustatic sea-level variation ( Figure 17). Independent uplift constraints across the continent are honored by a single set of erosional parameters, which implies that erosion is not sensitive to climate variation.

Implications
Inverting for uplift rate as a function of time, U(t), yields insights into the regional uplift of Australia. Drainage of the Eastern Highlands is consistent with a two phase uplift history ( Figure 15). The older phase of uplift occurs between 120 and 80 Ma and has a maximum uplift rate of 0.05 mm per year. It roughly coincides with the age of rifting. A second phase of uplift occurs between 70 and 20 Ma and has a maximum uplift rate of 0.04 mm per year. It coincides with the age of intraplate Cenozoic volcanism. In areas where we have constraints on uplift from basalt flows, uplift immediately precedes volcanism. A similar relationship is observed for the Eucla Basin (Figures 16e and 16h). Gippsland Basin [Rahmanian et al., 1990;Lowry and Longley, 1991;Norvick and Smith, 2001]. Hashed polygon = minimum uplift constraints from basalt flows [Wellman, 1979b]. Inverted triangle with error bars = palaeoelevation of creeks [Holdgate et al., 2008].  ; stippled boxes = siliciclastic sediment flux in the Gippsland Basin [Weber et al., 2004]. Inset shows geological units: white = terrestrial; gray = marine; black = unconformity; v letters = volcaniclastic sediment; stippled = siliciclastic sediment; rectangles = carbonate sediment [Bernecker and Partridge, 2005]. (g-l) Shoalhaven River. Basalt flows in Figures 15g and 15k from Nott et al. [1996] and McDougall [1982, 1985]; thermochronology in Figure 15k from Persano et al. [2005] (sample 99-OZ-12). (m-r) MacLeay River. Basalts and youngest marine deposits in Figure 15q from Wellman [1987] and Langford et al. [1995] (Figure 1a). (s-x) Herbert River. Age of Coral Sea in Figure 15u from Gaina et al. [1999]. In Figure 15w A multiple uplift history of the southeastern highlands is reflected in the stratigraphic record of the Gippsland Basin (Figure 15f ). Both phases of uplift are associated with basin-wide unconformities [Bernecker and Partridge, 2005]. The second phase of uplift is associated with development of onshore and offshore low-angle unconformities across southeastern Australia and with major channel incision in the Gippsland Basin between 44 and 40 Ma [Johnstone et al., 2001;Holdgate et al., 2003].
Sedimentary flux can be predicted from river profile evolution (e.g., Figure 15f ). In the Gippsland Basin, onset of carbonate sedimentation corresponds with a decrease in sediment flux at ∼30 Ma. The maximum clastic flux estimated from isopach maps roughly coincides with peak uplift rate but does not agree with the predicted sediment flux. There is some correspondence between the maximum rate of uplift and denudation predicted from thermochronologic analyses (Figures 15q-15r)  .

Uplift as Function of Time and Space
Having established the values of erosional parameters from one-dimensional river profile inversion, we can now address the more general problem and attempt to constrain uplift rate history as a function of both time and space by jointly inverting 254 river profiles from across Australia [Roberts et al., 2012a] (Figure 6 and Table A1). First, we define an X × Y× T grid where X, Y, and T are spatial and temporal vertices, respectively (Table 2). We set the grid spacing, Δx = Δy = 500 km, and invert for a period lasting 120 Ma with a time step of Δt = 10 Myr. Optimal values of n = 1, m = 0.3, and v = 5.96 m 0.4 Myr −1 , determined from 1-D river profile analysis, are used. Figure 18 shows  along each river in the downstream direction (equation (7)). The result shows the maximum resolvable age of knickzones input at any point along the river. Oldest ages are at the river mouth and decrease toward the head. If the model is not damped, then each vertex is only sensitive to uplift ages younger than the oldest G within the four surrounding cells. Figures 18b-18e shows histograms of these ages for four vertices. Each vertex is sensitive to uplift ages back to at least 120 Ma. By damping the model, maximum resolvable age increases.
We set U(x, y, t) = 0 as the starting condition and seek the smoothest uplift history which minimizes the misfit between observed and calculated river profiles. An additional four terms are included in the trial function of equation (8), which minimize slope and curvature in both the X and Y directions. At each model iteration, uplift rate is bilinearly interpolated between vertices to obtain values along each river profile. During optimization, the residual misfit between observed and calculated river profiles, 2 , decreases from 8.9 to 2.4. An optimal smoothing value is obtained by running multiple inverse models as before (Figure 13b). Good fits were obtained for many different river profiles (e.g., Figures 19e and 19f ). The largest residual misfits occur for river profiles which drain the southeastern highlands. These discrepancies are a consequence of the coarse grid spacing, which cannot resolve short-frequency topographic variations.  Figures 19b-19q. (b-q) Gray lines = observed river profiles; numbers in brackets = river numbers in Figure 6 and cataloged in Appendix A; black lines = calculated river profiles. Notice excellent fits to rivers draining western plateau and poorer fits along the Eastern Highlands. This discrepancy arises due to the coarse vertex spacing (500 km) which cannot capture shorter wavelength elevation variations along the southeastern highlands.
The resultant cumulative uplift history, ∫ t 0 U(x, y, t)dt, is shown in Figure 20. It is important to emphasize that a priori assumptions about the spatial and temporal pattern of uplift were not included. Instead, the inversion algorithm produces the smoothest uplift history which best fits all 254 river profiles. A comparison between the final cumulative uplift and a surface envelope fitted to drainage divides shows that we can accurately resolve most major topographic elements, despite the coarse grid spacing. The resultant uplift history also honors the distribution of coastal and marine strata. The first time step is at 110 Ma and shows uplift commencing along the northeast margin. The minor discrepancy at 110 Ma between predicted topography and distribution of marine strata west of the Atherton Tablelands arises from a lack of spatial nodes on the drainage divide. In general, our results suggest that the Eastern Highlands were uplifted before central and western Australia.
In Figures 21-22, we compare the growth of topography along different swaths. Comparisons with independent uplift estimates from the western Plateau and from the northern Eucla Basin demonstrate that temporal and spatial variation in uplift rate histories produces a better fit than that obtained from U(t) models (compare Figure 21b with 17h, and Figure 21c with 16e and 16h). These results suggest that our estimates of uplift of the Western Plateau since ∼50 Ma are robust. The pattern of uplift shows diachroneity from north to south. In the Flinders Ranges, increased uplift rates since ∼30 Ma are older than the inferred 10-5 Ma onset of the Australian compressional stress field, which is considered responsible for fault-controlled uplift in these ranges [Sandiford et al., 2004]. This result suggests that there might have been dynamic uplift preceding fault-related uplift.
Regional uplift of the Eastern Highlands appears to vary through space and time. The northern highlands were rapidly uplifted between 120 and 100 Ma and were subsequently uplifted at a constant rate. The southeastern highlands, including New England, were uplifted in two phases. Uplift between 120 and 80 Ma was regionally distributed with decreasing amplitude toward the north. Uplift between 80 and 10 Ma becomes progressively younger from north to south and corresponds to the duration of volcanic activity ( Figure 5). Our results are consistent with many independent uplift constraints (Figures 22b-22e).

Discussion
Continent-wide analysis of drainage patterns suggests that present-day topography was primarily generated during the Cenozoic Era. Admittance studies and tomographic models suggest that the Eastern Highlands are dynamically supported by a thermal anomaly, which is probably located within or just beneath the thermal boundary layer of the lithosphere. This inference contrasts with other regional models, which assume that isostatic compensation occurs at the base of the crust [Wellman, 1979a]. An important question then concerns the age of this dynamic support. Correlation between the amplitude of dynamic support and the elevation of prominent knickzones on either side of the drainage divide imply that onset of dynamic support can be determined provided that the age of knickzone formation is constrained. Inversion of river profiles suggests that regional dynamic uplift occurred between 80 and 10 Ma when volcanism occurred along the crest of the highlands. Possible links between uplift and Cenozoic volcanism have previously been downplayed due to the inferred antiquity of the Australian landscape which was thought to be responding to isostatic rebound [e.g., Bishop and Goldrick, 1999]. This inference was based upon observed rates of river incision and on the minimum age on uplift constrained by basaltic flows. In our river profile modeling, we have matched both incision rates and minimum age uplift constraints while demonstrating a temporal link between uplift and Cenozoic volcanism (Figures 15). Persistence of dynamic support in Eastern Australia over the last 80 Ma indicates that the inferred mantle thermal anomaly is either aligned north-south (i.e., parallel to Australia's rapid northward trajectory) or coupled to the lithospheric plate. Finn et al. [2005] proposed that volcanism in eastern Australia and western Antarctica is contemporaneous and spatially linked by low shear wave velocity zones within the upper mantle which could be interpreted as support for the first hypothesis. However, dynamic support of eastern Australia has a pronounced right-angle bend at ∼27 • S. Immediately to the west, there is no evidence for successive waves of uplift and subsidence since a continuous succession of Cretaceous marine strata is preserved in the Darling CZARNOTA ET AL.
©2014. American Geophysical Union. All Rights Reserved.  [Lowry and Longley, 1991;Veevers et al., 1991]; dashed band = expected relationship for knickzones initiated at the coast between 100 and 80 Ma; knickzones initiated inland will plot below this box.
Basin (Figures 1 and 3c). These observations suggest that the mantle thermal anomaly is confined to the thermal boundary layer of the lithospheric plate, which is consistent with seismic tomographic evidence.
Uplift of the Eastern Highlands is not only associated with growth of present-day dynamic support. There is evidence for an early (120-80 Ma) phase of uplift along the length of the Eastern Highlands. This event correlates with rifting along the eastern margin and with thermochronologic apatite fission track ages with CZARNOTA ET AL.
©2014. American Geophysical Union. All Rights Reserved.  [Gurnis et al., 2012]. Plus and minus signs = uplift and subsidence colored by the isostatic state prior to onset of vertical motion: blue = drawn down; black = unperturbed; red = uplifted. Red triangles = volcanoes; black dots = marine sediments of labeled age or younger. long mean tracks located east of the Great Escarpment [Rahmanian et al., 1990;Lowry and Longley, 1991;Gleadow et al., 2002]. The closing temperature for apatite fission tracks is 110 ± 10 • C, which corresponds to a minimum depth range of 1.8-2.2 km, based on a geothermal gradient calculated from mantle xenoliths Cull et al., 1991]. Thus, fission track data should be insensitive to formation of the < 1 km Great Escarpment. Yet there is good agreement between 60 Ma denudation predicted using fission track data and uplift of New England (Figures 15q-15r and 22c). It is straightforward to envisage how exhumation along normal faults can result in fission track cooling ages which correspond to the age of rifting along the margin between 100 and 80 Ma. However, younger ages might record the upstream propagation of inland Journal of Geophysical Research: Solid Earth 10.1002/2013JB010436 knickzones established during rifting. There is a spatial association between elevation along the Eastern Highlands and < 80 Ma fission track ages along the coastal strip (Figure 23c). To investigate whether these ages are related to the propagation of knickzones, we compare them with the landscape response time, G . Ages of knickzones which originate at the coast should plot within the diagonal dashed band in Figure 23d whereas knickzones originating inland will fall below the band. A cluster of data following this trend suggests that the fission track data and our landscape evolution model are mutually consistent. Together, these results indicate that the initial phase of rift-related uplift and erosion locked in the observed apatite fission track ages; yet it was the second phase of uplift and erosion which exposed these samples beneath the Great Escarpment. The second uplift event is unlikely to generate a denudational signal large enough to be measured by thermochronometry. Our uplift history also corresponds to the thermal evolution of the coastal plain, determined using lower temperature apatite (U-Th)/He thermochronometry (Figures 15e  and 15k) [Persano et al., 2005].
The precise uplift mechanism of the first, rift-related, uplift event remains unclear. Lister and Etheridge [1989] suggested that syn-rift uplift and mafic underplating were generated by a subplate thermal anomaly at the time of rifting. However, mapping of oceanic crustal thickness variation abutting the eastern Australian continental margin reveals an average crustal thickness of 6.4 km, which is 0.7 km thinner than the global average. These values imply that asthenospheric mantle was cool at the inception of the Tasman Sea (84-64 Ma) [Czarnota et al., 2013]. Alternatively, there may have been a sharp gradient in asthenospheric temperature across the margin. Another explanation is that syn-rift uplift was caused by rebound from a previously drawdown state. Middle Cretaceous intracontinental flooding of eastern Australia has long been inferred to be the consequence of dynamic drawdown above a subducting slab [Gallagher, 1990;Gurnis et al., 2000;Matthews et al., 2011]. On removal of this negative mantle anomaly, the region should have rebounded to achieve isostatic equilibrium. This inference is supported by subsidence analysis combined with vitrinite and fission track data, which suggests that the central highlands underwent regional uplift with no change in upper crustal geothermal gradient prior to the onset of seafloor spreading [Raza et al., 2009]. In contrast, oceanic crustal thickness is almost double the global average within the Coral Sea, where spreading began at ∼ 62 Ma [Shor, 1967;Gaina et al., 1999]. This observation implies a thermal anomaly within the mantle grew beneath northern Australia immediately before this time, consistent with a second phase of uplift along the Eastern Highlands (Figures 15 and 22).
It is probable that rapid south-north plate motion has caused diachronous uplift of the Western Plateau and subsequent drawdown of the Northwest Shelf ( Figure 19) [Czarnota et al., 2013]. Regional uplift is generated by emergence of Australia from a region of drawdown which is consistent with offshore estimates of residual topography and with the gradient of long-wavelength free-air gravity anomalies (Figure 1a). On the Northwest Shelf, regional uplift attained an unperturbed elevation before drawdown commenced. The amplitude of this drawdown is significant, and so our estimates of the amplitude of regional uplift of an area encompassing the Kimberley is likely to be a minimum. Timing of uplift of the Western Plateau corresponds with the age of Manganese oxide formation throughout Western Australia [Dammer et al., 1999]. This correspondence suggests that development of low topographic gradients in a landscape may be critical for the growth of these supergene deposits. Kimberlitic intrusions in northwestern Australia coincide with the highest rates of uplift in the region. We attribute both processes to changes in the thermal regime at the base of the lithosphere.

Conclusions
The spatial and temporal pattern of Australian dynamic topography is difficult to constrain away from coastal strips where uplifted marine terraces and strata provide useful constraints [Sandiford, 2007]. We suggest that, when combined with local constraints for fluvial erosion, drainage patterns provide useful information. We have modeled an inventory of 254 longitudinal river profiles using a combination of inverse algorithms which permit the temporal and spatial history of uplift rate to be extracted. Our results suggest that the epeirogeny of Australia is largely a Cenozoic phenomenon. Uplift of the Eastern Highlands began at about 120 Ma and intensified during the Cenozoic Era associated with volcanism. Variations in vertical motion correlate with the sign of long-wavelength gravity anomalies which suggests that there are significant gradients in dynamic topography across the east coast of Australia. This spatial pattern of uplift does not match a simple north-south trend, although there is evidence for a north to south diachroneity. Immediately offshore, the oceanic floor is drawn down and basin sediments record anomalous Cenozoic subsidence Journal of Geophysical Research: Solid Earth 10.1002/2013JB010436 [Czarnota et al., 2013]. In central and western Australia, inversion of river profiles, calibrated using the uplift histories of marine deposits, suggests that widespread regional uplift is dominantly of Cenozoic age. A marked north to south diachroneity is consistent with Australia's rapid, northward translation over slower moving or stationary long-wavelength convective cells within the mantle (Figure 24) [Czarnota et al., 2013].

Appendix A: Drainage Inventory
Vital statistics of drainage networks used in our analysis are given in Table A1.