A Thin Elastic Plate Model for Thermally Contracting Young Oceanic Lithosphere: Insights From Comparison With Modern Seafloor Observations

We find at fast‐ and intermediate‐spreading seafloor that their ridge‐parallel bathymetric profiles between two neighboring fracture zones, excluding the part of the seafloor inward to fracture zone (FZ) valley, are predominantly upward concave. The temporal evolution of the bathymetric profiles from the lithosphere formed at the Chile Rise is characterized by (a) the rapid growth of the middle deflection to about 200 m relative to the ends for the first few millions of years and (b) a steady state afterward. We show that these characteristics and the upward‐concave sense of bending can be reasonably explained as the flexure of a thin elastic plate contracting thermally from the top while cooling. The best‐fitting model needs only about 10% of the thermal bending moment based on the half‐space cooling model and the free‐end assumption. Our model is consistent with the recent observations that oceanic lithosphere is cut open at a FZ valley, which disprove the previous assumption that ocean floor is bent down forming the valley walls.

Thermal contraction of oceanic lithosphere was invoked for explaining bathymetry data from global seafloor mapping effort unfolded in 1980s'. For instance, Parmentier and Haxby (1986, hereafter PH86) modeled the observed bathymetric profiles on the lithosphere sectored between oceanic fracture zones using a thin elastic plate model. Their decision to include the FZ valleys in the bathymetric profiles to be modeled (Figure 1a) led to the thin elastic plate model with both ends bending downwards (Figure 1a). The model implies that the internal layers of the oceanic lithosphere would conform to the FZ valley walls (Figure 1a) rather than be crosscut by them.
However, as more in situ observations and higher resolution global seafloor data have become available over the last few decades, further deviation between our understanding of oceanic lithosphere and the implications of PH86 model grew. The deviation would soon become non-negligible as our knowledge has recently been advanced about much higher resolution on interplay between tectonics and magmatism that characterizes FZ valley and  Parmentier & Haxby, 1986). The sense of bending is consistent with "growing plate" cooling model (e.g., Wessel, 1992). (b) The model from this study that excludes FZ valleys. The senses of bending moment are opposite to (a). Model A implies internal layering conforming to the FZ valleys while model B is consistent with the internal layering cut and exposed by them. The sense of bending is consistent with "static plate" cooling model (this study). the behavior of adjacent oceanic lithosphere than a few decades ago when the PH86 model was proposed (e.g., Gregory et al., 2021;Grevemeyer et al., 2021;Kolhi et al., 2021;Wang et al., 2022). With a simple view, the deviation can be described in a few key points. First, where the inner walls of a FZ expose a section of oceanic lithosphere, the inner wall lithology implies the "transverse ridge" along a FZ is a product of upward bending (Bonatti et al., 2005;Cannat et al., 1991;Gregory et al., 2021;Juteau et al., 1995;Marjanović et al., 2020;Ren et al., 2022). Then, seafloor observations reveal that geological processes about flexure evolution from axis to off-axis following the flowline of a spreading corridor is seemingly more complicated as (a) not only both sides of fracture zones are off-set of lithosphere/depths-one or both sides can be characterized by ridges and other tectonic inheritances from spreading-rates dependent mid-ocean ridge axes; and (b) incipient propagators, globally observed near-axes and off-axis seamount chains and volcanic addition make the topography of lithosphere sectored by fracture zones elevated from adjacent seafloor. With improved resolution of seafloor bathymetry since PH86, for example, the effort of multi-scale bathymetry compilation (Ryan et al., 2009) and new satellite altimetry-based global coverage (Sandwell et al., 2014), and advanced knowledge on the seafloor observations today, we revisit thermal contraction within young oceanic lithosphere, which is intrinsically related to the origin and nature of flexure, using such modern data.
In this study, we start with presenting the evolution of seafloor topography between transform faults over 20 Ma using a publicly accessible ETOPO 1-arc min global bathymetry compilation with satellite altimetry data ( Figure 2). We then calculate the flexure of an elastic thin plate driven by bending moment based on the half-space cooling model. The point of departure from the PH86 model is that we exclude FZ walls from the bathymetric profiles to be modeled (Figure 1b), motivated by the observed cross sections of oceanic lithosphere exposed on the FZ walls described above. We show that our thin elastic plate model well captures the main characteristics of bathymetric evolution within a case study area, the Chile Rise oceanic lithosphere, despite limitations stemming from the simplifying assumptions we made.

