A distinctive strata-bound vein array occurs in the basal Pierre Shale exposed along the shores of Lake Francis Case, a reservoir on the Missouri River in south-central South Dakota. Typically 2–4 meters in thickness, the array consistently outcrops over a >50-km distance, a significant areal footprint. Ash layers define the upper and lower bounds of the vein array. Two, suborthogonal, preferred directions of vertical veins (northeast and southeast strike) define a regional pattern. By volume, vertical veins comprise 1–2% of the rock. Thinner, more discontinuous, and irregular horizontal veins also occur. Comparisons between array orientations and the joint/vein pattern in the immediately underlying marls of the top of the Niobrara Chalk identify distinct differences. Traverse data suggest that the vein arrays are characterized by uniform horizontal extension. Vertical veins in the array are typically 1–2 centimeters thick and contain massive jarosite, selenite, and fibrous gypsum. The abundance of jarosite and fibrous gypsum distinctly correlates with position in the weathering profile, and these phases are interpreted as due to replacement of original selenite during modern weathering. However, for initial vein array formation, the following suggests that they are not related to modern weathering and formed at depth: (1) a lack of correlation of vein width/frequency with position in the weathering profile; (2) the regional extent; (3) the consistent preferred orientations; (4) the uniform horizontal extension; and (5) the coarse-grained character of the selenite. The consistent strike pattern suggests influence of a regional stress field. The mechanism/timing of vein array formation is unclear. Formation due to diagenetic processes, which are especially significant in mud rocks, would explain the strata-bound character and isotropic horizontal strain and is considered most likely. Formation during glacial loading is one intriguing possibility. Localization of the vein array may be due to the organic-rich character of the host Burning Brule Member of the Sharon Springs Formation.


While quite common, veins found in the Cretaceous Pierre Shale and underlying Niobrara Chalk in the Great Plains have not been studied in any significant way in terms of their orientation patterns, age of formation, and genesis. Because of tight shale plays in the petroleum industry, there is renewed interest in the fracture mechanics and history of mud rocks in general (Gale et al., 2014), and major advances have been made in understanding how fine-grained mud rocks, in particular, are prone to diagenetically driven deformation (e.g., Laubach et al., 2010; Cartwright, 2011). Examples of such deformation include the Cretaceous strata in the Great Plains. To the west of the study area, in the Denver-Julesburg Basin, the Niobrara Chalk produces significant oil and gas and is described as containing polygonal fault systems that formed during early burial diagenesis (Sonnenberg and Underwood, 2013). Maher (2014) also argues that widely distributed normal faults in the Niobrara Chalk of South Dakota, Nebraska, and Kansas are part of polygonal fault systems.

In this context, this study: (1) describes a distinctive vein array with a particularly large footprint that is well exposed for more than 50 km along the shores of Lake Francis Case, an approximate 170-km-long reservoir behind Fort Randall Dam on the Missouri River in south-central South Dakota (Fig. 1); (2) argues that it is related to diagenetically driven deformation; and (3) explores ideas on when and why the vein array developed. The shorelines provide excellent exposures of rocks that are otherwise very poorly exposed.


Exposures along Lake Francis Case shorelines include the top of the Niobrara Chalk and overlying Pierre Shale (Nichols et al., 1994; Shurr et al., 1994; Bertog, 2012), both deposited in the Cretaceous Western Interior Seaway. The contact is a disconformity. The base of the overlying Pierre Shale is composed of the Burning Brule Member within the Sharon Springs Formation, the former of which is up to 24 m thick (Bertog, 2012). The Burning Brule Member is a particularly organic-rich unit that has produced underground fires, likely due to the oxidation of pyrite, and in the study area, it varies from <1 up to 6 m thick. Within the basal portion of the Burning Brule Member are a series of volcanic ash layers, known as the Ardmore series. The vein array is distinctly strata-bound and consistently occurs within the Burning Brule Member, with the Ardmore ash series at its base, and another thinner ash series at its top (Fig. 2). Herein, we will refer to this structural association as the Basal Pierre vein array (BPVA). The overlying unit is the Boyer Bay Member within the Sharon Springs Formation (Bertog, 2012). The underlying Niobrara Chalk is actually a marl with a significant clay component. Veins are also locally distributed within it, but these have a distinctly different character, as described below. Overlying units in the area consist of Pleistocene tills and outwash deposits as well as younger loesses.


