Giant polygonal anhydrite ridges with diameters of 2–4 km have been identified in the NE margin of the Southern Permian Basin in Poland. The primary ridges reach heights of 80–120 m, and secondary ridges, <40 m in height, subdivide them into smaller cells ~1 km in diameter. The ridges were up to ~40% higher before the alteration of gypsum to anhydrite. Both seismic and gravity surveys indicate the presence of similar structures over an area of thousands of square kilometers. Two alternative hypotheses, gypsum diapirism and free water convection, are proposed to explain their formation.

We report a previously unknown type of evaporitic structure formed by a kilometer-scale polygonal pattern of ridges composed of the Zechstein cyclothem Lower Werra Anhydrite in the NE part of the Southern Permian Basin (SPB). In the traditional view, the development of basal Zechstein evaporites was influenced by the inherited basin relief, with rapid gypsum accumulation on shallow shelves capable of accentuating bathymetric variations (Sonnenfeld, 1984; Richter-Bernburg, 1985). Syndepositional fault reactivation (Rockel and Ziegenhardt, 1979; Paul, 1993) served as an additional factor determining the location of gypsum platforms in the marginal part of the Zechstein basin.

Elongated anhydrite platforms, spanning tens of square kilometers, are a common feature in the Fore-Sudetic homocline at the SE margin of the SPB (Dyjaczyński and Peryt, 2014; Słotwiński and Burliga, 2023). Similar local elevations of the Zechstein anhydrites in the form of domes, swells, diapirs, and ridges have also been documented around the Harz Mountains of Germany (Paul, 2014). Various mechanisms have been proposed to account for their origin, including deformation related to salt doming (Borchert, 1959), primary gypsum diapirism (Paul, 2014), brine flow in porous convection cells (Fulda, 1929), undulations due to gypsum-to-anhydrite conversion (Hemmann, 1972), and tectonic processes such as large-scale folding during basin inversion (Richter-Bernburg, 1985). Our recent discovery may provide new insights into gypsum sedimentation and the formation of anhydrite elevations.

The study area is located in the westernmost part of the East European craton (Fig. 1A). Beneath the Permian complex, there lies a 2-km-thick Ordovician to Silurian shale sequence (Poprawa, 2019). Although the region was affected by Caledonian tectonics (Mazur et al., 2016), the study area hosts mostly flat-lying beds dissected by scarce faults, which terminate within the Silurian section without visibly affecting the Zechstein strata. The Zechstein Sea transgressed over an almost level surface comprised of Upper Silurian shale and, locally, a few meters of Rotliegend conglomerates. The basal Zechstein strata, the Kupferschiefer and the Zechstein Limestone, are 10 m thick, predominantly comprising subtidal deposits (Peryt and Piątkowski, 1977). Subsequently, the area witnessed extensive sulfate and chloride deposition in a peripheral basin (Peryt, 1989), exhibiting a variety of mainly subaqueous facies (Czapowski, 1987; Peryt, 1994). Numerous boreholes in the Puck Bay area allowed the distinction of the Lower Werra Anhydrite (A1d) highs called platforms, their slopes, and basins. A1d platforms, predominantly composed of massive anhydrite, exhibit pseudomorphs after upright-growth gypsum crystals and reach thicknesses up to 150 m (Peryt, 1994). Slope sequences, typically 50–150 m thick, are built of anhydrite displaying graded bedding and brecciation, as well as faulting and folding. Basinal sequences, <50 m thick, comprise brecciated and folded laminated anhydrite.

After the sea-level fall at the end of Zechstein Limestone deposition (Peryt and Piątkowski, 1977), the brine body became very shallow in the Puck Bay area, with a nearly flat depositional surface at the initiation of A1d deposition. However, the notable variation in the A1d facies and thickness may suggest that the relief during deposition of its lower part was significant. The influence of synsedimentary faults on the basin floor and A1d accumulation was previously postulated to account for the observed pattern of sulfate platforms and their substantial thickness variations (Peryt, 1994). The A1d sequences exhibit a deepening-upward trend in both the basinal and slope zones. At the onset of Oldest Halite (Na1) deposition, the basin was differentiated into deeper and shallower segments. The Na1 deposits leveled the relief formed by A1d deposition (Czapowski, 1987), and 20-m-thick uniform Upper Anhydrite (A1g) deposits indicate negligible bathymetric variations at this stage (Peryt, 1994).

The Zechstein deposits, with a total thickness of ~200 m, are covered by an almost undisturbed Mesozoic–Cenozoic sequence, and the base of the A1d deposits lies at depths ranging from 800 m in the north to 1000 m in the south of the study area.

The O-L three-dimensional (3-D) seismic survey revealed a distinctive polygonal pattern of ridges within the A1d deposits in the Puck Bay area (Figs. 2C and 2D). The ridges rise to 80–120 m above the surrounding pans, in which the anhydrite thickness reaches 25–35 m. The ridges display a primary pattern of 2–4-km-wide polygons, locally dissected by minor ridges up to 40 m in height. The slope angles of the primary ridges are up to 20°. The ridges are surrounded and covered by Na1, with a maximum thickness of 210 m. The basement of the A1d ridges appears to be free of elevations or intercepting faults (Figs. 2B and 2C).

