The Subhercynian Basin: An example of an intraplate foreland basin due to a broken plate

The Late Cretaceous, intraplate shortening event in Central Western Europe is associated with a number of marine basins of relatively high amplitude and short wavelength (2-3 km depth and 20-100 km width). In particular, the Harz Mountains, a basement uplift on a single, relatively steeply dipping, basement thrust, have filled the adjacent Subhercynian Cretaceous Basin with their erosive product, proving that the two were related and synchronous. The problem of generating subsidence of this general style and geometry in an intraplate setting is dealt with here, by using an elastic flexural model con5 ditioned to take account of basement thrusts as weak zones in the lithosphere. Using a relatively simple configuration of this kind, we reproduce many of the basic features of the Subhercynian Cretaceous Basin and related basement thrusts. As a result, we suggest that overall, these basins share many characteristics with larger scale, foreland basins associated with collisional orogens on plate boundaries.

basins in compressional settings on earth. One of the first people to recognise the distinct character of these Late Cretaceous basins (Voigt, 1962) also used the term "Vortiefe" (foredeep), although this does not refer to the much later concept of foreland basins, and Voigt (1962) discussion of subsidence mechanisms is limited to an analogy with much larger, pre-plate tectonic "geosynclines". 15 Although more recently, a number of different mechanisms (Voigt et al., 2008;von Eynatten et al., 2008) have been associated with SCB formation and to a lesser degree, subsidence, an explicit analysis of the mechanics of elastic flexural bending has not been presented.
2 Foreland Basin character of the SCB Much recent work on the SCB has focused on the geometry and unconformities of the Late Cretaceous basin fill. These have 20 been explained by progressive development of syn-sedimentary, partly detached cover folding above an incipient, Harz basement thrust (Voigt et al., 2004(Voigt et al., , 2006von Eynatten et al., 2008). The SCB is underlain by a Mesozoic cover of Triassic, Jurassic and some Early Cretaceous material, which in turn overlie Zechstein evaporites (Fig. 3). The latter forms a detachment to the underlying Hercynian basement. This cover was also present across the Harz Mountains, prior to uplift, but has since been removed by erosion, estimated at as much as ∼6km (von Eynatten et al., 2019). The Mesozoic cover beneath the SCB forms 25 a passive structural marker, recording the net effect of all Late Cretaceous deformation, and shows a strongly asymmetric, synclinal geometry, with a near vertical, sometimes overturned limb along the HNBF (see also appendix), and a very gently dipping, monocline, locally disturbed by the previously mentioned internal basin structures, towards the north of the Harz. As a result of this deformation, Late Cretaceous deposits of the SCB are sometimes strongly affected by syn-sedimentary deformation. Permian Rotliegend Basins lying unconformably on the Hercynian Basement but below the Zechstein detachment, in 30 turn form a passive marker of post-Hercynian deformation but excluding Late Cretaceous detachment folding of the Mesozoic cover (Fig. 3a).
As is typical of many foreland basins, the infill of the SCB shows a progressive syn-sedimentary deformation synchronous with subsidence. In the case of the SCB however, traces of the characteristic sedimentary geometries associated with classical foreland basins such as progressive onlap further towards the foreland, or tilting of older units and thickness changes as the basin grows in magnitude with increasing load, Homewood et al., 1986;Burkhard and Sommaruga, 1998) may be masked by a combination of erosion and internal deformation. It is however interesting that Turonian and Lower 5 Coniacian marly limestones and micrites, thicken somewhat towards the southern margin of the SCB, perhaps indicating a slightly earlier, probably submarine, initiation of subsidence and thrusting, than is generally assumed based on the first arrival of eroded material from the uplifting Harz Mountains (von Eynatten et al., 2008). A broader insight into the distribution of subsidence can, however, be deduced from the isopachs of Late Cretaceous sediments, which show a pronounced depocentre up to 2500 m deep, directly adjacent to the central portion of the HNBF in the area of Wenigerode (Fig. 1), part of the ∼ 40 10 km long main SCB.
Systematic change in sediment type throughout a basin's history is also an often cited characteristic of foreland basins. It has been argued that early sedimentation as orogenic wedge development begins adjacent to the basin is characterised by relatively deep water sediments, often turbiditic in character, with little or no clastic input due to the underdeveloped initial thrust wedge and lack of significant erosional product . Such basins are sometimes referred to as "underfilled" (Sinclair,15 1997), as for example the earliest stages of the Swiss Molasse sensu lato in the form of the Helvetic flysch, and sensu stricto in the form of the Lower Marine Molasse (Burkhard and Sommaruga, 1998;Homewood et al., 1986). Later, as the orogenic load develops, sedimentation becomes increasingly clastic dominated and fills with eroded product of the orogenic wedge.
Eventually, sedimentation may move from marine to lacustrine as the basin becomes overfilled, and sedimentary bypassing may begin. In the SCB, a relatively deep water to shallow water transition during basin development could be argued for if 20 the Turonian (deeper water, pelagic limestones, no clastic input) and Coniacian (marly limestones, with some clastics) (Voigt et al., 2004) are considered to be related to the very earliest stages of basin subsidence and flexure. Nevertheless, such an interpretation is speculative and the main thickness of the SCB fill, when erosive product from the emergent Harz is clearly present is of Santonian age and later, and consists of shallow water facies, generally interpreted as tidal plain and estuarine deposits (Voigt et al., 2008). 25 3 Flexure of plates, flexural amplitudes, wavelengths and loads Gunn (1943a, b) developed the first models for elastic bending of the lithosphere as part of wider studies of gravity anomalies (Watts, 2001). Quite remarkably in the pre-plate tectonic era, Gunn created end member models of elastic bending of "plates" both continuous and broken with both point and distributed loads (Fig. 4 a and b). For instance Gunn (1943a) modelled the Hawaiian island chain as a distributed load sitting upon an unbroken portion of the Pacific plate. Gunn (1947) also modelled 30 flexure of oceanic lithosphere in subduction zones, using an end-loaded plate with a "free" end, again without realising the plate tectonic implications of his work (Fig. 4 c). Later, as the idea of lithospheric plates developed, studies of flexural bending continued, (Walcott, 1970a, b, c) with the aim of better determining the magnitude and temporal evolution of elastic paramters of the lithosphere like flexural rigidity and elastic thickness.
In collisional belts and orogens such as the Alps, the elastic bending causing foreland basin development on their margins, due to the vertical loading of the underthrust continental lithosphere by an overriding stack of nappes or thrust sheets  has been often studied by using some kind of edge-loaded plate with the load set to correspond to the force exerted 5 by the weight of the orogenic wedge (Fig. 4 c). This type of model makes intuitive sense in these cases, since there are two plates interacting and the edge of the underthrust plate is depressed and bent into an elastic curve which is thought to be the main mechanism of foreland basin subsidence (Burkhard and Sommaruga, 1998;Homewood et al., 1986).
These elastic flexural end member models also neatly illustrate the problem posed by an intraplate thrust and associated basin. In the case of the Harz and the SCB for instance, the Harz basement uplift both represents the load for, and at the same 10 time is part of the same plate as the SCB (Fig. 4 d and e). Effectively, this means the plate is loading itself. The fundamentally different response of edge-loaded (or broken) plate models compared to intraplate-loaded, continuous plate models is also important in this context. Elastic flexural models assume that the lithosphere has a finite elastic thickness (h o ) which supports loads applied locally over a wider area by transferring the resistance to the load into more widespread elastic bending of the plate. In this context, an unbroken plate loaded in its middle by a relatively small load would usually be characterised by a 15 broad, long wavelength, and shallow, low amplitude flexural response, forming two similarly proportioned sedimentary basins either side of the load due to flexure. By contrast, an edge-loaded plate supports the load in a cantilever type of configuration, loaded at one end, and clamped at the other, thus maximising the moment of the load and lithospheric bending (cf Gunn (1943b)).
In this respect, the subsidence profile of the SCB is skewed towards the end load or broken plate situation. The question then naturally arises, can we imagine an end loaded, intraplate basement uplift to be a natural analogue for the HNBF and SCB?