Seafloor Observations
As being inspired by the approach PH86 had taken in assessing their model (e.g., Figure 4 in Parmentier & Haxby, 1986), we made a first-order assessment on the maturation/evolution of oceanic lithosphere along the age flow line between two neighboring fracture zones. Although "long-lived" fracture zones are ubiquitous in the world's ocean, there are only a few that are suitable to use as a basis for our modeling. The lithosphere we should focus on are: (a) fast-to-intermediate spreading lithosphere bounded by two continuous transform faults/ fracture zones observed within the ETOPO bathymetry compilation grids; (b) crustal ages and spreading rates (i.e., fast-∼74 mm/yr half rate and intermediate 28-34 mm/yr half rate, e.g., Buck et al., 2005) were determined by unambiguous interpretation of observed marine magnetic anomalies in previously published literature; and (c) the extent of the ridge to off-axis reaches up to 20 Myr in age and is as little tectonically and volcanically disturbed/overprinted as possible. Many of these lithosphere segments include propagator wakes (particularly pronounced in the Mid-Atlantic Ridge segments and their transform faults in the South Atlantic), and volcanic ridges and seamounts, being found unideal to capture our numerical model parameterization and results. Consequently, there is a very limited number of seafloor segments available for our study. In addition to Chile Rise, both east and west flanks of fast spreading East Pacific Rise 9N (110 mm/yr full rate) and eastern flank of the Udintsev FZ (76 mm/yr full rate, and the profiles examined in PH86) are assessed for reference (Figures S1 and S2 in Supporting Information S1). To confirm the segments are mostly free of major structural overprint, we have examined vertical gravity gradient data as our check point (e.g., Sandwell et al., 2014) as well.
Ridge-parallel profiles (X-X′ in Figure 1) were extracted in between the edge of the lithosphere dissected by Valdivia and Guano fracture zones (Figures 2c and 2d). Ridge-parallel profiles of depth and free-air gravity anomaly were obtained with 0.2 km data sampling frequency along the profile and with 2 km spacing along the spreading flow-line from the 0 (ridge) to ∼20 Ma age crust estimated from their spreading rate (Tebbens & Cande, 1997). To extract representative profiles of the bathymetry and gravity data to compare with numerical results, we: (a) group the sampled depth and gravity profiles into the five time periods, 0-3, 3-6, 6-10, 10-15, and 15-20 Ma; (b) smoothed each profile by convolving it with a linear kernel (Text S1 in Supporting Information S1); and (c) fitted the smoothed profiles to a polynomial ( Figures S3 and S4 in Supporting Information S1). We did not explicitly consider the effects of sediments (see Text S2 in Supporting Information S1).