The BPVA consists of a 2- to 4-m-thick interval with multiple directions of strata-bound, sub-vertical, intersecting veins that range from millimeters to several centimeters thick (Figs. 2–3). Vein spacing is often tens of centimeters to a meter. Thinner and more irregular horizontal veins are also common. Together, they produce a box-work fabric (Fig. 3). Shore exposures of the BPVA occur from the area where South Dakota Highway 44 (SD 44) crosses the lake in the south, northward to just south of Chamberlain, a linear distance of minimally 50 km (Fig. 1). At all the investigated localities within this stretch, the BPVA array has been present, and we interpret this structural feature to be continuous over this distance. Its extent away from the river is unconstrained because of a lack of outcrops at an appropriate stratigraphic level. Despite excellent outcrops at the appropriate level, the BPVA is absent along the lake shore north of Chamberlain. The Burning Brule Member is condensed here (Bertog, 2012), which suggests a distinct stratigraphic control on its development.

Abundant, massive, fine-grained jarosite dominates many of the veins, but other veins also contain clear selenite, and mixtures of the two are common (Fig. 3). Due to recent wave erosion associated with lake impoundment, topographic spurs are truncated by the cliffs so that the material in the cliff core is less weathered, and weathering increases toward the spur margins. The amount of jarosite distinctly increases with proximity to a weathering surface. The veins dominated by coarse selenite have very planar and matching walls, consistent with a mode 1 tensile fracture (Fig. 3C). Vein spacing is independent of the degree of weathering and the proximity to weathering surfaces. Veins dominated by jarosite also have more irregular and knobby walls. In the weathering transition zone, mixtures of jarosite and selenite were observed, usually with the jarosite on the vein margins and the selenite in the interior. These relationships clearly indicate to us that the jarosite is a weathering replacement of original selenite veins. Jarosite is known to be a relatively common weathering product of black shales under oxidizing, wet conditions (e.g., Schultz et al., 1980; Joeckel et al., 2007), and it is of renewed interest because of its discovery on Mars (Elwood Madden et al., 2004).

The vertical veins typically traverse the entirety of the shale between constraining ash layers, which they truncate against. Some of the veins locally traverse some of the thinner ash layers and continue into the next shale bed. Using the fracture-height classification of Hooker et al. (2013), as described in Gale et al. (2014), they form a hierarchal pattern where the thicker veins transgress some of the thinner ashes of the Ardmore series (Fig. 3). The veins appear to retain their thickness throughout the shale interval, although at one outcrop some tapering from the bounding ash beds into the shale was observed. The horizontal veins are thinner, much more discontinuous, and often wavy, and they both merge with and cross cut the vertical veins, which suggests that they are at least partially coeval.

In map view, one vertical fracture direction can be suborthogonal to and truncate against the other, forming a longitudinal and cross-fracture ‘ladder’ pattern (Fig. 3); however, the orientation of the longitudinal set varies (preferred orientations of the veins are discussed below). More complex polygonal patterns also exist in map view in close proximity (within meters) to orthogonal patterns. Mostly due to the jarosite overprint, any cross-cutting relationships for the vertical veins were obscured.

Topographically, the BPVA appears to be resistant to erosion, forming distinct bench cliffs along the shoreline, although the ash layers may also play a role in this topographic expression, and the bench is not present to the north where the BPVA is absent. In a number of localities, the topographic benches above the cliffs are capped by Quaternary deposits and soils. The BPVA can have a distinct yellowish color from a distance due to the jarosite replacement (Fig. 2).


Structural analysis was conducted along the shoreline in two general areas separated by about 36 km (Figs. 1 and 4): (1) N and south of where SD 44 crosses the lake (sites 1–3, Fig. 1); and, (2) an area near Boyer GPA (Game Production Area) south of Chamberlain (site 5, Fig. 1). Given that the dominant variance was in strike and not dip (the veins were sub-vertical) versions of circular histograms were constructed where, in increments of 1 degree each, ray length is a function of the number of readings in the data set 15 degrees to either side of that ray (Fig. 4). This approach has several advantages, including a better definition of the preferred direction, a lack of dependence on bin-bound positions, and the ability to better model polymodal statistical distributions. A statistical model composed of up to four normally distributed components and a uniform component was constructed by reiterative best-fit matching by eye between the observed plot and the model. The average model versus observed difference was 4% and 7% for the two study areas, and 13% for joints in the underlying Niobrara (described below). The two study areas show well-developed preferred directions at 135, 140, 43, and 47 degrees azimuth, which are interpreted as reflecting two regional, orthogonal directions. Suborthogonal geometries were also seen at the scale of individual outcrops (Fig. 4).