Broken Plates
At the present day, as shown by the DEKORP deep seismic profile (Group, 1999) as well as global and local GPS and seismic 25 data (Sella et al., 2002;Tesauro et al., 2006), the Central Western European lithosphere, in the region of the Late Cretaceous intraplate shortening structures, is an integral, stable part of the European portion of the Eurasian plate. Viewed from today's perspective, there is no reason to assume anything other than a continuous, elastic, lithospheric underpinning for the region, including that of the Harz and SCB.
In the Late Cretaceous however, the situation was quite different with active, basement thrusting affecting a portion of the 30 crust and lithosphere. The geometry and magnitude of shortening on the HNBF are difficult to constrain due to a variety of factors, mostly related to erosional loss of hanging wall marker beds, and hence geometric constraints on the hanging wall itself. One of the few attempts to do so (Tanner and Krawczyk, 2017) used kinematic forward modelling with an inclined shear 4 https://doi.org/10.5194/se-2020-185 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License. algorithm and obtained a ∼ 55 • dipping fault with a listric geometry reaching horizontal at a depth of 23 km. We present a new kinematic, forward model of the HNBF (see appendix) using a fault-parallel flow deformation algorithm because it most closely matches the hypothesis that basement uplifts represent rigid blocks rotating on faults with circular arc geometry (Erslev, 1986). We consider this more realistic than the simple shear algorithm used by Tanner and Krawczyk (2017) because the simple shear assumption requires pervasive shearing of the Variscan basement at high angles to pre-existing fabrics such as bedding 5 folded about SW-NE-trending axes or SW-NE-striking cleavage planes. We obtain an estimated initial dip of 50 • , with a listric geometry, reaching horizontal at 17 km depth. Our model suggests that during a ∼ 10 Ma pulse of tectonic activity the HNBF moved at an average rate of ∼ 0.7 mmyr −1 . The fault geometry also implies that the HNBF transects most or all of the elastic portion of the crust within a horizontal distance of ∼ 15 km. Principles of conservation of line length, area or volume also mean that shortening above the basement detachment level must be accommodated somewhere through the rest of the crust 10 and lithosphere. This may happen locally on a lower ramp joining the main basement detachment from below and inducing a finite offset of the lithosphere, or in a more diffuse, ductile, pervasive shortening of the lower crust and mantle lithosphere by an amount equal to thrust displacement (Fig. 4d). Interestingly the DEKORP deep seismic profile shows the Moho deepening by several kilometres in the vicinity of the HNBF (Group, 1999), which could at least partially be attributed to lower crustal ductile compensation of basement uplift induced shortening. 15 Given these new results it seems reasonable to consider the HNBF in the Late Cretaceous, to have formed an elastic break in a broken lithosphere, effectively meaning that intraplate shortening on the HNBF resulted in a broken plate, which could for instance be modelled as an edge-load, cantilever type system. In fact in a different context, Marotta et al. (2000) used an end load configuration to model the large scale elastic bending of the lithosphere as the cause of Late Cretaceous to present evolution of the North German Basin (NGB). This model implies a broken plate configuration (at the point of the HNBF) but 20 was used principally to explain the present day Moho geometry to the north of the Harz as a product of present day elastic stresses.
Here, we propose a simple modification of the classic, elastic flexure model to allow for large and sudden discontinuities in elastic thickness within a plate. We use this to simulate Late Cretaceous basement thrusts and associated basins as discrete, short wavelength elastic flexures of a broken plate. Breaks in the plate are simulated in the flexure model by sudden reductions 25 of the elastic thickness to very low values (<1 km). We use a specific finite difference formulation referred to as a "half-station" method (Cyrus and Fulton, 1968) and apply it to the 4th order differential equation for elastic flexure of the lithosphere including a coefficient of flexural rigidity, D(x), itself a function of a elastic thickness h(x) varying with position. The numerical method is described in detail in the appendix.
Using this scheme, we can apply arbitrary breaks at any point in a "continuous" plate and thus simulate the possible effects 30 of basement faults on the flexural response of the plate to loading. Broadly speaking, this type of plate behaves as follows.
If at some point in the plate there is a narrow zone of close to zero elastic thickness, a load applied entirely one side or the other of this weak zone, will make the loaded segment of the plate behave like an end-loaded plate. The segment of the plate on the opposite side of the elastic break, to which no load is applied, remains unaffected by the load. In detail it should be pointed out that in standard flexure models, end-loading of plates is usually achieved by a so-called "line load" applied as a shear force boundary condition which concentrates the entire downward force of an orogenic load at a single point. In reality, this load is distributed over a finite width of 10's or 100's of kilometres (Fig. 4c). In our model however, such loads are applied to a plate interior with an elastic break in it (Fig. 4d), and are numerically represented as distributed loads of a width that can be set to match any load geometry. The loads correspond to hanging wall material of a basement thrust. We assume that hanging wall material above the thrust lacks any elastic strength due to plastic deformation during thrust emplacement which is 5 suggested by the overturned forelimb of the Harz hanging wall material. Hence, the additional material acts like a superposed load such as is imagined for the volcanic Hawaiian island chain, and the hanging wall material is not able to cause an additional, elastic bending-related force (Fig. 4d). If there are a succession of basement thrusts in a plate interior, it is possible to some degree to simulate these as multiple discrete elastic breaks separated by zones of normal elastic thickness with multiple loads corresponding to associated basement thrusts. 10

