The effect of stress distribution on the shape and direction of hydraulic fractures in layered media-DTU Orbit (02/11/2019)

media DTU Orbit (02/11/2019) The effect of stress distribution on the shape and direction of hydraulic fractures in layered media In this study, the propagation of a hydraulic fracture in a three-layer medium is investigated using a robust finite element method. The top and bottom layers have similar properties, while the soft middle layer has variable thickness and stiffness. Both two-dimensional plane-strain and three-dimensional models are used to evaluate the vertical and lateral growth of the fracture, initiated from the bottom layer. In 2D model, the competition between the gravity forces and the stress concentration at the layer boundary dictates whether the fracture grows upward or downward, while in the 3D model, the lateral growth of the fracture is also captured. Results from the 2D model show that in majority of cases (6 out of 9) the fracture reaches to the middle layer, and due to the lower minimum principal stress in this layer, an extensive propagation occurs in this layer. When the lateral growth is considered utilising the 3D model, depending on the stress concentration and the position of the fracture in the bottom layer, the downward or lateral growth may become dominant. Thus, the geometry of the fracture can vary between a long and narrow shape towards a round shape. The lateral growth of fractures in the 3D model undermines the suitability of 2D plane-strain model results.


Introduction
Hydraulic fracturing technique has been extensively used to improve reservoir performance in hydrocarbon [1] as well as geothermal reservoirs [2], and for preconditioning the rock masses in underground mines [3]. Accurate modelling of hydraulic fracture geometry in the presence of stress gradients and sudden stress jumps is a fundamental requirement for treatment design. The targeted reservoirs or rock masses are commonly made of several layers which makes a reliable prediction of fracture geometry, based on the conventional single-layer solutions, a formidable task.
Hydraulic fracturing is a complicated process to model as it requires at least three components: (i) deformation of the matrix due to the hydraulic loading of the fluid, (ii) the highly nonlinear flow through the fracture, and (iii) the propagation of the fracture. Additional complexities can be envisaged by taking into account additional (but sometimes important) factors such as presence of layers with different properties, changes in in situ stresses, the leak-off of the injected fluid into the rock matrix, and the presence of a natural fracture or fracture network [4]. Among these factors, the presence of layers and in situ stress variation are studied in this work.
Numerous of studies have shown the effects of layer stiffness, fracture toughness and in situ stress on hydraulic fracture height containment in layered formations. Simonson et al. [5] studied the stress intensity factor as the fracture approached the interface of a layer with different Young's modulus. They found that the stress intensity factor increased to infinity when a fracture approached a low-modulus layer from a high-modulus layer, thus promoting growth toward the low-modulus layer. When the fracture approached a high-modulus layer from a low-modulus layer, the stress intensity factor decreased towards zero, thus making it harder for the fracture to grow towards the high-modulus layer. They therefore stated that a layer with a higher Young's modulus acts as a barrier of fracture propagation from a lower modulus formation layer. Similarly, Brenner and Gudmundsson [6] concluded that the aperture was greater in soft layers than in stiff layers. Gu et al. [7] found that both types of modulus contrast, low-to-high and high-to-low, hinder fracture height growth, but caused by different mechanisms. However, laboratory experiments by Daneshy et al. [8] and in situ experiments carried out by Warpinski et al. [9] showed that contrast in elastic moduli alone was not sufficient to stop the growth of a fracture across an interface, later Smith et al. [10] came to a similar conclusion. Simonson et al. [5] also argued that a difference in in situ stress between two layers could serve as a barrier to vertical extension of the fracture. Teufel and Clark [11] demonstrated that the vertical growth of a hydraulic fracture could be contained by a weak shear strength along the interface of the layers and more importantly the increase in the minimum principal (horizontal) in situ stress in the bounding layer. Laboratory experiments made by Warpinski et al. [12] showed that the in situ stress contrast was sufficient to control the fracture height containment. Although, the gravitational gradient of the minimal principal in situ stress was not considered. Thus, it is generally recognized that the in situ stress contrast is the dominant parameter controlling fracture height growth and that Young's modulus contrast is less important. Daneshy et al. [13] mentioned that blunting of the fracture tip or sliding at the interface also were a likely mechanism for fracture containment.
The conventional single-layer solutions such as KGD [14,15] and PKN [16,17] models were not capable of predicting the fracture footprint in layered media, due to the sensitivity of fracture footprint to the changes in in situ stress across layer interfaces. Thus, the development of Pseudo-3D (P3D) models arose in the 1980s [18][19][20] attempting to model the vertical growth of the fracture across layers in layered formations. P3D models can handle the fracture height growth while maintaining the lateral coupling for fluid pressure, so the initial 3D problem is reduced to a 2D elastic problem and a 1D fluid problem, which reduces the computational time [21]. The model builds on the work of Simonson et al. [5] and is an extension of the PKN model, thus applicable in cases where the height of the fracture remains small compared to its length. P3D assume that the reservoirs elastic properties are homogeneous and averaged over all layers containing the fracture height. This of course is a simplification of a real reservoir which can consist of multiple layers with heterogeneous properties on all relevant scales. To remove these deficiencies, a planar-3D model (PL3D) was first proposed by Clifton et al. [22]. The PL3D models describe the fracture geometry by a 2D mesh of cells, either on a moving triangular mesh or a fixed rectangular mesh. The fluid flow is accounted for in full 2D and coupled to the 3D elastic response of the rock. The fracture propagation is represented in planar direction and oriented perpendicular to the far-field minimum principal in situ stress. Planar-3D models are more accurate than 2D-and pseudo-3D models, but can be computationally expensive. 2D models are less expensive, but the results are often simplistic. The pseudo-3D models provide a compromise and are the most used model in the industry. Other studies have tried to create better and more accurate models, for instance, Zhang et al. [23] presented a new P3D model that aimed to include the effect of elastic modulus contrast. Baykin and Golovin [24] studied the effect of permeability contrast on the hydraulic fracture propagation using a PL3D model with poroelasticity. A more extensive review of the P3D and PL3D models can be found in [1,21,25].
The recent developments in 3D hydraulic fracture modelling have made it possible to investigate previously unknown 3D interaction mechanisms on fracture height growth in a multi-layered medium. Zhang et al. [26] developed a model using a finite element approach where they studied the effect of in situ stress contrast, modulus contrast, tensile strength contrast and viscosity of the fracturing fluid. Often it is assumed that the layers in a layered medium are perfectly bonded together which enables the fracture to penetrate the bedding directly without any deflection. Chang et al. [27] investigated these interfaces and how it affected the behaviour of a propagating fracture. Li et al. [28] used a model based on the finite element method to capture the fracture growth in a layered rocks. Four types of interaction between a fracture and an interface were found: penetration, arrest, T-shaped branching and offset. Khanna et al. [29] studied the situation where a hydraulic fracture was required to be contained in the pay zone, and how the fracturing fluid was able to control the height. They also investigated the interaction between multiple fractures and the fracture spacing. Huang et al. [30] investigated the height growth under the influence of stress barriers and natural fractures using the finite element method for the mechanical response and finite volume method for the fluid flow.
In this study, the problem of hydraulic fracturing in a three-layer media, in which the middle layer is softer than two identical top and bottom layers, is investigated. A robust finite element model is utilised to simulate the hydraulic fracturing process within the Linear Elastic Fracture Mechanics (LEFM) framework [31,32]. Different properties and thicknesses for the middle layer are considered. Both 2D plane-strain and 3D fracture models are utilised for the simulations and the results are presented and discussed.