The asymmetric development of the preferred orientations in the Boyer Game Production Area (site 5, Fig. 1) is distinctly a function of outcrop type (mainly vertical cliffs) and outcrop face orientation. Cliffs trending northwest–southeast caused the northeast–southwest trending veins to be encountered and measured more often, resulting in the apparent difference in contribution to the distribution of the two directions. In contrast, the outcrops at study sites farther south were more three dimensional, and included sub-horizontal benches that allowed more of a plan view and less of an outcrop bias. The presented plots do not weight by vein width, length, or orientation; however, data from traverses described below indicate that both vein sets are equally developed in terms of extension magnitude.

While an orthogonal regional set is interpreted to dominate the data set, statistical modeling shows that there is a significant additional uniform component (modeled as 31 and 32%), and locally other minor modes can be modeled (Fig. 4), which conforms to field observations. The modeled standard deviations are also larger than standard deviations for joints in well-lithified sedimentary rocks (Engelder and Delteil, 2004).

Orientations were also taken from underlying barren fractures (joints) in the top of Niobrara Chalk meters to tens of meters below the vein network from nine sites along a 50-km stretch of river that includes the sites where measurements on the BPVA were taken (Fig. 4). In aggregate, our data defined two regional preferred directions at 15 and 105 degrees azimuth, with distinctly smaller standard deviations than seen in the BPVA despite having come from more widely separated sites. Observations indicate that the 105-degree direction is primarily a longitudinal joint set (Fig. 4), while the 15-degree direction represents cross joints. The two were clearly asymmetrically developed in outcrops that provided a map view. Shurr et al. (1994) discuss joint sets in the area with a more regional perspective and identify some of the same general directions of preferred orientations, but without the stratigraphic distinctions made herein.

In relatively unweathered outcrops, the joints in the uppermost Niobrara Chalk lack gypsum mineralization except at one locality where selenite was observed. In the more weathered outcrops, fibrous gypsum is ubiquitous within these sub-vertical joint sets (discussed at greater length below). Less weathered surfaces also showed plumose structures, typically with a horizontal propagation axis. The height pattern was hierarchal, with common truncation at layers richer in clay or organics. Typical joint spacings were an order of magnitude greater than the BPVA spacings, but joints also were clustered, and spacing was also a function of bedding thickness.

A close look at the subordinate modeled preferred orientations in the BPVA of the southern area shows a common orientation with one of the two dominant sets in the underlying Niobrara Chalk (e.g., at 9°, Fig. 4). A second subordinate BPVA set differs by 13 degrees from the other dominant joint set in the Niobrara Chalk. A simple interpretation is that the strike distributions at the two levels are fundamentally different, but with some local subordinate common orientations that may represent a temporal overlap of fracture formation or reactivation.


Where possible, traverses were made with traverse orientation, vein position on traverse, vein orientation, vein thickness, and descriptive attributes recorded. An attempt was made to make traverses in different directions. Assuming that the veins are extensional (and not involving host rock replacement or alteration), elongation and shear strains were computed for 37 linear traverses from sites 1–3 and 5 (Fig. 1). Figure 5 shows a composite plot for 20 traverses that had a length greater than 5 m. Plots show significantly greater variability for traverses shorter than 5 m, which is interpreted as a sampling effect, where more variability exists for shorter traverses that sample fewer veins and are less representative.

Horizontal elongation as a function of traverse direction is interpreted to be uniform, and it ranges up to 5% with most values between 1 and 2% (Fig. 5). Given that the veins with a greater amount of jarosite have more irregular margins, in addition to replacing the earlier selenite fill, some replacement of the host wall material and associated vein thickening likely occurred during subsequent weathering. Considering these field observations, a 1% linear strain is considered a conservative estimate. If the strain was due to volume changes (discussed below), the amount would be ∼2%. This analysis ignores any contribution from the thinner horizontal veins.