Flexure models of the broken Central European lithosphere
Our numerical models allow two principle parameters to be varied. The first and most important is the load, which can be defined at every grid point along the model plate. The term "load" refers to a force acting at any point along the length of the plate. These forces are often implicitly understood to be due to things like the sedimentary infill of a basin created by flexure of the plate, or due to additional material placed on part of a plate by overthrusting and stacking of the crust or by wholesale 15 overriding of one plate by another. In these cases, the local force will be equal to the product of the gravitational constant, density and a "load" thickness, be it sedimentary thickness in a basin or the thickness of an orogenic wedge. However, we can also imagine loads due to elastic stress across faults as was actually discussed by Gunn (1943b). Nevertheless, for a flexure model, any force, regardless of its origin, can be expressed as an equivalent crustal or sedimentary load which can even be negative in the case of erosional removal of material.

20
The second fundamental parameter that may vary in our models is the effective elastic thickness, which is usually indirectly determined by a combination of flexure modelling and gravity data (Walcott, 1970a). Our model defines a value of effective elastic thickness at each node along a plate, with abrupt variations (from "normal" lithospheric values to close to zero in the space of a few nodes), considered to represent plate weakening by faults transecting at least the elastic crust.
These are the two variables we use to fit our models to the HNBF / SCB scenario. All models assume a horizontal, plate wide 25 compressive "force" of 1e11 N m −1 , approximately in line with generally assumed plate wide forces from modelling studies of the present day in Central Europe (Jarosiński et al., 2006). We model a 1000 km finite segment of plate with fixed ends. The default elastic thickness used is 20 km.
The regional cross section extending from approximately the Main to the Elbe (Fig. 3) shows the tight, asymmetrical, geometry of the SCB. It is noticeable how the SCB (Late Cretaceous) is mirrored in the Rotliegend, which lies below the 30 Zechstein detachment that accommodates shortening in the Mesozoic cover below the SCB, and is effectively a passive marker for flexure and subsidence of the SCB itself. The Rotliegend is unaffected by detachment folding in the younger units above the Zechstein.
The cross section also shows a number of significant, Late Cretaceous basement thrusts besides the HNBF. These are, from south to north, the Thuringian Forest, the HNBF itself, a possible blind basement thrust below the Huy anticline, and the Allertal Fault, Haldensleben Fault and the Gardelegen Fault which partly bound and dissect the Flechtingen High. All of these structures appear to transect most of the upper crust even if their precise geometry at depth is uncertain (Group, 1999;Scheck et al., 2002;Kossow and Krawczyk, 2002). Hence, a number of elastic breaks in the European lithosphere besides the HNBF 5 may be appropriate for elastic flexure models.
Basement thickness and Moho depth also appears to vary along the south-north transect, at least in the region covered by the DEKORP deep seismic line Group (1999). Basement thickness is significantly reduced beneath the main part of the North German Basin (NGB) relative to its northern and southern terminations. It has been suggested that the shallower Moho and thinner crust beneath the NGB may explain its unique, Mesozoic and Cenozoic subsidence history (Group, 1999). However, 10 the precise mechanism of Mesozoic and in particular Cenozoic subsidence of the NGB remains a subject of discussion (Voigt et al., 2008;Malz et al., 2020). Perhaps most importantly for the question of the SCB and other Late Cretaceous basins, it seems reasonably clear that the NGB is mechanically distinct and hence behaves differently from its margins throughout the Mesozoic and Cenozoic periods (Malz et al., 2020;Voigt et al., 2008;Nielsen et al., 2005).
In the region from the SCB to the Flechtingen High, the Rotliegend, a relatively robust passive marker of Late Cretaceous leads to a south dipping, relatively straight segment which again, matches the Rotliegend geometry well. We make no attempt in this model to account for the NGB behaviour.
A second model is shown in Fig. 6. Here, we add an additional upward force whose origin is unspecified, to try to match the second Rotliegend kink south of the Flechtingen High. We find that an upward force equivalent to 5000 m of crustal material spread over 3 km distance, is sufficient to quite well match the Rotliegend geometry. As was already discussed, any force 5 can be applied as an equivalent crustal load, in this case a negative one implying removal of crustal material. Clearly, there is no evidence for any such erosional event, meaning that if this additional load is the explanation for the Rotliegend geometry then it must have a different origin. Possible candidates could include fault related elastic stresses, or deeper crustal anomalies, although the latter would be very difficult to imagine in such a way that it acted on a zone of only 3 km width at the surface. It is possible that other, more distributed load geometries could have a similar effect, and these might be more compatible with 10 deeper anomalies. However, we have no clear explanation for what such a load represents, but note that applying it allows a very good match of the Rotliegend geometry.
6 Implications of the model Our models suggest a broken European lithosphere played a key role in allowing intraplate, compressional basins to form during Late Cretaceous, intraplate shortening in Central Western Europe. In order to examine how this lithosphere behaves 15 elastically, we can compare it to a "normal" unbroken lithosphere subjected to the same loads. Figure flexure-unbroken shows an unbroken lithosphere with the same Harz Mountains and Flechtingen High loads as were used for the preceding models of the SCB. The most obvious feature of the continuous lithosphere is the width of the basin, which is ∼ 300 km. Amplitude of the flexure is by contrast strongly reduced (∼ 500 m). As can be seen from the "bending" curve, the second derivative of the flexural curve u(x), which is also closely related to both the radius of curvature of the plate and the local bending moment and 20 is plotted in the lower part of Fig. 7, the plate is continuously curved. Bending, moment and also elastic stresses in the plate are approximately 2 orders of magnitude higher than for a broken plate configuration. The reason for this is that a continuous plate is able to spread any applied load over large distances (up to its total length in fact). This reduces maximum amplitude increases flexural wavelength and leads to overall higher bending stresses in the plate. In a broken lithosphere, a load's compensation is mostly restricted to the length available between any two breaks. Hence, short segments of lithosphere, such as that we suggest 25 underlies the SCB are compensating loads with relatively little bending since both ends of the segment are relatively "free", in contrast to an unbroken plate where the ends are fixed in some way, for instance a fixed value of u(0) and u(L), where L is the far end of the plate, hence providing a stable point around which a bending moment is generated. Thus, segments in a broken plate behave like semi-isostatic levers, with maximised subsidence and steep and almost constant dips into the load. and resistance to tilting of the plate mostly provided by isostatic forces from mantle material moving in to the uplifted end of 30 the segment. We stress that the plate segments are still redistributing load but behaving almost like rigid bodies undergoing tilting in doing so. boundaries, but also some systems that fall under the label of intraplate shortening used here. The latter are actually present across many parts of the earth, including the Andean Cordillera (Kley et al., 1999), and many parts of Central Asia linked to the ongoing India-Eurasia convergence (Cobbold et al., 1996), as well as Central and Western Europe (Krzywiec, 2006;Kley 5 and Voigt, 2008). Our model emphasises how important the role of the elastic structure of the lithosphere is in controlling basin formation during intraplate shortening events. In particular, the frequency of basement thrusts will determine segment length, which, as shown, is one of the key factors influencing the geometry of subsidence and hence basin depth and width.
It thus seems likely that basin geometries can vary considerably as basement thrust systems themselves vary, even if they are genetically strongly related. However, the general character of these intraplate basins is likely to be similar to the usually larger 10 scale, plate boundary related foreland basins.
Hence, we suggest calling the SCB an intraplate foreland basin, since this conveys the compressional nature of the subsidence, the tectonic setting and the basin's likely sedimentary evolution, and also the fact that it will contain large amounts of the eroded product of adjacent basement uplifts and their sedimentary cover, if any is present.