Similar polygonal anhydrite ridges are observed in the neighboring areas. South of the study area, the gravity map, enhanced to highlight short-wavelength anomalies anticipated from sources at depths of 500–1000 m (see Fig. 3A), exhibits a polygonal pattern of positive anomalies, which can be associated with triangular anhydrite ridges at 1 km depth, measuring 100–200 m in height and 550–720 m in width (Fig. S10 in the Supplemental Material1). Prominent A1d ridges were also identified on a seismic profile located 80 km east of the study area (see Fig. 3C). Their estimated height is up to 100 m, and the distance between them varies between 1 and 3 km. The variable spacing and the presence of apparent anhydrite platforms could be due to cross-section effects. The polygonal A1d ridges may be widespread over an area of at least several thousand square kilometers within the surveyed portion of the SPB.

The A1d was deposited as gypsum (Peryt, 1994). The gypsum to anhydrite transition is accompanied by a volume loss of up to 39% (Azam, 2007), which manifests as increased subsidence above the ridges in the analyzed seismic sections (Figs. 2A and 2B). The present thickness variations between the originally nearly flat A1d floor and Na1 roof are up to 30–50 m and reflect gypsum thickness reduction, with minima correlating well with the ridge locations (see Fig. 2E). For observed anhydrite ridge heights in the range of 80–120 m, the original height of the gypsum ridges was thus up to 170 m, and it was reduced by 28% on average. Local depressions above the primary ridges began to form after A1g deposition (Fig. 2A), and gentle deflections are evident in the upper part of the Zechstein strata and finally vanish in the lowest Triassic section. Thus, the original gypsum ridges were evidently formed no later than in the latest Permian. The dehydration process commenced when the base of the ridges was at a depth of ~200 m, and it took several million years to cease at 500 m of burial, i.e., slightly shallower than the typical maximum depth of 450–500 m (Jowett et al., 1993).

The formation of polygonal sulfate ridges in the SPB could be attributed to a variety of geologic processes. First, any potential influence of halite diapirism can be ruled out for this part of the SPB (Krzywiec, 2012). Based on the overall symmetric pattern and noncrescentic shape of the ridges, we have rejected the hypothesis of their formation as reticulate gypsum sand dunes (e.g., Szynkiewicz et al., 2010). The lack of clay-rich soil underlying and embedded in the ridges excludes the potential action of diverse processes linked to megapolygonal desiccation cracks forming in playa environments (El-Maarry et al., 2012; Warren, 2016). The underlying impermeable Lower Paleozoic shale is at odds with the possible hypothesis of patterned precipitation above an active fluid-flow system in the basement (Hovland et al., 2006; Zhu et al., 2023). In our discussion here, we consider gypsum diapirism and free water convection, both of which we deem viable, albeit alternative explanations for the formation of the large-scale polygonal ridges discovered in the SPB (Figs. 4A and 4B).

The mechanisms governing gypsum diapirism have been explored by Williams-Stroud and Paul (1997) and further elucidated with supporting evidence by Paul (2014). Following an early diagenetic stage, extending up to a maximum burial of 100–200 m, pure halite has a near-constant density of 2.17 g/cm3 (Warren, 2016), although it could be elevated by impurities. Solid gypsum exhibits a generally higher density of 2.30 g/cm3 (Williams-Stroud and Paul, 1997); however, even consolidated gypsum aggregates retain significant porosity (Sonnenfeld, 1984), which can be further increased due to water influx from underlying sections undergoing compaction and successive transition to anhydrite (Paul, 2014). A 10% porosity level is needed for a saturated gypsum aggregate to reach the density of pure halite. Inverted density stratification drives the Rayleigh-Taylor (R-T) instability, a process that gives rise to a diverse set of sedimentary and tectonic structures, including load casts, ridges, and diapirs (Ramberg, 1972). The growth rate of these structures depends on the density differences, as well as the viscosities and thicknesses of the involved layers (Turcotte and Schubert, 2002). For instance, in the case of a 100-m-thick layer characterized by a 0.1 g/cm3 density contrast and a viscosity on the order of 1017 Pa·s, characteristic of rock salt, the growth rate is sufficient to produce diapiric structures within 1 m.y. While the viscosity of gypsum is higher than that of rock salt in the power-law regime (Williams-Stroud and Paul, 1997), little is known about its rheologic behavior in the pressure solution–dominated, Newtonian regime.