Governing equations
In the present approach, fractures are modelled as surface discontinuities in the three-dimensional matrix. The deformation model is expressed satisfying the condition of equilibrium on a representative elementary volume (REV) of the porous medium. Fracture surfaces are not traction-free in the present model, and hydraulic loading are applied on the fracture walls. Assuming negligible shear tractions exerted from the fluid on the fracture surfaces, the fluid pressure is applied only in the normal direction to the fracture wall. The tractions on the fracture boundary Γ c are where p f is the fracture pressure, and n c is the outward unit normal to the fracture wall (on both sides of the fracture). The differential equation describing the deformation field for a representative elementary volume (REV) of fractured rock is given by where, D is the drained stiffness matrix, = ∇ + ∇ ε u u ( ) is the strain tensor in the porous medium, u denotes the displacement vector in the porous medium, F is the body force per unit volume, δ is Dirac delta, and x c represents the position of the fracture. Flow through the fracture is defined using the cubic law. The laminar flow through fracture as parallel plates can be written as where, a f is the fracture aperture, μ f is the fluid viscosity and c f is the fluid compressibility. Note that the term provides direct coupling between the displacement field and the fracture flow field, which is symmetric to the fracture pressure loading term, p n f C in Eq. (2). + u and − u are the displacements of the two opposing faces of the fracture. For the case of one-dimensional incompressible flow, the above equation reduces to the lubrication equation [33].