Shear strains also did not consistently vary with traverse strike, were an order of magnitude lower than the elongations, and have an average very close to zero. Shear strain magnitudes diminished with traverse length, and the longest traverses (about 14, 16, and 19 m) had essentially zero shear strains.


The weathering gradient in large cliff outcrops of the Niobrara Chalk along the shores of Francis Case is much more pronounced than in the Pierre Shale (which weathers more easily and deeply). Truncated topographic spurs expose relatively unweathered material in the core, and more highly weathered material at the flanks (Fig. 6). The fresh Niobrara Chalk marls are dark gray and, when broken, smell of organics; where weathered, they are very light gray to white, with common orange staining. In a traverse from the less to more weathered material, fibrous and porous gypsum veins are absent in the relatively unweathered rocks, and these veins increase in intensity as the degree of weathering increases, clearly indicating they are a weathering product.

These veins come in one of three orientations: (1) sub-vertical and parallel to the prevailing joint directions in the Niobrara Chalk strata (and interpreted to be utilizing pre-existing fracture sets); (2) sub-parallel to the slope and the weathering zone; and (3) sub-parallel to bedding. It is fairly common to see weathering/alteration zones along the sub-vertical veins in the weathering transition zone. These zones grade from grayer to light orange to whiter in coloration with perpendicular approach to the vein (Fig. 6). The veins that are sub-parallel to bedding are often dish-shaped, typically concave downward, and taper upward. Spacing of slope-parallel veins increases with the degree of weathering, locally reaching the centimeters to tens of centimeters range. Similar sets of veins have been seen in the Niobrara Chalk of Kansas and are interpreted as a common weathering phenomenon of these rocks. They have a distinctly different character than veins in the BPVA, even though they can be in close proximity.


At several localities, small fold-thrust structures were seen in association with the Ardmore ash series (Fig. 7). These were typically associated with bedding-parallel slip and appear to deform the jarosite veins, indicating they postdate vein formation. Not enough of these structures were observed to make orientation analysis meaningful. In one locality, they clearly deform pre-existing gypsum/jarosite veins in the BPVA and are considered to be related to modern weathering and swelling associated with drying and wetting of the smectite-rich ash layers.


The following discussion is organized around distinct questions we believe the above results raise.

Is the BPVA Weathering-related, or Did It Form at Depth?

Taken in aggregate, the following arguments lead us to conclude that the veins are not weathering related, but instead formed at some depth. Well-developed regional and sub-vertical preferred orientations are more consistent with a deeper setting than a shallow setting, where the influence of local topography would produce more variable orientations. The large footprint is also consistent with a deeper and more consistent formation environment, whereas an appropriate weathering environment would be localized. There is also the lack of a BPVA equivalent north of Chamberlain despite exposures of the appropriate stratigraphic interval and well-weathered profiles, suggesting some other controlling factor. The degree of jarosite replacement was a function of position in the weathering profile, but vein attributes (spacing, frequency, or strain) are independent of position in landscape or weathering profile. The coarse selenite vein fill instead of fibrous gypsum is consistent with more stable growth conditions that would be expected at depth versus the more fluctuating conditions found within the weathering profile. The relatively uniform horizontal strain also contrasts with a slope-related, non-uniform strain expected if vein formation was weathering related. Gypsum and jarosite that form during weathering are associated with significant expansion and present an engineering challenge (Hoover et al., 2004). Uni-directional expansion toward the slope would thus be expected and not uniform horizontal strain. The contrast between the BPVA and the common fibrous gypsum veins in the weathered portion of the Niobrara Chalk is particularly instructive and is attributed to their very different formation mechanisms.

Is Veining Tectonically or Diagenetically Driven?

The strata-bound character and stratigraphic control of the BPVA are crucial attributes to consider. Fracture systems argued to be generated by diagenetic processes are typically strata-bound and hosted in fine-grained sediment, and the field of structural diagenesis (Laubach et al., 2010) has established the not uncommon interplay between diagenesis and fracture development driven by volume changes and compaction. Best known amongst these fracture systems are perhaps polygonal faults (e.g., Cartwright, 2011), but they can also include vein systems (e.g., Maher and Shuster, 2012) and coal cleats (e.g., Laubach et al., 1998). Polygonal faults occur in the Niobrara Chalk and Pierre Shale and are interpreted as early burial features (Sonnenberg and Underwood, 2013; Maher, 2014). Potential diagenetic mechanisms for generation of these fracture systems include syneresis of colloidal material (silica, clay, and/or organic colloids), changes in smectite mineralogy, dewatering and compaction, and maturation of organics (Cartwright, 2011). Often, the specific diagenetic mechanism has not been clearly identified for a given occurrence.