The spacing of diapiric structures primarily depends on layer thicknesses and their viscosities (Turcotte and Schubert, 2002), with large values indicating a high viscosity ratio between the overburden and the buoyant layer. Numerical models of the R-T instability in 3-D (Fernandez and Kaus, 2015) typically demonstrate the formation of a hexagonal pattern of ridges converging at columnar triple junctions, resembling the discovered structures. However, the observed spacing in the SPB falls within the range of 10–20 and 20–40 times the initial thickness of the gypsum layer for the secondary and primary ridges, respectively. The latter spacing is considerably larger than the analytical and numerical predictions for the R-T–driven ridge structures (Turcotte and Schubert, 2002; Fernandez and Kaus, 2015). Hernandez et al. (2018) documented ridge spacing to layer thickness ratios ranging from 10 to 20 for the salt structures in the North Sea. While large diapir spacings in salt tectonics may be linked to the plastic behavior of the clastic overburden (Ismail-Zadeh et al., 2002), this is not a viable explanation for viscously deforming gypsum and halite. The spacing could also be influenced by the initial perturbations of the layer interface (Schmeling, 1987), potentially playing a role in gypsum diapirism. In addition, significant pressurization may arise in a gypsum layer transitioning gradually to anhydrite, which may lead to internal brecciation and the selective reworking of initial diapiric ridges through hydraulic fracturing into piercement structures (Paul, 2014), with spacing surpassing that between the initial diapiric ridges (Fig. 4A). Once transformed into highly viscous anhydrite, the ridges remain intact despite their high density relative to the surroundings.

We also considered the possibility that sulfate ridges may have originated through stable convection in the evaporitic basin. In a quiescent saline basin, a distinct separation occurs between the lighter surface layer, heated to the depth of solar absorption, and the colder, denser water beneath, delineated by a thermocline. The process of evaporation above the thermocline results in an increase in salinity, contributing to a negative water buoyancy. Arnon et al. (2016) demonstrated that the buoyancy resulting from temperature differences outweighs the impact of salinity saturation, establishing a stable thermohaline stratification. In contrast, the Dead Sea experiences seasonal density instability during winter due to surface-water cooling. This seasonal water mixing is unfavorable for the formation of extensive gypsum ridges, which requires stable convection cells. The location of the SPB, positioned near the equator during the late Permian, was conducive to such stability due to intense solar radiation. In this setting, convection may have been initiated by the precipitation of gypsum crystals above the thermocline, causing negative buoyancy in the brine. This is reminiscent of a process observed in estuaries, where rivers carry muddy water denser than the seawater below (Burns and Meiburg, 2015). In numerical simulations, stable convection cells form wherein descending polygonal currents deposit the mud on the seafloor.

Critical questions persist regarding the size and sustainability of the convective cells. To maintain isometric cells, their vertical dimension should be at least half the size of a polygon (see Fig. 4B). For the secondary polygons ~1 km across, the depth of the basin during the convective stage is then estimated to have been around 500 m, which appears to be a very high and hence controversial value. To estimate the timing of gypsum ridge deposition, we considered the evaporation rate, determined to be 1.13 m/yr for the Dead Sea (Hamdani et al., 2018), and the saturated gypsum solution, which is 2.6 kg/m3 at 30 °C. Assuming that 5% of the basin was covered by gypsum ridges with an average height of 50 m, we calculated that the evaporation of an ~2200 m brine column would have been necessary to precipitate this amount of gypsum. This mass evaporation was probably balanced with the inflow of brine into the saline basin. By using a doubled evaporation rate for the SPB, we arrived at an estimate of 1000 yr for gypsum precipitation and stable convection. However, if only the initiation of the polygonal ridge pattern (up to 10% of their final height and width) was controlled by convection, it would have required ~10 yr of stable convection to form. Subsequent ridge growth could be attributed to either gypsum diapirism or conventional gypsum sedimentation over the convection-initiated polygonal elevations of the basin floor.

The elevations of the Lower Werra Anhydrite located in NE Poland in the peripheral part of the SPB, previously documented in boreholes, are, in fact, narrow ridges forming polygons. The primary ridges, measuring 80–120 m in height, and the secondary ridges, ranging from 20 to 40 m in height, collectively create a system of giant polygons with diameters of 2–4 km and 1 km, respectively. The ridges developed on the flat bottom of the evaporitic basin, with no seismic-scale faults observed in the basement. The gypsum to anhydrite transition reduced the original height of the ridges by about one third. Using seismic and gravity data, we identified analogous structures in the neighboring areas, indicating their widespread occurrence in the studied part of the SPB.

Among several considered hypotheses, we suggest two that could potentially explain the formation of sulfate polygons of this size:

  1. Porous gypsum diapirs may potentially grow under consolidated halite. This process could be enhanced by water influx from the basal parts of the growing gypsum ridges successively transforming into anhydrite. The initial diapirs could undergo selective hydraulic reworking into piercement structures, giving rise to high-amplitude ridges.

  2. Free water convection may develop in the evaporite basin due to the negative buoyancy of a suspension of gypsum crystals precipitated in surface water above the thermocline. The crystals could be deposited at the bottom of the basin in a polygonal pattern. Assuming that the pattern is only initiated by convection, only a few years of stable convective cells are required.

1Supplemental Material. Description of geophysical data and additional figures with results of their interpretation. Please visit https://doi.org/10.1130/GEOL.S.27278178 to access the supplemental material; contact [email protected] with any questions.

Piotr Krzywiec and two anonymous reviewers are thanked for constructive comments. We thank ORLEN for providing the geophysical data. The research was supported and the publication fee was financed by The National Fund for Environmental Protection and Water Management. SLB is thanked for providing academic licenses for Petrel.