Finite element model
The governing equations are solved numerically using the finite element method. Spatial and temporal discretisation are accomplished using the Galerkin method and finite difference techniques, respectively. Displacements (three displacements for three dimensions) and fracture fluid pressures are defined as the primary variables. Quadratic unstructured elements are used for spatial discretisation of surfaces (quadratic triangles) and volumes (quadratic tetrahedra). The triangles on two opposite surfaces of a fracture are matched with each other, but do not share nodes, and duplicate nodes are defined for two sides of a fracture. The triangles are matched with faces of the tetrahedra connected to the fractures, and they share the same nodes. Using the standard Galerkin method, the primary variable  = p u { , } f within an element is approximated from its nodal values as where N is the vector of shape functions and   is the vector of nodal values. Using the finite difference technique, the time derivative of  is defined as where  + t dt and  t are the values of  at time t + dt and t, respectively. The set of discretised equations can be written in matrix form as   = , in which  is the element's general stiffness matrix, and  is the vector of right-hand-side loadings. Flow through deforming fracture is a highly nonlinear problem, thus a Picard iteration procedure is adopted to reach the correct solution within acceptable tolerance. For the current iteration + s 1 in current timestep + n 1, the solution-dependant coefficient matrices in the stiffness matrix are updated using weighted average solution vector  + The tolerance in the present simulations is set to 0.01, the number of iterations per step decreases from 7 iterations to 3 iterations per step as the time increases. The discretised equations are implemented in the Complex Systems Modelling Platform (CSMP, also known as CSP, [34]), an object-oriented application programme interface (API), for the simulation of complex geological processes and their interactions. The set of linear algebraic equations are solved with the algebraic multigrid method for systems, SAMG [35].

Fracture growth model
Within the framework of Linear Elastic Fracture Mechanics (LEFM), the stress intensity factors (SIFs) for three modes of fracture opening are computed using both the energy-based interaction integral (disk method, [36]) and the displacement correlation method (DC method, [37]). The disk method is less prone to numerical error and yield better approximations with significantly coarser meshes, while the DC method is computationally cheap and is able to yield very good approximations of the SIFs [38]. The accuracy of the SIFs computation can be further improved by using the so-called quarter-point elements [39]. The three SIFs are mode I (K I ) for opening due to tensile loading, mode II (K II ) for in-plane shearing due to sliding, and mode III (K III ) for out-of-plane shearing due to tearing. The crack grows if and when the equivalent stress intensity factor K eq overcomes the material toughness (k ic ). The equivalent SIF in the direction of propagation (θ p ) is calculated as [40] ⎜ ⎟ , and θ p is the propagation angle. The in-plane propagation angle (θ p ) and out-of-plane deflection angle (ψ p ) are determined using a modified maximum circumferential stress method that takes into account modal stress intensity factors under mixed loading [40]. The in-plane propagation angle (θ p ) is assumed to be perpendicular to the maximum circumferential stress σ 1 , thus, θ p can be calculated by [37] ∂ ∂ = ∂ ∂ < σ θ σ θ 0, and 0 The out-of-plane deflection angle ψ p is defined by the orientation of σ 1 as Although the hydraulic fracturing occurs mainly under tensile mode, the three SIFs are computed and used for the propagation angle computation. At every timestep, the SIFs and growth computations are performed at sixty locations along the fracture front, i.e. fracture tips, and if at least at one tip the equivalent SIF reaches the material toughness, the propagation triggers and the fracture tips are advanced according to the equivalent SIF value at each tip. Then, the fracture geometry is updated, the mesh is regenerated and the step is recalculated for the new fracture geometry.

