Integrated geological and gravity modelling to improve 3-D model harmonization—Methods and beneﬁts for the Saxony-Anhalt/Brandenburg cross-border region (North German Basin)

of the derived model, we discuss the new geophysical ﬁndings on the sedimentary and crustal structures of the cross-border region in the context of the regional geological setting. The cross-border region is dominated by an NW–SE oriented fault system that coincides with the Elbe Fault System. We interpret a low-density zone within the basement of the Mid-German Crystalline Rise as a northward continuation of the Pretzsch– Prettin Crystalline Complex into the basement of the North German Basin. Additionally, we observe two types of anticlines within the basin, which we link to provinces of contrasting basement rigidity. Our gravity modelling implies that the Zechstein salt has mostly migrated into the deeper parts of the basin west of the Seyda Fault. Finally, we identify a pronounced syncline that accommodates a narrow and up to 800 m deep Cenozoic basin.

. (a) Exemplary illustration of misaligned stratigraphic horizons at the boundary of independently prepared geological models. The models' calculated gravimetric effect (+basement) is not consistent with measured gravity data. (b) Structural geological input models of Brandenburg and Saxony-Anhalt without Cenozoic cover (case study). Before harmonization, layer interfaces are vertically displaced by up to 500 m in the cross-border region. salt structures), underexplored regions, or regions with limited data access or low-quality data. For these scenarios, the incorporation of additional geophysical data such as gravimetric, magnetic, and/or electromagnetic data has been proved to reduce model uncertainty (e.g. Jacoby & Smilde 2009).
Typically, 3-D models are made for local studies and projects that deal with particular scientific questions or economic problems (e.g. exploration of raw materials or the storage of greenhouse gases). A great number of such models have been developed over the years as demand for geological insight has grown alongside computing power. Present-day, these small-scale models are often merged together to form larger ones, which may have country-wide coverage (e.g. Anikiev et al. 2019;Deckers et al. 2019). One of the biggest obstacles to this work is that neighbouring geological models often do not fit together at their common or overlapping boundaries, and so require alignment or harmonization at them (Fig. 1a). Two approaches are applicable to achieve harmonization: (1) Merging and matching by constraining and harmonizing raw data (e.g. seismic or borehole data) in the vicinity of the problematic border, and then re-modelling the cross-border region. This procedure involves extensive and time-consuming modelling work and is often hampered by overlapping and contrasting jurisdictions and data protection laws.
(2) If no modelling constraints are available in the cross-border region, the neighbouring models can be combined by simple point connections or interpolation. However, both possibilities are either very inefficient or may produce unrealistic pronounced dips in the transition area.
Another possibility is the harmonization of individual geological models in the framework of comprehensive gravity modelling. Wherever stratigraphic and structural interfaces involve density contrasts, the incorporation of gravity data is a reliable approach to unravel the 3-D distribution of geological bodies in the subsurface (Nabighian et al. 2005). What is more, because the gravity acceleration is expected to be inherently smooth in the cross-border region and will not feature sharp artefactual boundaries. Adequate gravity data are usually available and, if not, their acquisition is less expensive than obtaining additional seismic data or drilling. Furthermore, the merging and harmonization of subsurface layers into a single geological model in the framework of gravity modelling also allows for simultaneous cross-checking of the pre-existing structural input models, as they have to satisfy the observed gravity data.
The test case for our study uses two separately developed statewide structural geological models of parts of the North German Basin (NGB) in the states of Saxony-Anhalt (Malz et al. 2020) and Brandenburg (Schilling et al. 2018;Fig. 1b). Both models comprise interfaces of sedimentary layers overlying crystalline basement, starting from the base of the Late Permian Zechstein group up to the surface. Only a few line drawings of vintage seismic data (Kobylkin 1962;Wohlfeil & Teicher 1969;Mitjagin et al. 1981) are available in the two states' southern cross-border area, which lies in the southeastern part of the NGB in the area around the town of Wittenberg (Fig. 2a). The interpreters of those line drawings describe several general uncertainties in the trends and depths of the sedimentary layer boundaries, which cannot be reduced on the basis of the small number of deep boreholes in the area. Consequently, the general geometry of the stratigraphic layers in the southern cross-border region of Saxony-Anhalt and Brandenburg is poorly constrained, the two state-wide models are misaligned at the state border and there is a particular need for data-driven cross-border harmonization.
In this study, we will present an efficient approach to the integration of gravity data into the geological modelling process that allows to refine and harmonize existing, independently built 3-D structural models. The approach is readily applicable to other regions and with models at other scales. In the first step, we combine and homogenize several gravity surveys to build a consistent database for gravity field interpretation. Initially, the interpretation is based on wavelength filtering, gravity gradient calculation, and 3-D Euler deconvolution, and focuses on constraining the locations of local faults and estimating the depths of gravity anomaly sources in the cross-border area of uncertainty. These new geological findings act as a set of constraints for the subsequent merging process. For merging, we combine all available a priori information (i.e. seismic data, petrophysical data and drilling information) with the gravity field interpretation in order to set up an initial density model for the entire compilation region. Subsequently, we refine this model by 3-D gravity modelling (forward and inverse) that aims to match and so explain the observed smooth regional gravity variation. This process harmonizes the geologic input models, and allows for alternative structural scenarios to be tested. Afterwards, the final model's sensitivity to geometrical changes is evaluated to assess its significance.