Fractures and veins can be localized in more competent layers and have unique strata-bound patterns, a phenomenon captured in the concept of mechanical stratigraphy (e.g., Gale et al., 2014). Gale et al. (2014) reviews fractures in shales in general and conclude that higher silica, carbonate, and feldspar content promote fracture development, and that both higher organic and clay content promotes more ductile behavior and inhibits fracture formation with an important exception, and this is if catagenesis drives fracture growth. The Burning Brule Member, which hosts the BPVA, is characterized as “the most organic-rich facies of the lower Pierre Shale Group” (Bertog, 2012, page 9), a trait likely contributing to its natural ignition. The unit is distinctly less cemented than the underlying Niobrara Chalk. The overlying Boyer Bay Member, defined in the study area, is both organic rich and siliceous (Bertog, 2012). Thus, under a regional loading (tectonic or otherwise) the BPVA would not be expected to concentrate veining and fractures. A localized layer-specific loading due to catagenesis or other diagenetically related volume change remains as a reasonable explanation

The indication of a relatively uniform horizontal elongation with very low shear strains for the BPVA is also consistent with a diagenetic origin driven by volume change/shrinkage. In this case, the vein extension is balanced by an unmeasured contraction in the host sediment and the actual layer elongation at a large scale is zero. In contrast, a uniform horizontal tectonic strain would be a special case, and the localization of tectonic strain into this particular, thin stratigraphic interval would be unexplained. Evidence for equivalent strain in adjacent strata is missing, and the magnitude of strain is well beyond that typically associated with regional cratonic fracture systems, but within the range of that expected by volume change-driven diagenetic deformation.

Fracture systems generated by diagenetic processes can have an unorganized pattern, or they can—and commonly do—exhibit local preferred directions. The preferred directions can be attributed to regional slope stresses or tectonic stresses that organize the fracture pattern (e.g., Laubach et al., 1998, for coal cleats; Møller Hansen et al., 2004, for polygonal faults; Maher and Shuster, 2012, for chalcedony veins). In such cases, while the fractures can be organized by and track the stress field at the time of fracturing and diagenesis, the loading that produces the fractures is interpreted to be due to local diagenetic processes and, specifically, volume contraction (shrinkage). The larger modeled standard deviations and uniform components for the BPVA strike distribution may also be characteristic of a diagenetic origin. Strata-bound chalcedony vein arrays in smectite-rich Tertiary strata in western Nebraska and South Dakota that are argued to be of diagenetic origin (Maher and Shuster, 2012) show significantly greater strike variability and a greater uniform component than is typical of fracture systems created by tectonic loading.

Taken together, an origin by diagenetic processes could explain the strata-bound character, the stratigraphic position in an incompetent unit, the large footprint, the uniform horizontal strain, the magnitude of strain, and the strike variability.

When and Why Did the BPVA form?

There is presently little direct constraint on the timing of formation of the BPVA in the Lake Francis Case area, and fractures, in general, are difficult to date (Gale et al., 2014). Two indirect arguments that consider the geochemical conditions that produce gypsum and the directions of the preferred orientations support a more speculative possibility: that of formation during glaciation and glacial loading. Gypsum development in organic-rich shales such as the Pierre Shale is often attributed to the oxidation of pyrite and production of sulfate, coupled with available Ca, for which the immediately underlying Niobrara Chalk is an obvious source. Sulfate in proglacial sediments in Svalbard has been attributed to sub-ice weathering of the bedrock (Szynkiewicz et al., 2013). Bain (1990) describes diagenetic gypsum development in glacio-lacustrine sediments of Ohio. Continental ice sheets are also known to drive large groundwater flow systems, generally outward and from beneath the glacier margin (Boulton et al., 1995; Provost et al., 2012). Establishment of an oxidative groundwater flow system could be expected to induce shallow diagenetic changes. While the Pierre Shale is typically considered impermeable, fluid access along fractures may have occurred, and the permeability in the immediately underlying chalks is higher. Increased pore pressures due to the overlying glacier may have also both kept fractures open and driven flow. In places, gypsum veins transgress the ash layers, indicating hydrologic breaching of these impermeable layers.