Model validation
The numerical finite element model has been validated extensively against analytical, experimental and numerical results in the literature [31,[41][42][43]. In this section, the laboratory experiment performed by Wu et al. [44] on fracking a three-layer media made of Polymethyl methacrylate (PMMA) blocks is simulated. The experiment was performed to study the propagation of a hydraulic fracture from an intermediate stress layer into higher and lower bounding layers. The PMMA layers were under distinct minimum principal stress ranging from 11.2 MPa in the top layer to 7 MPa in the middle layer and 5 MPa in the bottom layer. The distinct stresses were developed using the method proposed by Jeffrey and Bunger [45]. The layers thicknesses were 40, 50 and 180 mm for top, middle and bottom layers, respectively. The fracture was initiated in the centre of middle layer as shown in Fig. 1. The fluid is injected at a constant rate of 0.0023 ml/s through a 5 mm radius borehole. However, due to the compressibility of the injection system, the actual flow rate entering the fracture was different and was calculated at three stages as 0.0009 ml/s for 0 < t < 31 s, 0.0065 ml/s for 31 s < t < 151 s, and 0.0023 ml/s for 151 s < t [44]. The Young's modulus and Poisson's ratio of the PMMA is measured as 3.3 GPa and 0.4, respectively. The fracture geometry from present simulation at t = 144 s is shown in Fig. 1, together with the measured fracture footprint from the experiment. The induced fracture stops growing into the layer with higher stress, and instead it grows into the lower stress layer. Good agreement is found between both numerical and experiment results that further verifies the present model for hydraulic fracturing in layered media.

Simulation results
In this section, the propagation of both plane-strain (2D) and, radial (penny-shaped) and elliptic (3D) hydraulic fractures in threelayer media is investigated. The effect of stress gradient due to the gravity, as well as, the stress concentrations and jumps at the layer boundaries due to the stiffness contrast between layers, on the geometry and extension of hydraulic fractures are studied.

Model setup
A three-layer model is considered in all simulations with bottom and top layers to be identical and with large thicknesses, while the middle layer is softer and with finite thickness (Fig. 2a) layers. Since the middle layer have lower stiffness, its fracture toughness is set to a larger value (2 MPa m 0.5 ), compared to that of other two layers (1 MPa m 0.5 ), to account for the ductile behaviour [46]. In all cases, the initial fracture is positioned in the bottom layer. The initial fracture in plane-strain and radial cases is at a depth of 10 m below the interface between bottom and middle layers. In elliptic fractures, the depth of the initial fracture has been increased to 20 m and 30 m below the interface. The injection is performed through constant injection rate of Q = 0.001 m 3 /s of a fluid with viscosity μ = 0.0001 Pa s. Two types of far-field minimum principal stresses (horizontal) are applied to the model: (i) uniform stress of 15 MPa, (ii) a linearly increasing stress distribution with a gradient of 15 KPa/m to account for the gravity effects. The other two stresses (intermediate and maximum principal stresses) are given relatively large value (20 MPa) in order to induce a planar fracture. The stress gradient due to the gravity is calculated based on the difference between densities of saturated rock and formation water, assuming density of saturated rock 2500 kg/m 3 , and density of formation water 1000 kg/m 3 . Thus, a stress gradient of 15 KPa/m is applied to the model. The far-field minimum principal stress at the depth of injection (centre of the initial fracture) is set to 15 MPa for all cases. The middle layer is assumed to have a finite thickness and lower stiffness, thus under the far-field principal stresses, a stress concentration is developed in top and bottom layers at the interface with the middle layer, while the developed stress in the middle layer is lower than the applied boundary stress, as shown in Fig. 2b. This is due to higher compliance of the middle layer which results in the stress redistribution in the system. Higher thickness and lower stiffness of the middle layer create higher stress concentrations in the surrounding layers. The maximum excess stress at the boundary of stiff layer next to the middle layer varies between 0.1 and 3 MPa for different cases (Fig. 2b). Lower stiffness of the middle layer reduces the stress in that layer, while higher thickness of the layer slightly increases the stress in the middle layer. The gravity increases the stresses linearly in the layers with increasing depth, thus a local minimum in the stress profile is observed in the lower layer at a depth from the middle layer. Thus, moving upward from the local minimum point results in an increase in the stress due to the presence of a softer middle layer, while the stress also increases below the local minimum point due to the gravity.

Plane-Strain (KGD) fracture
In this section, the propagation of a plane-strain hydraulic fracture, also known as KGD fracture, in a three-layer media as shown in Fig. 2a, is investigated. In Fig. 3, the propagations of the KGD hydraulic fracture, initiated from 10 m below the interface, for different cases of middle layer's thickness and stiffness are shown. The fracture geometries for cases with the thickness of the middle layer equal to 10 m are shown on Fig. 3a. Included in this figure is a case with a homogeneous stiffness for reference. For the homogeneous case with no gravity, the fracture geometry is symmetric with respect to the initiation point as all other model parameters are homogeneous. Gravity pushes the fractures to grow upward (shown by dashed lines), while the stress concentration at the interface prevents or reduces the upward propagation of fractures. Thus for the homogeneous case, the fracture stops advancing downward after some steps due to the gravity effects. In three-layer model, all cases with thickness of the middle layer D = 10 m and gravity effects grow into the middle layer. However, the downward growth of the fracture in these cases is also active and no sign of growth stoppage is observed. Without the gravity effects, the cases with Young's modulus E = 5 and 10 GPa stop growing upward and continue growing downward. Results show that when the thickness of the middle layer is 5 m (Fig. 3b), all the cases reach to the middle layer at some point during the simulation, indicating that the stress concentration is not large enough to stop the upward propagation of the hydraulic fractures. When the hydraulic fracture reaches the middle layer, an extensive propagation of the fracture is observed due to lower stress, even though the fracture toughness is higher in the middle layer.
When the thickness of the middle layer is 20 m (Fig. 3c), majority of the fractures stop growing upward and they continue to grow downward. All cases without gravity effects grow downward. This is due to higher stress concentration near the interface. For instance in the case with lowest stiffness of the middle layer (E = 5 GPa), the fracture hardly grow upward and it mainly advances downward, however, in the case with highest stiffness (E = 15 GPa, i.e., lowest stiffness contrast), under the gravity effects, the fracture reaches to the middle layer, but it stops growing in that layer as the stress increases towards the centre of the layer. The upward height of the hydraulic fractures from the injection point at time t = 6 s for all cases with gravity effects is summarised in Fig. 4a. Majority of cases (6 out of 9 cases) manage to reach to the interface which is located 10 m above the injection point. This shows that the stress concentration at the boundary due to the layer stiffness contrast is not sufficiently large to overcome the gravity force and to prevent upward propagation of hydraulic fractures, unless the stiffness contrast and the thickness of the soft layer is large. For instance, if two reference points, each 10 m below (point A) and above (point B) the injection point is considered, the gravity driver induces 0.3 MPa extra stress at point A compared to point B (located at the interface). By inspecting the vertical distribution of the minimum principal stress in a three-layer media given in Fig. 2b, only in few cases, the stress at point B exceeds the stress at point A by 0.3 MPa. These are the cases that grow mainly downward.
In Fig. 4b, the ratio between downward to upward height advancement of the fractures are shown. When the height ratio is < 1 it means the fracture grows more upward. Thus, for middle layer thickness D = 5 m, the fracture grows more upward than downward as the height ratio < 1, while for D = 10 m, the hydraulic fractures grow more downward than upward as the ratio > 1. Finally, for D = 20 m, the hydraulic fractures grow mainly downward as the height ratio > > 1 and their upward growth is blocked by the high stress concentration at the interface.

Radial (Penny-Shaped) fracture
In this section, the effect of stress distribution on the shape of the radial (planar) fracture is investigated. In the 2D plane-strain (KGD) fracture, the hydraulic fracture has only two tips to grow (upward and downward in the present study), while in the real cases, the fracture can also propagate laterally. All the fractures are initiated including the gravity effects from a horizontal well located in the bottom layer, 10 m below the interface with the middle layer, with an initial radius of 1 m. In Fig. 5, the maximum vertical and lateral sizes of the induced fractures are compared for different values of middle layer thickness and Young's modulus. When the thickness of the middle layer is D = 5 m (Fig. 5a), the fractures show similar advancement in lateral and vertical directions, while in D = 10 m (Fig. 5b) Fig. 6. The dominant downward growth of the fracture can be observed for the case with E = 5 GPa and D = 20 m. In this case, highest concentration of the minimum principal stress is developed near the interface, thus the favourable direction of growth is downward and towards lower minimum principal stresses. As the Young's modulus in the middle layer increases to 10 and 20 GPa, the upward and lateral growth of the fracture also increases. This is due to lower concentration of the stresses near the interface (above the injection well) as a result of lower contrast in the layers stiffness. For the case with E = 5 GPa, when the thickness of the middle layer reduces to 10 and 5 m, again the upward and lateral growth of the fracture increases. This is due to lower stress concentration above the injection well for lower thickness of the middle layer as can be seen in Fig. 2b. Interestingly, in the intermediate cases shown in Fig. 6: (i) E = 10 GPa, D = 20 m, and (ii) E = 5 GPa, D = 10 m, the maximum lateral growth occurs at depths below the injection point where the local minimum stress is the lowest. This means that the both lateral and downward growth are dominant. In addition, when the fracture reaches to the middle layer, an extensive upward growth occurs in that layer due to lower minimum principal stress. For the case with highest stress concentration (E = 5 GPa and D = 20 m) the aperture distribution together with the geometry of the fracture at every propagation step are shown in Fig. 7. The maximum aperture occurs at higher depths where the local minimum principal stress is lower. The upward growth stops in early stages, while the lateral growth continues with lower speed than the downward growth. The lateral growth stops eventually and fracture continue to grow downward towards lower minimum principal stresses. This results in the shape of the fracture becoming long and relatively narrow.

Elliptic fracture
In many cases in real reservoirs, the hydraulic fracture grows along the well. Fractures induced from vertical wells as well as fractures in water injection horizontal wells are some examples. These fractures are used for water injection, and sweeping the oil towards the nearby production wells [47,48]. In these cases, several perforations along the well may be used for fracture initiation, thus the initial fracture does not have a circular disk shape. Commonly, these fractures are modelled in 2D assuming plane-strain (KGD) conditions. However, these fractures may not grow according to those assumptions. In this section the initial fracture is assumed to be an ellipse with major and minor radii equal to 10 m (horizontal) and 1 m (vertical), respectively. The injection is GPa and z = 20 m, the fracture grows in all directions, however, the downward growth is higher than the other directions. The maximum lateral growth occurs at depths greater than the injection well where the local minimum stress is the lowest. When the fracture depth changes to z = 30 m, the lateral growth becomes higher than the vertical growth as the size of the fracture increases. This is due to the fact that the minimum principal stress, shown in Fig. 2b, occurs at depth z = 30 m, and as the fracture grows vertically, the stress increases in both upward and downward directions. Thus the lateral direction becomes the favourable direction for growth. Again back to Fig. 2b, when the Young's modulus of the middle layer reduces to 10 GPa, the stress concentration in the upper part of the bottom layer (above the injection well) increases, so the elliptic fracture grows less upward as shown in Fig. 8. The depth of z = 20 m is not the favourable depth in terms of the local minimum stress, so the lateral growth of the fracture is restricted and the fracture propagates mainly downward. When the depth of the initial fracture increases to z = 30 m in this case, the stress concentration due to the layers stiffness contrast reduces at this depth, and the lateral growth of the fracture increases. However, the  downward growth is still more favourable as the location of the least local minimum principal stress is at a depth higher than 40 m. The lateral growth of the elliptic fractures undermines the plane-strain (KGD) assumption commonly used to model these fractures.

Conclusions
The propagation of hydraulic fractures in a three-layer medium in which the middle layer has lower stiffness is investigated in this study. The geometry of hydraulic fractures under plane-strain conditions (2D), and, radial and elliptic geometries (3D) in several cases of varying thickness and stiffness of the middle layer is studied using a robust finite element model. The 2D plane-strain results show that the gravity effects force the fracture to grow upward, while the stress concentration at the interface of middle layer prevents or reduces the upward growth of the fracture. Since the fracture under plane-strain conditions has only two directions to grow (upward and downward in this study), a competition between the gravity forces and stress concentration occurs. Results show that in majority of cases (6 out of 9) the fracture reaches to the middle layer and due to the lower stress in the middle layer, an extensive growth in the layer is observed. In the 3D model, the fracture can also grow laterally, thus the competition is between the lateral, upward and downward growth. Since the stress distribution curves, given in Fig. 2b, have a minimum at a certain depth below the interface, the fracture tend to grow towards that depth and in that depth the lateral growth will be favourable. The lateral growth of the elliptic fractures undermines the plane-strain assumption commonly used to model these fractures. Finally, the simulation results show that depending on the stress distribution, the shape of the fracture can be far from a simple plane-strain or radial shape.