Geological setting
The study area is located in Central Germany towards the southeastern end of the border of the federal states of Saxony-Anhalt and Brandenburg, close to the city of Wittenberg (Fig. 2, TK50-4142). Information on the structure of the crust and upper mantle base on interpretation of seismic refraction (Figs 2a and 3c) and seismological data, and gravity modelling in the KMGV (Lange 1973), SEMLJA (Pomeranzewa et al. 1971) and ZENTROSEIS research programs (Bormann et al. 1986). The main tectonic feature in the research area is the Wittenberg Fault (Fig. 3), a system of several major faults with varying dips. Some of these faults are assumed to displace the whole crust, down to the Moho at about 31-32 km depth (Bormann et al. 1986).
The thickness and depth range of the lower and middle crust (thought to have formed in Lower to Middle Proterozoic times) are still partly under debate. Seismic refraction profiles show a continuous reflection off the top of the lower crust at about 20-22 km (Bormann et al. 1986;Frischbutter & Lück 1997). The top of the middle crust is downthrown to the north along the Wittenberg  Martiklos et al. (2001). Naming convention of faults is adopted from Becker et al. (1989). (b) Map of the regional basement based on the work of Reinhold (2005). The outline of the Elbe Fault System is based on the work of Franke & Hoffmann (1999a,b) and Scheck et al. (2002). (c) Basement setting in the research area with available seismic refraction profiles (TS III profile-KMGV survey (Lange 1973); GSS-11 & GSS-6 profiles-SEMLJA survey (Pomeranzewa et al. 1971); THUBRA, FLELAU & GRIMBU II profiles-ZENTROSEIS survey (Bormann et al. 1986)). Naming and outlines of the tectonic blocks are taken from Bachmann et al. (2008)  Fault, deepening from about 4-6 km depth to as much as 11 km. In contrast, gravity modelling in the northwestern part of the study area indicates dense material of the lower crust lies at about 14-18 km depth and indicate the middle crust in the area to be thin (Bormann et al. 1986;Maystrenko & Scheck-Wenderoth 2013).
During the Ordovician to Early Carboniferous, the research area was part of the collision zone between Gondwana and Laurussia, whose action led to a complex division of the upper crust into zones of contrasting lithologies and densities. The Middle Devonian initiation of subduction between the microcontinents of East Avalonia and Armorica led to the rise of a volcanic arc (Henningsen & Katzung 2006). Today, the remnants of this arc are represented by dense magmatic rocks embedded in the Mid-German Crystalline Rise (MGCR; Oncken 1997), which comprises large parts of the study area's basement (Figs 3b and c). Ongoing plate convergence culminated in the Early Carboniferous Variscan orogeny, which resulted in the formation of the Northern Phyllite Zone and the uplift and subsequent erosion of the MGCR (Bachmann et al. 2008). At greater depths, high pressures and temperatures in the lower crust accompanied the ascent of syn-collisional, comparably light granitic to granodioritic magma (Henningsen & Katzung 2006), which were encountered by several boreholes in the Pretzsch-Prettin Crystalline Complex (PPCC; Fig. 3a; Kopp & Röllig 1996). A syn-collisional magmatic origin is also assumed for dense diorites and granodiorites that contribute to the prominent gravity high close to Wittenberg (Fig. 4b, GH-W;Conrad 1996;Röber et al. 1998). Finally, the collision of North Gondwana with Armorica in Early to Late Carboniferous times, led to the metamorphic overprinting of sediments and volcanic rocks in the Southern Phyllite Zone between the MGCR and the Saxothuringian structural unit (Figs 3b and c;Bachmann et al. 2008). The regional extents of these basement units and the zones of different rock composition/density are still a matter of debate (e.g. Böhnert 1988;Becker et al. 1989;Eidam et al. 1998).
In the Late Carboniferous, major changes in tectonic plate movements led to the development of the mainly SE-NW striking Elbe Fault System (Franke & Hoffmann 1999a;Scheck et al. 2002), which subdivided the basement into several small tectonic blocks and provided pathways for the extrusion of intermediate to acid magmatic rocks. On the basis of gravity anomaly fitting, Conrad (1995) postulated the presence of such post-collisional granitoidsmost likely light, acidic granites (Lehmann 1969)-as the source of the southeastern part of the Elster-Schweinitz gravity low (Fig. 4b, GL-ES).
The initial formation of the NGB has been attributed to additional heating, thermal destabilization of the crust, and increasing sedimentary loading north of the Wittenberg Fault (Scheck et al. 1999). Thermal subsidence of the basin continued until the end of the Early Permian with the deposition of 500-1000 m-thick aeolian and fluvial playa sediments (Rappsilber 2003;Bachmann et al. 2008) in the research area. In the Late Permian, the epicontinental transgression of the Zechstein sea led to the deposition of a thick pile of evaporites of contrasting densities (salt, gypsum and anhydrite), carbonates (limestone and dolomite), shales and a minor portion of sandstone (Ziegler 1989). Vintage seismic data and borehole information provide local evidence for the presence of up to 900 m of low-density Permian salt deposits of the Staßfurt subgroup in the study area (Kobylkin 1962;Wohlfeil & Teicher 1969;Mitjagin et al. 1981;Kalusa & Meyer 1983). A combined interpretation of seismic and potential field data provides some evidence for a possible thickened formation of dense Permian anhydrite (>300m) in the basin's marginal area around Jessen (Fig. 3a;Lehmann 1969), but the data resolution is too low to reveal its entire extent.
Uniform basin subsidence and phases of E-W-(Triassic) and NE-SW-directed (Jurassic to Early Cretaceous) extension (Ziegler 1990) led to the deposition of up to 3 km of German Triassic (i.e. Buntsandstein, Muschelkalk, Keuper) to Cretaceous sediments in the study area (Wohlfeil & Teicher 1969;Kalusa & Meyer 1983). These sediments, of various compositions and densities, were deposited in phases of continental, epicontinental and marine conditions as sea-level changes fluctuated. In the early Late Cretaceous, the study area experienced N-S-to NE-SW-directed far-field compressive stress related to plate convergence between Africa and Europe (Kley & Voigt 2008). This led to the uplift and erosion of basement blocks as existing normal faults were reactivated in a reverse sense (Fig. 3a).
During Cenozoic times, the cross-border region between Saxony-Anhalt and Brandenburg experienced alternating marine transgressions and regressions (Bachmann et al. 2008) and the deposition of marine and continental sediments of about 100 m thickness south of the Wittenberg Fault and up to 600 m thickness in the basin north of it (Wohlfeil & Teicher 1969). Locally, lignite of considerable thickness is intercalated in the Cenozoic cover.