The study area is near the terminus of the North America Pleistocene ice sheets, which the Missouri River broadly follows. Some early Pleistocene glacial deposits occur west of the Missouri River (Fig. 1), and, thus, it is reasonable to infer associated glacial loading for the study area during earlier glaciations. The vein network pattern is aligned with the ice sheet front, with one preferred vein orientation sub-parallel to the ice sheet margin, and the other sub-perpendicular. Whether glacial loading can produce fractures in underlying bedrock, as suggested by Mörner (1978), is in some dispute; however, hydrofractures have been described from deformed glaciogenic sediments in Ireland that were overridden in an advance (Rijsdijk et al., 2010). With their weak character, the Pierre shale is considered more like these sediments than bedrock. Clark (1982) also argues that crustal flexure associated with continental glaciation exceeded the strength of Devonian shales in the northeast United States, producing fractures that now play a role in shale gas production.

Pyrite oxidation and associated gypsum development are typically associated with a volume increase and expansion (e.g., Hoover et al., 2004). This trait is seemingly at odds with a shrinkage origin and needs to be considered. The gypsum formation may reflect diagenetic processes in adjacent strata from the underlying Niobrara Chalk and overlying Pierre Shale with fluid migration along fractures into the BPVA, which could have undergone a different suite of diagenetic processes with net shrinkage. The particularly organic-rich nature of the Burning Brule Member is perhaps significant in localization of the vein array, and a possible diagenetic candidate to consider is sub-surface oxidation of the organics. Petsch et al. (2001) document an aerobic microbial role for kerogen degradation in a modern weathering profile that may play an important role in the carbon cycle. A similar degradation at greater depths seems a reasonable possibility. In a positive feedback loop, initial fracturing of the Burning Brule could induce diagenetic changes that drive further fracturing and fluid flow, which allows additional oxidation of the organics. In addition, a glacially related, more vigorous groundwater flow regime could aid mass transfer aiding volume reduction.

Budai et al. (2002) argue that calcite fill fractures in the Antrim Shale, a gas-producing formation in the Michigan Basin, are directly related to glacially related meteoric recharge along pre-existing fracture networks opened by the glacial loading (or increased pore pressures). Based on the isotopic evidence from the calcite fill, Budai et al. argue that the meteoric recharge initiated methanogenesis. While this represents a different chemical pathway in vein mineralization to what is proposed here for the BPVA, it helps demonstrate the possibility of an interplay between glacial-related fracture and microbial activity.

In summary, during glaciation both loading and a new groundwater flow regime that could induce diagenetic changes are expected, and sulfate and gypsum production is associated with sub-ice glacial waters. Fracturing associated with glacial loading (e.g., Clark, 1982) could provide initial fluid access, inducing shallow diagenesis, which then triggers stratigraphically localized, diagenetically driven fracturing. While poorly constrained, sub-ice formation during glaciation is an initial working hypothesis for the BPVA.


A previously undescribed, extensive vein network with a large footprint exists in the basal Pierre Shale of the Lake Francis Case area of South Dakota. The strata-bound character, a uniform horizontal strain pattern, the magnitude of the strain, and the fine-grained mud host rock are all traits consistent with an origin by diagenetically driven deformation in the subsurface. Jarosite replacement seen in the veins is due to subsequent weathering. The common preferred orientations over substantial distances suggest influence by a regional stress field. The orientation distribution differs significantly from that in the underlying Niobrara Chalk strata, likely reflecting different times of formation. The stratigraphic localization is not accidental, but likely due to the organic-rich character of the Burning Brule Member. The age of vein formation is unconstrained, but on a speculative note, formation during loading by glaciation could produce both the flux of more oxygen-rich waters triggering a phase of diagenesis and the stress field that would align the veins.

The GDL Foundation is thanked for crucial funding that supported students and fieldwork. In addition, Dave Loope, Steve Sonnenberg, and an anonymous reviewer are thanked for review feedback, which led to substantial improvements in this manuscript. The significant editorial guidance of RMG co-editor Art Snoke and RMG managing Brendon Orr is also much appreciated. Any errors are those of the authors.

Gold Open Access: This paper is published under the terms of the CC-BY license.