Application of Thin Elastic Plate Model
Toward building our new model, we adopt a few simplifying assumptions for approximating lithosphere as a thin elastic plate that were made first in PH86 and then used in subsequent studies (e.g., Haxby & Parmentier, 1988;Parmentier & Haxby, 1986;Wessel, 1992;Wessel & Haxby, 1990). First, newly-formed oceanic lithosphere is assumed to freely contract in the ridge-parallel direction (x) as it cools down. Second, the cooling and freely-contracting lithosphere can still bend in a consistent way with the deviatoric strain, − , where is the depth average of . If bending is prohibited due to kinematic constraints on the lithosphere, the corresponding bending stress would arise. Third, the lithosphere is assumed to be elastic when its temperature (T) is at or lower than a brittle-ductile transition temperature, T l . Thermal stress, if any, would dissipate completely and sufficiently fast by viscous relaxation when T > T l . Finally, depth-dependent of temperature within the elastic portion of lithosphere is assumed to be linear, an approximation well-justified in the light of the half-space or a plate cooling model. With elastic thickness at time t denoted as h(t), the depth (z) distribution of temperature within the elastic portion is given as T(z, t) = T l z/h(t).
The above assumptions yield an expression for ridge-parallel thermal strain ( ) : where is the linear thermal expansion coefficient. This expression holds for 0 ≤ ≤ ℎ( ) for any ( ) and is zero or negative (i.e., contractional) for all depths, 0 ≤ ≤ ℎ( ) . The mean horizontal strain, , is defined as 10.1029/2023GL103511 5 of 10 For the linear top-down cooling considered here, becomes a constant, Then, we get deviatoric horizontal strain, ′ from Equations 1 and 2: The bending moment (M) required for a flat, stress-free, thin elastic plate to generate the flexure compatible with the deviatoric thermal strain given by Equation 3 can be expressed in terms of the deviatoric stress, ′ : where ′ is ∕(1 − ) and E and are the Young's modulus and the Poisson's ratio (Turcotte & Schubert, 2014). From Equations 3 and 4, we get the following expression: The moment in Equation 5 generates upward-concave bending of the plate (Figure 1b). Upward concavity is the sense of bending expected for the "static plate" cooling from the top. The "static plate" model is one of the end-member models for how thermal strain accumulates in oceanic lithosphere over time while cooling (e.g., PH86; Wessel, 1992). Another end-member is the "growing plate" cooling model (Wessel, 1992), in which the brittle portion is assumed to have contracted freely before a new layer that has just become brittle is added at the bottom and starts contracting. The sense of bending in the growing plate model is depicted in Figure 1a. Although the way thermal strain actually accumulates in cooling oceanic lithosphere must fall between these end-member modes, we adopt the "static plate" model in this study because it is consistent with the sense of bending seen in the bathymetric profiles ( Figure 2).
The equation for a thin elastic plate with the isostatic restoring force and without any horizontal boundary force is where ( ) is the flexural rigidity, ( ) 3∕ ( 12 ( 1 − 2 )) (e.g., Turcotte & Schubert, 2014;Watts, 2001), ρ m and ρ w are the mantle and water densities, and g is the gravitational acceleration.
The bending moment given by Equation 5 is an overestimation because the deviatoric thermal contraction cannot occur freely in oceanic lithosphere because of the water pressure and the resistance to bending at the FZ boundaries. Thus, we introduce effective bending moment and flexural rigidity, ( ) and ( ) , using them for solving Equation 6. With a misfit defined as ‖ obs − ‖∕‖ ‖ with ‖ ‖ being the Euclidean norm of a vector , the optimal values of ( ) and ( ) are those that minimizes the sum of the misfits of the Chile Rise profiles for the five time periods; and they are 0.09 ( ) and 0.9 ( ) , respectively (Text S3 in Supporting Information S1). The coefficients, 0.09 and 0.9 are constant and the cooling effects is introduced through ( ) and ( ) .
We apply the boundary conditions for the pinned ends: (0, ) = 0 and ′′ (0, ) = (t)∕ ( ) on the left (i.e., x = 0); and ( ) = 0 and ′′ ( ) = (t)∕ ( ) on the right (i.e., x = L) end, where L is the length of the elastic plate. The details of acquiring approximate solutions to (6) under these boundary conditions are provided in Supporting Information S1 (Text S4). Another possible set of boundary conditions is to set ′′′ instead of to be zero. The qualitative characteristics of the modeled flexure are retained with either choice of boundary conditions.
In summary, the best-fitting deflections are those of a 200 km-long thin elastic plate with pinned ends for = 500°C, ν = 0.25, E = 100 GPa, ρ m = 3,300 kg/m 3 , ρ w = 1,030 kg/m 3 , g = 9.8 m/s 2 , ( ) = 0.09 ( ) , and ( ) = 0.9 ( ) . The modeled deflection profiles are averaged over each of the time intervals used for the bathymetric profiles (see Section 2.1). To facilitate the comparison of the thin plate deflections with the bathymetric profiles, we vertically adjusted the bathymetric profiles such that their left-end depth becomes to 0 m (Figure 3).