Previous gravity modelling
Previous gravity studies in the research area focused on different targets. The early gravimetric investigations in the region mainly studied the normal anomaly field and concentrated on the general tectonic setting (Closs 1939;Siemens 1948;Küttner & Bein 1963;. Detailed models of the basement's lithologies and structures were obtained by 2-D density modelling along coincident seismic refraction profiles (Pomeranzewa et al. 1971;Bormann et al. 1986) as well as by filtering, gradient calculations and modelling along a single profile in the scope of combined potential field studies (Becker et al. 1989;Conrad 1995). Some of these results were more recently incorporated into large-scale 3-D gravity and geothermal models of the Central European Basin System, with the aim to understand the regional crustal structure and geothermal conditions (e.g. Scheck et al. 1999;Kuder 2002;Yegorova et al. 2007;Kaban et al. 2010;Maystrenko & Scheck-Wenderoth 2013;Anikiev et al. 2019). Other 2-D gravimetric studies have focused on the sedimentary column, for lignite and ore prospecting (Eichner 1981;Böhnert 1988) or for tracing a hydrocarbon trap associated with a regional anhydrite platform (Lehmann 1969). Except for the study of shallow lignite exploration, the interpretation of gravity gradients as well as filtered anomaly maps revealed limited stratigraphic information, and led Lehmann (1969) to suggest an integrated approach and modelling in future work.

DATA
Our merging and harmonization strategy of 3-D geological models is applied to an area at the southeastern part of the NGB in the southern cross-border region between the German states Brandenburg and Saxony-Anhalt (Wittenberg, TK50-4142, Fig. 2). The main study area is extended with a buffer of 12.5 km to allow for more confident detection of large-scale tectonic features and faults that continue outside it, as well as to fully account for gravity signals from close-by geological features during the subsequent modelling process. In total, the gravity data set covers an area of about 50×50 km. The standard coordinate system we use is the ETRS89/UTM zone 32N (EPSG code: 25832). All gravity values are stated in the frequently used unit mGal, where 1 mGal = 10 -5 m s −2 (SI unit). All density values are stated in the widely used cgs unit g cm −3 , where 1 g cm −3 = 1000 kg m −3 (SI unit).

Structural models
The geological 3-D SKUA-GOCAD models of the NGB in Saxony-Anhalt (Malz et al. 2020) and Brandenburg (Schilling et al. 2018) serve as structural input models (Fig. 1b). These models were separately developed during the TUNB (Deep underground of the NGB) project, which is a collaboration of several states of Northern Germany to prepare a large-scale 3-D subsurface model. The models include 13 main horizons, from the surface down to the base of the Permian Zechstein formation, not all of which are present in our study area. Throughout most of the research area, these models are constrained by borehole data, regional depth maps of major seismic horizons (Reinhardt & Gruppe Regionales Kartenwerk 1968-1991 and an extensive set of interpreted depth-converted seismic reflection sections for structural interpretations, fault sticks, narrow folds and outlines of diapirs. However, the southern cross-border region between Saxony-Anhalt and Brandenburg is poorly sampled, leading to large uncertainties in the modelled surfaces and misalignments between the models at the border.

Gravity data
The incorporated gravity data were acquired between the years 1934 and 1969 over the course of five campaigns (Table 1). In total, 5464 readings ( Fig. 4a) are available at an average point spacing of 680 m (1.5 gravity stations per square kilometre). This is better than the resolution required to discern gravity anomalies with wavelengths suitable for regional studies (e.g. Hinze et al. 2013). The data were digitized during the last decades and are today available from different institutions. Originally, the gravity surveys were tied to different gravimetric nets and different processing techniques were applied to the recorded data. A common datum and workflow for reprocessing are described in section 4.1.

Seismic data
Vintage seismic reflection data in the study area provide local constraints on the setting of the Cenozoic sediments down to the base of the Permian Zechstein formation (Kobylkin 1962;Wohlfeil & Teicher 1969;Mitjagin et al. 1981;Kalusa & Meyer 1983;Fig. 2a). Mostly, these constraints come in the form of line drawings together with general uncertainties concerning the depths and trends of the interpreted interfaces. These data are uninterpreted across areas of geological complexity, such as halokinetic or intensely faulted sequences.
Information on the structure of the crust and upper mantle down to about 40 km depth comes from combined studies of seismic refraction (Fig. 3c) and seismological data of the KMGV (Lange 1973), SEMLJA (Pomeranzewa et al. 1971) and ZENTROSEIS research programs (Bormann et al. 1986). The main findings are a lateral varying composition in the upper crust and a change of With this in mind, we incorporated three main crustal layers as well as the upper mantle into our modelling process, in order to account for long wavelength gravity anomalies. The depths and trends of the crustal interfaces are based on main reflections in the deep seismic profiles. Neither the original interpretations, nor newer reinterpretations, of the seismic refraction data in the study area were able to discern distinct velocity gradients in the individual crustal layers (Schulze & Lück 1992;Lück 1996;Frischbutter & Lück 1997).

Borehole data and petrographic information
Stratigraphic data were incorporated from 89 wells (Fig. 2a, triangles). These data are available from the borehole database of the Geological Survey of Saxony-Anhalt. Mesozoic and Palaeozoic sediments in the basin area were encountered in 27 boreholes to depths of 500-2350 m, whereas boreholes outside the NGB reach to depths of about 100-500 m. Downhole geophysical logs (especially formation density) are available for 13 of the boreholes (Fig. 2a, yellow triangles). Additional data on the petrographic compositions of crystalline rocks in the study area are provided by thin sections from 9 boreholes (Fig. 2a, blue triangles;Erzberger et al. 1967a,b;Röber et al. 1996Röber et al. , 1997Eidam et al. 1998) and density measurements of rock samples from 5 boreholes (Fig. 2a, light blue triangles; Carl 1988). Kopf (1967) prepared a first compilation of density gradients in the Mesozoic and Palaeozoic sediments in the northeastern part of the NGB, which is based on different geophysical and petrophysical approaches. Subsequently, Köhler & Eichner (1973) extended the study of density gradients to the entire eastern part of the NGB by determining saturation densities of Mesozoic rock samples from 94 boreholes in laboratory experiments. In the framework of our study, the resulting density-depth relation charts of Köhler & Eichner (1973) were digitized and approximated with polynomial equations (Fig. 5a). The obtained densities and gradients were cross-checked against a collection of local average densities from Becker et al. (1989) and are noted in Appendix A.

Densities
The density characteristics of the Permian Zechstein formation are complex. In the basin area of the study region, the Zechstein formation is dominated by thick salt deposits of the Staßfurt subgroup, with densities typically constant at between 2.18-2.22 g cm −3 . Substantially different densities are identified for the younger and older Zechstein cycles, however, with their locally distinct layers of limestone, dolomite, gypsum and anhydrite along the basin margin. To account for this, we inserted additional layers at the top and bottom of the Zechstein layer. Based on borehole information, we adopted an average density of 2.85 g cm −3 for these parts of the Zechstein formation that differ lithologically from the less dense Staßfurt halite.
Eight boreholes with gamma-gamma borehole measurements (formation density) provided additional density information in cpm (counts per minute) for the Mesozoic and Palaeozoic sediments in the basin area ( Fig. 2a, yellow triangles). Sadly, formulas necessary for calibration of these data were not documented in the logging reports, and so only a rough conversion to g cm −3 units was possible. Accordingly, these logs are only qualitatively used to evaluate density gradients and prominent density contrasts in the subsurface at the borehole positions ( Fig. 5b).
Basement densities in the southern part of the model are assigned based on laboratory measurements at rock samples of five boreholes in the transition area between the PPCC and its surroundings (Fig. 2a, light blue triangles;Carl 1988). The small database of 27 measurements reveals slightly lower densities in the PPCC (Fig. 6). The only information on basement density beneath the basin comes from a single gamma-gamma (formation density) log at the basin margin.
The only petrophysical data for the deeper (Palaeozoic and Proterozoic) parts of the crust come in the form of P-wave velocities from seismic refraction studies (without velocity gradients). The velocities were converted to densities according to the relationship of Brocher (2005). No density-depth or velocity-depth relationships are available for these layers and, consequently, we assigned them constant densities.

M E T H O D S
Once the basic parameters of the study are fixed (data boundaries, projection system, relevant horizons, etc.) and the input data are gathered, our modelling workflow (Fig. 7) consists of the following steps: (i) compilation and homogenization of a consistent gravity  (Köhler & Eichner 1973). The relationship was adopted for the density model voxels. (b) Comparison between measured (gamma-gamma-log, black line) and modelled (red line) densities in the voxels coincident with the borehole. database, (ii) fault detection and gravity interpretation throughout the entire model domain with a special focus on the cross-border area, (iii) forward modelling, including model harmonization and (iv) inverse modelling. In the following sections, we give detailed descriptions of the preparation and incorporation of the input data as well as of the modelling steps. The described approach focuses on the extraction of enough new geological information in the sediment column and upper crust to allow merging and harmonization of the structural input models (Schilling et al. 2018;Malz et al. 2020) in the cross-border region. Application of these methods to other regions or for models of different scales might require choices of different wavelengths or gradients for gravity interpretation and incorporation of additional modelling constraints.

Compilation of gravity database
The first step is to compile a consistent and uniformly processed gravity database (Fig. 4a). For this purpose, we selected the approach of Skiba (2011), which was used in processing of a standardized gravity map of Germany (Skiba & Gabriel 2010). In this approach, the gravity data are tied to the International Gravity Standardization Network 1971 (IGSN71) and heights refer to the German First Order Levelling Network 1985 (DHHN85). In our study, calculation of the free-air anomaly encompasses latitude correction with the formula of Somilgiana (Moritz 2000), free-air correction according to Heiskanen & Moritz (1967), and an atmospheric correction based on Wenzel (1985). Derivation of the complete Bouguer anomaly (CBA) required the additional calculation of a spherical Bouguer correction with first term according to Lafehr (1991) and second term based on Cogbill (1979). DEM5 data (GeoBasis-DE/LVermGeo LSA 2013) for Saxony-Anhalt and DEM25 data (GeoBasis-DE/LGB 2018) for Brandenburg were available for the near-field topographic correction and SRTM data (Earth Resources Observation And Science (EROS) Center 2017) for the far-field topographic correction within a radius of 166.7 km.
Ten boreholes with gamma-gamma (formation density) and calliper logs in the Cenozoic cover were available for the  determination of a suitable reduction density (Fig. 2a, yellow triangles). The combined investigation of the logs with the petrographic details of the boreholes reveal average densities of 1.8-1.9 g cm −3 for the sediments from the surface down to the reduction level at sea level. This comparatively low reduction density is attributable to the widespread presence of porous lignite in the Cenozoic cover in the study area. This reduction density appears to be relatively constant throughout the research area, and so we applied a uniform value of 1.85 g cm −3 .
The consistency of the data was assessed by leave-one-out crossvalidation, which removes one data value at a time and predicts the associated value at the removed station by means of geostatistical analysis of the surrounding values. Stations with differences of greater than 0.5 mGal between the digitized and predicted values were evaluated and, where necessary, rejected. Finally, measurement points with present-day elevations differing by more than 5 m from their protocolled heights and/or with inadequate positioning were also discarded. Remaining minor differences between the processed surveys were levelled by a linear approximation.

Gravity field transformations for first-order structural analysis
A detailed study of the gravity field delivers interpretations that we use as a priori inputs to the subsequent modelling process. The wavelengths of gravity anomalies vary according to the depths, sizes and density contrasts of their sources (e.g. Militzer & Weber 1984). Several filtering and processing techniques allow for decomposition of the gravity field to emphasize the signals of interest (Table 2; e.g. Blakely 1995). Standard interpretation of the gravity anomaly field comprises filtering to separate the long wavelength signals of regional or deep-seated geological features from short wavelength signals of local or shallow structures (Fig. 8). At this stage, the depths and edges of causative structures can be only roughly estimated by rules of thumb, because of the superposition of gravity signals from unrelated nearby sources. More precise estimates can be obtained by the calculation of derivatives of the gravity field and the application of so-called depth-to-source techniques. Even then, however, interpretation of the gravity field remains subject to the nonuniqueness problem, which means that a continuum of solutions Estimation of source depths x x x Figure 8. Separation of the regional (a) and local (b) gravity anomaly field by wavelength filtering of the CBA (Fig. 4b). Long wavelength anomalies are caused by varying rock compositions in the basement and regional differences in sediment thickness and sedimentary facies (Lehmann 1969;Conrad 1996). exists to explain an observed gravity signal (Skeels 1947;Saltus & Blakely 2011). Furthermore, superposition of anomalies leaves open the possibility of destructive interference, which may leave some sources uninterpretable. Consequently, it is crucial to consider all available constraints prior to interpretation to limit the number of possible solutions in gravity modelling and interpretation.

Application of gravity gradients
Gravity field gradient products are suitable for highlighting the edges of buried density contrasts, and thus for increasing the resolution of source boundaries (e.g. Militzer & Weber 1984;Blakely 1995). In the framework of this study, we incorporate the horizontal gradient (a combination of the Gxz and Gyz component of the Eötvös tensor, HDR; Cordell 1979) to highlight density contrasts at local fault systems (Fig. 9a). The HDR maximum of gravity over a nearly vertical and isolated fault plane is situated above the plane's top edge. In contrast, at dipping faults the maximum is offset from the top of the fault plane. Its true location can be estimated with the approach by Grauch & Cordell (1987). Where several faults contribute to an HDR maximum, or large terrain effects occur, the true boundary location can only be assumed (Grauch & Cordell 1987).
To gain rough insights into fault extents at depth, we bandpass filtered the CBA by to isolate wavelengths like those caused by a synthetic vertical fault in a certain depth range, for example faults at 5-10 km depth can be expected to raise anomalies with wavelengths of about 25-50 km. Subsequently, we calculated and plotted the HDR for each filtered grid, ready for fault interpretations (Fig. 9b). Preservation of an HDR maximum through multiple plots can be interpreted in terms of a fault's depth extent. Here, the problem of non-uniqueness remains, for example if a shallow, wide, sediment body with undulating surfaces raises a long wavelength anomaly, whose edges might be misinterpreted in the HDR in terms of a deep fault. Furthermore, only faults that dip by more than 45 • and with a distinct density contrast at the fault plane will produce a significant change in the HDR.
The vertical gravity gradient (Gzz component of the Eötvös tensor) represents the rate of change of the vertical component of gravity with height (Evjen 1936). Its signal is positive over a source body, is negative outside it, and crosses zero over vertical density contrasts (Fig. 10a). In comparison to the standard gravity field, the vertical gravity gradient is characterized by narrower anomaly shapes and emphasizes shallower sources.
In contrast, the gradient tilt angle (Miller & Singh 1994) responds equally well to shallow and deep sources (Fig. 10b). For this reason, it is used to reveal subtle deeper sources, whose responses are usually obscured by those of shallower sources. Additionally, the maximum value of the real part of the hyperbolic tangent function of the tilt angle (Cooper & Cowan 2006) is calculated to emphasize the boundaries of sources detected by the tilt angle (Fig. 10c).  (Kobylkin 1962;Wohlfeil & Teicher 1969) are also plotted. The inset in the upper right corner is a rose diagram showing the strike directions of the main faults in the study area. (b) HDRs of the bandpass filtered CBA products, which allow for rough tracing of the main faults with depth (see the text for further explanation).

Estimation of source depths
A quantitative interpretation of the gravity data focuses on the extraction of source body characteristics, for example depth, size or shape. First evidence for the approximate depth of an anomaly's source can be estimated by the application of so-called depth-tosource techniques.
In the framework of this study, a 3-D Euler deconvolution (Thompson 1982) was applied to the CBA using the code of Pašteka et al. (2009;Fig. 7a), to obtain new insights into the main geological features in the upper 5 km (the depth range of the input models).
The deconvolution requires the prior calculation of gravity field gradients and the selection of a suitable structural index (SI), a parameter that indicates the degree of gravity field homogeneity (Reid et al. 1990). Originally, the SI is defined as an integer and varies for different source types . To locate the depths to the tops or centres of the main anomaly sources, such as evaporite-bearing layers of the Zechstein formation, basement and deep-seated intrusions, we selected SIs of 1 (appropriate for horizontal lines/cylinders, faulted thin beds, or vertical line ends) and 2 (appropriate for spherical sources; Reid & Thurston 2014). In contrast, a SI of 0 is applied for sources that are expected to resemble the edges of thin sheets and faults.
We use a moving window with a size of seven grid points (point spacing of 1000 m) during the deconvolution. A set of linear equations is solved at each window position to locate the source . Selection of the obtained solutions is based on evaluation of their standard deviations and condition numbers (Pašteka et al. 2009;Fig. 11). Only clusters of depth solutions are considered during the subsequent gravity modelling.

Model set up and incorporation of constraints
3-D density modelling was performed with the interactive gravity and magnetic modelling software IGMAS+ (Götze & Lahmeyer 1988;Schmidt et al. 2020). Set-up of the initial density model is based on the geological 3-D SKUA-GOCAD models of the NGB in Saxony-Anhalt (Malz et al. 2020) and Brandenburg (Schilling et al. 2018;Fig. 1b). For this, the horizons of the input models are transferred to 26 2-D serial model sections in IGMAS+ (Fig. 12b), all oriented perpendicular to the NW-SE strike direction of the main tectonic structures. We chose a 2 km spacing between the sections in the main study area, and 4 km in the extended study area surrounding it. All model sections were laterally extended by 250 km to minimize edge effects. The selected reference density (outside the extended model sections) is 2.67 g cm −3 . To account for long wavelength anomalies of deep-seated sources (not included in the 3-D SKUA-GOCAD models), we introduced additional layers for the Permo-Carboniferous sediments, three crustal layers (upper, middle, lower), and the upper mantle (Fig. 12a) based on borehole information and interpretations of the nearby seismic refraction profiles. In total, the density model comprises 17 layers from the surface down to 40 km depth.
We use free-air gravity anomaly data and its horizontal gradient for density modelling. The surface of the model's Quaternary layer (topmost layer) is taken at the topography of the DEM5 (GeoBasis-DE/LVermGeo LSA 2013) for Saxony-Anhalt and the DEM25 (GeoBasis-DE/LGB 2018) for Brandenburg. Structures interpreted from the gravity anomaly field are imported as point information, to allow the correct positioning and integration of faults and main anomaly sources in the cross-border region. Line drawings of depthconverted 2-D reflection seismic profiles (Kobylkin 1962;Wohlfeil & Teicher 1969;Mitjagin et al. 1981) are included to permit consideration of their depth information. Additionally, constraints were included from boreholes that penetrate Mesozoic or older strata.
The density model was gridded in 125 m steps in the horizontal directions and in 25 m depth steps, in order to implement vertical density gradients in the stratigraphic horizons (Fig. 13). In total, about 50 million voxel cubes were generated, the maximum number possible with IGMAS+ . Subsequently, based on the polynomial approximations to the density-depth relation charts of Köhler & Eichner (1973; see section 3.5), density gradients are assigned to most layers (Tables 3 and A1). As they are assigned individually to each layer, density gradients can change at the layer interfaces and at fault planes.

Modelling
In a first step, we applied a forward modelling approach (Fig. 7b, orange box), whereby iterative changes were made to the model geometry or density parameters with the aim to achieve a close fit of the calculated gravitational response to the observed free-air gravity anomaly. Forward modelling is done with a bottom-to-top approach, starting with the long wavelength anomalies. Due to its strong variability in lithology and density, the upper crust was subdivided into several parts. The layer boundaries in the cross-border region were initially harmonized by observing the effect of continuing each of the bounding input models across the state boundary (Fig. 12c). Here, the states' individual models and a model with a different Zechstein facies are tested in the cross-border region. In this example, a northeastward extrapolation of the model for state 2 (across the border) results in the best fit to measured gravity. The structural setting and comparison to gravity data prior harmonization and gravity modelling (simple merging) is shown in Fig. 1(a).
Where neither of the 3-D input models produced a plausible match to the observed gravity, either an interpolation of the two or an independent solution was used to produce an acceptable fit. In addition to the free-air gravity anomaly, gradient modelling (HDR) was used to model the positions and dip angles of major faults. Since the Mesozoic layers are locally well-constrained by the existing structural models, the layer properties were first adjusted by varying densities within the range of their standard deviations (Köhler & Eichner 1973). If this did not produce an acceptable fit, then the layer's geometry was iteratively changed. Various geological scenarios were tested during forward modelling, to keep full control on the density model's geometry and geological plausibility. The residual anomaly field was inspected to evaluate the consistency of the tested geological scenario to the measured gravity data (Fig. 14).
A physical parameter inversion (Saether 1997) for layer densities was used to further improve the model's fit. For this, the error of the gravity anomaly was selected to be about ± 0.05 mGal, taking into consideration the standard deviations of the individual measurements in the gravity surveys (Table 1) and potential inaccuracies in the processing routines. The standard deviations of densities in the Mesozoic layers with their incorporated density gradients were assigned based on Köhler & Eichner (1973). The standard deviations adopted for all other layers are listed in Table 3. An automatic geometry inversion was applied to subvolumes in areas of higher misfit between neighbouring model sections to improve the final density model (Fig. 12b). In a last step, the harmonized and depth-adjusted layers of the IGMAS + model were incorporated in a volumetric 3-D SKUA-GOCAD model, which was parametrized with the density values.
We tested the gravity anomaly's sensitivity to changes in individual layer geometries to assess the significance of the respective layer for the gravity modelling process. This test aims to highlight those layer boundaries whose depth changes result in obvious variations of the calculated gravity anomaly and are consequently most significant to the success of the modelling process. By identifying less influential boundaries, the test also highlights layers that may be subject to larger geometrical uncertainty. For each layer, the test involves perturbing its base vertically at the midpoint of a 4 kmwide window by enough to change the modelled gravity anomaly by 0.1 mGal. The perturbations are listed in Table 3. To account for the decreasing gravity signal of a mass with increasing distance from the gravity stations (surface), the uppermost and lowermost positions of the base of the respective layer were altered. The sensitivity also depends on the density contrast at the interface.

R E S U LT S A N D I N T E R P R E TAT I O N
Our gravity modelling workflow (Fig. 7) allows us to successfully merge the structural input models and set up a harmonized geological model that is consistent with both the observed gravity anomaly and all incorporated geological constraints. The residual anomaly of the input models prior to our process (Fig. 15c) shows long wavelength misfits with amplitudes of up to 10 mGal (average 3.5 mGal), indicating areas that require structural changes. The subsequent gravity modelling and adjustments to the geological model allow us to reduce the maximum residual anomaly to 3 mGal, and its average to 0.7 mGal (Fig. 15e). The residual anomalies are small enough to be attributable to local shallow sources, which are not object of this study.
Gravity-based harmonization of 3-D geological models unavoidably requires changes to the input model geometries. To assess the significance of the depth information and shape of the modelled structures obtained by gravity modelling, the calculated gravity anomaly's sensitivity to geometrical changes of its individual layers (section 4.3.2). The tests showed the greatest sensitivities to the near-surface Cenozoic layers, and to the Zechstein and Muschelkalk layers, with their distinct density contrasts to adjacent layers. As a consequence, interfaces of the Jurassic and inner boundaries of the Buntsandstein are only modelled to follow the general trend of the surrounding layers. The contribution of the top basement surface with its density contrast of > 0.03 g cm −3 to the overlying Permo-Carboniferous, is significant in areas, which have a distinct density contrast (> 0.03 g cm −3 ) to the overlying Permo-Carboniferous. Modelling of the intracrustal interfaces revealed that density and/or topographic variations deeper than about 21 km give rise to very long wavelength changes that are indistinguishable from static shifts to the modelled anomaly within the study area (Table 3, last column). Consequently, this study delivers no new constraints on the structure of the deeper crust.
Effective merging and harmonization of 3-D geological models in underexplored or unconstrained areas by means of gravity modelling involves revealing new information on the geological setting in the study area. These must be described and discussed to demonstrate their plausibility and contribution to the modelling process. In our case study of the NGB, new information is obtained for the sedimentary column and the upper crust. In the following, we will highlight those that are relevant for the merging process in the cross-border region, as well as describing their implications for the general geological setting by comparison to previous largescale density models. A more detailed description of the geological results is provided in Appendix B.
The investigation of gravity gradients and Euler depth solutions (Figs 9-11) provided a detailed depiction of the fault system, tracing of geological structures and the assignment of discrete upper crustal Table 3. Layer parameters and their respective uncertainties. The incorporated density gradients are noted in Table A1. The last column states values of the sensitivity analysis. These values indicate perturbations to the base of the respective layer necessary to result in a variation of 0.1 mGal in the model gravity anomaly. The first number indicates the necessary perturbation at the highest position of the layer, and the second number at its lowest position. Layer  Figure 14. Left side: tests of different structural scenarios for an anticline associated with the Elster Fault (1-basement and sediment cover similarly folded; 2-detachment fold with mildly faulted basement; 3-fault steps in basement and cover). Right side: the resulting modelled free-air gravity anomaly and residual anomalies, after subtraction of the observed field. The average residual anomalies suggest that scenarios 1 and 3 are plausible, whereas scenario 2 is not. blocks in the cross-border region, which served as constraints in the modelling and merging process. We observe indications for mainly NW-SE and NE-SW striking faults in the gravity gradient maps. These coincide with interpretations of faults in vintage seismic data (Wohlfeil & Teicher 1969;Bormann et al. 1986; Fig. 9a) and allow us to map the continuity and orientation of faults throughout the study area. The highest gradients are observed along the Wittenberg Fault, where gradient modelling of the HDR (Figs 16  and 18) indicates a steepening of the fault from northwest (45 • ) to southeast (70 • ), and the displacement of upper and middle crustal blocks down to at least 20 km depth (Fig. 9b). Minor faults (e.g. Elster Fault, Zahna Fault) are interpreted to cut the Mesozoic and Palaeozoic sediments and the upper crust. Interpretation of the tilt angle (Fig. 10b) and gravity modelling reveal pronounced basement topography and an inhomogeneous density structure in the upper crust. This is in agreement with local studies (Bormann et al. 1986;Becker et al. 1989), but was not resolved in previous regional density models (e.g. Scheck et al. 1999;Kuder 2002;Yegorova et al. 2007;Kaban et al. 2010;Maystrenko & Scheck-Wenderoth 2013, Anikiev et al. 2019. Modelling of the long wavelength gravity anomalies provides new indications of the boundaries of the low-density PPCC and its extension into the marginal area of the NGB (Figs 16 and 17). Densities in the surrounding upper crust are about 0.07-0.14 g cm −3 higher than in the PPCC (Carl 1988;Conrad 1995;Figs 7, 16 and 17).
Generally, the sedimentary cover is strongly folded (Fig. 18). Tests of different structural settings of anticlines in the study area (e.g. along the Elster Fault; Fig. 14) suggest anticlines in the area of the Northern Phyllite Zone sole out on a flat Zechstein base. In contrast, a folded basement and thick uniform Zechstein layer are observed in association with the anticline along the Elster Fault within the MGCR. Combined with new evidence that the    (Figs 18 and 19), our findings on the topography of the Zechstein formation are of great importance for the merging and harmonization process of the structural input models because the Mesozoic sediment layers follow the topography of the Zechstein. East of the Seyda Fault, gravity modelling and borehole data indicate the Staßfurt Salt to be mostly substituted by dense anhydrite and dolomite (up to 350 m thickness; Figs 18 and 19). Finally, gravity modelling revealed the extent of a pronounced Cenozoic basin, locally drilled in the past, of up to 550 m thickness in the cross-border region (Figs 18 and C1). The smooth synclinal structure of this basin causes, together with the Cenozoic basin and the low upper crustal density of the PPCC the Gravity Low of Elster-Schweinitz. In general, these stratigraphic and lithological details were neither resolved by in previous regional density models (e.g. Maystrenko & Scheck-Wenderoth 2013) nor predicted by the structural input models (Schilling et al. 2018;Malz et al. 2020).

D I S C U S S I O N
The geological significance of new findings obtained by gravity modelling always require a discussion in the context of the regional geological setting and with regard to model accuracy, which is provided in Appendix C. The discussion here is more appropriate to our study's focus on the modelling strategy for merging and harmonizing 3-D geological input models in the framework of 3-D gravity modelling. It concentrates on the assessment of uncertainties in the input data and modelling approach.

Uncertainty estimation of input data and constraints
Generally, uncertainties in 3-D density models originate from the input data, processing and modelling procedures, and the problem of non-uniqueness in potential field data interpretation. Malz et al. (2020) discuss sources of uncertainty in the input geological 3-D models used in this study. These include objective uncertainties such as depth uncertainties of several metres related to the quality and limitations of underlying seismic data and depth maps, or smaller depth uncertainties in the incorporated borehole data and model resolution controlling the model's level of detail. Additionally, the authors highlight more subjective uncertainties related to the interpreter's experience of structural interpretation, borehole-seismic correlation, fault interpretation, and other underlying geological concepts (e.g. Bond 2015).
The additional input data for gravity modelling similarly subject to uncertainties. The reported standard deviation of the gravity data readings (Table 1) and potential inaccuracies of the processing routines indicate the input gravity data may be reliable at the level of about ± 0.05 mGal. This can be considered to be of minor importance for the interpretation and modelling procedure.
Depth estimates from 3-D Euler deconvolution require assumptions to be made about the gravity anomaly's source geometry, but also that a suitable window length is chosen for the sampling density and grid spacing of the input data, and the geological complexity they portray . This all depends strongly on the interpreter's knowledge of the local geology and level of experience in application of the depth-to-source technique. Different SI  values will alter the resulting solutions in depth and, to an extent, produce clusters of solutions in different areas, complicating the interpretation of main anomaly sources (Fig. 11). The fault interpretation is mainly based on the study of gravity gradients. Interference between the fault's gravity signal and those of adjacent faults or other sources causes a blurring of the anomaly of interest and hinders precise tracing and dip estimation. For an isolated non-vertical fault, the dip-related shift of the HDR signal can be corrected by using the scheme of Grauch & Cordell (1987). Significant HDR signals can only be expected from faults that dip at ≥ 45 • , and whose fault planes presenting a distinct density contrast. Consequently, low angle detachment faults or nearly horizontal interfaces will not be observed in the HDR (or tilt angle). The incorporation of the HDR to the modelling of main faults in the research area allowed a dip estimate of about 15-20 • (for faults that dip at ≥ 45 • ), hampered by the complex sediment setting at these interfaces.
The densities and density gradients in our modelling are mostly based on compilations of laboratory measurements, direct readings of formation density (gamma-gamma) logs, converted densities from acoustic log readings, and estimated densities based on petrographic descriptions. Laboratory studies partly investigated different types of densities (e.g. dry bulk density, grain density) with a precision of about ± 0.03 g cm −3 (Schubert 1977). These densities are subsequently converted to saturation densities, since these are considered to represent in situ conditions (Kopf 1967). Density gradients in any specific stratigraphic horizon represent an average gradient, whose deviation from the entire database is estimated by the standard deviation [about ± 0.05 g cm −3 for the compilation of Köhler & Eichner (1973) ; Fig. 5a]. Densities converted from acoustic logs and gamma-gamma logs in Mesozoic and Palaeozoic sediments in the research area are assumed to be precise to within about ± 0.05 g cm −3 (Schubert 1977). The reliability of densities converted from P-wave velocities, as for the lowermost basement layers (using seismic refraction data), strongly depend on the precision of the velocity determination and could be only estimated by a manual perturbation test of the model to be in the region of ± 0.1 g cm −3 . Generally, the model densities can only represent an approximation to the real subsurface density distribution. To prevent an unnecessarily complex model density distribution for the gravity modelling, the model layers are mostly assigned constant density gradients, or even uniform average densities (Fig. 5b). Local deviations from these conditions are certain to exist in the research area, for example due to differences in void filling, diagenesis, facies variability, and differential compaction. These will be not resolved by our approach. Where it is necessary to account for regional inhomogeneity of a stratigraphic horizon, a layer subdivided to allow for more realistic gravity modelling.

Uncertainty introduced by the modelling approach
In general, gravity forward modelling combined with the nonuniqueness problem in interpreting potential field data introduce an unquantifiable subjective uncertainty. In the scope of this study, we have tried to account for some of this uncertainty by testing multiple geological scenarios that appear, based on our experience and understanding of the literature, to be plausible for the region's geology and to agree with the incorporated a priori constraints. However, by this approach it cannot be excluded that a model based on another plausible geological concept may explain the observed gravity data with a similar accuracy.
In general, we have focused on the region's main gravity anomalies, given that our aims are to produce model harmonization and to generate new insights into the structural setting in the research area. A geometric inversion was applied to improve the model in areas of significant local misfit between the model sections. Beyond this, however, detailed modelling of all of the region's small-scale anomalies and local geometry changes is beyond the scope of this study.
Finally, gravity modelling alone cannot determine absolute density values. Gravitational effects are caused by density contrasts at the layer boundaries. Absolute density values can only be determined on the basis of laboratory measurements or geophysical borehole logs. Consequently, densities obtained by gravity modelling can be only as precise as the assigned absolute densities, which are based on other methods.
The accuracy of depth information extracted from the density model is hard to estimate numerically, due to the non-uniqueness problem and the resulting large number of plausible solutions. However, to assess at least the significance of a layer to the gravity modelling process, the sensitivity of the calculated gravity anomaly to depth changes of individual layers can be tested. We found that layer interfaces with distinct density contrasts and near-surface layers are highly relevant for fitting the gravity signal and so contribute most strongly to the modelled geological setting. Accordingly, in our study, these are (i) the Cenozoic layers, (ii) the Muschelkalk and Zechstein and (iii) much of the top of crystalline basement (Table 3). Layer interfaces with a smaller density contrast, like the Jurassic and inner boundaries of the Buntsandstein are only modelled to follow the general trend of the adjacent layers.

C O N C L U S I O N S
Our study presents a strategy for the harmonization of 3-D geological models in the framework of 3-D gravity modelling, which is exemplarily applied to a case study at the southern boundary of the NGB in the cross-border region between Saxony-Anhalt and Brandenburg. Based on a consistently compiled database of gravity, seismic, borehole and petrophysical data, we find that gravity interpretations (filtering, gradient calculation and 3-D Euler deconvolution) provide robust new insights into the cross-country-wide fault system and the spatial extent of main structural features. Incorporation of these insights as new geological constraints, together with the assignment of density-depth gradients for the sediment layers to existing 3-D geological model geometries, form the basis for 3-D forward and inverse gravity modelling of the cross-border region. This model forms the context for testing, improvement, and harmonization of the structural input models and, where necessary, competing geological scenarios related to them. The resulting parametrized 3-D geological model is consistent with the measured gravity data and reveals new cross-country-wide insights into structures of the sedimentary cover and the crystalline basement in a complex regional tectonic setting. Our main geological findings are: (1) Gradient calculation and filtering of the gravity data allow us to continue fault interpretations from vintage seismic reflection throughout the study region, and so to give a detailed impression of its dominating fault systems, which mainly comprise NW-SE and NE-SW striking faults. Gravity interpretation indicates that some of these faults reach into deeper crustal layers, consistent with published models of seismic refraction data.
(2) Modelling of the long wavelength gravity anomalies, together with interpretation of gravity gradients, provides a map of the regional basement topography and allows for delineation of basement zones of contrasting density. Combined interpretation with petrophysical information from the monzogranites and granodiorites of the PPCC indicates that a low-density zone extends into the NGB's marginal area. Rocks of higher mafic mineral contents and higher densities are modelled in the surrounding area.
(3) The maximum thickness of Zechstein salt is located in the basin's distal parts, where it accumulates in the cores of anticlines and at the Gravity Low of Elster-Schweinitz. East of the Seyda Fault, the Zechstein formation is mainly represented by massive anhydrites and dolomites, which substitute the halite-bearing formations and thus represent either the margin of the ancient Zechstein basin or comprise evidence for considerable salt mobility in the basin.
(4) Anticlines in the northern part of the research area represent salt pillows over a flat Zechstein base, similar to anticlines in the Altmark region further northwest. In contrast, an elevated faulted basement is observed in the core of an anticline at the Elster Fault. We suggest that the different style of this can be related to focusing of Late Cretaceous compressional stress at the NGB's boundary or to the lower rigidity of the MGCR basement in comparison to that of the Northern Phyllite Zone.
(5) A pronounced syncline with an up to 800 m deep Cenozoic basin is observed between the Elster and Zahna Fault. The formation of this basin is likely related to halokinesis or subrosion of Zechstein evaporites.

A C K N O W L E D G E M E N T S
Research was performed in the scope of the TUNB (Deep underground of the North German Basin) project, which is funded by the BGR (Federal Institute for Geosciences and Natural Resources, Germany). Additional financial support is provided by the European Commission under the ERANET Cofound action 'Establishing the European Geological Surveys Research Area to deliver a Geological Service for Europe ' (GeoERA;Grant No.: 731166). Neptune Energy (formerly ENGIE) is thanked for data supply and permission for publishing results based on the provided data. Special thanks to Hans-Jürgen Götze, Sabine Schmidt and the sediment basin modelling group of the GFZ (German Research Centre for Geosciences) for providing the interactive gravity and magnetic modelling software IGMAS+. Furthermore, we thank Emerson-Paradigm, whose SKUA-GOCAD software was used for academic research under the non-profit organizations license agreement. The TUNB working group of the LBGR (Geological Survey of Brandenburg, Germany) is very much thanked for the supply of gravity data as well as relevant parts of their state-wide 3-D geological model. We are grateful to the editor, two anonymous reviewers and Magdalena Scheck-Wenderoth for their valuable comments and time they spent on reviewing the manuscript. Finally, we thank Graeme Eagles for improving the manuscript in content and grammar.

D E C L A R AT I O N O F C O M P E T I N G I N T E R E S T
The authors declare that they have no competing financial interests or personal relationships that could have influenced the work reported in this paper.

DATA AVA I L A B I L I T Y
The data underlying this paper will be shared on reasonable request to the corresponding author. Table A1. Equations for the digitized best-fitting curves of the density-depth relation charts for sedimentary rocks in the eastern part of the NGB by Köhler & Eichner (1973). The equations were assigned to the voxels of the density model to account for increasing densities within a layer, due to increasing compaction. Variable z is the depth in meters (e.g. −2340 m) and ρ is the density in g cm −3 . The depth interval indicates the depth range to which the equations can be applied for density calculation. In general, the last summand of the formulae needs to be deleted to obtain only the density gradient.

B1. Results on the local fault system
We interpreted the fault system in the study area based mainly on the horizontal and vertical gradients of the CBA and inspection of clusters of Euler depth solutions (Figs 9-11). The calculated gradients highlight several mainly NW-SE and NE-SW trending features of differing extents and amplitudes. Local interpretations of faults in seismic data (Wohlfeil & Teicher 1969;Fig. 9a, red lines) coincide with the observed gradient maxima, which encourages us to interpret the maxima as gravimetric expressions of faults that offset layers of contrasting densities. This combined analysis of seismic and gravity data allows for profound interpretation of fault orientation and extent in 3-D space and reveals substantial new information on the local fault system.
A wide, high-amplitude gradient high is observed for the NW-SE trending Wittenberg Fault (Fig. 9a), which we interpret as the main fault in the study area. The HDR maximum peaks in its northwestern part, where the 45 ± 10 • dipping major reverse fault juxtaposes Mesozoic strata against Palaeozoic basement over a height of up to 3 km. A second HDR maximum in the wavelength filtered CBA (Fig. 9b, left light grey plane with question mark) implies the existence of an additional nearly vertical fault, confirming the interpretation of such a fault in a previous seismic refraction study along the nearby THUBRA profile (Bormann et al. 1986;Fig. 3c).
The prominent HDR signature of the Wittenberg Fault decreases markedly southeast of Wittenberg (Fig. 9a). This arises from (i) a reduced fault offset, (ii) a change of basement density, and (iii) the fault crosscutting the low-density zone of the PPCC (Figs 16  and 17b). Additionally, narrowing of the HDR signal, indicates steepening of the fault to a dip of 65 ± 10 • , a conclusion that can be supported by interpretation of a seismic reflection profile close to Wittenberg (Wohlfeil & Teicher 1969;Fig. 2a). This trend continues towards the Annaburg Fault (Fig. 9a) in the southeastern part of the study area. Here, the Wittenberg Fault probably displaces crustal blocks down to about 20 km depth (Fig. 9b), which is consistent with the interpretation of the seismic refraction profile GRIMBU II (Frischbutter & Lück 1997) further southeast of the study area (Fig. 3c).
The Wittenberg Fault is accompanied by NW-SE trending subordinated faults, which offset the Mesozoic to Late Palaeozoic stratigraphic pile (Fig. 9a). Interpreted seismic reflection profiles (Wohlfeil & Teicher 1969) and the results of the seismic refraction profiles THUBRA and GRIMBU II (Bormann et al. 1986;Frischbutter & Lück 1997;Fig. 3c) imply these faults to be vertical to steeply SW-dipping features (e.g. Elster Fault). However, our gravity modelling revealed even gently NE-dipping faults in the more distal parts of the basin in the study region (e.g. the Zahna and Schweinitz faults).
Additionally, we observe a cluster of NE-SW oriented faults in the HDR. Beside several minor NE-SW oriented gradient maxima, a continuous narrow maximum traces the Seyda Fault (Fig. 9a). This reverse fault offsets the Mesozoic strata by up to 200 m in its northeastern part but by less and less towards the south. Here, gravity modelling reveals a continuous increase of the basement surface and a facies boundary of the Zechstein formation rather than fault offsets.

B2. Results on main geological structures and gravity anomalies
Modelling of the regional gravity field's long wavelength anomalies provides new estimates of the density distribution and crustal architecture in the study area (Figs 16 and 17). This required the implementation of zones with different densities in the upper crust. A density of 2.63 ± 0.03 g cm −3 is assigned to model the slightly negative gravity anomaly in the area of the PPCC (Figs 4b, 16 and 17b). This zone is bordered to the south by the Annaburg and Kemberg faults, in agreement with previous geological studies (Eichner 1981;Becker et al. 1989;Fig. 17b). Although only tentatively interpretable, owing to its location at the study area's outer edge, the low-density zone is interpreted to be slightly thinner in its southeastern part than previously postulated on the basis of borehole interpretations (Becker et al. 1989;Röllig et al. 1990). Interfering gravity signals of different depths and orientations hamper the identification of the western boundary of the low-density zone (Fig. 17a). The density model however clearly indicates that the zone continues towards the north, where it is bordered by the Zahna Fault. We interpret this as the northward extension of the PPCC into the marginal area of the NGB. Gravity modelling suggests a density of 2.70 ± 0.03 g cm −3 both west and east of the low-density zone (Figs 16 and 17b). Higher densities, of 2.73 ± 0.05 and 2.77 ± 0.03 g cm −3 , are modelled in the southern and northern parts of the research area. This highly variable upper crust shows a distinct top (crystalline) basement topography (Fig. 16) that is significantly displaced along the Wittenberg Fault and at the PPCC (up to 3 km offset). A stepwise drop of the basement is observed in the eastern part of the study area (0.5-1 km offset).
The general structural setting of the sedimentary cover is strongly folded. A prominent anticline is observed in the density model along the Elster Fault, where seismic reflection profiles indicate several subvertical faults to uplift the crystalline basement to the Mesozoic strata (Wohlfeil & Teicher 1969). In the centre of the anticline, the seismic sections only provide evidence for the Cenozoic strata. Here, the interpretation of slightly increased gravity gradients and a clustering of southwest dipping Euler depth solutions (SI = 0, vertical contact; Fig. 18B, km 10) help to unravel the subsurface structure. To do so, contrasting geological scenarios were tested to investigate the structure of the anticline (Fig. 14). Scenario 1 assumes a uniformly thick salt layer across the anticline; all subsalt and supra-salt layers show a similar layering. Scenario 2 implies a significantly thickened Staßfurt Salt layer, similar to salt pillows or salt-cored detachment folds as known from the northern foreland of the Gardelegen Fault further to the northwest (Malz et al. 2020). In scenario 3, the anticlinal structure is cut by a set of subvertical faults, as indicated by vintage reflection seismic profiles (Wohlfeil & Teicher 1969).
Residual anomalies between observed and calculated free-air gravity anomaly of the different scenarios reveal the highest misfit for scenario 2 (Fig. 14), which we therefore consider unlikely. Scenarios 1 and 3 are more likely since each result in acceptable fits to the observed gravity (Fig. 14). We discuss these scenarios further in the context of the regional geological framework in Appendix C2. Modelling of anticlines in the northern part of the research area (Fig. 19), reveal a similar setting to that known from the Altmark region, with a flat Zechstein base and thick salt accumulations in the anticlinal cores (scenario 2, salt pillow).
In general, the salt is thickest (> 500 m) close to the boundary of the Gravity Low of Elster-Schweinitz and in the anticlines of the northern part of the study area. Minor salt accumulation (< 250 m) is observed directly west of the Seyda Fault (Fig. 19). In contrast, east of the Seyda Fault, where the basement is elevated between the Wittenberg Fault and the Schönewalde Fault, gravity modelling indicates very thin salt (< 40 m), with its place taken instead by anhydrite and dolomite (mostly of the Werra formation, z1), which is locally proved by drilling information (Figs 18E and 19).
The prominent Gravity Low of Elster-Schweinitz is clearly caused by the combination of several anomaly sources: (i) the low basement density (2.63 g cm −3 ), (ii) the smooth synclinal structure northeast of the anticline at the Elster Fault, and (iii) the fill of Cenozoic sediments in the syncline core (Fig. 18A). Gravity modelling suggests increasing thicknesses of the Eocene and Palaeocene sediments within the Cenozoic basin, which is supported by two boreholes. These strata reach a maximum thickness of 550 m in the centre of the gravity low. Densities of 2.1-2.25 g cm −3 in the thickened Eocene and Palaeocene layers are low, compared to the surrounding strata. The Cenozoic basin continues towards the Gravity Low of Belzig. Here, the final density model indicates a clear deepening of the basin, which is filled with Jurassic and younger strata of low densities (Fig. 18A).

C1. Discussion on basement setting
Gravity modelling reveals areas of contrasting densities in the upper crust, which may indicate differences of composition and origin (Figs 16 and 17). A remarkable continuous zone of relatively low density of 2.63 ± 0.03 g cm −3 extends along much of the PPCC and further north into the marginal area of the NGB (Figs 16  and 17). Density measurements of PPCC rock samples (average density 2.63 g cm −3 (Fig. 6; Carl 1988)) and a gamma-gamma-log (2.60-2.65 g cm −3 ; Figs 5b and 17b) in the NGB's marginal area support the modelled densities. Petrographic studies of rock samples from this area (Fig. 17b) indicate quartz contents of 30-40 per cent, more plagioclase than alkali feldspar, biotite content of 1-10 per cent, and a porosity of 1-2 per cent (Erzberger et al. 1967a;Erzberger et al. 1967b;Carl 1988;Röber et al. 1997). Based on the modelled density and petrographic composition, we suggest that the central part of the PPCC consists of a biotite granodiorite, which extends northwards into the marginal area of the NGB.
East and west of the PPCC, a modelled density of 2.70 ± 0.03 g cm −3 is consistent with laboratory density measurements in a range 2.66-2.74 g cm −3 (Fig. 6) for described monzogranites, quartz monzonites and granodiorites in the research area (Carl 1988). Petrographic studies indicate a slightly lower portion of quartz, a balanced ratio of feldspars and a higher portion of dense mafic minerals like biotite (5-19 per cent) and amphibole (up to 26 per cent; Röber et al. 1996Röber et al. , 1997Eidam et al. 1998), which contribute to increased densities compared to the PPCC. South of the Annaburg Fault, drilled diorites and granodiorites with a considerable proportion of amphibole and a modelled density of 2.70 ± 0.03 g cm −3 represent the transition area to the Southern Phyllite Zone, with its phyllites and mica schists (2.73 ± 0.05g cm −3 ; Fig. 17b). Basement densities are more tentatively assigned at the outer border of the modelling area and associated with the Gravity High of Wittenberg (to 2.79 ± 0.05 g cm −3 ) and the area east of the Schönewalde Fault (to 2.74 ± 0.05 g cm −3 ). A better understanding of the origin of these structures (Fig. 17b) would require an extension of the model area, which is beyond the scope of this study.
The transition between the MGCR and Northern Phyllite Zone is expected to lie in the northern part of the study area (Figs 16 and 17;e.g. Oncken 1997;Reinhold 2005;Franke 2014). Gravity modelling of the long wavelength anomalies of the free-air gravity anomaly in this region is possible with a single density of 2.77 ± 0.03 g cm −3 and does not require additional zoning. However, interpretation of the tilt angle (Fig. 17a) and HDR (Fig. 9a) indicate the existence of separate basement blocks in this area. Such high densities of plutonic rocks in the MGCR are known for quartz monzonites, quartz diorites, and tonalities enriched with > 15 per cent of mafic minerals, as observed south of the PPCC. The subdued magnetic field over the area of the Northern Phyllite Zone excludes the presence of normal diorites like those known from the adjacent outer rim of the Prettin and Dessau complexes further west, which are characterized by high susceptibilities and distinct magnetic anomaly highs. Alternatively, densities ranging between 2.72-2.79 g cm −3 are known for phyllites and various schists from the Southern Phyllite Zone Figure C1. Comparison between (a) Local field of CBA (highpass-filtered for wavelengths shorter than 13 km) with marked prominent gravity highs and (b) Schematic illustration of modelled main structural features. The encircled gravity highs of the local field mark the tip of tilted basement blocks (H1), encountered amphibolitic diorites and granodiorites (H2), drilled clay shales (H3) as well as areas of thin Cenozoic sediment coverage or increased content of mafic minerals (H4). (Becker et al. 1989), which therefore may also be present in the transition area between the PPCC and the Northern Phyllite Zone.

C2. Discussion on fault system and style of main structural features
The main faults in the central and southern part of the study strike NW-SE (Figs 9 and 10), parallel to the Elbe Fault System (Fig. 3b). Although known to offset Mesozoic strata, some of these faults are also assumed to reach deeper, having been established during east-west extension/transtension and initial formation of the NGB in Upper Carboniferous times (Franke & Hoffmann 1999a;Scheck et al. 2002).
Minor faults observed in the Mesozoic stratigraphic column show irregular trends (Fig. C1a), which may indicate a variable polyphase deformation history (Franke & Hoffmann 1999a;Scheck et al. 2002). In this regard, N-S and NE-SW trending anticlines in the northern part of the study area are assumed to be located north of the MGCR in the area of the Northern Phyllite Zone (Fig. C1b, yellow features). Similar to the Altmark region, thick salt accumulations are observable in these anticline's cores; that is they are formed as salt pillows or detachment folds. Based on regional structural and stratigraphic correlations, Malz et al. (2020) provide timing constraints for the initiation of salt movements in the adjacent Altmark region, where salt pillow formation may be associated with Late Jurassic to Early Cretaceous NE-SW directed extension. Alternatively, N-S directed compression in the Late Cretaceous led to folding of the sedimentary cover and maybe even initiated halokinesis (Ziegler 1990;Malz et al. 2020).
In the central and eastern part of the research area, our results indicate the presence of relatively small (a few kilometres wide) spurs of uplifted basement rocks, which dip towards the NE (Fig. C1). Similar small basement uplifts can be observed at several locations along the northern boundary of the MGCR (e.g. the Paschleben Spur in the Subhercynian Basin or the Hornburg Anticline and the Kyffhäuser Mountains south and southeast of the Harz; Fig. 3b) and probably indicate contrasting mechanical properties in the basement rocks. Folding and basement uplift in Central Europe is well known to be caused by Late Cretaceous compressional forces, due to the convergence of the African, Iberian and European plates (Ziegler 1990;Henningsen & Katzung 2006;Kley & Voigt 2008). Hence, we interpret reverse faults in our research area to have formed during late Cretaceous times.
In this context, a further striking feature in our research area is a narrow, short wavelength anticline associated with the Elster Fault. A test of three different geological settings revealed two to sufficiently explain the observed free-air gravity anomaly (Fig. 14). In scenario 3, the crest of the anticline is cut by closely spaced, steeply dipping faults. These faults, if present, must date from Mesozoic or later times, and their steepness suggests they acted in strikeslip, which seem unrealistic with the regional geologic framework in mind. In contrast, scenario 1, in which the cover rocks are all similarly folded over a reverse fault affecting the basement, shows the smallest residual anomaly and most consistency with the dips of clustered Euler solutions (Fig. 14). This scenario differs strongly from those of the Altmark region and the northern part of our research area, which typically formed due to thin-skinned contraction or salt pillow formation without basement involvement (Malz et al. 2020). Nevertheless, modelling results of the Elster Fault raise the question to the particular mechanical and kinematic conditions resulting in similar structures along the MGCR. Probably the contrasting anticline structures are responses to the varying rigidity of basement rocks. On one hand, these contrasts promoted the development of a distributed fault pattern in the MGCR areas that permitted folding down to the basement. On the other hand, larger fault spacing is observed in the Northern Phyllite Zone and Rhenohercynian crust, where Cretaceous compression mostly affected the supra-salt sediments and relatively few large-offset basement faults evolved.
A pronounced syncline with a remarkable Cenozoic basin is observed in the area between the anticline along the Elster Fault and the Zahna Fault (Fig. C1b). Similar oval shaped depressions filled by Cenozoic sediments adjacent to salt-cored anticlines are also reported in the Altmark region (Malz et al. 2020). The depression's shape and setting may be either referred to ongoing halokinetic movements in the Cenozoic (Malz et al. 2020) or to subrosion of the Zechstein evaporites.