15
The Subhercynian Cretaceous Basin is part of a wider system of intraplate, compressional basins across Western Central Europe whose subsidence and geometry can be well explained by a combination of elastic flexure and rigid tilting of lithosphere that has been tectonically segmented by basement thrusting. A relatively simple modification of the elastic flexure equations gives a reasonable, first order model of this process. These basins are likely to be relatively narrow, and deep relative to load size, although their dimensions will be strongly influenced by the frequency of individual basement thrusts, and hence the length of 20 plate segments.
Code and data availability. Both code and data for model runs are freely available from the author. 9 https://doi.org/10.5194/se-2020-185 Preprint. Discussion started: 27 November 2020 c Author(s) 2020. CC BY 4.0 License.

Appendix A: A broken plate formulation for intraplate loading
The general problem of variable thickness elastic lithosphere in elastic flexure equations has been discussed in several papers (Van Wees and Cloetingh, 1994;Manríquez et al., 2013;Garcia et al., 2014). It is not generally amenable to analytical solution, which has been the preference of most flexure studies. Our approach starts with a general formulation of the problem as follows u(x) is the deflection of the plate at any position x, along its length. D (the flexural rigidity) varies in space and is implicitly a function of x. The value of D is given by Eh 3 /(12 * (1 − ν 2 )), where E is the elastic modulus of the lithosphere, h is the effective elastic thickness of the lithosphere and is actually the parameter in D that varies in space, and ν is Poisson's ratio.
P is a constant representing a plate wide, horizontal stress. k = ρ mantle * g and represents a restoring force due to displaced 10 mantle, and q(u) is the load term, which consists of fixed, imposed loads generating subsidence q load * ρ load * g and also infill loads of the resulting basins q(u) inf ill * ρ inf ill * g. For instance, as shown in Fig. A1, a load due to an orogenic wedge on the edge of a plate would constitute the "fixed" load q o in our model, and is given directly at each node. The subsidence this load causes generates basins, which are then "filled" by q(u) inf ill , with material of density corresponding usually to water or sediments. This formulation of the load term and the restoring force results in a non-linear equation where the infill and 15 resulting elastic equilibrium has to be calculated iteratively. The model results in this paper show this process by plotting the initial calculated subsidence due to static loads as a black line, with the final, iteratively derived subsidence, including the infill plotted as a red curve. Static loads are shown in grey, whilst infill is shown with vertical red lines. The most common, analytical formulation of the flexure equation actually combines the restoring force due to displaced mantle with the infill as a load as part of the constant k, thus making the equation apparently linear. In fact, the equation should only be treated as linear under 20 certain special conditions, although this is rarely stated (cf. Gunn (1943a)).
We apply a finite difference operator for the second derivative, directly to the term in brackets, which is itself, discretised as a second derivative finite difference scheme, something referred to as the half station method (Cyrus and Fulton, 1968).
collecting terms, we obtain discretising the remaining parts of the equation then gives where the two load terms, (q o ) i and q(u i ) represent the static, fixed load and the iteratively calculated infill load respectively.
If we gather all coefficients into a matrix A and form a matrix equation, the resulting system is of the form a non-linear series of equations in u. We reformulate this as a recursive matrix fixed point problem of the form where r is the iteration step, and solve.
An alternative finite difference scheme can be derived by instead applying the product rule of differentiation to the problem before discretisation. For instance, Buiter (2000) shows Crucially however, this form of the equation assumes a smoothly varying elastic thickness function. Hence, for our problem, with abrupt variations in elastic thickness, it is not applicable. It can easily be verified that a finite difference scheme derived from this form of the equation will give an identical set of equations to the half-station method if D is constant everywhere.
Appendix B: An estimate of Harz shortening, and fault geometry from kinematic modelling foreland (Fig. 3). The southern end of the backlimb was thus tentatively placed at the transition from the very gentle dip of the Thuringian syncline to the slightly steeper one of the southern Harz.
In a second step we modelled the footwall syncline in the Subhercynian Basin using the Trishear algorithm (Erslev, 1991;Allmendinger, 1998). The fit shown in Fig. A2 was obtained with a trishear triangle having its apex initially 2 km below the top basement surface, an opening angle of 55 • , and upward movement slower than the fault tip (propagation-to-slip ratio 0.6).