Results
Ridge-parallel bathymetry profiles from the fracture-zones bounded lithosphere formed at Chile Rise show temporal changes in their concavity for crustal ages less than 10 Ma and then an apparently steady state over time older than the 10 Ma crustal age. To best capture the first-order lithosphere curvature evolution over time, we have divided timeseries of lithosphere evolution to 0-3, 3-6, 6-10, 10-15, and 15-20 Ma spans of crustal age group (Figure 3). During the first 0-3 Ma, we observe slight upward "bulge" in the middle of the averaged profile, seemingly superposed on a broadly upward-concave profile (Figure 3a). The older (3-6 and 6-10 Ma) age seafloor exhibits reduced magnitude of the middle bulge (Figures 3c and 3e). The middle part becomes about 200-220 m deeper than the ends. The 10 Ma or older seafloor shows upward concave profiles without visible perturbations. The depth difference between the middle and the ends of the profiles stays near 200 m (Figures 3g Figure 3. The best-fitting thin plate model's deflections (black, left column), the polynomial-fitted bathymetric profiles (red, left column), and the polynomial-fitted free-air gravity anomaly (ΔFAA) profiles (blue, right column) from the Chile Rise for (a, b) 0.1-3, (c, d) 3-6, (e, f) 6-10, (g, h) 10-15, and (i, j) 15-20 Ma. Note that the center of flowline has been shifted over time due to the curvature of the long-lived fracture zones. Each of the Chile Rise profiles are smoothed with filtering processes described in Section 2.1. 7 of 10 and 3i; Figure S2 in Supporting Information S1). The averaged and smoothed profiles of free air gravity anomaly (FAA) also show clear upward concavity with the difference between maximum and minimum (ΔFAA) of 10-15 mGal for the first 10 Ma (Figures 3b, 3d, and 3f). The upward concavity is discernible but less pronounced in the gravity anomaly profiles for later time periods (Figures 3h and 3j) with the end-to-center difference decreasing to 5-10 mGal.
The maximum downward deflection of the best-fitting thin elastic plate is about 120 m in the 0-3 Ma profile (Figure 3a) but increases to about 240 m in the next 3 Ma interval (Figure 3c). The deflection profile barely changes for the following 4 Ma (Figure 3e) but after 10 Ma, the downward deflection decreases to about 200 m in the 10 to 15 Ma interval (Figure 3g), and to 180 m in the 15-20 Ma (Figure 3i). The origin of the relatively large misfits of 0.2-0.4 (defined in Sec. 2.2) for the first three time periods is a middle bulge seen in Figures 3a and 3c that is most likely attributed to underlying on-and off-axes magmatism (Figures 2a and 2b). Also identifiable in the gravity anomaly profiles (Figures 3b and 3d), this short-wavelength feature is not accounted for in our model. The thin plate deflections show much better fitting to the bathymetric profiles for 10-15 and 15-20 Ma. It is notable that the misfit between the 10-15 Ma profiles is only about 0.02.

Discussion
The temporal evolution of the bathymetric profiles is characterized by the rapid initial deepening to about 200 m (Figures 3a and 3c) and the profile reaches a more or less steady state after 6 Ma (Figures 3e, 3g, and 3i). This observation is seemingly ubiquitous in lithosphere segments bounded between fracture zones ( Figure S2 in Supporting Information S1). Assessing this temporal evolution of lithosphere deflection provides an important insight on the validity of our model and guides us to future effort. The magnitude of the deflection observed in other fast-intermediate spreading lithosphere changes similarly with time ( Figure S2 in Supporting Information S1), suggesting that the flexure is "frozen in" at the crustal age no older than 6 My where thermal gradient within lithosphere (e.g., 1,250°C isotherm) becomes steepest to moderate (e.g., Audhkhasi & Singh, 2022). This steady state seen in the 10 Ma or older bathymetric profiles partially is also found in the temporal evolution of modeled deflection. The central deflection low reaches ∼220 m in subsequent 3-6 Ma time period after it reaches ∼100 m, about half of the full deflection amount, in the 0-3 Ma time period. Afterward, the magnitude of the maximum deflection slightly decreases with time but stays around 200 m (6 Ma∼) (Figure 4a). The brittle thickness of the lithosphere, h(t), increases as the square root of time according to the half-space cooling model adopted in this study. The magnitude of bending moment and the flexural rigidity are proportional to h 2 (t) and h 3 (t) (Figure 4b); and to t and t 3/2 (Figure 4c). As a result, the to ratio decreases as t −1/2 ( Figure 4d). As this ratio determines the deflection magnitude in our thin plate model, the deflection magnitude also decreases in time as seen in the periods after 6 Ma in our model.
One implication of the to ratio changing as t −1/2 is that the ratio eventually approaches zero and thus the deflection magnitude also approaches zero although slowly. However, such behavior is not necessarily expected for the topographic expression in older lithosphere because the origin of the upward concave shape is the product of the balance between deviatoric thermal contraction and resistance to it, not bending moment as we assumed for convenience in this study.
The upward concavity in the gravity anomaly profiles corroborates the validity of our elastic flexure model. Since not isostatically compensated, the upward-concave flexure of an elastic plate by thermal contraction must produce free-air gravity anomaly of a self-similar pattern. When applied to the central low of ∼200 m in the bathymetric profiles, the Bouguer plate correction formula, 2 GΔρΔh, yields ΔFAA of about 15 mGal for Δρ equal to the density difference between basalt (2,800 kg/m 3 ) and seawater (1,030 kg/m 3 ) and Δh equal to 200 m. ΔFAA is about 10 mgal for 0-3 and 6-10 Ma (Figures 3b and 3f) and about 15 mgal for 3-6 Ma (Figure 3d). While the density structure needs to be better understood, these similar magnitudes of ΔFAA are consistent with our view that the 200-m amplitude of the upward concave bathymetric profiles is a product of thermal contraction of oceanic lithosphere as depicted in Figure 1b.
The diminished magnitudes of ΔFAA for 10-20 Ma is apparently puzzling in the light of the deflection magnitude of 200 m maintained for these ages. However, it can be understood in terms of the upward continuation effects on the gravity anomaly profiles. The sea-level gravity anomaly profiles would be the upward continuation of those due to the uncompensated seafloor topography. The low-passing effect of upward continuation also reduces the overall amplitude, and the reduction magnitude is proportional to the seafloor depth. Because the seafloor depth increases with age, so does the amount of reduction in the gravity anomaly amplitude. Another possible explanation for the diminished gravity anomaly amplitudes would be geological processes such as the density-lowering alteration of oceanic crust due to hydrothermal circulation. The alteration near the FZ valleys is expected to be more extensive than in the middle of the lithospheric section. Accumulated over a sufficiently long time (i.e., >10 Ma), the density-lowering alteration occurring more extensively near the fracture zones might counteract the gravity anomaly produced by the upward concave bathymetry.
The middle bulge seen in the Chile Rise bathymetry profiles for ages younger than 10 Ma (Figures 3a, 3c, and 3e) can be attributed to the combination of relatively thin crust over melt-rich region within relatively young, warm lithosphere, which is close to the spreading axis generating high isotherm gradient toward off axis (e.g., Cochran & Buck, 2001). Comparing with other lithosphere portions that are formed at fast to intermediate spreading rates, the curvature in young (0-3 Ma) lithosphere can be unique to their underlying magmatism. The lithosphere curvature is clearly influenced by the nature and distribution of many off-axis seafloor products, including seamount chains and incipient rifts although magnitudes and upward concavity of the curvature seem to be very similar (Figures S1 and S2 in Supporting Information S1).

Conclusions
The ridge-parallel bathymetric profiles from the 0-20 Ma fracture zones-bounded lithosphere formed at the Chile Rise have an upward concave shape of an amplitude of about 200 m. A thin elastic plate going through thermal contraction as illustrated in Figure 1b can explain this observation. The free-air gravity anomaly is lower in the middle of the profiles for the Chile Rise lithosphere younger than 10 Ma by about 10-15 mgal. The magnitude is consistent with the Bouguer correction magnitude for the 200 m bathymetric low and the density difference between basalt and seawater. The consistency between the bathymetric and gravity anomaly profiles supports our view that oceanic lithosphere bounded by fracture zones can develop uncompensated upward concave flexure by thermal contraction.
A logical next step upon our study would be to construct a comprehensive thermomechanical numerical model to release the major simplifying assumptions made in this study and originally in PH86: for example, the completely free shortening by mean thermal strain at all ages. This effort of refining currently widely accepted half-space cooling model and more accurately understanding crustal heat extraction via hydrothermal cooling (e.g., Prigent et al., 2020), coupled with recent observations that highlights complicated thermal conditions within lithosphere (e.g., Audhkhasi & Singh, 2022), will advance our knowledge on thermal properties of lithosphere.