Mixed-dimensional poromechanical models of fractured porous media

Abstract

We combine classical continuum mechanics with the recently developed calculus for mixed-dimensional problems to obtain governing equations for flow in, and deformation of, fractured materials. We present models in both the context of finite and infinitesimal strain, and discuss nonlinear (and non-differentiable) constitutive laws such as friction models and contact mechanics in the fracture. Using the theory of well-posedness for evolutionary equations with maximal monotone operators, we show well-posedness of the model in the case of infinitesimal strain and under certain assumptions on the model parameters.

1 Introduction

The general topic of fractured porous media is of great importance in applications from biomedicine, to industrial materials, to subsurface geophysics. Its successful mathematical treatment requires a combination of classical elasticity [1] with contact mechanics [2], and poromechanics [3] all of which must be extended to allow for complex geometric descriptions. Much progress has been made recently on the understanding of fluid flow in fractured porous media utilizing the conceptual framework of mixed-dimensional geometries [4, 5], which allows for lower-dimensional representations of fractures and their intersections.

Despite the importance and recent attention, the mathematical modeling and analysis of flow and deformation in fractured porous media is still far behind the needs of numerical analysts and practitioners. As a response to this, the current paper has two main aims: First, to provide the first consistent and frame-invariant mathematical model for fractured porous media on mixed-dimensional geometries. Second, to provide a well-posedness theory covering a broad class of problems of relevance to applications.

1.1 Introduction to modeling and analysis of fractured porous media

A realistic model of flow in fractured porous media necessarily requires a mathematical description of both the fluid flow as well as the mechanical response. This combination includes important nonlinearities stemming from the finite strain theory itself, combined with the contact-mechanical problem in the fracture and finally the nonlinear dependence of fluid flow on the fracture aperture. These nonlinearities appear in the context of a problem that essentially has a saddle-point structure due to the coupling of flow and deformation. To date, and to the best of the knowledge of the authors, this problem has not been successfully analyzed in full. In this contribution, we will exploit recent developments in the form of mixed-dimensional calculus, together with abstract results from the theory of nonlinear monotone evolutionary equations, to provide a frame-invariant and self-consistent theory that, together with well-posedness results for poromechanics of fractured media, extends well beyond existing analysis.

The present work needs to be seen in context of three independent developments over the last two decades. Firstly, in terms of modeling, it has been recognized since the work of Martin et al. [5], that fluid flow in fractured materials can be successfully modeled using a co-dimension one representation of the fracture. We will refer to such models, where the underlying geometry is composed of domains with different topological dimension, as mixed-dimensional. By now, mixed-dimensional models for fluid flow in fractured media are well established, both from the perspective of well-posedness [6], as well as their approximation properties relative to the underlying equidimensional problem [7, 8]. Secondly, the present authors have developed a general framework for considering mixed-dimensional models of this type, where using the language of exterior calculus and differential forms, basic concepts of mixed-dimensional functions and operators are established [4]. This leads to a mixed-dimensional functional analysis, which has been shown to inherit many of the tools associated with standard functional analysis. Thirdly, Picard and collaborators have developed the existence theory for evolutionary equations in the setting of maximally monotone operators [9, 10]. In particular, they have shown how poroelasticity can be analyzed in this framework [11] and that the setting is well suited to handle nonlinearities such as arise for contact problems [12]. The combination of these three developments is the foundation that allows us to consider the poromechanical contact problem which lies at the heart of poromechanics for fractured media. However, a key missing ingredient in the above is the representation of poromechanics as a mixed-dimensional model.

Earlier works have considered this problem using more standard approaches. Girault et al. have considered coupled poromechanics for fracture with a mixed-dimensional formulation for flow in the sense of a lower-dimensional flow representation within the fracture [13]; however, in their analysis they have disregarded the nonlinearities associated with changes in fracture aperture, both as it pertains to the flow problem, but also the contact mechanics. Furthermore, geometric complexity is ignored as only a single fracture is considered. A different perspective was taken by Yotov et al., who considered the problem in an equidimensional sense using Stokes’ equation for the flow in the fracture, but again considering only infinitesimal aperture changes such that contact was disregarded [14]. Bonaldi et al. show well-posedness for the case where nonlinearities arise due to multiphase flow, but consider only linearized mechanics [15]. This expands on similar results for the single-phase case in [16]. Finally, we mention also the work of Cusini et al., which address numerical method for this coupled problem [17]. While they consider geometric complexity, they limit their discussion to quasi-static, small-strain kinematics. This limitation, in particular, implies that only small slip-lengths are allowed in the resulting problem. None of the works discussed above considered finite strain modeling. Numerical and other modeling contributions have been summarized in two recent reviews [14, 18], where important contributions relevant for this paper include the work of Jha and Juanes [19], Garipov et al. [20], Norbeck et al. [21], Berge et al. [22] and Stefansson et al. [23].

With the above background in mind, we here summarize the main contributions of this paper:

• 1. A frame-invariant formulation of finite strain suitable for fractured media within the context of mixed-dimensional calculus, allowing for a large class of complex fracture networks, and its correspondence to classical finite strain theory.

• 2. Governing equations for finite strain poromechanics of fractured media expressed in terms of mixed-dimensional variables and operators for infinitesimal strain, while allowing for contact mechanics, frictional sliding and lubrication theory for flow in fracture.

• 3. Well-posedness theory for a linearized strain model, in the presence of contact mechanics and friction, under certain constraints on the constitutive laws.

1.2 Outline

The remainder of this paper is structured as follows. Section 2 introduces the fundamental definitions used in the formulation and analysis of mixed-dimensional models. We discuss the mixed-dimensional continuum assumption that is central in handling the different length scales inherent to these models. The admissible geometry is then introduced and we keep track of the connectivity between subdomains using directed acyclic graphs (DAGs). These DAGs allow us to create function spaces containing scalar and vector-valued functions that are relevant to modeling poromechanics in fractured media. All functions are defined on smooth reference domains and we use concepts from exterior calculus to appropriately map these to physical space.

Section 3 derives invariant strain measures for the mixed-dimensional setting. The definitions follow a “top-down" approach in which we a strain measure is formed as the linearization of a rotationally invariant finite strain. Additional attention is given to the volumetric strain as it forms the key term that couples the flow and mechanics equations in Biot poroelasticity models.

The mixed-dimensional poromechanics model is presented in Sect. 4. The model consists of the physical conservation principles of mass and momentum supplemented by appropriate constitutive laws. Two models are derived, based on the finite and linearized strain measures of Sect. 3, respectively. This section concludes with a discussion relating our model to classic models of (poro)elasticity.

Section 5 focuses on the well-posedness analysis of our model. We introduce a set of simplifying assumptions on the constitutive laws that ensures that the system can be analyzed as an evolutionary equation. Using the assumed monotonicity of our relations, we obtain well-posedness of our model in temporally weighted spaces. To close the section, example models are presented that contain conventional choices for the constitutive laws and satisfy our assumptions.

We conclude the paper in Sect. 6, highlighting the necessary aspects of our model that ensure physical relevance and well-posedness. The paper is supplemented by Appendix , which gives the background on evolutionary equations necessary for the well-posedness analysis.

2 Preliminaries: mixed-dimensional modeling and analysis

In this section, we make precise the problem setting, its geometry and the operators adapted to the mixed-dimensional problems. The first subsection is more general and introductory in nature providing a continuum mechanical perspective on mixed-dimensional modeling, while the remaining sections lay the mathematical foundation for the exposition that follows.

2.1 Problem setting and motivation

Classical continuum mechanics (which we will refer to as fixed-dimensional whenever needed to avoid confusion) is the modeling tool which allows for the derivation and statement of the classical field equations [1, 24]. In particular, it leads to the development of conservation laws and constitutive laws satisfying suitable notions of frame invariance. A key building block for continuum mechanics is the assumption that the notion of a continuum is a reasonable modeling choice. We choose to formulate this as follows (the precise statement of the continuum assumption is not essential, see, e.g., the thorough discussion in [24]):

Definition 2.1

(Fixed-dimensional continuum assumption) For a domain $(Y\subset {\mathbb {R}}^{n})$, there exists a scale of consideration , such that for any quantity of interest , and a (n)-dimensional ball centered on x and with radius , the integral below is well-defined, and the approximation is sufficiently accurate for the applications of interest:

(2.1)

In other words, our formulation of the fixed-dimensional continuum assumption states that a point evaluation of a quantity (say, porosity of a porous material) can be approximated by a (say, volume) integral of characteristic size , and that this approximation is accurate enough that the precise size (and indeed shape) of the integral is immaterial. As a classical example, one notes that for porosity, it is typically taken as a modeling assumption that a scale of consideration exists, the so-called “Representative Elementary Volume" (on the order of 10 to 100 times the mean grain size), wherein the porosity is well defined [25, 26]. At lower scales, the integral in (2.1) will be strongly affected by the precise number of grains in the integration volume. In the continuation, we will only be interested in continuum scales, and omit the overbar on the continuum quantity.

The classical continuum assumption is suitable for a wide range of applications, and underlies the vast majority of real-world industrial computations in applied engineering.

Fig. 1

Left: Illustration of a domain Y, with a high aspect-ratio inclusion $\Psi _i$, with maximal inscribed and minimum covering balls shown in red and blue, respectively. Right: Illustration of the same domain Y, where the high aspect-ratio inclusion is now modeled by a lower-dimensional manifold $\Omega _i$

Our interest herein is in problems for which the geometry in consideration contains high-aspect-ratio inclusions $\Psi _{i}$, indexed by i, that in some sense interferes with the continuum assumption. To be concrete (although the argument is more general), we consider thin fractures and their intersections. We characterize these high-aspect ratio inclusions by two length-scales, and , corresponding to the diameter of the maximal inscribed and minimum covering ball, respectively (as illustrated for a single manifold in the left part of Fig. 1). We are furthermore interested in the case where the small length scale violates the continuum assumption, i.e., where the following ordering holds:

(2.2)

For such problems, the standard fixed-dimensional continuum assumption cannot be applied, since the high-aspect feature (and variables within it) cannot be appropriately defined. Depending on the needs of the application at hand, we may nevertheless be inclined to consider as the appropriate modeling scale, in which case we have no other choice than to represent the thin inclusions as manifolds of lower dimension. We note that the intersection of such manifolds will have even lower dimension yet.

The above provides the motivation for considering a mixed-dimensional continuum assumption, wherein we are still interested in a domain $Y\subset {\mathbb {R}}^{n}$, but where we allow this domain to contain a set of manifolds of topological dimension $d<n$ (as illustrated for a single manifold in the right part of Fig. 1). We formulate this extension of Definition 2.1 as follows:

Definition 2.2

(Mixed-dimensional continuum assumption) Any inclusion $\Psi _{i}\subset Y$ which satisfies (2.2) can be well-represented by a $d_{i}$-dimensional manifold $\Omega _{i}$. Moreover on the scale of consideration , then for any quantity of interest , and a $d_{i}$-dimensional ball centered on ${\textrm{x}}\in \Omega _{i}$, the integrals below are well defined, and the approximation is sufficiently accurate for the applications of interest:

(2.3)

Here, we use the notation $\Psi _{i}^{\perp }$ to indicate the cross section of $\Psi _{i}$ orthogonal to $\Omega _{i}$, and we denote the measure of integration perpendicular and parallel to $\Omega _{i}$ by $\textrm{d}V_{\perp }$ and $\textrm{d}V_{\parallel }$, respectively.

We remark that our definitions of continuum assumptions suffer from the usual weaknesses [24], in that they are poorly adapted to quantities near boundaries and for variables which have macroscopic discontinuities. These issues can be resolved by appealing to more technical definitions. However, as we will not deal with the issue of upscaling in the following, but merely use the result that continuum variables can be assumed to be sufficiently accurate for applications of interest, we will not elaborate these details further.

In the following, we assume that we are always within the setting of Definition 2.2, and proceed to make the notion of mixed-dimensional continuum variables precise, and apply the framework of continuum mechanics to derive the governing equations for the poroelastic response in fractured porous media.

2.2 Geometry

To initialize our description of a mixed-dimensional problem, we first consider an admissible mixed-dimensional partitioning. The partitioning of Fig. 1 is, in a sense, too simple since it does not keep track of the boundaries between domains of different dimension. To achieve this, we herein introduce structured partitions of the domain along with the corresponding graph representations. These graphs first provide a canonical way to describe the connectivity between subdomains. Additionally, these definitions give the structure that allows us to define spaces of mixed-dimensional functions and the associated semi-discrete operators in Sects. 2.3 and 2.4. For a detailed exposition of these concepts in the scalar setting, we refer to [4].

We will only consider problems embedded in a $n$-dimensional Cartesian ambient domain, and we are primarily concerned with the case $n = 3$. Thus, let $Y\subset {\mathbb {R}}^{n}$ be given, and let it be decomposed into non-overlapping, oriented manifolds $\Omega _{i}$ of topological dimension $d_{i}$ such that $Y = \cup _{i\in I}^{}\Omega _{i}$ with $I$ the index set. In order to distinguish the domain and the partition, we will refer to the partition as $\Omega$ without a subscript.

We will not allow for arbitrary partitions, and therefore introduce a concept of admissible partitions. This requires some preliminaries. We first give each manifold $\Omega _{i}$ some additional hierarchical structure [4]. Each $\Omega _{i}$ is $C^{1}$-diffeomorphic to a smooth reference domain denoted by $X_{i}$ and we denote the corresponding mapping by $\phi _{0,i}:X_{i}\rightarrow \Omega _{i}$. We then endow each manifold with a directed acyclic graph, defined as follows.

Definition 2.3

A rooted directed acyclic graph (DAG) ${\mathfrak {S}}_{i}$ with $i \in I$, is conforming to $\Omega _{i}$ if for all nodes $j \in {\mathfrak {S}}_{i}$:

• There exists a root $s_{j}\in I$ such that $\phi _{0,j}(X_{j}) = \Omega _{s_{j}}$. Moreover, we assume $s_{i} = i$ for each root $i$ for convenience.

• For each descendant $l \in I_j$, where $I_{j}$ is the set containing the descendants of a node $j\in {\mathfrak {F}}$, a differentiable map $\phi _{j,l}:X_{l}\rightarrow {\overline{X}}_{j}$ exists with bounded derivative. We denote its range by $\partial _{l}X_{j}$. Compound maps telescope in the sense that

  $$\begin{aligned} \phi _{k,l} = \phi _{k,j}\circ \phi _{j,l} \end{aligned}$$

  for each ancestor $k\in {\mathfrak {S}}_{i}$ of $j$.

• The descendants uniquely cover the parent node in the sense that

  $$\begin{aligned} \bigcup _{j\in {\mathfrak {S}}_{i}}^{}\phi _{i,j}(X_{j}) ={\overline{X}}_{i} {\setminus } \phi _{0,i}^{-1}(\partial Y\cap \partial \Omega _{i}). \end{aligned}$$

  In other words, each point $x_{i}$ in reference domain $X_{i}$ and on its boundary is uniquely associated to a node $j\in {\mathfrak {S}}_{i}$ and a point $x_{j}\in X_{j}$ such that $x_{i} = \phi _{i,j}(x_{j})$. For $x_{i}$ on the boundary $\partial X_{i}$, we have $j$ a descendant of $i$ whereas for $x_{i}$ in the interior of $X_{i}$, we have $j = i$. All points that are mapped to the physical boundary $\partial Y$ by $\phi _{0, i}$ are exempt from this rule.

From this we see that each $\Omega _{i}$ is indeed a manifold, and furthermore, we have access to a partition of its boundary through the DAG ${\mathfrak {S}}_{i}$. This partitioning is illustrated in Fig. 2.

Fig. 2

Left: Illustration of reference domains and mappings to a physical domain with a single slit. The domains that map to $\Omega _{2}$ and their respective images under the mappings $\phi _{i,j}$ are highlighted in red. To comply with Definition 2.4, each $X_{j}$ with equal $s_{j}$ is considered equal. Right: The structure of the forest ${\mathfrak {F}}$ and its component DAGs. For each node $j$, the value of $s_{j}$ denotes index of the domain with which $\phi _{0,j}(X_{j})$ coincides in the physical domain. The $k$-forests ${\mathfrak {F}}^{k}$, depicted in blue, are introduced in Sect. 2.4

Based on the structure given in Definition 2.3, we can now provide a global structure to partitions $\Omega$ of $Y$ as follows [4].

Definition 2.4

A forest ${\mathfrak {F}} :=\bigcup _{i\in I}^{}{\mathfrak {S}}_{i}$ is conforming to $\Omega$ if the DAGs ${\mathfrak {S}}_{i}$ are conforming to $\Omega _{i}$ for all $i\in I$ in the sense of Definition 2.3, and if for any $j_{1},j_{2}$ such that $s_{j_{1}} = s_{j_{1}}$, it holds that $X_{j_{1}} = X_{j_{2}}$ and $\phi _{0,j_{1}} = \phi _{0,j_{2}}$.

Our main concern is partitions with conforming forests, and for clarity, we encode this in the following definition:

Definition 2.5

A partition $\Omega$ of $Y$ is admissible if a conforming forest ${\mathfrak {F}}$ exists. For any admissible partition, we denote the product space of reference domains as ${\mathfrak {X}}_{i}:=\prod _{j\in {\mathfrak {S}}_{i}}^{}X_{j}$ and ${\mathfrak {X}} :=\prod _{i\in I}^{}{\mathfrak {X}}_{i}$.

Definition 2.5 allows for a rather large generality of domains, including curved and self-intersecting domains, and multiple examples are provided in the cited reference [4]. In the present context, a relevant illustration for the case of a single fracture is provided below.

Example 2.1

In the case of a single fracture, as was discussed in Fig. 1, it creates a geometry as illustrated in Fig. 2. Note that in this example, only the domain $\Omega _{4}$ has a non-contractible reference domain $X_{4}$ associated with it, which comes in part from the fact that it has two boundaries (from “the top" and from “the bottom") neighboring the fracture $\Omega _{3}$.

Remark 2.1

This structure naturally allows for lower-dimensional domains to terminate at, or even exist entirely on, the global boundary $\partial Y$. In turn, boundary conditions can naturally be inherited on the fractures through the mappings $\phi _{0,j}$.

As a convention for indexing, we will use as above $i\in I$ for the roots of the DAGs, $j\in {\mathfrak {S}}_{i}$ for the components of DAG $i$, and finally we will write $j\in {\mathfrak {F}}$ for the components of the full forest. We moreover define the following index set

$$\begin{aligned} I^{d} :=\left\{ i\in I \mid d_{i} = d\right\} , \end{aligned}$$

(2.4)

Additionally, let $I_{j}^{d} :=\left\{ l\in I_{j} \mid d_{j} = d\right\} .$ We will moreover use $I^{d<n}$ to denote the set $\left\{ i\in I\mid d_{i}<n\right\}$.

For any domain $X_{i}$ with $\Omega _{i} = \phi _{0,i}(X_{i})$, we denote by ${\textbf{F}}_{i} :=D\phi _{0,i}$ the Fréchet derivative of the $C^{1}$ mapping $\phi _{0,i}$, defined as the linear operator such that for any vector $v\in {\mathbb {R}}^{d_{i}}$ and any point $x\in X_{i}$ then

$$\begin{aligned} {\textbf{F}}_{i}(x)v = \lim _{\epsilon \rightarrow 0}\frac{\phi _{0,i}(x+\epsilon v)-\phi _{0,i}(x)}{\epsilon }. \end{aligned}$$

(2.5)

Note that in terms of vector–matrix notation, which we conform to herein, we represent ${\textbf{F}}$ by a matrix whose rows correspond to gradients of the components of $\phi _{0, i}$.

We will need appropriate extensions of the mappings $\phi _{0, i}$, so that we can transform vectors in ${\mathbb {R}}^{n}$. For root nodes $i\in I$, i.e., the manifolds of dimension $d_i$, we define these in the following way:

Definition 2.6

Let $i\in I$. Let ${\hat{X}}_{i}$ be an open domain such that $\dim ({\hat{X}}_{i}) = n$ and ${\overline{X}}_{i}\times \left\{ 0\right\} ^{n-d_i} \subset {\hat{X}}_{i}$. Then, the extended mapping ${\hat{\phi }}_{0,i}:{\hat{X}}_{i}\rightarrow \hat{\Omega }_{i}\subset Y$ is defined such that

• ${\hat{\phi }}_{0,i} = \phi _{0,i}$ in ${\overline{X}}_{i}\times \left\{ 0\right\} ^{n-d_i}$.

• Orthogonality with respect to $X_{i}$ and $\Omega _{i}$ is preserved, i.e., the standard basis vector(s) ${\textbf{e}}_{d}$ for $d>d_{i}$ is/are mapped to $D{\hat{\phi }}_{0,i}{\textbf{e}}_{n}\perp T\Omega _{i}$ with $T\Omega _{i}$ the tangent bundle of $\Omega _{i}$.

• ${\hat{\phi }}_{0,i}$ has a fixed scaling with respect to the orthogonal complement of $X_{i}$, i.e., it holds that in ${\overline{X}}_{i}\times \left\{ 0\right\} ^{n-d_i}$.

We recall that the $d_{i}$ dimensional volume spanned by the derivative of a map $D\phi _{0,i}$ is given by

$$\begin{aligned} \textrm{vol}(D\phi _{0,i}) :=\sqrt{\det ((D\phi _{0,i})^{T}D\phi _{0,i})}. \end{aligned}$$

(2.6)

The definition of the extended mappings implies that for $i\in I^n$, then simply ${\hat{\phi }}_{0,i} = \phi _{0,i}$.

In Definition 2.6, the extension is given a length-scale designated by a parameter . This is an important point with respect to modeling, because this implies that we choose a conceptually arbitrary length-scale transversely to the fractures. Indeed, this is a necessary consequence of the mixed-dimensional continuum assumption: Since the transverse opening of a fracture is negligible (relative to the metric of the problem), if we are to nevertheless measure this opening, it must be measured in a different metric. This situation is analogous to multi-scale expansions encountered in homogenization, where the model depends in a non-negligible way on a (fine-scale) coordinate, which has negligible extent relative to the main (coarse-scale) coordinate of the problem (see, e.g., [27]).

For the leaves of the DAGs, corresponding to, e.g., boundaries of the solid matrix and the tips of the fractures, we have more freedom to choose an extended mapping, resulting in an arbitrary extended coordinate system around the boundaries of domain. We specify the definition of the extended mappings on the boundaries of domains as follows:

Definition 2.7

Let $i\in I$ and let $j\in I_{i}$ be a descendant. Let ${\hat{X}}_{j}\subset {\mathbb {R}}^{n}$ with ${\overline{X}}_{j}\times \left\{ 0\right\} ^{n-d_{j}}\subset {\hat{X}}_{j}$. Then, an extended mapping ${\hat{\phi }}_{i,j}:{\hat{X}}_{j}\rightarrow {\hat{X}}_i$ satisfies

• ${\hat{\phi }}_{i,j} = \phi _{i,j}$ in ${\overline{X}}_{j}\times \left\{ 0\right\} ^{n-d_{j}}$.

• Orthogonality with respect to $X_{j}$ and $X_i$ is preserved, i.e., the standard basis vector(s) ${\textbf{e}}_{d}$ for $d>d_{i}$ is/are mapped to $D{\hat{\phi }}_{i,j}{\textbf{e}}_{d}\perp TX_i$.

• ${\hat{\phi }}_{i,j}$ has a fixed scaling with respect to orthogonal complement of $X_{j}$, i.e., it holds that $\textrm{vol}(D{\hat{\phi }}_{i,j}) = \textrm{vol}(D\phi _{i,j})(x)$ in ${\overline{X}}_{j}\times \left\{ 0\right\} ^{n-d_{j}}$.

The extended mapping to the physical domain is given by composition, ${\hat{\phi }}_{0,j} = {\hat{\phi }}_{0,i} \circ {\hat{\phi }}_{i,j}$.

The above definition does not uniquely specify an extended mapping whenever $d_i-d_j \ge 2$; however, the precise choice of extensions has no impact on the following derivations, and we therefore omit a further specification.

For mechanics, we will be interested in deformation. Thus, we will allow for the domain $Y$, the partition $\Omega _{i}$, and the mappings $\phi _{0,i}$ to be time-dependent. However, we will not allow for structure of the partition to change, i.e., the DAGs ${\mathfrak {S}}_{i}$, the topological dimensions $d_{i}$ and the identification of root nodes $s_{i}$ are not variable and neither are the mappings $\phi _{i,j}$ for $j\in I_{i}$. Nor do we allow for any macroscopic opening of the fractures, that is to say, we a priori assume that the dynamics stay within the range of validity of the mixed-dimensional continuum assumption. Nevertheless, we allow for sliding of the fractures, as well as fracture opening on the scale of . In order to capture this, we extend the definition of coordinate mappings:

Definition 2.8

For $j, k\in {\mathfrak {F}}$ with $s_{j} = s_{k}$, we define the coordinate mapping

$$\begin{aligned} \phi _{j, k} :=\phi _{0,j}^{-1}\circ \pi _{j,k} \circ \phi _{0,k}, \end{aligned}$$

with $\pi _{j, k}(x_{k})$ identifying the point on $\Omega _{j}$ closest to $x_{k}$, i.e.,

$$\begin{aligned} \pi _{j, k}(x_{k}) :=\text {argmin}_{x_{j}\in \Omega _{j}}\left| x_{k}-x_{j}\right| . \end{aligned}$$

Note that if $j, k$ are members of the same conforming DAG, then $\pi _{j,k}$ is the identity operator. In turn, Definition 2.8 generalizes the definition of $\phi _{j,k}$ to nodes of different DAGs. Moreover, the mixed-dimensional continuum assumption assures that the expression $\left| x_{k}-x_{j}\right| = {\mathcal {O}}( \ell _\epsilon )$, and as such we infer that $\pi _{j, k}(x_{k})$ should be uniquely defined in our context. We will interpret opening of fractures that are sufficiently large to lead to non-uniqueness of $\pi _{j, k}(x_{k})$ as outside the range of validity of the modeling scales.

Remark 2.2

We recall that since the mappings $\phi _{0,i}$ are time-dependent, the coordinate map $\pi _{j, k}$ will also be. Thus mappings $\phi _{j,k}$, when $j$ and $k$ belong to different DAGs ${\mathfrak {S}}_{i_{1}}$ and ${\mathfrak {S}}_{i_{2}}$, will also in general be time-dependent.

Concluding this section, we provide an overview of the most important definitions related to mixed-dimensional geometries (Table 1).

Table 1 Summary of geometrical concepts

2.3 Mixed-dimensional function spaces

As a basis for mixed-dimensional modeling of the poromechanical system, we start by defining the appropriate variables on the geometry from Sect. 2.2. A critical aspect of our model is that the variables will be defined on reference domains, whereas the equations describing the physical model relate to the physical domain. It is therefore important to correctly transform functions between these domains.

We approach these transformations systematically by using concepts from the field of exterior calculus, in particular the equivalent representation of functions as differential forms (for a concise introduction, see, e.g., [28]). The pullback operator then provides the appropriate transformation mappings and, additionally, we obtain a canonical definition for trace operators. However, in order to make this presentation accessible to a broader audience, we only briefly exploit the calculus of differential forms, and then translate the definitions in terms of the representation by “standard" functions. We refer the interested reader to [4] for more details on the mixed-dimensional exterior calculus framework. Readers not familiar with exterior calculus are encouraged to skip ahead to Examples 2.3 and 2.4 and use these as a guide to the exposition.

We start with the following key building block, which is illustrated in the right part of Fig. 2:

Definition 2.9

The $k$-forest ${\mathfrak {F}}^{k}\mathfrak {\subseteq F}$ is defined for $0\le k\le n$ as the subgraph induced by the nodes

$$\begin{aligned} \bigcup _{ \begin{array}{c} i\in I \\ d_{i}\ge n-k \\ \end{array} }\left\{ j\in {\mathfrak {S}}_{i} \mid d_{i}-d_{j}\le n-k\right\} . \end{aligned}$$

As is apparent here, we keep track of the codimension between $X_{j}$ and $X_{i}$ and the difference between $n$ and $k.$ Later in this subsection, these determine on which boundary segments of given codimension we define function traces. For ease of reference, we summarize the integer values that play important roles in this section in Table 2.

Table 2 Summary of integers describing mixed-dimensional functions

We continue with our brief exposition of mixed-dimensional differential forms (a full account is given in [4]). The following five definitions suffice for the purposes of this work (confer, e.g., [28] for a concise introduction to fixed-dimensional differential forms).

Definition 2.10

For any $0\le k\le n$, let the mixed-dimensional reference domain ${\mathfrak {X}}^{k}$ be denoted ${\mathfrak {X}}^{k} :=\coprod _{j\in {\mathfrak {F}}^{k}}^{}X_{j}$ with $\coprod$ denoting the disjoint union. Each ${\mathfrak {X}}^{k}$ is thus a collection of subdomains that corresponds to a $k$-forest ${\mathfrak {F}}^{k}$, exemplified in Fig. 2.

We continue by defining the linear forms that are continuous on each $X_{j}\subseteq {\mathfrak {X}}^{k}$.

Definition 2.11

For any $0\le k\le n$, let the mixed-dimensional locally continuous $k$-forms with values in ${\mathbb {R}}^{p}$, for $p\in \{1,n\}$, be denoted ${\tilde{C}}{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$, and defined as a product space of alternating differential $k_{j}$-linear forms

$$\begin{aligned} {\tilde{C}}{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p}) :=\prod _{j\in {\mathfrak {F}}^{k}}^{}C^{1}\Lambda ^{k_{j}}(X_{j},{\mathbb {R}}^{p}), \end{aligned}$$

where the local order is given by $k_{j} = d_{i}-(n-k)$ for $j\in {\mathfrak {S}}_{i}$, and where $C^{1}\Lambda ^{k_{j}}(X_{j},{\mathbb {R}}^{p})$ are bounded $C^{1}$-continuous forms on $X_{j}$.

We refer to ${\mathfrak {a}}\in {\tilde{C}}{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$ as a mixed-dimensional differential form since it contains elements defined on manifolds of different dimensionalities. An element of such a function space will be denoted using the Gothic font. To extract a local form on, say, $X_{j}$ from ${\mathfrak {a}}$, we will use the notation $\iota _{j}{\mathfrak {a}}\in C^{1}\Lambda ^{k_{j}}(X_{j},{\mathbb {R}}^{p})$. This allows us to generalize the normal algebraic operators by insisting that they commute with $\iota$, thus if also ${\mathfrak {b}}\in {\tilde{C}}{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$, then $\iota _j({\mathfrak {a}}+{\mathfrak {b}})= \iota _j({\mathfrak {a}})+\iota _j({\mathfrak {b}})$, and similarly for subtraction, multiplication and division.

Example 2.2

We recall that a differential $k$-linear form $a\in C^{1}\Lambda ^{k}(X_{j}, {\mathbb {R}})$, takes as argument $k$ vectors from the tangent space $TX_{j} ={\mathbb {R}}^{d_{j}}$. Thus, $a(x):({\mathbb {R}}^{d_{j}})^{k_{i}}\mathbb {\rightarrow R}$. Moreover, the skew-symmetric (alternating) properties of $C^{1}\Lambda ^{k}(X_{j}, {\mathbb {R}})$ ensure that permutations of vectors alternate signs, e.g., for $v_{1},v_{2}\in TX_{j}$, and $k = 2$ then $a(x)(v_{1},v_{2}) = -a(x)(v_{2},v_{1})$.

Coordinate transformations are naturally handled in the context of exterior calculus through the pullback operator, defined next.

Definition 2.12

For any differentiable mapping $\phi : X\rightarrow \phi (X)$, the pullback $\phi ^{*}$ of the differential form $a\in C^{1}\Lambda ^{k}(\phi (X),{\mathbb {R}}^{p})$ is the unique operator such that $\phi ^{*}a$ satisfies for $v_{1},\ldots ,v_{k}\in TX$

$$\begin{aligned} (\phi ^{*}a)(v_{1},\ldots ,v_{k}) = a((D\phi ) v_{1},\ldots ,(D\phi ) v_{k}) \end{aligned}$$

Additionally, we have the trace operator that maps continuous differential forms to forms defined on the boundary.

Definition 2.13

For the boundary $\partial X$, the trace of the differential form $a\in C^{1}\Lambda ^{k}(X, {\mathbb {R}}^{p})$ is denoted $Tr_{\partial X} a\in C^{1}\Lambda ^{k}(\partial X,{\mathbb {R}}^{p})$ and is defined as the restriction of $a$ to the manifold $\partial X$.

We combine the locally continuous differential forms in order to establish a notion of globally continuous forms by exploiting the pullback and trace operators.

Definition 2.14

For any $0\le k\le n$, let the mixed-dimensional continuous $k$-forms be denoted $C{\mathfrak {L}}^{k} ({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$, and defined as

$$\begin{aligned} C{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p}) :=&\left\{ {\mathfrak {a}}\in {\tilde{C}}{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p}) \mid \phi _{i,j}^{*}\textrm{Tr}_{\partial _{j}X_{i}}\iota _{i}{\mathfrak {a}} = \varepsilon _{i,j}\iota _{j}{\mathfrak {a}},\right. \\&\left. \forall i,j\in {\mathfrak {F}}^{k} \text { with } i\in I \text { and } j\in I_{i}\right\} . \end{aligned}$$

Here, $\varepsilon$ is the orientation indicator that takes the value $\varepsilon _{i,j} = 1$ if $\partial _{j}X_{i}$ and $\phi _{i,j}(X_{j})$ have the same orientation, and $\varepsilon _{i,j} = -1$ otherwise.

Remark 2.3

The orientation $\varepsilon _{i,j}$ can directly be calculated by verifying whether the composition ${\hat{\phi }}_{j}^{-1}{\hat{\phi }}_{i}$ of extended coordinate maps preserves orientation in ${\mathbb {R}}^{n}$.

An important detail is that the space $C{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$ is not globally continuous, as the continuity is only imposed within each DAG ${\mathfrak {S}}_{i}$. Pre-empting later developments, we note that, e.g., deformations in $C{\mathfrak {L}}^{0}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ are therefore allowed to be discontinuous across fractures in physical space.

The above definitions of mixed-dimensional differential forms, as well as their pullback and trace operators, allow us now to consider representations as mixed-dimensional functions as used in the remainder of this paper. We first start by identifying the standard representation of differential forms in terms of classical functions from multivariate calculus. This discussion will exclusively consider the case of $n = 3$. Similar representations are used for lower dimensions (see, e.g., [29]).

Definition 2.15

At any given point $x\in X_{j}$, the space of differential forms $C^{1}\Lambda ^{k}(X_{j},{\mathbb {R}}^{p})$ has $p \left( {\begin{array}{c}d_{j}\\ k\end{array}}\right)$ degrees of freedom. The standard representation $\grave{a}=\mathbb {r}\acute{a}$ of a differential form $\acute{a}\in C^{1}\Lambda ^{k}(X_{j},{\mathbb {R}}^{p})$ is given for $p = 1$ with respect to the standard basis for ${\mathbb {R}}^{n}$ as follows (when we need to distinguish between the form and its representation, we denote the form by an accent aigu, and the representation by an accent grave):

• 1. For $k = 0$, the differential forms coincide with functions $\grave{a}\in C^{1}(X_{j})$, thus $\grave{a} =\acute{a}$.

• 2. For $k = 1$, the differential one-forms $\Lambda ^{k}(X_{j})$ are represented by vector functions $\grave{a}\in C^{1}(X_{j},TX_{j})$. This representation is the Riesz representation, which satisfies for vector fields $v_{1}\in TX_{j}$ that $\acute{a}(v_{1}) = v_{1} \cdot \grave{a}$.

• 3. For $k = n-1$, the differential forms $\Lambda ^{k}(X_{j})$ are represented by “flux" functions $\grave{a}(x)\in C^{1}(X_{j},TX_{j})$. This representation satisfies for vector fields $v_{1},v_{2}\in TX_{j}$ that $\acute{a}(v_{1},v_{2}) = \textrm{vol}(\grave{a},v_{1},v_{2})$, where $\textrm{vol}$ is the volume of the parallelepiped spanned by its arguments.

• 4. For $k=n$, the differential forms $a \in \Lambda ^{k} (X_j)$ coincide with “density” functions $\grave{a}(x) \in C^1 (X_j )$. This representation satisfies for vector fields $v_1,v_2,v_3 \in TX_j$ that $\acute{a}(v_1, v_2, v_3) = \grave{a} \textrm{vol} (v_1, v_2, v_3)$

Since both pullback and trace act differently depending on the order of the form, we use the vernacular “flux" to distinguish representations of $n-1$ forms from representations of $1$-forms, and similarly “density" to distinguish representations of $n$-forms from representations of $0$-forms.

Remark 2.4

From Definition 2.15, it is apparent that for $n\le 2$ (and importantly for our context, domains $X_{j}$ for $d_{j}\le 2$), the choice of representation is not unique, since, e.g., the possibility $1 = k = d_{j}-1$ exists. This is a classical observation, and is resolved in the current context by the following conventions: 1) For $p = n$, representations as (vector) functions, i.e., 1. and 2. in Definition 2.15, are preferred over 3. and 4. When needed, we emphasize this representation by the subscript $\mathbb {r}_{-}$. 2) For $p = 1$ representations as fluxes and densities, i.e., 3. and 4. in Definition 2.15, are preferred over 1. and 2. When needed, we emphasize this representation by the subscript $\mathbb {r}_{+}$.

Once a choice of representations has been established, we now have a one-to-one correspondence between mixed-dimensional differential forms and their function counterparts, and the inverse representation $\mathbb {r}^{-1}$ is thus well defined. The definitions of mixed-dimensional function spaces are now implied by the previous developments.

Definition 2.16

For any mixed-dimensional form $\acute{{\mathfrak {a}}}\in C{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$ we denote its standard mixed-dimensional representation as $\grave{{\mathfrak {a}}} =\mathbb {r}_{\pm }\acute{{\mathfrak {a}}}\in C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$, iff $\mathbb {r}_{\pm }\iota _{i}\acute{{\mathfrak {a}}} = \iota _{i}\grave{{\mathfrak {a}}}$ for all $i\in I$. For $p = n$ the spaces of continuous mixed-dimensional functions on ${\mathfrak {X}}^{k}$ that are relevant for this paper are given for $k = 0, 1$ by the choice:

$$\begin{aligned} C({\mathfrak {X}}^{k},{\mathbb {R}}^{n}) :=\mathbb {r}_{-}C{\mathfrak {L}}^{k}({\mathfrak {X}}^{k},{\mathbb {R}}^{n}). \end{aligned}$$

These are referred to as mixed-dimensional vector functions ($k = 0$) and matrix functions ($k = 1$).

For $p = 1$, the spaces of continuous mixed-dimensional functions on ${\mathfrak {X}}^{k}$ that are relevant for this paper are given for $k = n-1,n$ by the choice:

$$\begin{aligned} C({\mathfrak {X}}^{k}, {\mathbb {R}}) :=\mathbb {r}_{+}C{\mathfrak {L}}^{k}({\mathfrak {X}}^{k}, {\mathbb {R}}). \end{aligned}$$

We refer to the latter spaces as mixed dimensional fluxes ($k = n-1$) and densities ($k = n$).

Definition 2.17

The space of continuous mixed-dimensional functions with vanishing trace is given by

$$\begin{aligned} \mathring{C}({\mathfrak {X}}^{k},{\mathbb {R}}^{p}) :=\{\grave{{\mathfrak {a}}}\in C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})\mid \textrm{Tr}_{\partial _{Y}X_{j}}\iota _{j}\acute{{\mathfrak {a}}} = 0, \forall j\in {\mathfrak {F}}^{k}\}, \end{aligned}$$

with $\partial _{Y} X_{j} :=\phi _{0,j}^{-1}(\partial Y\cap \partial \Omega _{j}).$

On $C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$, we introduce a component-wise inner product as follows (we use angled brackets to denote inner products, and reserve parenthesis for tuples):

$$\begin{aligned} \left\langle \mathfrak {a, b}\right\rangle _{{\mathfrak {X}}^{k}}&:=\sum _{j\in {\mathfrak {F}}^{k}}^{}\left\langle \iota _{j}\mathfrak {a,}\iota _{j}{\mathfrak {b}}\right\rangle _{X_{j}}, \qquad \quad \forall \mathfrak {a, b} \in C({\mathfrak {X}}^{k},{\mathbb {R}}^{p}). \end{aligned}$$

(2.7)

This naturally induces an $L^{2}$-norm

$$\begin{aligned} \Vert {\mathfrak {a}}\Vert _{{\mathfrak {X}}^{k}}&:=\sqrt{\left\langle \mathfrak {a, b}\right\rangle _{{\mathfrak {X}}^{k}}}, \qquad \quad \forall {\mathfrak {a}} \in C({\mathfrak {X}}^{k},{\mathbb {R}}^{p}). \end{aligned}$$

(2.8)

The space of $L^{2}$ integrable functions can now be defined as the closure of the continuous functions with respect to this norm [4]:

Definition 2.18

For $0\le k\le n$ and $p \in \{1, n\}$, let the space of mixed-dimensional square integrable functions on ${\mathfrak {X}}^{k}$ be defined as

$$\begin{aligned} L^{2}({\mathfrak {X}}^{k},{\mathbb {R}}^{p}) :=\overline{C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})}. \end{aligned}$$

We emphasize that this definition is a direct result of the representation of differential in terms of conventional function spaces. As is clear from the definition (and motivated by Definition 2.15), even when the space is generated with $p = 1$, some function components may be vector-valued, as we will see in the examples below.

Example 2.3

The following two cases exemplify the spaces for $p = 1$ that are relevant to our model.

• Let $k = n$, then ${\mathfrak {F}}^{n}$ is given by the roots $I$. The number of degrees of freedom $p \left( {\begin{array}{c}d_{j}\\ k\end{array}}\right) = \left( {\begin{array}{c}n\\ n\end{array}}\right) = 1$ in this case, so we have $L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}}) = \prod _{i\in I}^{}L^{2}(X_{i}, {\mathbb {R}})$. This is the space in which we will define scalar density functions such as fluid pressures. In Fig. 2, these roots have indices $i$ with $1\le i\le 4$.

• Let $k = n-1$ and ${\mathfrak {a}}\in L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}})$. Then for each root $i\in I$ with $d_{i}\ge 1$, $\iota _{i}{\mathfrak {a}}$ is given by a vector function in $L^{2}(X_{i},{\mathbb {R}}^{d_{i}})$. Moreover, for each $j\in I_{i}^{d_{i}-1}$, we have that $\iota _{j}{\mathfrak {a}}\in L^{2}(X_{j}, {\mathbb {R}})$, i.e., a scalar distribution on each boundary that models an interface between manifolds of codimension one. This will form our space in which we define the fluid flux, both internal to each subdomain and across interfaces. In Fig. 2, these functions are then vectors on the root nodes $i = \{ 3,4\}$ and scalars on the nodes $j$ with $5\le j\le 8$. Note that for $\Omega _{3}$, there are three nodes $j\in \left\{ 3,7,8\right\}$ with $s_{j} = 3$. This allows us to separate the tangential flux inside the fracture ($j = 3$) and the normal flux entering from the two sides ($j \in \left\{ 7,8\right\}$).

Example 2.4

Similarly, we give examples of the two relevant spaces for $p = n$.

• Let $k = 0$. We have that ${\mathfrak {F}}^{0}$ consists of the roots $i\in I^{n}$ and all their descendants. Thus, for ${\mathfrak {a}}\in L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ and $i\in I^{n}$, we have that $\iota _{j}{\mathfrak {a}}\in L^{2}(X_{j},{\mathbb {R}}^{n})$ for all $j\in {\mathfrak {S}}_{i}$. We will use this space of ${\mathbb {R}}^n$-valued (vector) distributions to model the displacement of the solid. In Fig. 2, these functions are defined on the root with index $4$ and the descendant nodes $j$ with $7\le j\le 10$.

• Let $k = 1$ and $n = 3$. The 1-forest ${\mathfrak {F}}^{1}$ then contains all roots $i\in I^{2}\cup I^{3}$ and their descendants $j$ with $d_{j}\ge d_{i}-2$. For the roots $i\in I^{3}$, we have that $\iota _{i}{\mathfrak {a}}\in L^{2}(X_{i},{\mathbb {R}}^{3\times d_{i}})$ and for $j\in I_{i}^{1}\cup I_{i}^{2}$, it follows that $\iota _{j}{\mathfrak {a}}\in L^{2}(X_{j},{\mathbb {R}}^{3\times d_{j}})$. On the other hand, for $i\in I^{2}$, we obtain $\iota _{i}{\mathfrak {a}}\in L^{2}(X_{i},{\mathbb {R}}^{3})$ and $\iota _{j}{\mathfrak {a}}\in L^{2}(X_{j},{\mathbb {R}}^{3})$ for $j\in I_{i}^{0}\cup I_{i}^{1}$. This space of ${\mathbb {R}}^{3\times d_{j}}$-valued (tensor)-distributions will be used to model displacement gradients, allowing us to define stresses and strains in the bulk matrix and its boundaries, as well as across fractures. In Fig. 2, these functions are defined on the same domains as the flux functions discussed in Example 2.3. However, for $n = 3$, the flux and stresses will have different domains of definition.

The final spaces required are those that are defined on all reference domains, which becomes particularly useful when we consider volumetric strains in Sect. 3.3 and 3.4. These are defined for $p\in \{1,n\}$ as, e.g.:

$$\begin{aligned} L^{2}({\mathfrak {X}}, {\mathbb {R}}^p) :=\prod _{j\in {\mathfrak {F}}}^{}L^{2}(X_{j}, {\mathbb {R}}^p). \end{aligned}$$

(2.9)

Analogous definitions extend to $L^{\infty }({\mathfrak {X}}, {\mathbb {R}}^p)$ and $C^m({\mathfrak {X}}, {\mathbb {R}}^p)$.

In summary, the mixed-dimensional function spaces are defined using their equivalent representations as differential forms of order $k$. This gives us access to pullback and trace operators, as is illustrated in Fig. 3. In turn, function spaces in the physical domain are defined in the next subsection such that their pullback onto reference domains ${\mathfrak {X}}^{k}$ have certain regularity properties.

Fig. 3

The pullback $\grave{\phi }^{*} :=\mathbb {r}\acute{\phi }^*\mathbb {r}^{-1}$ of a function is defined using its representation as a differential form. Illustrated here is the case of $p = 1$ and $k = n-1$, i.e., the flux functions, for which the operator $\grave{\phi }^{*}$ is known as the Piola transform. The trace operator on functions is defined analogously as $\grave{\textrm{Tr}}:=\mathbb {r} \acute{\textrm{Tr}} \mathbb {r}^{-1}$. Due to the commutativity of this diagram, we omit the accents when denoting these operators

2.4 Differential operators

The standard differential operators on manifolds can be extended to the setting of mixed-dimensional geometries. However, in order to achieve the proper coupling between the domains, as required from physical relevance (i.e., the use of the differential operators in conservation laws), manifolds of adjacent dimensionality must be coupled via so-called jump operators.

This section presents the mixed-dimensional gradient and divergence operators, assuming continuous functions of sufficient regularity. For a rigorous exposition of all mixed-dimensional differential operators (including the curl), we refer again to [4], which follows a classical construction of Čech and de Rham cohomology.

Let us start by defining the jump operator by $\mathbb {d}:C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})\rightarrow C({\mathfrak {X}}^{k+1},{\mathbb {R}}^{p})$ that maps between subdomains of codimension one. We define this mapping by introducing a key set of indices.

Definition 2.19

For any root $i\in I$, let the index set $J_{i}$ be given by

$$\begin{aligned} J_{i} :=\left\{ j\in {\mathfrak {F}} \mid s_{j} = i \text { and } j\in I_{{\hat{\jmath }}} \text { for some } {\hat{\jmath }}\in I^{d_{i}+1}\right\} . \end{aligned}$$

In other words, this is the set of nodes that coincide in the physical domain with $i$ and have a root of dimension $d_{j} = d_{i}+1$. Then, the jump $\mathbb {d}_{\Phi }{\mathfrak {a}}\in C({\mathfrak {X}}^{k+1},{\mathbb {R}}^{p})$ of a continuous mixed-dimensional function ${\mathfrak {a}}\in C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$ is defined on a subdomain $X_{i}$ with $i\in I$ by the signed sum

$$\begin{aligned} \iota _{i}(\mathbb {d}_{\Phi }{\mathfrak {a}}) =(-1)^{n-k}\left( \sum _{l\in J_{i}}^{}\varepsilon _{i,l}\phi _{l,i}^{*}\iota _{l}{\mathfrak {a}}\right) ,\,\,\forall i\in {\mathfrak {F}}^{k+1}\cap I \end{aligned}$$

Where the pullback is used to map the function to the appropriate subdomain $X_{i}$. It is defined through the representation of functions as differential forms, cf. Fig. 3.

The jump $\mathbb {d}_{\Phi }$ is extended to descendants $j\in I_{i}$ by imposing that $\mathbb {d}_{\Phi }{\mathfrak {a}}$ is in $C({\mathfrak {X}}^{k+1},{\mathbb {R}}^{p})$. The remaining values $\iota _{j}\mathbb {d}_{\Phi }{\mathfrak {a}}$ with $j\in {\mathfrak {F}}^{k+1}{\setminus } I$ are thus determined by trace values, cf. Definition 2.14.

For an illustration of the domain and range of the jump operator, we refer to Fig. 2, recalling that ${\mathfrak {X}}^{k}$ is the set of subdomains corresponding to the $k$-forest ${\mathfrak {F}}^{k}$.

We note, as emphasized in Remark 2.2, that when $\Phi = \Phi (t)$ is time-dependent, then so are the operators $\pi _{l,i}$ and the mappings $\phi _{l,i}$, and hence also the definition of the operator $\mathbb {d}_{\Phi (t)}$. As a notational shorthand, we denote this time-dependent jump in reference space as $\mathbb {d}_{t} =\mathbb {d}_{\Phi (t)}$ when emphasis is needed, and otherwise simply write $\mathbb {d} = \mathbb {d}_{\Phi }$ also for the jump operator on reference space to declutter the presentation. This dependence of $\mathbb {d}_{\Phi }$ on the mapping $\Phi$ has the (intended) consequence that the jump term $\mathbb {d}_{t}$ remains local in physical space for two points in contact, even when the domains they belong to are sliding relative to each other.

Remark 2.5

By the mixed-dimensional continuum assumption, . This definition allows for geometries that are slightly more general than a conforming forest. Indeed, for geometries with a conforming forest $\pi _{l,i}(x) = x$. Furthermore, the closest point projection is Lipschitz continuous in the limit of infinitesimal smooth deformations. As a consequence, for a configuration $\Phi$ and a smooth perturbation $\Psi$ tangential to all boundaries, the derivatives considered from the “left" and “right" limits coincide, such that for all $i\in I^{2}$:

$$\begin{aligned} \lim _{\epsilon \rightarrow 0} \epsilon ^{-1}\iota _{i}(\mathbb {d}_{\Phi +\epsilon \Psi }\Phi ) = \lim _{\epsilon \rightarrow 0} \epsilon ^{-1}\iota _{i}(\mathbb {d}_{(\Phi +\epsilon \Psi )-\epsilon \Psi }(\Phi -\epsilon \Psi )) = \lim _{\epsilon \rightarrow 0} \epsilon ^{-1}\iota _{i}(\mathbb {d}_{\Phi }(\Phi -\epsilon \Psi )). \end{aligned}$$

This limit will be useful when considering the linearized theories later.

Next, we will construct the mixed-dimensional differential operators. Formally, these can be defined based on the differential forms and the exterior derivative [4], and then defining the differential operators on representations by requiring that commutation holds. However, as this introduces more formalisms than what is needed in the current exposition, we will here present the differential operators directly on the mixed-dimensional functions, with an understanding that mixed-dimensional exterior derivatives can similarly be defined on the mixed-dimensional forms.

We consider first the gradient, for which we set $k = 0$, and define the local gradient operator $\nabla :C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ as the standard gradient on each reference space, i.e., for ${\mathfrak {a}}\in C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ let $\nabla {\mathfrak {a}}\in L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ be such that

$$\begin{aligned} \iota _{j}(\nabla {\mathfrak {a}}) ={\left\{ \begin{array}{ll} \nabla (\iota _{j}{\mathfrak {a}}), &{} \forall j\in {\mathfrak {F}}^{1}\cap {\mathfrak {F}}^{0}, \\ 0, &{} \forall j\in {\mathfrak {F}}^{1}{\setminus }{\mathfrak {F}}^{0}. \\ \end{array}\right. } \end{aligned}$$

We emphasize that this operator takes the gradient not only of the components defined on the roots $j\in I^{n}$, but also on its boundaries $j\in I_{i}$.

Definition 2.20

The mixed-dimensional gradient on vector functions ${\mathfrak {D}}_{\Phi }:C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n}),$ is defined as

$$\begin{aligned} {\mathfrak {D}}_{\Phi }{\mathfrak {a}}&:=\nabla {\mathfrak {a}} + \mathbb {d}_{\Phi }{\mathfrak {a}},&\forall {\mathfrak {a}}&\in C({\mathfrak {X}}^{0},{\mathbb {R}}). \end{aligned}$$

(2.10)

The mappings $\Phi$ are usually implied from the context, we will then omit the subscript and write ${\mathfrak {D}}$.

Similarly, for the mixed-dimensional divergence, we first define $(\nabla \cdot ):C({\mathfrak {X}}^{n-1}, {\mathbb {R}})\rightarrow L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})$ such that

$$\begin{aligned} \iota _{j}(\nabla \cdot {\mathfrak {b}}) ={\left\{ \begin{array}{ll} \nabla \cdot (\iota _{j}{\mathfrak {b}}), &{} \forall j\in {\mathfrak {F}}^{n}\cap {\mathfrak {F}}^{n-1}, \\ 0, &{} \forall j\in {\mathfrak {F}}^{n}{\setminus }{\mathfrak {F}}^{n-1}. \\ \end{array}\right. } \end{aligned}$$

Definition 2.21

The mixed-dimensional divergence $({\mathfrak {D}}_{\Phi } \cdot ):C({\mathfrak {X}}^{0}, {\mathbb {R}})\rightarrow L^{2}({\mathfrak {X}}^{1}, {\mathbb {R}})$ is defined as

$$\begin{aligned} {\mathfrak {D}}_{\Phi }\cdot {\mathfrak {a}}&:=\nabla \cdot {\mathfrak {a}}+\mathbb {d}_{\Phi } {\mathfrak {a}},&\forall {\mathfrak {a}}&\in C({\mathfrak {X}}^{0}, {\mathbb {R}}). \end{aligned}$$

(2.11)

The mappings $\Phi$ are usually implied from the context, we will then omit the subscript and write $({\mathfrak {D}} \cdot )$.

It is important to note that these operators do not possess the same adjointness properties as the conventional gradient and divergence operators since they are defined on different $k$-forests. However, the adjoints of mixed-dimensional operators do play a vital role in our model and, to properly define these, we consider the differential operators as instances of densely defined unbounded linear operators (see, e.g., [30]) on $L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{p})$. Taking the adjoint (see Definition A.6) of the mixed-dimensional gradient and divergence then leads us to the co-gradient $({\mathbb {D}}_{\Phi } \cdot )$ and co-divergence $({\mathbb {D}}_{\Phi })$, respectively.

Definition 2.22

Let the mixed-dimensional co-gradient be denoted $({\mathbb {D}}_{\Phi } \cdot ): \textrm{dom}({\mathbb {D}}_{\Phi } \cdot )\subseteq L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ and the mixed-dimensional co-divergence be denoted ${\mathbb {D}}_{\Phi }:\textrm{dom}({\mathbb {D}}_{\Phi })\subseteq L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})\rightarrow L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}})$, defined such that for ${\mathfrak {b}} \in L^2 ({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ and ${\mathfrak {c}} \in L^2 ({\mathfrak {X}}^{n},{\mathbb {R}})$,

$$\begin{aligned} \left\langle {\mathbb {D}}_{\Phi }\cdot \mathfrak { b,a}\right\rangle _{{\mathfrak {X}}^{0}}&= -\left\langle {\mathfrak {D}}_{\Phi }\mathfrak {a,b}\right\rangle _{{\mathfrak {X}}^{1}},&\forall {\mathfrak {a}}&\in \mathring{C} ({\mathfrak {X}}^{0},{\mathbb {R}}^{n}), \end{aligned}$$

(2.12a)

$$\begin{aligned} \left\langle {\mathbb {D}}_{\Phi }\mathfrak {c,a}\right\rangle _{{\mathfrak {X}}^{n-1}}&= -\left\langle {\mathfrak {D}}_{\Phi }\cdot \mathfrak { a,c}\right\rangle _{{\mathfrak {X}}^{n}},&\forall {\mathfrak {a}}&\in \mathring{C}({\mathfrak {X}}^{n-1}, {\mathbb {R}}). \end{aligned}$$

(2.12b)

As with the gradient and divergence, we will in later sections mostly omit the subscript $\Phi$.

The differential operators ${\mathbb {D}}_{\Phi } \cdot$ and ${\mathbb {D}}_{\Phi }$ coincide with the conventional divergence and gradient on the roots, complemented by so-called half-jump operators on the boundaries $\partial _{j}X_{i}$ that relate $\iota _{i}{\mathfrak {a}}$ and $\iota _{s_{j}}{\mathfrak {a}}$. We refer the interested reader to [4] for explicit representations. On the other hand, Definition 2.18 defines the mixed-dimensional gradient and divergence on the more regular spaces $C({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$. We do not wish to require such regularity in the weak formulation of our model and we therefore expand the definition.

Definition 2.23

Let the mixed-dimensional gradient and divergence with boundary conditions be given by

$$\begin{aligned} \mathfrak {\mathring{D}}_{\Phi }&:\textrm{dom}(\mathfrak {\mathring{D}}_{\Phi })\subseteq L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n}),&\mathfrak {\mathring{D}}_{\Phi }&:=(-{\mathbb {D}}_{\Phi } \cdot )', \end{aligned}$$

(2.13a)

$$\begin{aligned} (\mathfrak {\mathring{D}}_{\Phi } \cdot )&:\textrm{dom}(\mathfrak {\mathring{D}}_{\Phi } \cdot )\subseteq L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}})\rightarrow L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}}),&(\mathfrak {\mathring{D}}_{\Phi } \cdot )&:=(-{\mathbb {D}}_{\Phi })'. \end{aligned}$$

(2.13b)

Remark 2.6

The circular accent on these operators indicates that the functions in the respective domains have vanishing trace on $\partial _{Y}{\mathfrak {X}}^{k}$. This is a direct consequence of using test functions ${\mathfrak {a}}\in \mathring{C}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})$ in (2.12). Moreover, the role of boundary conditions could have been reversed (and indeed generalized), but we will retain the choice implied above for simplicity of exposition.

To conclude this section, we emphasize that all differential operators (and the co-differentials) introduced herein are densely defined, unbounded linear operators mapping as $L^{2}({\mathfrak {X}}^{k},{\mathbb {R}}^{p})\rightarrow L^{2}({\mathfrak {X}}^{k+1},{\mathbb {R}}^{p})$, cf. Appendix . In particular, the density of $C({\mathfrak {X}}^{k}, {\mathbb {R}})$ in $L^{2}({\mathfrak {X}}^{k}, {\mathbb {R}})$ was shown in Theorem 3.1 of [4]. A consequence of this statement is that the spaces $\textrm{dom}({\mathbb {D}}_{\Phi } \cdot )$, $\textrm{dom}({\mathbb {D}}_{\Phi })$, $\textrm{dom}(\mathring{{\mathfrak {D}}}_{\Phi })$ and $\textrm{dom}(\mathring{{\mathfrak {D}}}_{\Phi } \cdot )$ are all Hilbert spaces with respect to their respective graph norms.

3 Mixed-dimensional strain measures

Scalar elliptic mixed-dimensional equations are well understood [4], and the case of $p = 1$ and $k\in \left\{ n-1, n\right\}$ leads to the standard equations used for mixed-dimensional models of flow in fractured porous media [4, 6, 31]. The main outstanding challenge in constitutive modeling is thus the correct treatment of the mechanical deformation in the mixed-dimensional setting. This is the topic of this section.

Our approach in this development is to follow the “top-down" modeling associated with classical continuum mechanics, in the tradition of, e.g., [1, 3, 24, 32], adapted to the mixed-dimensional geometry and spaces presented in Sect. 2. Thus, we obtain a mixed-dimensional finite strain theory directly for the geometric representation ${\mathfrak {F}}$. The converse approach, which we will not pursue in this work, would be to take the standard theory of mechanics as applied to the domain $Y$ with its high-aspect inclusions $\Psi _{i}$, and derive a finite strain theory for ${\mathfrak {F}}$ through an upscaling based on the limit process of . We will discuss the relationship between the results obtained in this work and classical theory as posed on $Y$ in Sect. 4.5.

3.1 Recollection of fixed-dimensional finite strain theory

To provide context for the mixed-dimensional strain measure introduced later, we briefly recall the standard setting of finite strain theory. We recall from (2.5) in Sect. 2.2 that for domains $X$ and $\Omega = \phi (X)$, we denote by ${\textbf{F}} = D\phi$ the derivative of the $C^{1}$ mapping $\phi$. Then,

Definition 3.1

The right Cauchy–Green deformation tensor ${\textbf{C}}:T\Omega \rightarrow T\Omega$ is defined for a configuration $\phi$ as

$$\begin{aligned} {\textbf{C}}(\phi ) :={\textbf{F}}^{T}{\textbf{F}}. \end{aligned}$$

(3.1)

Since we are only concerned with problems embedded in ${\mathbb {R}}^{n}$ with Cartesian coordinates, we will in the following use the same notation for all associated tensors. However, we will not have need for the full generality of tensor calculus as all variables are defined on subsets of ${\mathbb {R}}^{d_{j}}$. We will therefore not distinguish notationally between “raising and lowering indexes".

In our geometric setting, the reference domain $X$ is without physical meaning, and we will be concerned with a time-dependent physical configurations, represented by $\phi (t)$ and where the initial state is denoted ${\underline{\phi }} :=\phi (t = 0)$. These are naturally compared on the reference domain $X$, since the deformation tensor is rotationally invariant here. Thus, we have

Definition 3.2

Green–Lagrange strain tensor with respect to the configurations $\phi$ and ${\underline{\phi }}$ is defined by the 2-tensor

$$\begin{aligned} {\textbf{E}}(t) :=\frac{1}{2}({\textbf{C}}(\phi (t))\mathbf {-C}({\underline{\phi }})). \end{aligned}$$

(3.2)

The normalization factor $\frac{1}{2}$ is somehow arbitrary, but is typically included to ensure that the linearized strain becomes dual to the divergence operator on symmetric tensor functions. If furthermore the deformation is infinitesimal from the baseline configuration ${\underline{\phi }}$, i.e., that $\phi (t) ={\underline{\phi }}+u(t)$, and ${\textbf{F}} ={{{{\underline{\varvec{F}}}}}}+Du$, with $\left| Du\right| \ll 1$, in the case of $d = n$ leads to

$$\begin{aligned} {\textbf{E}}(t)&=\frac{1}{2}(({{{{\underline{\varvec{F}}}}}} + Du)^{T}({{{{\underline{\varvec{F}}}}}} + Du)-{{{{\underline{\varvec{F}}}}}}^{T}{{{{\underline{\varvec{F}}}}}})\\&=\frac{1}{2}({{{{\underline{\varvec{F}}}}}}^{T}Du+(Du)^{T}{{{{\underline{\varvec{F}}}}}})+(Du)^{T}Du. \end{aligned}$$

The linearized strain tensor is obtained by retaining the first-order terms in $\left| Du\right|$, as summarized below.

Definition 3.3

The linearized strain tensor with respect to the deformation $u(t) = \phi (t)-{\underline{\phi }}$ is defined by the 2-tensor

$$\begin{aligned} {\textbf{e}}(t) :=\frac{1}{2}({{{{\underline{\varvec{F}}}}}}^{T}Du(t)+(Du(t))^{T}{{{{\underline{\varvec{F}}}}}}). \end{aligned}$$

(3.3)

When expressed as a linear operator on $u(t)$, we refer to this operator as the symmetric gradient, and write

$$\begin{aligned} D_{s}u(t):={\textbf{e}}(t). \end{aligned}$$

(3.4)

We remark that this definition simplifies whenever the reference configuration $X$ is equal to the baseline physical configuration $\Omega$, since in this case ${\underline{\phi }}(x) = x$, and ${{{{\underline{\varvec{F}}}}}} = {{{{\underline{\varvec{F}}}}}}^{T} ={\textbf{I}}$. However, due to the nature of the mixed-dimensional geometries of interest herein, this will in general not be the case in our context. For example, fractures are not restricted to be located on the $xy$-plane, but the corresponding reference domains are.

3.2 Mixed-dimensional finite strain

We follow the same approach to derive a mixed-dimensional finite strain theory. To proceed, we first make precise the needed extensions of fixed-dimensional calculus to the mixed-dimensional setting. In particular, we have already defined the mixed-dimensional differential operators in Sect. 2.4. While these are in principle sufficient to obtain a mixed-dimensional strain, a richer strain notion can be obtained by also considering the derivative of the mixed-dimensional extended coordinate mappings, defined in Definitions 2.6 and 2.7. In the same manner as above, we therefore introduce

Definition 3.4

The derivative of the mixed-dimensional extended deformation ${\hat{\Phi }}$ is denoted $\varvec{{\mathfrak {F}}} :=D{\hat{\Phi }}$, and satisfies $\iota _{j}\varvec{{\mathfrak {F}}} = D{\hat{\phi }}_{j}$ for all $j\in {\mathfrak {F}}$.

Note that the forest ${\mathfrak {F}}$ and the derivative $\varvec{{\mathfrak {F}}}$ should not be confused.

It is an important point that the mixed-dimensional setting now deviates from the fixed-dimensional case, in that ${\mathfrak {D}}\Phi \ne D\Phi \ne {\mathbb {D}}\Phi$. These notions of a derivative of the configuration $\Phi$ (all in a sense “gradients") have important distinctions. The derivative of the deformation $\varvec{{\mathfrak {F}}}$, defined on ${\mathfrak {X}}$, contains information of the deformation of each $\Omega _{i}$, but has no information about the relative placements of domains. On the other hand, the mixed-dimensional gradient ${\mathfrak {D}}\Phi$, defined on ${\mathfrak {X}}^{1}$, contains information on relative placements (due to the jump operator $\mathbb {d}$), but only contains information regarding the deformation of the top-dimensional domains, i.e., those domains $\Omega _{i}$ where $d_{i} = n$. It is therefore clear that ${\mathfrak {D}}\Phi$ contains the desired physical information (since the lower-dimensional domains are fractures–voids–and their precise deformation is immaterial). Conversely, the deformation $\varvec{{\mathfrak {F}}}$ is required for coordinate transformations, and can be thought of as a “fabric" onto which to project vectors.

The above discussion suggests:

Definition 3.5

The mixed-dimensional right Cauchy–Green deformation tensor is defined for a configuration $\Phi$ as

$$\begin{aligned} \varvec{{\mathfrak {C}}} :=(\Pi ^{1}\varvec{{\mathfrak {F}}}^{T}){\mathfrak {D}}(\Pi ^{0}\Phi ), \end{aligned}$$

(3.5)

where $\Pi ^{k}$ is the restriction from ${\mathfrak {X}}$ to ${\mathfrak {X}}^{k}$.

Remark 3.1

The deformation $\varvec{{\mathfrak {F}}}$ based on the extended mappings ${\hat{\phi }}$ allows for transforming vectors (and thus forces) in ${\mathbb {R}}^{n}$ appropriately. Alternative suggestions for a “symmetric" deformation tensor, such as expressions of the type $\varvec{{\mathfrak {F}}}^{T}\varvec{{\mathfrak {F}}}$ or $({\mathfrak {D}}(\Pi ^{0}\Phi ))^{T}{\mathfrak {D}}(\Pi ^{0}\Phi )$, can be seen to be unsuitable, as the former contains no information of relative placements of domains, while the latter only retains the magnitude of displacements across a fracture, without orientation information.

Proceeding as in Sect. 3.1, we will use the mixed-dimensional right Cauchy–Green deformation tensor as the basis for defining a strain measure on the reference domain. We therefore consider a time-dependent mapping $\Phi = \Phi (t)$ from which we obtain a time-dependent deformation tensor $\varvec{{\mathfrak {C}}}(t)$. By again identifying time $t = 0$ as the reference time with ${\underline{\Phi }} :=\Phi (t = 0)$, then

Definition 3.6

The mixed-dimensional Green–Lagrange strain tensor with respect to the configurations $\Phi$ and ${\underline{\Phi }}$ is defined by

$$\begin{aligned} {\mathfrak {E}}(t) :=\varrho (\varvec{{\mathfrak {C}}}(t)-\underline{\varvec{{\mathfrak {C}}}}), \end{aligned}$$

(3.6)

where the mixed-dimensional gradients ${\mathfrak {D}}_{t}$ are evaluated based on the configuration at time $t$, such that in particular, $\underline{\varvec{{\mathfrak {C}}}} = \underline{\varvec{{\mathfrak {F}}}}^T{\mathfrak {D}}_{\Phi (t)}\underline{\Phi }$. Moreover, the normalization factor $\varrho$ is assigned the value $\iota _{j}\varrho =\frac{1}{2}$ for $j\in {\mathfrak {S}}_{i}$ and $i\in I^{n}$ and the value $\iota _{j}\varrho = 1$ otherwise. The justification for this choice will become apparent in Lemma 5.1.

Example 3.1

We consider the interpretation of the mixed-dimensional Green–Lagrange strain tensor on domains of various dimensionality:

• 1. For top-dimensional domains, $i\in I^{n}$, then as in the fixed-dimensional case,

  $$\begin{aligned} \iota _{i}{\mathfrak {E}}(t) ={\textbf{E}}_{i}(t). \end{aligned}$$

  (3.7)

• 2. On the boundaries of the top-dimensional domains $i\in {\mathfrak {S}}_{j}$, where $j\in I^{n}$, the deformation tensor is given by $\iota _{i}\varvec{{\mathfrak {C}}} = (\hat{{\textbf{F}}}_{i})^{T}{\textbf{F}}_{i}$ with $\hat{{\textbf{F}}}_{i} :=D{\hat{\phi }}_{i}$. It is thus represented by a ${\mathbb {R}}^{n}\times {\mathbb {R}}^{d_{i}}$ matrix. Then, the strain takes the form

  $$\begin{aligned} \iota _{i}{\mathfrak {E}}(t) =\frac{1}{2}(\hat{{\textbf{F}}}_{i}^{T}(t){\textbf{F}}_{i}(t)-\hat{{{{{\underline{\varvec{F}}}}}}}_{i}^{T}{{{{\underline{\varvec{F}}}}}}_{i}) =\frac{1}{2}\left( {\begin{array}{c}{\textbf{F}}_{i}^{T}(t){\textbf{F}}_{i}(t)-{{{{\underline{\varvec{F}}}}}}_{i}^{T}{{{{\underline{\varvec{F}}}}}}_{i}\\ {\textbf{0}}\end{array}}\right) . \end{aligned}$$

  (3.8)

  Note that due to Definition 2.7, the “extended" components of $\hat{{\textbf{F}}}_{i}$ are orthogonal to ${\textbf{F}}_{i}$ (whose columns are vectors in the tangent space $T\Omega _{i}$), thus the last $n-d_{i}$ rows of $\iota _{i}{\mathfrak {E}}(t)$ are identically zero, justifying the claim that the precise choice of extensions in Definition 2.7 is immaterial for the developments.

• 3. For domains $i\in I^{n-1}$ (the fractures), the mixed-dimensional gradient of the deformation is simply $\iota _{i}{\mathfrak {D}}\Phi = \iota _{i}\mathbb {d}_{t}\Phi$, i.e., the jump in $\phi _{j}$ between the two $n$-dimensional neighbors to $\Omega _{i}$. Thus,

  $$\begin{aligned} \iota _{i}{\mathfrak {E}}(t) =\hat{{{\textbf {F}}}}_{i}^{T}(t)(\iota _{i}\mathbb {d}_{t}\Phi (t))-\hat{{{{{\underline{\varvec{F}}}}}}}_{i}^{T}(\iota _{i}\mathbb {d}_{t}\underline{\Phi }). \end{aligned}$$

  As above, it is natural to decompose it into its parallel and normal components, denoted by subscripts $\Vert$ and $\perp$, respectively, which takes the form

  

  (3.9)

  Here, we denote the $n$^th row of $\hat{{\textbf{F}}}_{i}^{T}(t)$ by $\hat{{\textbf{F}}}_{i,n}^{T}(t)$, which we note equals , where ${\textbf{n}}_{i}$ is the normal vector orthogonal to $\Omega _{i}$ preserving the orientation of ${\hat{\phi }}_{i}$. Therefore, the expression for the strain can be simplified to

  

  (3.10)

  Here, we have decomposed the displacement jump into its orthogonal and parallel components,

  $$\begin{aligned} (\iota _{i}\mathbb {d}_{t}{\underline{\Phi }})_{\perp }&:={{{{\underline{\varvec{n}}}}}}_{i}^{T}(\iota _{i}\mathbb {d}_{t}\underline{\Phi }){} & {} \text {and}&(\iota _{i}\mathbb {d}\Phi _{0})_{\parallel }&:=\iota _{i}\mathbb {d}\Phi _{0}-{{{{\underline{\varvec{n}}}}}}_{i}(\iota _{i}\mathbb {d}\Phi _{0})_{\perp }. \end{aligned}$$

  Moreover, by the definition of the jump operator, the jump in the direction parallel to the fracture is identically zero, $(\iota _{i}\mathbb {d}_{t}\Phi (t))_{\parallel } = 0$, and this term can be omitted from (3.10). We furthermore note that by the mixed-dimensional continuum assumption the jump in the direction perpendicular to the fracture is of order , thus . In contrast, sliding is measured as $(\iota _{i}\mathbb {d}_{t}{\underline{\Phi }})_{\parallel }$, which measures the slip of the two fracture surfaces from the initial state until the current configuration. We thus arrive at the final expression for the strain in fractures,

  

  (3.11)

• 4. Since ${\mathfrak {D}}\Phi$ is void on domains $\Omega _{j}$ with $j\in I^{d<n-1}$, so is $\iota _{i}{\mathfrak {E}}(t)$ for all $i\in {\mathfrak {S}}_{j}$.

Again, we emphasize that the measure of opening of a fracture has arbitrary scale, depending on the choice of . This implies that ${\mathfrak {E}}(t)$ is a multi-scale strain measure, which we will return to in Sect. 4 (Example 4.2). We close this section by verifying that the mixed-dimensional finite strain ${\mathfrak {E}}(t)$ is rotationally and translationally invariant.

Lemma 3.1

Let $\Phi (t)$ be a rigid body motion relative to ${\underline{\Phi }}$. Then, ${\mathfrak {E}}(t) = 0$.

Proof

A rigid body motion can be described by a rotation matrix $R(t)$ and a vector $V(t)$, both independent of space and the rotation satisfying $R^{-1}(t) = R^{T}(t)$. Then, $\Phi (t) = R(t){\underline{\Phi }}+V(t)$, i.e., for all $i\in {\mathfrak {F}}$ the local mapping is given by

$$\begin{aligned} \iota _{i}\Phi (t) = R(t)\phi _{0, i}(0)+V(t). \end{aligned}$$

Then, since differentiation is a linear operator with constants in its null-space, we have both ${\textbf{F}}_{i}(t) = R(t){{{{\underline{\varvec{F}}}}}}_{i}$ and $\hat{{\textbf{F}}}_{i}(t) = R(t)\underline{\hat{{\textbf{F}}}}_{i}$, while by the same argument the jump operator satisfies $\iota _{i}\mathbb {d}_{t}\Phi (t) = \iota _{i}\mathbb {d}_{t}(R(t){\underline{\Phi }}) = \iota _{i}R(t)\mathbb {d}_{t}{\underline{\Phi }}$.

Now a direct substitution gives

$$\begin{aligned} \hat{{\textbf{F}}}_{i}^{T}(t){\textbf{F}}_{i}(t) =(R(t)\underline{\hat{{\textbf{F}}}}_{i})^{T}R(t){{{{\underline{\varvec{F}}}}}}_{i} =\underline{\hat{{\textbf{F}}}}_{i}R^{T}(t)R(t){{{{\underline{\varvec{F}}}}}}_{i} =\underline{\hat{{\textbf{F}}}}_{i}{{{{\underline{\varvec{F}}}}}}_{i} \end{aligned}$$

and

Comparison with the local expressions for $\iota _{i}{\mathfrak {E}}(t)$ provided in Example 3.1 verifies the lemma. $\square$

3.3 Mixed-dimensional linearized strain

When considering a constitutive theory, our primary variable will be the displacement of the top-dimensional domains ${\mathfrak {u}}\in L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$. To be able to calculate the linearized strain in the remaining, lower-dimensional subdomains, we require an extension operator onto the domain of $\Phi$, which we define as:

Definition 3.7

A bounded linear operator $\Xi : L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}},{\mathbb {R}}^{n})$ is an ${\mathfrak {X}}^{0}$-extension operator if it is a right inverse of the restriction $\Pi ^{0}$, i.e.,

$$\begin{aligned} \Pi ^{0}\Xi {\mathfrak {u}}&:={\mathfrak {u}},&\forall {\mathfrak {u}}&\in L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n}). \end{aligned}$$

(3.12)

Remark 3.2

We note that the most natural choice for $\Xi$ is an averaging operator, such that for a fracture $i\in I^{2}$, its displacement $\iota _{i}(\Xi {\mathfrak {u}})$ is defined as the average displacement of the rock on the two sides. Such operators allow us to consider the representation as the mean of neighboring (rock) positions in $I^{n}$, such that for any $i\in I^{d<n}$

$$\begin{aligned} \iota _{i}(\Xi {\mathfrak {u}}) =\frac{1}{\left| {\hat{I}}_{i}^{n}\right| }\sum _{j\in {\hat{I}}_{i}^{n}}^{}\phi _{j,i}^{*}(\iota _{j}{\mathfrak {u}}). \end{aligned}$$

We study the role of extension operators in more detail in Sect. 3.4.

We obtain a linearized strain by considering deformations $\Phi (t)$ such that $\Phi (t) = \Xi {\mathfrak {u}}(t)+\underline{\Phi }$, and where the mixed-dimensional gradients in $\mathfrak {Du}$ are small in the sense that for all $i\in {\mathfrak {F}}$, and for all $x\in X_{i}$, it holds that

$$\begin{aligned} \Vert (\iota _{i}{\mathfrak {D}}{\mathfrak {u}})(x)\Vert \ll \Vert {{{{\underline{\varvec{F}}}}}}_{i}(x)\Vert . \end{aligned}$$

(3.13)

Using this we define the linearized strain as ${\mathfrak {e}}(t)$, obtained by omitting “small" terms. More precisely, we define the linearized strain as the Fréchet derivative of the finite strain in the following sense:

Definition 3.8

For some initial mapping ${\underline{\Phi }}$, and some deformation ${\mathfrak {u}}\in C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$, let the mixed-dimensional linearized strain be defined as

$$\begin{aligned} {\mathfrak {e}}({\mathfrak {u}}) :=D{\mathfrak {E}}(\underline{\Phi })(\Xi {\mathfrak {u}}), \end{aligned}$$

(3.14)

where $D{\mathfrak {E}}({\underline{\Phi }})(\Psi )$ is the Fréchet derivative of ${\mathfrak {E}}$ at ${\underline{\Phi }}$ acting on the perturbation $\Psi$. When expressed as a linear operator on ${\mathfrak {u}}(t)$, we refer to this operator as the symmetric gradient. Thus, ${\mathfrak {D}}_{s}: C^{1}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ is defined as

$$\begin{aligned} {\mathfrak {D}}_{s}{\mathfrak {u}} :={\mathfrak {e}}({\mathfrak {u}}). \end{aligned}$$

(3.15)

Example 3.2

Continuing from Example 3.1, we consider the interpretation of the mixed-dimensional linearized strain tensor on domains of various dimensionality, keeping in mind that ${\underline{\Phi }}$ is the unstrained state, i.e., ${\mathfrak {E}}({\underline{\Phi }}) = 0$.

• 1. For top-dimensional domains, $i\in I^{n}$, then as in the fixed-dimensional case, the linearized strain tensor takes the form

  $$\begin{aligned} \iota _{i}{\mathfrak {e}}({\mathfrak {u}}) =\frac{1}{2}({{{{\underline{\varvec{F}}}}}}_{i}^{T}Du_{i}(t)+(Du_{i}(t))^{T}{{{{\underline{\varvec{F}}}}}}_{i}). \end{aligned}$$

  (3.16)

• 2. On the boundaries of the top-dimensional domains $i\in {\mathfrak {S}}_{j}$, where $j\in I^{n}$, the linearized strain tensor is represented by a ${\mathbb {R}}^{n}\times {\mathbb {R}}^{d_{i}}$ matrix. It has the explicit form

  $$\begin{aligned} \iota _{i}{\mathfrak {e}}(t) =\frac{1}{2}(\hat{{{{{\underline{\varvec{F}}}}}}}_{i}^{T}Du_{i}(t)+(D{\hat{u}}_{i}(t))^{T}{{{{\underline{\varvec{F}}}}}}_{i}) =\frac{1}{2}\begin{pmatrix} {{{{\underline{\varvec{F}}}}}}_{i}^{T}Du_{i}(t)+(Du_{i}(t))^{T}{{{{\underline{\varvec{F}}}}}}_{i} \\ {\textbf {0}} \end{pmatrix}. \end{aligned}$$

  (3.17)

  As in the case of the finite strain, the last $n-d_{i}$ rows of $\iota _{i}{\mathfrak {e}}(t)$ are identically zero.

• 3. For domains $i\in I^{n-1}$ (the fractures), we calculate the Fréchet derivative as

  

  (3.18)

  In the final line, we have used the continuity of the jump operator, as elaborated in Remark 2.5. Substituting in the extended deformation $\Xi {\mathfrak {u}}$, we now obtain

  

  (3.19)

  Here, the extension operators vanish since they are identity operators on the top-level domains (by Definition 3.7).

• 4. As in the finite deformation case, the linearized strain is void for $i\in I^{d<n-1}$.

Remark 3.3

Example 3.2 illustrates that the extension $\Xi$ plays no role in the final expressions for the linearized strain. However, this situation changes when considering the linearization of volumetric strain in the next section. Secondly, we emphasize the trivial (but sometimes forgotten) fact that while the finite strain is rotationally invariant, its linearization is not. The importance in deriving a linearized strain from a rotationally invariant quantity is thus to ensure consistency in the limit of small deformations.

Remark 3.4

It is an important detail that while in the finite strain case the differential operators are time-dependent via their dependence on the jump operator $\mathbb {d}_{\Phi (t)}$, it is clear from Definition 3.8 and Example 3.2 (see, e.g., (3.19)), that the differential operator in the linearized strain is evaluated at ${\underline{\Phi }}$, and thus not time-dependent.

Fig. 4

Geometry considered in Example 3.3, where all mappings can be chosen such that $\phi _{0,j}(0)=0$ and that $(F_j)_{\Vert } = I$ near the origin

We give a second example to be more concrete.

Example 3.3

Consider a circular fracture defined by the unit disk in the plane as illustrated in Fig. 3. We choose an initial mapping that is the identity mapping near the origin, i.e., $X_{1} = B^{2}(x)\subset {\mathbb {R}}^{2}$, such that ${\underline{\phi }}_{0,1}(x) =\left[ x, 0\right] \in {\mathbb {R}}^{3}$ and $\Omega _{1} ={\underline{\phi }}_{0,1}(X_{1})$. Let the domains “above" and “below" be enumerated $2, 3$ with mappings $\phi _{0, 2}(x) = x$ and $\phi _{0, 3} = x$ on their respective domains $X_{2}$ and $X_{3}$, and let the extended mapping of $X_{1}$ be defined such that for $(x,y)\in X_{1}\times [-\epsilon ,\epsilon ]$ for some $\epsilon >0$. Then on the fracture, $x\in X_{1}$, the fully linearized strain is simply

(3.20)

while on the lower boundary of the fracture, indexed by say $j = 5$ such that $s_{5} = 1$ and with $\underline{\hat{\phi }}_{0,5} =\underline{{\hat{\phi }}}_{0,1}$, the strain is

$$\begin{aligned} \iota _{5}{\mathfrak {e}}(t;x) =\frac{1}{2}((D\iota _{5}{\mathfrak {u}})^{T} + D\iota _{5}{\mathfrak {u}}) = \frac{1}{2}\begin{bmatrix} (D_{\Vert }\iota _{5}{\mathfrak {u}}_{\Vert })^{T}+D_{\Vert }\iota _{5}{\mathfrak {u}}_{\Vert } \\ 0 \end{bmatrix}. \end{aligned}$$

(3.21)

Thus, $\iota _{5}{\mathfrak {e}}(t)$ is represented by a $3\times 2$ matrix similar to the horizontal components of the linearized strain. We note that when seen together, the fracture strain and (either of) the boundary strains can be combined and considered as a representation of the full strain.

For the interior of the matrix, we recover the standard linearized strain, such as in, e.g., $X_{3}$

$$\begin{aligned} \iota _{3}{\mathfrak {e}}(t;x) =\frac{1}{2}((D\iota _{3}{\mathfrak {u}})^{T} + D\iota _{3}{\mathfrak {u}}). \end{aligned}$$

(3.22)

3.4 Mixed-dimensional volume measure and the matrix trace operator

We will see in the continuation that flow is naturally formulated with pressures in $L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})$ and fluxes in $L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}})$, whereas the mechanics is naturally formulated with displacements in $L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ and strains in $L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$, cf. Examples 2.3 and 2.4. In order to develop appropriate coupling between flow and mechanics, we will also need operators mapping between these spaces, particularly to capture the effect of volume changes in the lower-dimensional domains (and conversely, the impact of fluid pressure on the total stress). As in the preceding sections, we present the volume measure in the setting of finite deformation first, and subsequently its linearization.

3.4.1 Finite deformation volumetric strain

By continuity, the determinant of the derivative of the transformation $J :=\det (D{\hat{\Phi }})$ contains the volume of the physical configuration relative to the volume of the reference domains. This is sufficient for the top-dimensional domains, $i\in I^{n}$; however, for lower-dimensional domains $i\in I^{d<n}$, it is of interest to include not only the static weight , but also the change associated with the jump in displacement. We therefore define the mixed-dimensional volume ${\mathfrak {J}}$ as follows.

Definition 3.9

Let $\omega _{i}^{j}\in C^{1}(X_{i})$ be a set of nonnegative weights for $i\in I^{d<n}$ and $j\in {\mathfrak {F}}$ with $s_{j} = i$. The mixed-dimensional volume density ${\mathfrak {J}}(\Phi )$ is defined such that:

• 1. For $i\in I^{n}$, then $\iota _{i}{\mathfrak {J}}(\Phi ) = \textrm{vol}({\textbf{F}}_{i})$, defined in (2.6).

• 2. For $i\in I^{d<n-1}$, then

  

  (3.23)

  where $J_{i}^{n-1}$ is the set of indexes $j\in \bigcup _{l\in I^{n-1}}^{}{\mathfrak {S}}_{l}$ such that $s_{j} = i$.

• 3. For $i\in I$, and $j\in I_{i}$, then $\iota _{j}{\mathfrak {J}}(\Phi ) = 0$.

The above definition is motivated as follows (confer also Fig. 5): As stated in the mixed-dimensional continuum assumption, Definition 2.2, it is natural to consider that the idealization of the fracture has some effective opening , and that in general its volume per unit area changes linearly with perpendicular opening (or closing) $\iota _{i}(\mathbb {d}\Phi )_{\perp }$ according to a proportionality constant $\omega _{i}^{i}$. Similarly, the cross-sectional area associated with an intersection will have some lower limit , and change proportionally to the opening of nearby fractures meeting at that intersection, according to weights $\omega _{i}^{j}$. We recognize that the definition stated in point 1 above can be written as special case of point 2 (by introducing the convention that for $i\in I^{n}$), but retain separate definitions for pedagogical clarity. Finally, it is typically not relevant to consider volume changes of the surfaces (point 3).

Fig. 5

Illustration of multi-scale contact mechanics, adapted from Oden and Martins [33]

3.4.2 Linearized volume change

We briefly recall under the small deformation assumption in fixed-dimensional continuum mechanics, the linearization of the volumetric change of a deformation becomes the matrix trace. To be precise, let us as in Sect. 3.3 consider a small deformation $\Phi (t)$ such that $\Phi (t) = \Xi {\mathfrak {u}}(t)+{\underline{\Phi }}$, where (3.13) holds. Then, since the determinant commutes with the matrix product, for any $i\in I^{n}$ we have the relationship (recall that in the case when $d_{i} = n$, then $\textrm{vol} = \det$, and for $d_{i}<n$, we have $\textrm{vol}({\textbf{F}}_{i})^{2} = \det ({\textbf{F}}_{i}^{T}{\textbf{F}}_{i})$):

$$\begin{aligned} \begin{aligned} \frac{J_{i}}{{\underline{J}}_{i}}&=\frac{\det ({\textbf{F}}_{i})}{\det ({{{{\underline{\varvec{F}}}}}}_{i})} = \det ({{{{\underline{\varvec{F}}}}}}_{i}^{-1}{\textbf{F}}_{i}) = \det ({{{{\underline{\varvec{F}}}}}}_{i}^{-1}D({\underline{\phi }}_{i}+u_{i})) = \det (I+{{{{\underline{\varvec{F}}}}}}_{i}^{-1}Du_{i})\\&= \det (I+{{{{\underline{\varvec{C}}}}}}_{i}^{-1}{{{{\underline{\varvec{F}}}}}}_{i}^{T}Du_{i})\\&= 1+\text { Trace }({{{{\underline{\varvec{C}}}}}}_{i}^{-1}{{{{\underline{\varvec{F}}}}}}_{i}^{T}Du_{i})+{\mathcal {O}}(\left| Du\right| ^{2})\\&= 1+\text { Trace }({{{{\underline{\varvec{C}}}}}}_{i}^{-1}{\textbf{e}}(u_{i}))+{\mathcal {O}}(\left| Du\right| ^{2}). \end{aligned} \end{aligned}$$

(3.24)

where we see that the linear term is the trace of the linearized strain, scaled by the deformation tensor of the undeformed state. Thus,

$$\begin{aligned} \frac{(DJ_{i})(u_{i})}{\det ({{{{\underline{\varvec{F}}}}}}_{i})} = \text { Trace }({{{{\underline{\varvec{C}}}}}}_{i}^{-1}{\textbf{e}}(u_{i})). \end{aligned}$$

(3.25)

In the same spirit, we first consider the linearization of the mixed-dimensional volume change, and secondly identify its interpretation as a matrix trace.

Definition 3.10

For some initial mapping ${\underline{\Phi }}$, and some deformation ${\mathfrak {u}}\in C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$, the linearized mixed-dimensional volume change ${\mathfrak {j}}\in L^{2}({\mathfrak {X}}, {\mathbb {R}})$ is defined as

$$\begin{aligned} {\mathfrak {j}}({\mathfrak {u}}) :=D{\mathfrak {J}}(\underline{\Phi })(\Xi {\mathfrak {u}}), \end{aligned}$$

(3.26)

where $D{\mathfrak {J}}({\underline{\Phi }})(\Psi )$ is the derivative of ${\mathfrak {J}}$ at ${\underline{\Phi }}$ for the perturbation $\Psi$.

Example 3.4

Continuing from Example 3.2, we consider the interpretation of the mixed-dimensional linearized volumetric strain on domains of various dimensionality.

• 1. For $i\in I^{n}$, then we obtain as in the fixed-dimensional case

  $$\begin{aligned} \frac{\iota _{i}{\mathfrak {j}}({\mathfrak {u}})}{\det ({{{{\underline{\varvec{F}}}}}}_{i})} = \text { Trace }({{{{\underline{\varvec{C}}}}}}_{i}^{-1}\iota _{i}{\mathfrak {e}}({\mathfrak {u}})). \end{aligned}$$

  (3.27)

• 2. For $i\in I^{n-1}$, then we first calculate the Fréchet derivative based on its action on $\Psi$ as (using that for a conforming forest, $(\iota _{i}\mathbb {d}_{{\underline{\Phi }}}{\underline{\Phi }})_{\perp } = 0$ by definition, and introducing the shorthand notation $\mathbb {d}_{\Delta }({\underline{\Phi }}+\epsilon \Psi ) =\mathbb {d}_{{\underline{\Phi }}+\epsilon \Psi }({\underline{\Phi }}+\epsilon \Psi )-\mathbb {d}_{{\underline{\Phi }}}({\underline{\Phi }}+\epsilon \Psi )$):

  

  (3.28)

  Here, we have used the continuity of the closest point projection in the definition of the jump operator to conclude that $\lim _{\epsilon \rightarrow 0}(\iota _{i}(\mathbb {d}_{\Delta }\Phi ))_{\perp } = 0$. Now it follows that:

  

  (3.29)

  Here, we have used that for $i\in I^{n-1}$, it holds that $\iota _{i}\mathbb {d}\Xi {\mathfrak {u}} = \iota _{i}\mathbb {d}{\mathfrak {u}}$.

• 3. For $i\in I^{d<n-1}$, a similar calculation as above gives

  

  (3.30)

  Here, we have used that the initial state is a conforming forest, thus since $s_{i} = s_{j}$ it holds that $\underline{\phi }_{j,i}^{*}$ is the identity.

• 4. For $i\in I$, and $j\in I_{i}$, then $\iota _{j}{\mathfrak {j}}({\mathfrak {u}}) = 0$.

The above calculation shows that for fractures and intersections, $i\in I^{d<n}$, the linearized volume change has two components, which have the natural interpretations of transverse opening and longitudinal stretching.

Motivated by the fixed-dimensional case, we wish to express the linearized volume change in terms of the linearized strain. We begin by observing that with the representation of $\Xi$ as defined above, and $i\in I^{d<n}$, then

$$\begin{aligned} \begin{aligned}&\text { Trace }\left( {{{{\underline{\varvec{C}}}}}}_{i}^{-1}\frac{1}{2}({{{{\underline{\varvec{F}}}}}}_{i}^{T}D(\iota _{i}\Xi {\mathfrak {u}})+(D(\iota _{i}\Xi {\mathfrak {u}}))^{T}{{{{\underline{\varvec{F}}}}}}_{i})\right) \\&\qquad \qquad =\frac{1}{\left| {\hat{I}}_{i}^{n}\right| }\sum _{j\in {\hat{I}}_{i}^{n}}^{}\text { Trace }\left( {{{{\underline{\varvec{C}}}}}}_{i}^{-1}\frac{1}{2}({{{{\underline{\varvec{F}}}}}}_{i}^{T}D(\iota _{j}\hat{{\mathfrak {u}}})+(D(\iota _{j}\hat{{\mathfrak {u}}}))^{T}{{{{\underline{\varvec{F}}}}}}_{i})\right) \\&\qquad \qquad =\frac{1}{\left| {\hat{I}}_{i}^{n}\right| }\sum _{j\in {\hat{I}}_{i}^{n}}^{}\text { Trace }({{{{\underline{\varvec{C}}}}}}_{i}^{-1}(\iota _{j}{\mathfrak {e}})_{\Vert }). \end{aligned} \end{aligned}$$

(3.31)

This calculation again exploits that the initial state is a conforming forest, since for $s_{i} = s_{j}$ it holds that ${{{{\underline{\varvec{F}}}}}}_{i} ={{{{\underline{\varvec{F}}}}}}_{j}$. Moreover, we recall that for $i\in I^{n-1}$, it holds that . This allows us to introduce the following:

Definition 3.11

The mixed-dimensional matrix trace operator $\mathfrak {T'}:L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})\rightarrow L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})$ is defined as follows. For ${\mathfrak {u}}\in L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$, and ${\mathfrak {e}} = {\mathfrak {D}}_{s}{\mathfrak {u}}\in L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ the trace operator satisfies

$$\begin{aligned} {\mathfrak {j}}({\mathfrak {u}}):={\mathfrak {J}}(\underline{\Phi })\mathfrak {T'e}={\mathfrak {J}}(\underline{\Phi }){\mathfrak {T}}'({\mathfrak {D}}_{s}{\mathfrak {u}}). \end{aligned}$$

(3.32)

Example 3.5

Continuing from Example 3.4, we consider the interpretation of the mixed-dimensional matrix trace on domains of various dimensionality.

• 1. For $i\in I^{n}$, then as in the fixed-dimensional case

  $$\begin{aligned} \iota _{i}({\mathfrak {T}}'{\mathfrak {e}}) = \text { Trace}({{{{\underline{\varvec{C}}}}}}_{i}^{-1}\iota _{i}{\mathfrak {e}}). \end{aligned}$$

• 2. For $i\in I^{d<n}$, then

  $$\begin{aligned} \iota _{i}({\mathfrak {T}}'{\mathfrak {e}}) = \sum _{j\in {\hat{I}}_{i}^{n-1}}^{}\omega _{i}^{j}(\iota _{j}{\mathfrak {e}})_{\perp }+\frac{1}{\left| {\hat{I}}_{i}^{n}\right| }\sum _{j\in {\hat{I}}_{i}^{n}}^{}\text { Trace }({{{{\underline{\varvec{C}}}}}}_{i}^{-1}(\iota _{j}{\mathfrak {e}})_{\Vert }). \end{aligned}$$

• 3. For $i\in I$, and $j\in I_{i}$, then $\iota _{j}({\mathfrak {T}}'{\mathfrak {e}}) = 0$.

Fig. 6

The mixed-dimensional trace operator $\mathfrak {T'}$ maps strains defined in $L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ to volumetric changes defined in $L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})$ for $i\in I^{2}$ (black), $i\in I^{1}$ (blue), $i\in I^{0}$ (red)

Example 3.6

Let $n = 2$ and consider the geometry from Fig. 6.

• Let $i\in I^{2}$. The linearized volumetric strain is given, as in the fixed-dimensional case by the trace of the linear strain, depicted by black arrows in Fig. 6.

• Let $i\in I^{1}$. Volumetric changes of $\Omega _{i}$ are measured using two metrics: changes with respect to tangential dilation are captured using the strain on the adjacent skins whereas changes in the aperture are measured using the locally defined strain vector. This is illustrated by the blue arrows in Fig. 6.

• Let $i\in I^{0}$. In this case, volumetric changes of $\Omega _{i}$ are measured using the perpendicular components of the strain vectors in the adjacent fractures. This is illustrated by the red arrows in Fig. 6.

• Let $n = 3$ and consider an extrusion of the geometry from Fig. 6 in the third dimension. Then, for $i\in I^{1}$, $\Omega _{i}$ represents an intersection line between fractures. Here, the tangential stretching is captured using the trace of the three-dimensional strain in the four corner lines, whereas the opening is measured using the change in aperture of the adjacent fractures.

We also define a mapping of pressures to stresses, generalizing the identity tensor. The following definition ensures that the duality from the fixed-dimensional case is preserved in reference space.

Definition 3.12

Let the mixed-dimensional identity operator ${\mathfrak {T}}: L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})\rightarrow L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$ be such that

$$\begin{aligned} \langle \mathfrak {Ta, b}\rangle _{{\mathfrak {X}}^{1}}&:=\langle \mathfrak {a, T'b}\rangle _{{\mathfrak {X}}^{n}},&\forall {\mathfrak {a}}&\in L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}}) \text { and } {\mathfrak {b}}\in L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n}). \end{aligned}$$

(3.33)

We summarize the linearized operators constructed in Sect. 3.3 and 3.4.2 with the diagram shown in Fig. 7. This diagram also includes the co-symmetric-gradient, $({\mathbb {D}}_{s} \cdot )$, which will be defined in Sect. 4.5. Its duality with ${\mathfrak {D}}_{s}$ is discussed in Sect. 5.1.

Fig. 7

The canonical structure of coupled poromechanics

4 Mixed-dimensional poromechanics

The preceding sections have laid the geometric foundations for considering deformation in mixed-dimensional geometries. In this section, we use these foundations to develop a theory of poromechanics. Our presentation will first establish the spatial structure of the system, as inferred from the geometry and differential operators defined above. We then consider the physical modeling of fluid flow, thence mechanics, and finally summarize the complete model. We close this section by discussing the relationship between mixed-dimensional poromechanics and classical equidimensional models. The exposition follows the standard modeling approach for poromechanics, which in the fixed-dimensional case is carefully reviewed by, e.g., Coussy [3], and summarized recently in [34].

We emphasize that all modeling in this section is considered on the reference domains ${\mathfrak {X}}$. This implies that all variables and derivatives are with respect to these domains, such that, e.g., when considering a mass density $\rho _{i}$, it is understood that this is mass density per unit volume form on $X_{i}$. With density as a concrete example, and domains $i$ with $d_{i} = 3$, this implies that density is a function of thermodynamics state (pressure and temperature), as well as the volumetric deformation (pullback of $n$-forms). On domains $i$ with $d_{i}<3$, we immediately get the interpretation of mass per measure, i.e., that the physical units of mass density are not simply mass per volume, but rather mass per area, length, or point.

4.1 Primary variables and structure

We consider the quasi-static poroelastic system as the composition of a mass conservative description of fluid flow, together with a mechanical system at equilibrium, as formulated on the reference domains ${\mathfrak {X}}^{k}$. As such, we consider displacement ${\mathfrak {u}}\in C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$, total stress ${\mathfrak {s}}\in C({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$, fluid mass density $\rho ({\mathfrak {p}})$ (represented for simplicity by pressure ${\mathfrak {p}}\in C({\mathfrak {X}}^{n}, {\mathbb {R}})$ in absence of thermal effects) and fluid mass flux ${\mathfrak {q}}\in C({\mathfrak {X}}^{n-1}, {\mathbb {R}})$ as primary variables. For the sake of exposition, we will frequently also refer to the strain ${\mathfrak {E}}\in C({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$, which we consider as a secondary variable derived from the displacement. In this section, we will assume enough regularity for all the definitions below to be valid: a weak formulation where regularity is considered in more detail is presented in Sect. 5.

Before we detail the physical constitutive laws, we note that the structure of the function spaces and constitutive laws essentially dictate the underlying canonical structure of the spatial operators. This is already summarized in Fig. 7 above, which, with reference to the primary variables defined above, can be summarized in terms of the following structures.

Conservation structure for fluid mass is stated on the space $C({\mathfrak {X}}^{n}, {\mathbb {R}})$, and can thus (apart from a source term ${\mathfrak {r}}_{{\mathfrak {m}}}$) only involve the quantities pressure ${\mathfrak {p}}$, volumetric strain ${\mathfrak {J}}(\Xi {\mathfrak {u}}+{\underline{\Phi }})$ and divergence of flux ${\mathfrak {D}} \cdot {\mathfrak {q}}$.

Similarly, the balance of mechanical forces is stated on the space $C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$, and can thus (apart from external forces ${\mathfrak {r}}_{{\mathfrak {s}}}$) only involve the quantities displacement ${\mathfrak {u}}$ and co-gradient of stress ${\mathbb {D}} \cdot {\mathfrak {s}}$. There are no other operators that map to this space.

We will see that the system has two constitutive laws. For the fluid, this is a binary relationship on $C({\mathfrak {X}}^{n-1}, {\mathbb {R}})$, involving the flux ${\mathfrak {q}}$ and the gradient ${\mathbb {D}} {\mathfrak {p}}$.

For the mechanical forces, this is in principle a binary relationship on $C({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$, involving the total stress ${\mathfrak {s}}$ and the strain ${\mathfrak {E}}$. However, as noted in Sect. 3.4.2, the linearized volumetric strain is actually a trace map from the strain space to the mass density space, and thus its dual, the identity map, enters the force balance.

In the following, we will present the full temporal modeling for finite strain and nonlinear constitutive laws. As a guide to this development, we provide the following structure, which we consider the canonical “Laplacian" of mixed-dimensional poromechanics, obtained from linearly combining the admissible relationships (with unit weights):

$$\begin{aligned}&\text {Force balance:}&\partial _t^2{\mathfrak {u}} - {\mathbb {D}}_s \cdot \overbrace{\left( {\mathfrak {D}}_s {\mathfrak {u}} + \mathfrak {Tp} \right) }^{{\mathfrak {s}}}&= {\mathfrak {r}}_{{\mathfrak {m}}}, \end{aligned}$$

(4.1a)

$$\begin{aligned}&\text {Mass balance:}&\partial _t{\mathfrak {p}} + {\mathfrak {T}}^* {\mathfrak {D}}_s {\mathfrak {u}} + {\mathfrak {D}} \cdot \underbrace{\left( - {\mathbb {D}} {\mathfrak {p}} \right) }_{{\mathfrak {q}}}&= {\mathfrak {r}}_{{\mathfrak {s}}}. \end{aligned}$$

(4.1b)

4.2 Fluid flow

As surveyed in the previous section, the conservation law for fluid mass can be stated as:

$$\begin{aligned} \partial _{t}{\mathfrak {m}} + {\mathfrak {D}} \cdot {\mathfrak {q}} = {\mathfrak {r}}_{{\mathfrak {m}}}, \end{aligned}$$

(4.2)

where the fluid mass content satisfies a constitutive relation dependent the mechanical strain and fluid pressure,

$$\begin{aligned} \mathfrak {m = m}(\mathfrak {E, p})\in C({\mathfrak {X}}^{n}, {\mathbb {R}}). \end{aligned}$$

(4.3)

The fluid mass flux is assumed to satisfy a linear proportionality with a gradient, specifically the co-divergence ${\mathbb {D}}$, of fluid pressure, which corresponds to Darcy’s law for intact porous rock ($i\in I^{n}$) and laminar flow fractures ($i\in I^{n-1}$) [5, 6]:

$$\begin{aligned} {\mathfrak {q}} = \kappa (\mathfrak {E,q})(\mathbb {-D}{\mathfrak {p}}+{\mathfrak {r}}_{{\mathfrak {q}}}). \end{aligned}$$

(4.4)

With ${\mathfrak {r}}_{{\mathfrak {q}}}$, the contribution due to gravity. Note that we have allowed the material coefficient to depend on the strain, in order to accommodate that the fracture conductivity may be altered when there is slip along the fracture. Moreover, the dependency on the fluid flux accommodates nonlinear relationships such as the Darcy–Forchheimer law.

Remark 4.1

The permeability of a rock domain $\iota _{i}\kappa (\mathfrak {E,q})$ with $d_{i} = 3$ is the permeability on the reference space $X_{i}$. Given a permeability $K_{i}$ on the physical domain, these are related by the usual transformation rules [1], e.g.,

$$\begin{aligned} K_{i} = \det ({\textbf{F}}_{i})^{-1}{\textbf{F}}_{i}(\iota _{i}\kappa (\mathfrak {E,q})){\textbf{F}}_{i}^{T}. \end{aligned}$$

(4.5)

In particular, this implies that a scalar permeability $K_{i}$ in physical space may nevertheless be represented by an anisotropic tensor $\iota _{i}\kappa$ on $X_{i}$ (and opposite). On the other hand, it is important to note that (4.5) preserves symmetry properties, thus symmetry of $\iota _{i}\kappa$ corresponds to symmetry of $K_{i}$. Similar comments apply to all material properties introduced in this and the following sections.

4.3 Mechanical response

We consider balance of momentum as the basis for modeling the mechanical response. We state the equilibrium assumption recognizing that ${\mathbb {D}} \cdot$ (the co-gradient) is a divergence operator on the top-dimensional domains and their boundaries. Thus for a stress variable ${\mathfrak {s}}$ (with the interpretation of a Piola-Kirchhoff stress of the second kind), a change in momentum $\rho _{r}\partial _{t}^{2}{\mathfrak {u}}$ in view of a mixture mass density for fluid and rock $\rho _{r}$ leads us to the following balance of momentum:

$$\begin{aligned} \rho _{r}\partial _{t}^{2}{\mathfrak {u}}-{\mathbb {D}} \cdot (\varvec{{\mathfrak {F}}}{\mathfrak {s}}) ={\mathfrak {r}}_{{\mathfrak {s}}}, \end{aligned}$$

(4.6)

with ${\mathfrak {r}}_{{\mathfrak {s}}}$ describing body forces. This vector equation has components associated with the basis vectors of ${\mathbb {R}}^{n}$. The deformation $\varvec{{\mathfrak {F}}}$ enters the momentum balance to transform forces on $T{\hat{X}}_{i}\rightarrow {\mathbb {R}}^{n}$ [1] with ${\hat{X}}_{i}$ the extended domain defined in Definitions 2.6 and 2.7.

The mechanical stress state depends on both the strain and the fluid pressure. As we are primarily interested in moderate deformations, we restrict our attention to linear (Saint Venant–Kirchhoff) stress–strain response in the porous rock and its boundaries, $j\in {\mathfrak {S}}_{i}$, with $i\in I^{n}$ [32, 35]:

$$\begin{aligned} \iota _{j}{\mathfrak {s}}(\mathfrak {E,p}) = C_{j}:\iota _{j}{\mathfrak {E}}-\iota _{j}(\alpha \mathfrak {Tp}). \end{aligned}$$

(4.7)

With $\alpha$ a positive, symmetric linear operator generalizing the Biot–Willis constant [3] and $C_{j}$ the fourth-order stiffness tensor.

Physical reality demands a greater generality than linear mechanical response of the fracture, since both friction and contact mechanics will in general be given by nonlinear relationships. Moreover, for finite deformation, determination of contact itself may be a nonlinear problem. Therefore, we allow for a (possibly nonlinear and non-local in space) constitutive law ${\mathfrak {A}}_{i}$ between stress and strain to be defined on fractures $i\in I^{n-1}$ and their boundaries [2]. We impose this by introducing ${\mathfrak {A}}_{i}$ as a binary relation (see Appendix ) such that

$$\begin{aligned} ( \iota _{i}({\mathfrak {s}}+\alpha \mathfrak {Tp}), \iota _{i}(\hat{\gamma }+{\check{\gamma }}\partial _{t}){\mathfrak {E}})\in {\mathfrak {A}}_{i}. \end{aligned}$$

(4.8)

Here, ${\hat{\gamma }}$ and ${\check{\gamma }}$ are given parameters with ${\hat{\gamma }}+{\check{\gamma }} = 1$ and $\hat{\gamma }{\check{\gamma }} = 0$ on each $X_{i}$ that allow us to model relationships between either stress and strain or stress and strain rate, respectively. In Sect. 4.6, we will give examples of particular choices of friction and contact laws.

We combine the constitutive laws for the bulk and fracture subdomains in the notation of a mixed-dimensional binary relation:

$$\begin{aligned} ({\mathfrak {s}}+\alpha \mathfrak {Tp, }({\hat{\gamma }}+\check{\gamma }\partial _{t}){\mathfrak {E}}) \in {\mathfrak {A}}. \end{aligned}$$

(4.9)

Remark 4.2

The presence of the material law (4.7) on boundaries implies that the constitutive modeling is general enough to allow for materials with coated boundaries (say, covered by a thin membrane), or other disturbances of the material parameters associated with the boundaries of the domain. A perfectly homogeneous material with no disturbance in material parameter at its boundaries will then be a degenerate case of the model with ${\mathfrak {A}}_{j} = 0$ (for $j\in I_{i}$ with $i\in I^{n}$).

4.4 Governing equations for mixed-dimensional finite strain poromechanics

We summarize the above developments in the system of equations for mixed-dimensional poroelastic fractured media.

By taking the right-hand sides ${\mathfrak {r}}_{{\mathfrak {s}}}$ and ${\mathfrak {r}}_{{\mathfrak {m}}}$ together with the mixture density $\rho _{r}$ as given, as well as appropriate initial and boundary conditions, the above set of equations are formally closed. Superficially, we identify that we have 7 equations for the 7 (mixed-dimensional) unknowns $\left[ {\mathfrak {u}}, {\mathfrak {s}}, \varvec{{\mathfrak {F}}}, {\mathfrak {m}}, {\mathfrak {q}}, {\mathfrak {E}}, {\mathfrak {p}}\right]$. A careful counting of scalar and vector equations on domains of various dimensionality supports this claim. The field equations should be supplemented by appropriate boundary conditions, we will discuss one such choice in Sect. 5.

A well-posedness theory for the general finite deformation mixed-dimensional model is not within reach, since it would require us to simultaneously address open questions in both contact mechanics and poromechanics. On the other hand, in the context of infinitesimal deformation, fairly general results are nevertheless possible to establish. In the next sections, we will therefore restrict our attention to linearized strain and make specific assumptions on the constitutive laws that allow us to rigorously establish an example of a well-posed model for mixed-dimensional poromechanics.

4.5 Governing equations for mixed-dimensional linearized strain poromechanics

We continue by considering the model (4.10) in the context of the mixed-dimensional linearized strain from Definition 3.8. Thus, we assume that the strain is given by

$$\begin{aligned} {\mathfrak {E}}({\mathfrak {u}})\approx {\mathfrak {e}}({\mathfrak {u}}) ={\mathfrak {D}}_{s}{\mathfrak {u}}. \end{aligned}$$

Our exposition simplifies (and as we will see, symmetrizes) by including the transformation of the reference configuration in the definition of the divergence operators (in analogy to the fixed-dimensional case, see, e.g., [1]).

Definition 4.1

For an initial configuration $\underline{{\hat{\Phi }}}$ with derivative $\underline{\varvec{{\mathfrak {F}}}}$, let the mixed-dimensional co-symmetric-gradient $({\mathbb {D}}_{s} \cdot ):C^{1}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})\rightarrow C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ be defined such that

$$\begin{aligned} {\mathbb {D}}_{s}\cdot {\mathfrak {s}} :={\mathbb {D}} \cdot (\underline{\varvec{{\mathfrak {F}}}}{\mathfrak {s}}). \end{aligned}$$

Remark 4.3

Recall that in the fixed-dimensional case, the symmetric gradient is adjoint to the divergence applied to symmetric tensor fields. This forms a key tool in the well-posedness for fixed-dimensional elasticity models and we note that this property remains valid in the mixed-dimensional setting, as we show in Sect. 5.2.

Secondly, we linearize the fluid mass content relationship ${\mathfrak {m}}(\mathfrak {e,p})$. Consistent with the relationships discussed in [3], we obtain

$$\begin{aligned} \mathfrak { m}(\mathfrak {e,p})\approx {\mathfrak {m}}_{0}+\mathfrak {T'}(\alpha {\mathfrak {e}})+\beta {\mathfrak {p}}. \end{aligned}$$

Here, ${\mathfrak {m}}_{0}$ is the initial mass content, $\mathfrak {T'}$ is the mixed-dimensional identity operator from Definition 3.12, $\alpha$ is the generalized Biot–Willis constant from (4.7), and $\beta$ is the specific storativity.

We summarize these simplifications in the following system of equations

4.6 Connection to classical continuum mechanical formulations

In this section, we identify how fixed-dimensional formulations of contact mechanics and fluid flow in fractured porous rocks appear as special cases of the governing equations derived in the preceding sections and as summarized in Sect. 4.4. This identification will serve as support for our claim that our development is a consistent generalization of accepted mathematical descriptions, extended to our current setting of mixed-dimensional geometries.

4.6.1 Contact mechanics

We illustrate the connection to contact mechanics by omitting all fluid considerations, and restricting our attention to a single fracture in the context of linearized strain, such as described in Example 3.3, and illustrated by Fig. 4. The governing equations, assuming linearized strain, then simplify to (indexing of $i$ corresponds to indexing in Example 3.3).

Hooke’s law for solids:

$$\begin{aligned} \sigma _{i}&= C_{i}:\varepsilon (u_{i})&\text {for } i = 2,3. \end{aligned}$$

(4.12a)

Conservation of momentum for solids:

$$\begin{aligned} \rho _{r,i}\partial _{t}^{2}u_{i}-\nabla \cdot (\underline{{{\hat{{\textbf {F}}}}}}_{i}\sigma _{i})&= r_{s,i}&\text {for } i = 2,3. \end{aligned}$$

(4.12b)

Hooke’s law for surfaces:

$$\begin{aligned} \sigma _{j}&= C_{j}:\varepsilon (u_{j})&\text {for } j= 4,5. \end{aligned}$$

(4.12c)

Conservation of momentum on surfaces:

$$\begin{aligned} \rho _{r,j}\partial _{t}^{2}u_{j}-\nabla _{\Vert } \cdot (\underline{\hat{{\textbf{F}}}}_{j}\sigma _{j})+(\sigma _{j-2} \cdot n_j-\sigma _{1}) = r_{s,j} \text {for } j= 4,5. \end{aligned}$$

(4.12d)

Constitutive law for fracture:

$$\begin{aligned} F(\varepsilon _{i},{\dot{\varepsilon }}_{i},\sigma )&= 0&\text {for }i = 1. \end{aligned}$$

(4.12e)

These equations embody several classical formulations. Surface stress of elastic materials is an important phenomenon at some scales [35], as we will discuss in Sect. 4.6.4. In this subsection, we will disregard these terms, and consider a degenerate model with $C_j = 0$ and hence $\sigma _j = 0$ on domains $j = 4,5$. In the quasi-static case, $\rho _{r,i}\partial _{t}^{2}u_{i} = 0$, the above equations now exactly correspond to the classical equations of elasticity and material contact, as summarized in the book of Kikuchi and Oden (see, e.g., the exposition in chapters 2 and 13 of [2]). We highlight the importance of two particular choices for the constitutive law for fracture.

Example 4.1

(Signorini problem) Decomposing as in Sect. 4 the fracture strain into its tangential part (sliding) $\iota _{1}{\mathfrak {e}}_{\Vert } =\frac{1}{2}{{{{\underline{\varvec{F}}}}}}_{i}^{T}(\iota _{i}\mathbb {d}_{\underline{\Phi }}{\mathfrak {u}})_{\parallel }$ and perpendicular part (opening) , we first state the constitutive model for frictionless contact (see, e.g., equations (2.31) of [2]), where the absence of friction is given by:

$$\begin{aligned} \iota _{1}{\mathfrak {s}}_{\Vert } = 0, \end{aligned}$$

(4.13a)

and contact mechanics is given by the Karush–Kuhn–Tucker (KKT) triplet

$$\begin{aligned} \begin{aligned} (\iota _{1}{\mathfrak {s}}_{\perp }) (\iota _{1}{\mathfrak {e}}_{\perp })&= 0,\\ \iota _{1}{\mathfrak {s}}_{\perp }&\le 0,\\ \iota _{1}{\mathfrak {e}}_{\perp }&\ge 0. \end{aligned} \end{aligned}$$

(4.13b)

The two inequalities state that the materials cannot interpenetrate, and the contact cannot be tensile, while the equality states that one of the two conditions must hold as an equality. We recognize that equations (4.13) are of the form given by the binary relation (4.8), with $\iota _{1}{\check{\gamma }} = 0$.

Example 4.2

(Rough surfaces) The presence of in the perpendicular part of the strain measure, together with the macroscopic condition that the deformation preserves the mixed-dimensional nature of the problem, essentially provides a multi-scale representation of strain. Indeed with , we realize that fracture opening is measured relative to a finer scale than the rest of the deformation. This makes sense relative to the mixed-dimensional continuum assumption, Definition 2.2, since macroscopically, the fracture always has negligible transversal width, while the strain nevertheless measures perturbations of the fracture at the scale of .

At the length-scale transversal to a fracture opening, it is well known that the compression and crushing of micro-roughness can significantly impact the stress–strain response. This is illustrated in Fig. 5, which has been adapted from the classical presentation by Oden and Martins [33]. With access to such representations, it is suggested that an appropriate stress–strain model for a fracture (see Chapters 11 and 13 of [2]) depends on both rate of compression and effective opening (in these expressions the plus sign indicates that only positive values are considered, e.g., $(a)_{+} = \max (a,0)$, and $C_{1}^{l}$ correspond to material constants):

$$\begin{aligned} \sigma _{1,\perp } = -C_{1}^{1}(-\varepsilon _{1,\perp })_{+}^{C_{1}^{2}}-C_{1}^{3}(-{\dot{\varepsilon }}_{1,\perp })_{+}^{C_{1}^{4}}. \end{aligned}$$

(4.14)

Complementing the perpendicular stress–strain law for fracture is the Coulomb law of friction, expressed as the KKT triplet:

$$\begin{aligned} \lambda ^{2}(c_{1}^{5}\sigma _{1,\perp }-\left| \sigma _{1,\Vert }\right| )&= 0, \end{aligned}$$

(4.15a)

$$\begin{aligned} c_{1}^{5}\sigma _{1,\perp }-\left| \sigma _{1,\Vert }\right|&\ge 0, \end{aligned}$$

(4.15b)

$$\begin{aligned} \lambda ^{2}&\ge 0, \end{aligned}$$

(4.15c)

where $\lambda ^{2}$ is a Lagrange multiplier associated with sliding:

$$\begin{aligned} {\dot{\varepsilon }}_{1,\Vert } = \lambda ^{2}\sigma _{1,\Vert }. \end{aligned}$$

(4.16)

While our framework as stated does not allow for the full generality of (4.14), with both $C_{1}^{1}$ and $C_{1}^{3}$ nonzero, these contact laws are nevertheless admissible as a binary inclusion when either $C_{1}^{1}$ or $C_{1}^{3}$ are zero.

4.6.2 Fluid flow in rigid fractured porous media

In the absence of mechanical deformation, the mixed-dimensional equations (4.10) have previously been shown to be algebraically equivalent to reduced-dimensional formulations of flow in porous media [4]. These equations have become quite popular over the last decade, as illustrated by a recent literature review [18]. An important consideration is the validity of reduced-dimensional models for fracture flow. This has been extensively studied, and is by now well established in the single-phase regime considered herein [5, 8, 36].

4.6.3 Fluid flow in deformable porous media

Let us now validate the model equations in the case without fractures, wherein the mixed-dimensional geometry trivially reduces to the normal fixed-dimensional geometry, which is to say that $I = I^{n} = \{ 1\}$, and $\Omega _{1} = Y$.

A review of the definitions in Sect. 2.3 and 2.4 for this case of a single domain now verify that all the mixed-dimensional functions revert to their standard definitions from calculus. Thus, e.g., a deformation ${\mathfrak {u}}\in C({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$ is identically equal to $u\in C(Y,{\mathbb {R}}^{n})$, similarly, the differential operators also reduce to their fixed-dimensional counterparts ${\mathfrak {D}} \sim {\mathbb {D}} \sim \nabla$ while ${\mathfrak {D}} \cdot \sim {\mathbb {D}} \cdot \sim \nabla \cdot$.

Equations (4.10) are thus equivalent to the same equations written in “Latin letters", which correspond exactly to the standard model for poromechanics subjected to large deformations, as summarized in, e.g., Table 3.4 of [3].

4.6.4 Coated deformable solids

In this final example, we will consider the special case of elastic solids with surface coatings. The continuum theory for elastic material surfaces goes back to Gurtin and Murdoch [35], and our general mixed-dimensional model includes some aspects of their theory, notably the momentum balance and elastic constitutive law.

To illustrate this, we consider a single internal domain $\Omega _{2}$, together with a lower-dimensional domain contained on a part of its exterior boundary $\Omega _{4}\subset \partial Y$. The numbering is chosen so that the example can be considered as the upper domain of Fig. 4, and it is then the interpretation that the bottom part of this domain is coated. In terms of the reference configurations, the governing equations in the finite strain case are now given by equations (4.12a) and (4.12b) for the bulk material $i = 2$. For the surface $j = 4$, the surface stress is given by (4.12c), while the momentum balance simplifies to

$$\begin{aligned} \rho _{r,2}\partial _{t}^{2}u_{2}-\nabla _{\Vert } \cdot ({\textbf{F}}_{4}\sigma _{4})+\sigma _{2} \cdot n&= r_{s,j}&\text { for } j&= 4. \end{aligned}$$

(4.17)

We recognize our surface momentum balance (4.17) and our surface stress (4.12c) as equation (6.1) and the first equation of section 7 of reference [35], respectively. An important detail is the application to curved surfaces, in which case the surface divergence term becomes (omitting the subscript 4 for clarity):

(4.18)

For curved surfaces, this relationship expresses the fact that the gradient of the normal component of a surface enters the balance of momentum. This geometric relationship can also be expressed in terms of mean curvature, as in equation (2.9) of [35].

5 Well-posedness of the mixed-dimensional linearized strain model

In this section, we will analyze the linearized strain model (4.11). While strain is infinitesimal in the bulk, the model retains two important (and non-trivial) nonlinearities in the constitutive laws: First, the frictional contact law at the fractures, and secondly nonlinear relationship between pressure gradients and fluid flows. This implies that the well-posedness results presented in this section are a significant generalization of the analysis presented in previous work on flow and deformation in fractured media (see, e.g., the recent paper [13] and references therein).

However, our analysis below excludes one important nonlinear dependence: That of fracture permeability on the fracture opening. We will discuss this point in Sect. 6. Throughout this section, we will rely on the theory for maximally monotone binary relations, as well as the existence theory of evolutionary equations, both summarized in Appendix .

5.1 The space of symmetric tensor fields

By Definition 3.8, the linearized strain ${\mathfrak {e}}$ has symmetries in the sense that the representation of $\iota _{i}{\mathfrak {e}}$ with $i\in I^{n}$ is a symmetric tensor field in ${\mathbb {R}}^{n\times n}$. Moreover, for all boundaries $j\in I_{i}\cap {\mathfrak {F}}^{1}$, we recall that the tangential components $\iota _{j,\parallel }{\mathfrak {e}}$ form a symmetric tensor in ${\mathbb {R}}^{(n-1)\times (n-1)}$, whereas the normal components $\iota _{j,\perp }{\mathfrak {e}}$ are zero, cf. Example 3.2. These properties are captured in the definition of the following function space

$$\begin{aligned} {\mathfrak {G}}:=\left\{ {\mathfrak {e}}\in L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})\left| \begin{array}{ll} \text {asym}(\iota _{j,\parallel }{\mathfrak {e}}) = 0, &{} \quad \forall j\in {\mathfrak {F}}^{1}, \\ \iota _{j,\perp }{\mathfrak {e}} = 0, &{}\quad \forall j\in {\mathfrak {F}}^{1}\cap I_{i}, i\in I^{n} \end{array} \right. \right\} . \end{aligned}$$

By the usual Cauchy arguments, the stress–strain relationship ${\mathfrak {A}}$ retains these properties, i.e., we have that the stress ${\mathfrak {s}}$ belongs to ${\mathfrak {G}}$ as well.

Note that this space generalizes the space of symmetric tensors that is often used in fixed-dimensional linear elasticity. A strong tool in that setting is the adjointness between the divergence on symmetric tensors and the symmetric gradient on vector fields. We now show that this adjointness property is retained in the mixed-dimensional setting.

Lemma 5.1

The mixed-dimensional symmetric gradient and the co-symmetric gradient operators satisfy the following integration by parts formula for all ${\mathfrak {s}} \in \textrm{dom}({\mathbb {D}}_s \cdot ) \subset {\mathfrak {G}}$ and ${\mathfrak {u}} \in \textrm{dom}(\mathring{{\mathfrak {D}}}_s ) \subset L^2({\mathfrak {X}}^0, {\mathbb {R}}^n)$:

$$\begin{aligned} \left\langle {\mathbb {D}}_{s}\cdot {\mathfrak {s}}, {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{0}}+\left\langle {\mathfrak {s}}, \mathring{{\mathfrak {D}}}_s {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{1}} = 0. \end{aligned}$$

Proof

We first observe using Definition 4.1 and the duality from Definition 2.23 that

$$\begin{aligned} \left\langle {\mathbb {D}}_{s}\cdot {\mathfrak {s}}, {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{0}} :=\left\langle {\mathbb {D}} \cdot (\varvec{{\mathfrak {F}}}{\mathfrak {s}}), {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{0}} = -\left\langle {\mathfrak {s}}, \varvec{{\mathfrak {F}}}^{T}\mathring{{\mathfrak {D}}}{\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{1}}. \end{aligned}$$

Thus, it remains to show that $\left\langle {\mathfrak {s}}, \varvec{{\mathfrak {F}}}^{T}\mathring{{\mathfrak {D}}}{\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{1}} = \left\langle {\mathfrak {s}}, \mathring{{\mathfrak {D}}}_s {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{1}}$ for ${\mathfrak {s}} \in \textrm{dom}({\mathbb {D}}_s \cdot )$ and ${\mathfrak {u}} \in \textrm{dom}(\mathring{{\mathfrak {D}}}_s )$. We do this by dimension:

• 1. For $i\in I^{n}$, the symmetry of $\iota _{i}{\mathfrak {s}}$ gives us

  $$\begin{aligned} \left\langle \iota _{i}{\mathfrak {s}}, \iota _{i}\varvec{{\mathfrak {F}}}^{T}\mathring{{\mathfrak {D}}}({\mathfrak {u}})\right\rangle _{X_{i}}&=\left\langle \iota _{i}{\mathfrak {s}},{{{{\underline{\varvec{F}}}}}}_{i}^{T}D\iota _{i}{\mathfrak {u}}\right\rangle _{X_{i}}\\&=\left\langle \iota _{i}{\mathfrak {s}},\frac{1}{2}({{{{\underline{\varvec{F}}}}}}_{i}^{T}D\iota _{i}{\mathfrak {u}}+(D\iota _{i}{\mathfrak {u}})^{T}{{{{\underline{\varvec{F}}}}}}_{i})\right\rangle _{X_{i}}\\&=\left\langle \iota _{i}{\mathfrak {s}}, \iota _{i}\mathring{{\mathfrak {D}}}_{s}{\mathfrak {u}}\right\rangle _{X_{i}} \end{aligned}$$

• 2. For $j\in I_{i}$ with $i\in I^{n}$, the same argument as above applies for the tangential components. On the other hand, the normal components of $\iota _{j}{\mathfrak {s}}$ are zero, immediately giving us

  $$\begin{aligned} \left\langle \iota _{j}{\mathfrak {s}}, \iota _{j}\varvec{{\mathfrak {F}}}^{T}\mathring{{\mathfrak {D}}}({\mathfrak {u}})\right\rangle _{X_{j}} =\left\langle \iota _{j}{\mathfrak {s}}, \iota _{j}\mathring{{\mathfrak {D}}}_{s}{\mathfrak {u}}\right\rangle _{X_{j}} \end{aligned}$$

• 3. For $j\in {\mathfrak {S}}_{i}$ with $i\in I^{n-1}$, we obtain

  $$\begin{aligned} \iota _{j}\varvec{{\mathfrak {F}}}^{T}\mathring{{\mathfrak {D}}}({\mathfrak {u}}) =\hat{{{{{\underline{\varvec{F}}}}}}}_{j}^{T}\iota _{j}\mathbb {d}{\mathfrak {u}} = \iota _{j}\mathring{{\mathfrak {D}}}_{s}{\mathfrak {u}} \end{aligned}$$

Since this covers all $j\in {\mathfrak {F}}^{1}$, we have $\varvec{{\mathfrak {F}}}^{T}\mathring{{\mathfrak {D}}}({\mathfrak {u}}) =\mathring{{\mathfrak {D}}}_{s}$ on ${\mathfrak {X}}^{1}$ and thus $\left\langle {\mathbb {D}}_{s} \cdot {\mathfrak {s}}, {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{0}} = -\left\langle {\mathfrak {s}}, \mathring{{\mathfrak {D}}}_s {\mathfrak {u}}\right\rangle _{{\mathfrak {X}}^{1}}$. $\square$

We will require different treatment of the stress–strain relationships depending on whether the relationship concerns strains or strain rates. Recall that we have introduced the parameters ${\hat{\gamma }}$, respectively ${\check{\gamma }}$, in (4.8) to make this distinction. Using these parameters, we define the restricted identity operators as follows.

Definition 5.1

Let the restricted mixed-dimensional identity operators $\check{{\mathfrak {T}}},\hat{{\mathfrak {T}}}:L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})\rightarrow \tilde{{\mathfrak {S}}}$ be defined as

$$\begin{aligned} \check{{\mathfrak {T}}}&:={\check{\gamma }}{\mathfrak {T}}&\text {and}{} & {} \hat{{\mathfrak {T}}}&={\hat{\gamma }}\mathfrak {T.} \end{aligned}$$

Similarly, let the restricted mixed-dimensional trace operators be defined as $\check{{\mathfrak {T}}}' :={\check{\gamma }}{\mathfrak {T}}'$ and $\hat{{\mathfrak {T}}}' :={\hat{\gamma }}{\mathfrak {T}}'$.

5.2 Model equations for well-posedness analysis

Let us consider the system of equations from (4.11) with the goal of obtaining a system of four equations and four variables. In particular, we aim to rewrite the system in terms of fluid pressure ${\mathfrak {p}}$, fluid flux ${\mathfrak {q}}$, and two new variables; namely the bulk velocity ${\mathfrak {v}}$, and an augmented stress $\tilde{{\mathfrak {s}}}$:

$$\begin{aligned} {\mathfrak {v}}&:=\partial _{t}{\mathfrak {u}},&\tilde{{\mathfrak {s}}}&:={\mathfrak {s}}+\alpha \check{{\mathfrak {T}}}\mathfrak {p.} \end{aligned}$$

(5.1)

Note that $\tilde{{\mathfrak {s}}}$ corresponds to the mechanical stress in the fractures as this is the natural variable for which frictional contact laws are formulated. In the bulk, it equals the original, poroelastic stress ${\mathfrak {s}}$ since $\check{{\mathfrak {T}}}$ is zero there.

We proceed in four steps. First, substituting the definitions of ${\mathfrak {v}}$ and $\tilde{{\mathfrak {s}}}$ in (4.11a), the balance of forces becomes

$$\begin{aligned} \rho _{r}\partial _{t}{\mathfrak {v}}-{\mathbb {D}}_{s} \cdot (\tilde{{\mathfrak {s}}}-\alpha \check{{\mathfrak {T}}}{\mathfrak {p}}) ={\mathfrak {r}}_{{\mathfrak {s}}}. \end{aligned}$$

(5.2)

We make the following assumptions on $\rho _{r}$:

Assumption 5.1

$\rho _{r}$ is a coercive, linear operator with coercivity constant $c_{\rho }>0$.

Second, we consider the stress–strain relationships. For that, we first take the derivative of (4.11c) with respect to time and use the commutativity of $\partial _{t}$ and ${\mathfrak {D}}_{s}$ (cf. Remark 3.4):

$$\begin{aligned} \partial _{t}{\mathfrak {e}} = \partial _{t}{\mathfrak {D}}_{s}{\mathfrak {u}} = {\mathfrak {D}}_{s}{\mathfrak {v}}. \end{aligned}$$

Substituting this in (4.11d) together with the definition of $\tilde{{\mathfrak {s}}}$ from (5.1), we obtain

$$\begin{aligned} (\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}\mathfrak {p, }({\hat{\gamma }}\partial _{t}^{-1}+\check{\gamma }){\mathfrak {D}}_{s}{\mathfrak {v}})\in {\mathfrak {A}}. \end{aligned}$$

(5.3)

Note that $\partial _{t}^{-1}$ implies integration in time. In the surrounding bulk ($i\in I^{n}$), we assume that ${\hat{\gamma }} = 1$ (and ${\check{\gamma }} = 0$) giving us a stress–strain relationship $\hat{{\mathfrak {A}}}$. Conversely, we let ${\check{\gamma }} = 1$ in the fractures leading to a stress–strain rate relationship describing (frictional) contact $\check{{\mathfrak {A}}}$. More precisely, we define the restricted operators

$$\begin{aligned} \check{{\mathfrak {A}}}&:={\check{\gamma }}{\mathfrak {A}},&\text {and}{} & {} \hat{{\mathfrak {A}}}&:=\hat{\gamma }{\mathfrak {A}}, \end{aligned}$$

such that ${\mathfrak {A}} = \check{{\mathfrak {A}}} + \hat{{\mathfrak {A}}}$. We now continue by making the following assumptions:

Assumption 5.2

$\check{{\mathfrak {A}}}$ is bounded and $\check{c}$-maximal monotone relation for some $\check{c}>0$, c.f. Definition A.5. Moreover, $(0,0)\in \check{{\mathfrak {A}}}$. We emphasize that this means that for $({\mathfrak {s}}_{1},\dot{{\mathfrak {e}}}_{1}),({\mathfrak {s}}_{2},\dot{{\mathfrak {e}}}_{2})\in \check{{\mathfrak {A}}},$ we have

$$\begin{aligned} \left\langle \iota _{j}({\mathfrak {s}}_{1}-{\mathfrak {s}}_{2}),\iota _{j}(\dot{{\mathfrak {e}}}_{1}-\dot{{\mathfrak {e}}}_{2})\right\rangle _{X_{j}}&\ge \check{c}\left| \iota _{j}({\mathfrak {s}}_{1}-{\mathfrak {s}}_{2})\right| _{X_{j}}^{2},&\forall i&\in I^{n-1}, j\in {\mathfrak {F}}^{1}\cap {\mathfrak {S}}_{i}. \end{aligned}$$

Assumption 5.3

$\hat{{\mathfrak {A}}}$ is a coercive linear operator with coercivity constant ${\hat{c}}>0$. Thus, for $({\mathfrak {s}}_{1},{\mathfrak {e}}_{1}),({\mathfrak {s}}_{2},{\mathfrak {e}}_{2})\in \hat{{\mathfrak {A}}}$, it follows that

$$\begin{aligned} \left\langle \iota _{j}({\mathfrak {s}}_{1}-{\mathfrak {s}}_{2}),\iota _{j}({\mathfrak {e}}_{1}-{\mathfrak {e}}_{2})\right\rangle _{X_{j}}&\ge {\hat{c}}\left| \iota _{j}({\mathfrak {s}}_{1}-{\mathfrak {s}}_{2})\right| _{X_{j}}^{2},&\forall i&\in I^{n}, j\in {\mathfrak {F}}^{1}\cap {\mathfrak {S}}_{i}. \end{aligned}$$

With the assumed linearity of $\hat{{\mathfrak {A}}}$, we have

$$\begin{aligned} (\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}\mathfrak {p,}\hat{\gamma }\partial _{t}^{-1}{\mathfrak {D}}_{s}{\mathfrak {v}})\in \hat{{\mathfrak {A}}}\Leftrightarrow (\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}\mathfrak {p,}\hat{\gamma }{\mathfrak {D}}_{s}{\mathfrak {v}})\in \partial _{t}\hat{{\mathfrak {A}}}. \end{aligned}$$

Together with $(\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}\mathfrak {p,}\check{\gamma }{\mathfrak {D}}_{s}{\mathfrak {v}})\in \check{{\mathfrak {A}}}$, the stress–strain (rate) relationships become

$$\begin{aligned} (\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}\mathfrak {p, }{\mathfrak {D}}_{s}{\mathfrak {v}})\in \check{{\mathfrak {A}}}+\partial _{t}\hat{{\mathfrak {A}}}. \end{aligned}$$

(5.4)

Remark 5.1

More generally, we may assume 5.2 for ${\mathfrak {A}}_{i}$ if $\iota _{i}{\check{\gamma }} = 1$ and 5.3 if $\iota _{i}{\hat{\gamma }} = 1$. However, for ease of presentation, we herein consider the case where these parameters are determined by dimension and thus have 5.2 in the fractures and 5.3 in the bulk and its surfaces.

Our third equation concerns the mass balance. To capture volumetric change in terms of our four variables, we first use the decomposition induced by ${\check{\gamma }}$ and ${\hat{\gamma }}$ to rewrite

$$\begin{aligned} {\mathfrak {e}} = {\hat{\gamma }}{\mathfrak {e}}+{\check{\gamma }}{\mathfrak {e}} = \hat{{\mathfrak {A}}}(\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}{\mathfrak {p}})+\partial _{t}^{-1}\check{\gamma }{\mathfrak {D}}_{s}{\mathfrak {v}}. \end{aligned}$$

Taking the derivative of (4.11f) with respect to time and substituting this equality, we have

$$\begin{aligned} \partial _{t}{\mathfrak {m}} = \partial _{t}\hat{{\mathfrak {T}}}'\alpha \hat{{\mathfrak {A}}}(\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}{\mathfrak {p}})+\check{{\mathfrak {T}}}'\alpha {\mathfrak {D}}_{s}{\mathfrak {v}}+\partial _{t}\beta {\mathfrak {p}}. \end{aligned}$$

Inserting this in (4.11b), the mass balance equation becomes

$$\begin{aligned} \partial _{t}\hat{{\mathfrak {T}}}'\alpha \hat{{\mathfrak {A}}}(\tilde{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}{\mathfrak {p}})+\check{{\mathfrak {T}}}'\alpha {\mathfrak {D}}_{s}{\mathfrak {v}}+\partial _{t}\beta {\mathfrak {p}}+\mathfrak {D \cdot q} = {\mathfrak {r}}_{{\mathfrak {m}}}. \end{aligned}$$

(5.5)

Here, we assume that

Assumption 5.4

$\beta$ is a coercive, linear operator with coercivity constant $c_{\beta }>0$.

Finally, in (4.11e), i.e., Darcy’s law, we neglect the dependency of the permeability on the strain ${\mathfrak {e}}$:

$$\begin{aligned} {\mathfrak {q}} = \kappa ({\mathfrak {q}})(\mathbb {-D}{\mathfrak {p}}+{\mathfrak {r}}_{{\mathfrak {g}}}). \end{aligned}$$

This is equivalent to stating that $\kappa ^{-1}({\mathfrak {q}})\mathfrak {q = }-{\mathbb {D}}{\mathfrak {p}}+{\mathfrak {r}}_{{\mathfrak {g}}}$. This relation only contains two variables (due to the negligence of strain dependencies) and we can consider this law, which is possibly nonlinear, as a binary relation $\kappa ^{-1}$. Observing that the dependency on ${\mathfrak {q}}$ is implied in the notation, we arrive at the binary relation

$$\begin{aligned} (\mathfrak {q, }-{\mathbb {D}}{\mathfrak {p}}+{\mathfrak {r}}_{{\mathfrak {g}}})\in \kappa ^{-1}. \end{aligned}$$

(5.6)

Again, we make an assumption on the relation $\kappa ^{-1}$, namely that

Assumption 5.5

$\kappa ^{-1}$ is bounded and $c_{\kappa }$-maximal monotone for some $c_{\kappa }>0$. Moreover, $(0,0)\in \kappa ^{-1}$.

With these simplifications and assumptions in place, we have a system of four equations with four unknowns.

5.3 Weak formulation

To accommodate analysis of system (5.7), we next present the weak formulation of the poromechanics problem. The first step is to introduce the relevant function space on which to pose the problem. Following the observations from Sect. 2.3, we consider the following four function spaces for the variables:

$$\begin{aligned} {\mathfrak {v}} \in {\mathfrak {U}}&:=L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n}),&{\mathfrak {p}} \in {\mathfrak {P}}&:=L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}}), \end{aligned}$$

(5.8a)

$$\begin{aligned} \tilde{{\mathfrak {s}}}\in {\mathfrak {G}}&\subset L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n}),&{\mathfrak {q}} \in {\mathfrak {Q}}&:=L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}}). \end{aligned}$$

(5.8b)

Recall that ${\mathfrak {G}}$, as defined in Sect. 5.1, contains symmetry properties for stress and strain. Together, these spaces form the composite space $U$:

$$\begin{aligned} U :=\mathfrak {U\times P\times G\times Q}. \end{aligned}$$

(5.9)

The space $U$ is naturally endowed with a $L^{2}$-type inner product and norm, given by

$$\begin{aligned} \begin{aligned} \left\langle u_{1},u_{2}\right\rangle _{U}&=\left\langle \begin{bmatrix} {\mathfrak {v}}_{1}&{\mathfrak {p}}_{1}&\tilde{{\mathfrak {s}}}_{1}&{\mathfrak {q}}_{1} \end{bmatrix}^T,\begin{bmatrix} {\mathfrak {v}}_{2}&{\mathfrak {p}}_{2}&\tilde{{\mathfrak {s}}}_{2}&{\mathfrak {q}}_{2} \end{bmatrix}^T \right\rangle _{U} \\&:=\left\langle {\mathfrak {v}}_{1},{\mathfrak {v}}_{2}\right\rangle _{{\mathfrak {X}}^{0}}+\left\langle {\mathfrak {p}}_{1},{\mathfrak {p}}_{2}\right\rangle _{{\mathfrak {X}}^{n}}+\left\langle \tilde{{\mathfrak {s}}}_{1},\tilde{{\mathfrak {s}}}_{2}\right\rangle _{{\mathfrak {X}}^{n-1}}+\left\langle {\mathfrak {q}}_{1},{\mathfrak {q}}_{2}\right\rangle _{{\mathfrak {X}}^{n-1}}, \\ \left\| u\right\| _{U}&:=\sqrt{ \left\langle u,u\right\rangle _U }. \end{aligned} \end{aligned}$$

With the function spaces defined, we continue by considering all operators in the system as binary relations, including the linear operators. For example, we write

$$\begin{aligned} \rho _{r}\subseteq L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\times L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n}) \end{aligned}$$

also in the case that $\rho _{r}$ is simply multiplication by a constant (see Example A.2). To further emphasize this, all mixed-dimensional differential operators (see Sect. 2.4) are also interpreted as binary relations:

$$\begin{aligned} \begin{aligned} \mathring{{\mathfrak {D}}}_{s} \subset L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})\times L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n}),&~~~~~~({\mathbb {D}}_{s} \cdot ) \subset L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})\times L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n}), \\ (\mathring{{\mathfrak {D}}} \cdot ) \subset L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}})\times L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}}),&~~~~~~\mathbb {D \subset }L^{2}({\mathfrak {X}}^{n}, {\mathbb {R}})\times L^{2}({\mathfrak {X}}^{n-1}, {\mathbb {R}}). \\ \end{aligned} \end{aligned}$$

Recall that the domains of these differential operators are proper, dense subsets of the $L^{2}$ spaces, e.g., $\textrm{dom}(\mathring{{\mathfrak {D}}}_{s})\subset L^{2}({\mathfrak {X}}^{0},{\mathbb {R}}^{n})$. In turn, by searching the solution in these domains, we ensure that the solution has sufficient regularity for the corresponding differentials to be well-defined. This same argument enforces the boundary conditions on the variables. We note that our choice of enforcing boundary conditions on the $\mathring{{\mathfrak {D}}}$-type differential operators implies zero (clamped) conditions on the displacement, and similarly zero normal component (no-flow) conditions on the fluid flux. The development below would be equally valid with boundary conditions imposed via $\mathring{{\mathbb {D}}}$-type operators, which would correspond to zero normal stress (floating) conditions for mechanics, and zero fluid pressure (open) conditions for flow.

Finally, we incorporate the time dependency. Following [10], we introduce the exponentially weighted Bochner space $L_{\nu }^{2}({\mathbb {R}},U)$ as follows.

Definition 5.2

Given $\nu >0$, let $L_{\nu }^{2}({\mathbb {R}},U) :=\left\{ f: {\mathbb {R}} \rightarrow U \mid \int _{{\mathbb {R}}}^{}\left\| e^{-\nu t}f(t)\right\| _{U}^{2}\textrm{d}t<\infty \right\}$.

Let $\Vert u \Vert _{L_{\nu }^{2}({\mathbb {R}},U)}:= \int _{{\mathbb {R}}}^{}\left\| e^{-\nu t}f(t)\right\| _{U}^{2}\textrm{d}t$ be the endowed norm. The weight with positive $\nu$ ensures that causality is preserved. The time derivative is then introduced as an operator acting on this weighted space.

Definition 5.3

Given $\nu >0$, let $\partial _{0,\nu }:\textrm{dom}(\partial _{0,\nu })\subseteq L_{\nu }^{2}({\mathbb {R}},H)\rightarrow L_{\nu }^{2}({\mathbb {R}},H)$ be given by

$$\begin{aligned} \partial _{0,\nu } :=e^{\nu t}(\partial _{t}+\nu )e^{-\nu t}. \end{aligned}$$

The motivation behind this definition can be found in Definition A.14. With the function space and interpretation of operators in place, we arrive at the weak formulation (5.10) of the simplified hydromechanical problem (5.7).

Remark 5.2

We refer to (5.10) as the weak formulation since the problem is posed in a Hilbert space setting (the domains of the differential operators are Hilbert spaces [29, 30]). As a direct consequence, we recall that the solution, if it exists, is defined up to the equivalence classes of $L_{\nu }^{2}({\mathbb {R}},U)$.

Remark 5.3

The function space $L_{\nu }^{2}({\mathbb {R}},U)$ does not ensure more regularity than square integrability in space and (weighted) time. However, the presence of the differential operators ensures that the solution, if it exists, has sufficient regularity for these to be well-defined. For example, the term $\partial _{0,\nu }{\mathfrak {v}}$ ensures that ${\mathfrak {v}}\in \textrm{dom}(\partial _{0,\nu })$ and thus $\partial _{0,\nu }{\mathfrak {v}}\in L_{\nu }^{2}({\mathbb {R}},{\mathfrak {U}})$. Similarly, the solution has ${\mathfrak {v}}\in \textrm{dom}(\mathring{{\mathfrak {D}}}_{s})$ and thus $\mathring{{\mathfrak {D}}}_{s}{\mathfrak {v}}\in L^{2}({\mathfrak {X}}^{1},{\mathbb {R}}^{n})$.

Remark 5.4

This formulation does not allow for a constant effect of gravity throughout the past since ${\mathfrak {r}}_{{\mathfrak {g}}}$ is assumed to be in $L_{\nu }^{2}({\mathbb {R}},{\mathfrak {Q}})$. However, this effect can be properly incorporated by instead considering an initial-value problem on the real half-line ${\mathbb {R}}_{>0}$, see [12] and Remark 3.3 in [10].

5.4 Well-posedness of the weak formulation

In order to analyze problem (5.10) in the appropriate setting, we recognize that the binary relation has a favorable underlying structure. In particular, we recognize that problem (5.10) is an evolutionary equation [10] of the form

$$\begin{aligned} ( u, f)\in \partial _{0,\nu }M_{0}+M_{1}+A_{\nu }. \end{aligned}$$

Here, the first two components $M_{0}$ and $M_{1}$ are linear operators, given by

$$\begin{aligned} M_{0}&:=\Sigma '\begin{bmatrix} \rho _{r} &{}\quad &{}\quad &{}\quad \\ &{}\quad \beta &{}\quad &{}\quad \\ &{}\quad &{}\quad \hat{{\mathfrak {A}}} &{}\quad \\ &{}\quad &{}\quad &{}\quad 0 \\ \end{bmatrix} \Sigma ,&\quad \Sigma&\quad :=\begin{bmatrix} 1 &{}\quad &{}\quad &{}\quad \\ &{}\quad 1 &{}\quad &{}\quad \\ &{}\quad \alpha \hat{{\mathfrak {T}}} &{} \quad 1 &{}\quad \\ &{}\quad &{}\quad &{}\quad 0 \\ \end{bmatrix},&M_{1}&:=\begin{bmatrix} 0 &{}\quad &{} \quad &{}\quad \\ &{}\quad 0 &{}\quad &{}\quad \\ &{}\quad &{}\quad \check{c} &{} \quad \\ &{}\quad &{}\quad &{}\quad c_{\kappa } \\ \end{bmatrix}. \end{aligned}$$

(5.11)

Furthermore, $A_{\nu }\subseteq L_{\nu }^{2}({\mathbb {R}},U)\times L_{\nu }^{2}({\mathbb {R}},U)$ is a temporal extension (see Definition A.11) given by

$$\begin{aligned} A_{\nu } :=\left\{ u,v\in L_{\nu }^{2}({\mathbb {R}},U)\mid (u(t),v(t))\in A_{0}+A_{1}, \text { for a.e. } t\in {\mathbb {R}}\right\} . \end{aligned}$$

Where the spatial relations $A_{0}$ and $A_{1}$ are given by

$$\begin{aligned} A_{0}&:=\begin{bmatrix} 0 &{}\quad &{}\quad &{}\quad \\ &{}\quad 0 &{}\quad &{}\quad \\ &{}\quad &{}\quad \check{{\mathfrak {A}}}-\check{c} &{} \quad \\ &{} \quad &{} \quad &{} \quad \kappa ^{-1}-c_{k} \end{bmatrix},&A_{1}&:=\begin{bmatrix} 0 &{}\quad {\mathbb {D}}_{s} \cdot \alpha \check{{\mathfrak {T}}} &{}\quad -{\mathbb {D}}_{s} \cdot &{} \quad \\ \check{{\mathfrak {T}}}'\alpha \mathring{{\mathfrak {D}}}_{s} &{}\quad 0 &{} \quad &{}\quad \mathring{{\mathfrak {D}}} \cdot \\ -\mathring{{\mathfrak {D}}}_{s} &{}\quad &{} \quad 0 &{}\quad \\ &{}\quad {\mathbb {D}} &{} \quad &{}\quad 0 \end{bmatrix}. \end{aligned}$$

(5.12)

The decomposition of the evolutionary equation in terms of the binary relations $M_{0}$, $M_{1}$ and $A_{0}$, $A_{1}$ highlights the structure of the problem: $M_{0}$ and $M_{1}$ contain the weights of the time-derivative and diagonal terms, respectively, while $A_{0}$ and $A_{1}$ contain the nonlinearities and differential operators, respectively.

This identification of the problem allows us to use the solution theory of evolutionary equations, in particular we recall the following key theorem.

Theorem 5.1

(Well-posedness of autonomous evolutionary inclusions) Let $\nu >0$ and $r>\frac{1}{2\nu }$. Let $A_{\nu }\subseteq L_{\nu }^{2}({\mathbb {R}},U)\times L_{\nu }^{2}({\mathbb {R}},U)$ be a binary relation and $M_{0},M_{1}\subseteq U\times U$ linear, bounded mappings. Assume the following hypotheses:

• H1 $A_{\nu }$ is maximal monotone, time-translation invariant (autonomous), and satisfies

  $$\begin{aligned} \int _{-\infty }^{0}\text {Re}(\left\langle u_{1}(t)-u_{2}(t),v_{1}(t)-v_{2}(t)\right\rangle )e^{-2\nu t}\textrm{d}t&\ge 0,&\forall (u_{1},v_{1}),(u_{2},v_{2})&\in A_{\nu }. \end{aligned}$$

• H2 $\exists c>0$ such that $z^{-1}M_{0}+M_{1}-c$ is monotone for all $z\in {\mathcal {B}}_{{\mathbb {C}}}(r,r)$. Here, ${\mathcal {B}}_{{\mathbb {C}}}(r,r)$ denotes the open complex ball with radius $r$, centered at $r$.

Then for each $f\in L_{\nu }^{2}({\mathbb {R}},U)$, there exists a unique $u\in L_{\nu }^{2}({\mathbb {R}},U)$ such that

$$\begin{aligned} (u,f)\in \overline{\partial _{0,\nu }M_{0}+M_{1}+A_{\nu }}. \end{aligned}$$

Moreover, the solution operator $(\overline{\partial _{0,\nu }M_{0}+M_{1}+A_{\nu }})^{-1}$ is causal and Lipschitz-continuous with a Lipschitz constant bounded by $\frac{1}{c}$.

Proof

See Theorem 3.2 of [10]. $\square$

This theorem leads to the statement of our main well-posedness result.

Theorem 5.2

(The main result) If assumptions 5.1–5.5 are fulfilled, then the mixed-dimensional poromechanics problem (5.10) is well-posed. In particular, for any right-hand side $f\in L_{\nu }^{2}({\mathbb {R}},U)$, a solution $u\in L_{\nu }^{2}({\mathbb {R}},U)$ exists uniquely such that $(u,f)\in \overline{\partial _{0,\nu }M_{0}+M_{1}+A_{\nu }}$. Moreover, a $c>0$ exists such that the solution operator is bounded:

$$\begin{aligned} \left\| u\right\| _{L_{\nu }^{2}({\mathbb {R}},U)}\le \frac{1}{c}\left\| f\right\| _{L_{\nu }^{2}({\mathbb {R}},U)}. \end{aligned}$$

Proof

Lemmas 5.2 and 5.3, presented below, suffice to invoke Theorem 5.1. The bound is a direct consequence of the Lipschitz continuity of the solution operator. $\square$

The proof of Theorem 5.2 requires validating maximal monotonicity of several operators, for which we often need to take a difference between two elements $u_{1},u_{2}\in U$. As a short-hand notation, we denote

$$\begin{aligned} \delta u :=u_{1}-u_{2} \end{aligned}$$

and let $\left[ \delta \mathfrak {v,}\delta \mathfrak {p,}\delta \hat{{\mathfrak {s}}},\delta {\mathfrak {q}}\right] ^{T} :=\delta u$. Moreover, we omit the subscripts on inner products and norms for notational brevity.

Lemma 5.2

(H1) If Assumptions 5.2, 5.5 are satisfied, then [H1] is fulfilled.

Proof

First, $A_{0}$ is maximal monotone by 5.2 and 5.5. Second, $A_{1}$ is linear and skew self-adjoint and therefore also maximal monotone. Furthermore, $0$ is in the domain of both operators and $A_{0}$ is bounded due to 5.2, 5.5. We then invoke Lemma 7.1 from the appendix to conclude that the sum $A_{0}+A_{1}$ is maximal monotone.

Next, we note that $A_{\nu }$ is the temporal extension of $A_{0}+A_{1}$ and we have that $(0,0)\in A_{0}+A_{1}$. Proposition 2.5 of [12] then ensures that $A_{\nu }$ is maximal monotone.

Time-translation invariance follows directly from the definition of $A_{\nu }$. The positivity of the time integral is clear by the monotonicity of $A_{0}+A_{1}$. $\square$

Lemma 5.3

(H2) If Assumptions 5.1–5.5 are satisfied, then [H2] is fulfilled.

Proof

Let $(u_{1},v_{1}),(u_{2},v_{2})\in z^{-1}M_{0}+M_{1}$. The assumptions and Lemmas 5.4 and 5.5, presented below, imply that

$$\begin{aligned} \begin{aligned} \left\langle \delta u,\delta M_{0}u\right\rangle&\ge c_{0}(\left\| \delta {\mathfrak {v}}\right\| ^{2}+\left\| \delta {\mathfrak {p}}\right\| ^{2}+\left\| {\hat{\gamma }}\delta \tilde{{\mathfrak {s}}}\right\| ^{2}), \\ \left\langle \delta u,\delta M_{1}u\right\rangle&\ge c_{1}(\left\| \delta {\mathfrak {q}}\right\| ^{2}+\left\| {\check{\gamma }}\delta \tilde{{\mathfrak {s}}}\right\| ^{2}), \\ \end{aligned} \end{aligned}$$

for some $c_{0},c_{1}>0$. Using these bounds, we derive

$$\begin{aligned} \begin{aligned} \text {Re}(\left\langle \delta u,\delta v\right\rangle )&= \text {Re}(z^{-1}\left\langle \delta u, \delta M_{0}u\right\rangle +\left\langle \delta u,\delta M_{1}u\right\rangle ) \\&\ge \text {Re}(z^{-1})c_{0}(\left\| \delta {\mathfrak {v}}\right\| ^{2}+\left\| \delta {\mathfrak {p}}\right\| ^{2}+\left\| \hat{\gamma }\delta \tilde{{\mathfrak {s}}}\right\| ^{2})+c_{1}(\left\| \delta {\mathfrak {q}}\right\| ^{2}+\left\| \check{\gamma }\delta \tilde{{\mathfrak {s}}}\right\| ^{2}). \end{aligned} \end{aligned}$$

Since $z\in {\mathcal {B}}_{{\mathbb {C}}}(r,r)$, its real part satisfies $\text {Re}(z^{-1})\ge \frac{1}{2r}$ and thus we obtain:

$$\begin{aligned} \text {Re}(\left\langle \delta u,\delta v\right\rangle )\ge \frac{c_{0}}{2r}(\left\| \delta {\mathfrak {v}}\right\| ^{2}+\left\| \delta {\mathfrak {p}}\right\| ^{2}+\left\| \hat{\gamma }\delta \tilde{{\mathfrak {s}}}\right\| ^{2})+c_{1}(\left\| \delta {\mathfrak {q}}\right\| ^{2}+\left\| \check{\gamma }\delta \tilde{{\mathfrak {s}}}\right\| ^{2}). \end{aligned}$$

We then set $c = \min \left\{ \frac{c_{0}}{2r},c_{1}\right\}$ to conclude that

$$\begin{aligned} \text {Re}(\left\langle \delta u,\delta v\right\rangle )\ge c\left\| \delta u\right\| ^{2}. \end{aligned}$$

$\square$

Lemma 5.4

Given Assumptions 5.1, 5.3, and 5.4, then the relationship $M_{0}$ satisfies

$$\begin{aligned} \left\langle \delta u,\delta M_{0}u\right\rangle \ge c_{0}(\left\| \delta {\mathfrak {v}}\right\| ^{2}+\left\| \delta {\mathfrak {p}}\right\| ^{2}+\left\| \hat{\gamma }\delta \tilde{{\mathfrak {s}}}\right\| ^{2}), \end{aligned}$$

for some $c_0 >0$.

Proof

Assumptions 5.1,5.3 and 5.4 give us that

$$\begin{aligned} \left\langle \delta u,\delta M_{0}u\right\rangle \ge c_{\rho }\left\| \delta {\mathfrak {v}}\right\| ^{2}+c_{\beta }\left\| \delta {\mathfrak {p}}\right\| ^{2}+{\hat{c}}\left\| \hat{\gamma }(\delta \tilde{{\mathfrak {s}}}+\alpha {\mathfrak {T}}\delta {\mathfrak {p}})\right\| ^{2}. \end{aligned}$$

Next, we use the continuity of $\alpha$ to obtain

$$\begin{aligned} \begin{aligned} c_{\beta }\left\| \delta {\mathfrak {p}}\right\| ^{2}+{\hat{c}}\left\| {\hat{\gamma }}\delta (\tilde{{\mathfrak {s}}}+\alpha {\mathfrak {T}}\delta {\mathfrak {p}})\right\| ^{2}&\gtrsim \frac{1}{2}\left\| \delta {\mathfrak {p}}\right\| ^{2}+\frac{1}{2}\left\| -\delta \alpha {\mathfrak {p}}\right\| ^{2}+\left\| {\hat{\gamma }}\delta (\tilde{{\mathfrak {s}}}+\alpha {\mathfrak {T}}\delta {\mathfrak {p}})\right\| ^{2} \\&\gtrsim \left\| \delta {\mathfrak {p}}\right\| ^{2}+\left\| {\hat{\gamma }}\delta \tilde{{\mathfrak {s}}}\right\| ^{2}, \\ \end{aligned} \end{aligned}$$

in which $a\gtrsim b$ implies that a $c>0$ exists such that $a\ge cb$. The combination of these bounds now proves the result. $\square$

Lemma 5.5

Given Assumptions 5.2, 5.5, then $M_{1}$ satisfies

$$\begin{aligned} \left\langle \delta u,\delta M_{1}u\right\rangle \ge c_{1}(\left\| \delta {\mathfrak {q}}\right\| ^{2}+\left\| \check{\gamma }\delta \tilde{{\mathfrak {s}}}\right\| ^{2}), \end{aligned}$$

for some $c_{1}>0$.

Proof

Substituting the definitions and using assumptions 5.3, 5.5, we have

$$\begin{aligned} \left\langle \delta u,\delta M_{1}u\right\rangle = c_{\kappa }\left\| \delta {\mathfrak {q}}\right\| ^{2}+\check{c}\left\| \check{\gamma }\delta {\mathfrak {s}}\right\| ^{2}. \end{aligned}$$

The result therefore follows with $c_{1} = \min \left\{ c_{\kappa },\check{c}\right\}$. $\square$

5.5 Degeneracies: Maximal monotone contact relations with bounded inverse

For concrete applications, it may be desirable to relax some of assumptions 5.1–5.5. Such relaxations will often correspond to sending a physical parameter to zero, and can therefore be considered as degenerate limits of the base model (5.10) as analyzed in Sect. 5.4. As a general expectation, these limit models need to be analyzed on a case-by-case basis, as they will in principle imply that the Lipschitz constant in Theorem 5.2 is not bounded. To illustrate the implications of such degenerate limits, and to show how they can be treated within the theory as presented above, we include an analysis of maximal monotone contact relations with bounded inverse.

We focus on the assumption 5.2 regarding the frictional contact law $\check{{\mathfrak {A}}}$. This may be too restrictive for some conventional friction relations (we return to this in Sect. 5.6), in particular due to the bound on the constant $\check{c}>0$ and the assumed boundedness of the relation. In this section, we consider a more general class of models, namely those concerning maximal monotone contact relations with bounded inverse. To be concrete, we adopt the following relaxation of assumption 5.2:

Assumption 5.6

$\check{{\mathfrak {A}}}\mathfrak {\subseteq G\times G}$ is a maximal monotone relation with $(0,0)\in \check{{\mathfrak {A}}}$. Moreover, $\check{{\mathfrak {A}}}^{-1}$ is bounded.

A key generalization from 5.2 is that we now allow for the case $\check{c} = 0$, and as such the treatment below is a consideration of a degenerate limit of the main model equations from Sect. 5.2. The penalty we pay for allowing this degeneracy is that we lose control over the $L^{2}$-norm of the components $\iota _{j}{\mathfrak {s}}$ wherever $\iota _{j}{\check{\gamma }} = 1$. As a result, the problem now needs to be posed in a smaller space, which we define by

$$\begin{aligned} \hat{{\mathfrak {G}}} :={\hat{\gamma }}{\mathfrak {G}}. \end{aligned}$$

(5.13)

We remark that for $\hat{{\mathfrak {s}}}\in \hat{{\mathfrak {G}}}$, we have $\iota _{j}\hat{{\mathfrak {s}}} = 0$ for $\iota _{j}{\hat{\gamma }} = 0$. This space is therefore isomorphic to a restriction of ${\mathfrak {G}}$, but our definition (5.13) avoids the need for explicitly introducing restriction and inclusion operators. The endowed norm is

$$\begin{aligned} \left\| \hat{{\mathfrak {s}}}\right\| _{\hat{{\mathfrak {G}}}}^{2} :=\left\| \hat{{\mathfrak {s}}}\right\| _{{\mathfrak {X}}^{1}}^{2}&= \sum _{ \begin{array}{c} j\in {\mathfrak {F}}^{1}, \\ \iota _{j}{\hat{\gamma }} = 1 \end{array}} \left\| \iota _{j}\hat{{\mathfrak {s}}}\right\| _{X_{j}}^{2},&\forall \hat{{\mathfrak {s}}}&\in \hat{{\mathfrak {G}}}. \end{aligned}$$

(5.14)

Importantly, (5.14) forms a natural choice for this setting since it does not contain any terms from the fractures, i.e., on $X_{i}$ for $i\in I^{n-1}$ (and its descendants). The stress variable can be decomposed in a similar manner as the stress–strain relationship, and we therefore define

$$\begin{aligned} {\mathfrak {s}} = \hat{{\mathfrak {s}}}+\check{{\mathfrak {s}}} :={\hat{\gamma }}{\mathfrak {s}}+{\check{\gamma }}{\mathfrak {s}}. \end{aligned}$$

We proceed by eliminating $\check{{\mathfrak {s}}}$ from the system. For that, we note that since $\check{\gamma }{\mathfrak {D}}_{s}{\mathfrak {v}} = \mathbb {d}{\mathfrak {v}}$, the stress–strain relationship $(\check{{\mathfrak {s}}},\check{\gamma }{\mathfrak {D}}_{s}{\mathfrak {v}})\in \check{{\mathfrak {A}}}$ can be restated as

$$\begin{aligned} (\check{{\mathfrak {s}}},\mathbb {d}{\mathfrak {v}})\in \check{{\mathfrak {A}}}\Leftrightarrow (\mathbb {d}\mathfrak {v,}\check{{\mathfrak {s}}})\in \check{{\mathfrak {A}}}^{-1}\Leftrightarrow (\mathfrak {v,}\check{{\mathfrak {s}}})\in \check{{\mathfrak {A}}}^{-1}\mathbb {d}. \end{aligned}$$

The momentum balance equation, i.e., the first row of (5.10), is rewritten using this substitution as

$$\begin{aligned} \left( \begin{bmatrix} {\mathfrak {v}} \\ {\mathfrak {p}} \\ \hat{{\mathfrak {s}}} \\ {\mathfrak {q}} \\ \end{bmatrix},{\mathfrak {r}}_{{\mathfrak {s}}}\right) \in \begin{bmatrix} \rho _{r}\partial _{0,\nu }+\mathbb {d}'\check{{\mathfrak {A}}}^{-1}\mathbb {d}&{\mathbb {D}}_{s} \cdot \alpha \check{{\mathfrak {T}}}&-{\mathbb {D}}_{s} \cdot&0 \end{bmatrix} \end{aligned}$$

with $\mathbb {d'}$ the adjoint of the jump operator $\mathbb {d}$. On the other hand, the stress–strain relationship in the bulk is now given by

$$\begin{aligned} (\hat{{\mathfrak {s}}}+\alpha \hat{{\mathfrak {T}}}\mathfrak {p, }{\hat{\gamma }}{\mathfrak {D}}_{s}{\mathfrak {v}})\in \partial _{t}\hat{{\mathfrak {A}}}. \end{aligned}$$

Since the stress component $\check{{\mathfrak {s}}}$ does not influence the remaining equations, we are ready to pose the new problem (5.15) in the composite space

$$\begin{aligned} {\hat{U}}\mathfrak { :=U\times P\times }\hat{{\mathfrak {G}}}\mathfrak {\times Q,} \end{aligned}$$

endowed with the norm $\left\| \cdot \right\| _{{\hat{U}}} =\left\| \cdot \right\| _{U}$.

Again, we recognize the structure of problem (5.15) as an evolutionary equation and we note that the decomposition $\partial _{0,\nu }M_{0}+M_{1}+A_{\nu }$ now holds with $M_{0}$ unchanged and $M_{1}$ given by

$$\begin{aligned} M_{1} :=\begin{bmatrix} 0 &{}\quad &{} \quad &{}\quad \\ &{} \quad 0 &{}\quad &{}\quad \\ &{} \quad &{}\quad 0 &{}\quad \\ &{}\quad &{} \quad &{}\quad c_{\kappa } \\ \end{bmatrix}. \end{aligned}$$

(5.16)

On the other hand, $A_{\nu }$ is here the temporal extension of $A_{0}+A_{1}$ with

$$\begin{aligned} A_{0} :=\begin{bmatrix} \mathbb {d}'\check{{\mathfrak {A}}}^{-1}\mathbb {d} &{} \quad &{} \quad &{} \quad \\ &{}\quad 0 &{} \quad &{} \quad \\ &{} \quad &{} \quad 0 &{} \quad \\ &{} \quad &{} \quad &{} \quad \kappa ^{-1}-c_{k} \\ \end{bmatrix}, A_{1} :=\begin{bmatrix} 0 &{}\quad {\mathbb {D}}_{s} \cdot \alpha \check{{\mathfrak {T}}} &{}\quad -{\mathbb {D}}_{s} \cdot &{}\quad \\ \check{{\mathfrak {T}}}'\alpha {\mathfrak {D}}_{s} &{}\quad 0 &{} \quad &{}\quad {\mathfrak {D}} \cdot \\ -{\hat{\gamma }}{\mathfrak {D}}_{s} &{}\quad &{}\quad 0 &{} \quad \\ &{}\quad {\mathbb {D}} &{}\quad &{}\quad 0 \\ \end{bmatrix}. \end{aligned}$$

Theorem 5.3

If Assumptions 5.1 and 5.3–5.6 are fulfilled, then the mixed-dimensional poromechanics problem with maximal monotone contact relations (5.15) is well-posed.

Proof

Lemmas 5.6 and 5.7, presented below, show that hypotheses [H1]–[H2] are fulfilled for (5.15). Theorem 5.1 then provides the result. $\square$

Lemma 5.6

(H1) If Assumptions 5.5, 5.6 are satisfied, then [H1] is fulfilled.

Proof

The fact that $A_{0}$ is maximal monotone and bounded follows from 5.5, 5.6. The skew self-adjointness of $A_{1}$ is verified by

$$\begin{aligned} \left\langle -\hat{\gamma }{\mathfrak {D}}_{s}\mathfrak {v,}\hat{{\mathfrak {s}}}\right\rangle +\left\langle \mathfrak {v,}-{\mathbb {D}}_{s} \cdot \hat{{\mathfrak {s}}}\right\rangle =\left\langle -{\mathfrak {D}}_{s}\mathfrak {v,}\hat{{\mathfrak {s}}}\right\rangle +\left\langle {\mathfrak {D}}_{s}\mathfrak {v,}\hat{{\mathfrak {s}}}\right\rangle = 0. \end{aligned}$$

Now, the arguments from Lemma 5.2 provide the result. $\square$

Lemma 5.7

[H2] If Assumptions 5.1 and 5.3–5.5 are satisfied, then [H2] is fulfilled.

Proof

Since $M_{0}$ has remained unchanged, Lemma 5.4 gives us the bound

$$\begin{aligned} \left\langle \delta u,\delta M_{0}u\right\rangle \ge c_{0}(\left\| \delta {\mathfrak {v}}\right\| ^{2}+\left\| \delta {\mathfrak {p}}\right\| ^{2}+\left\| \delta \hat{{\mathfrak {s}}}\right\| ^{2}) \end{aligned}$$

for some $c_{0}>0.$ Continuing with $M_{1}$, its definition (5.16) directly gives us

$$\begin{aligned} \left\langle \delta u,\delta M_{1}u\right\rangle = c_{\kappa }\left\| \delta {\mathfrak {q}}\right\| ^{2}. \end{aligned}$$

Using these bounds, the arguments from Lemma 5.3 are followed to conclude the proof. $\square$

5.6 Exemplary models

We finalize this section by describing two example models. The two models differ on whether assumption 5.2 or 5.6 is satisfied. Since the only difference lies in the frictional contact law $\check{{\mathfrak {C}}}$, we first present the other components of the model.

Example 5.1

• Let the bulk density be given by a scalar $c_{\rho }>0$. Then, the corresponding binary relation, defined by multiplication with this scalar, is given by

  $$\begin{aligned} \rho _{r} :=\left\{ ({\mathfrak {v}}_{1},{\mathfrak {v}}_{2})\in L_{\nu }({\mathbb {R}},{\mathfrak {U}})\times L_{\nu }({\mathbb {R}},{\mathfrak {U}})\mid {\mathfrak {v}}_{2} = c_{\rho }{\mathfrak {v}}_{1}\right\} \\ \end{aligned}$$

  and satisfies 5.1.

• Let the bulk medium and its boundaries be isotropic materials. Then, the stress–strain relationship $\hat{{\mathfrak {A}}}$ can be described, using the Lamé parameters $\mu ,\lambda$, as

  $$\begin{aligned} \hat{{\mathfrak {A}}} :=\left\{ (\mathfrak {s,e})\in L_{\nu }({\mathbb {R}},{\mathfrak {G}})\times L_{\nu }({\mathbb {R}},{\mathfrak {G}})\mid \iota _{j}{\mathfrak {s}} = 2\mu \iota _{j}{\mathfrak {e}}+\lambda (I:\iota _{j}{\mathfrak {e}}_{\parallel })I, \forall j\in {\mathfrak {S}}_{i}, i\in I^{n}\right\} . \\ \end{aligned}$$

  Here, $I$ is the identity tensor in ${\mathbb {R}}^{d_j}$. It is easy to see that $\hat{{\mathfrak {A}}}$ satisfies 5.3 with ${\hat{c}} =(2\mu +n\lambda )^{-1}$.

• Similar to 5.1, we set $\beta$ as multiplication with $c_{\beta }>0$:

  $$\begin{aligned} \beta :=\left\{ ({\mathfrak {p}}_{1},{\mathfrak {p}}_{2})\in L_{\nu }({\mathbb {R}},{\mathfrak {P}})\times L_{\nu }({\mathbb {R}},{\mathfrak {P}})\mid {\mathfrak {p}}_{2} = c_{\beta }{\mathfrak {p}}_{1}\right\} \\ \end{aligned}$$

  which satisfies 5.4.

• Finally, the constitutive law relating flux and pressure is assumed to be given by Darcy(-Forchheimer) flow.

  $$\begin{aligned} \kappa ^{-1} :=\left\{ (\mathfrak {q,}-{\mathbb {D}}{\mathfrak {p}})\in L_{\nu }({\mathbb {R}},{\mathfrak {Q}})\times L_{\nu }({\mathbb {R}},{\mathfrak {Q}})\mid \mathbb {-D}{\mathfrak {p}} = (\kappa _{1}+\kappa _{2}\left| {\mathfrak {q}}\right| ){\mathfrak {q}}\right\} . \end{aligned}$$

  With the material parameters $\kappa _{1}$,$\kappa _{2}>0$, we have $c_{\kappa }>0$ in 5.5.

Next, we assume that the frictional contact law $\check{{\mathfrak {A}}}$ is defined as the direct sum of disjoint relations:

$$\begin{aligned} \check{{\mathfrak {A}}} :=\bigoplus _{ \begin{array}{c} i\in I^{n-1} \\ j\in {\mathfrak {S}}_{i}\cap {\mathfrak {F}}^{1} \end{array}} A_{j,\parallel }\oplus A_{j,\perp } \end{aligned}$$

(5.17)

Here, the subscript $\parallel$ indicates the tangential friction law and $\perp$ the perpendicular contact. It is clear that if each $A_{j.\parallel }$ and $A_{j.\perp }$ satisfies 5.2, respectively 5.6, then so does $\check{{\mathfrak {A}}}$. We dedicate the following two subsections to the description of two exemplary models that fulfill these respective assumptions.

5.6.1 Maximal monotone contact relations with bounded inverse

We start by considering relations that satisfy 5.6, since this assumption is easier to fulfil in practice. For the stress–strain relation in the parallel direction, let us consider Tresca friction. For a given threshold $\tau >0$, this law is given by

$$\begin{aligned} A_{j,\parallel } :=\left\{ (\sigma ,{\dot{\varepsilon }})\in \iota _{j,\parallel }\mathfrak {G\times }\iota _{j,\parallel }{\mathfrak {G}}\left| \begin{array}{c} \tau {\dot{\varepsilon }} =\left| {\dot{\varepsilon }}\right| \sigma \text { or }{\dot{\varepsilon }} = 0 \\ \left| \sigma \right| \le \tau \end{array} \text { a.e. on } X_{j} \right. \right\} \end{aligned}$$

(5.18)

Lemma 5.8

Tresca friction $A_{j,\parallel }$ of (5.18) satisfies 5.6, i.e., is maximal monotone, has $(0,0)\in A_{j,\parallel },$ and has bounded inverse.

Proof

The fact that $(0,0)\in A_{j,\parallel }$ is immediate and the boundedness of the inverse (cf. Definition A.8) follows from the inequality $\left| \sigma \right| \le \tau$. It remains to show maximal monotonicity. Let $(\sigma _{1},{\dot{\varepsilon }}_{1}),(\sigma _{2},\dot{\varepsilon }_{2})\in A_{j,\parallel }$. We consider the inner product of the differences $\left\langle \delta \sigma ,\delta \dot{\varepsilon }\right\rangle _{X_{j}}$ and distinguish three cases:

• 1. If $\tau {\dot{\varepsilon }}_{k} =\left| {\dot{\varepsilon }}_{k}\right| \sigma _{k}$ for both $k$, then

  $$\begin{aligned} \left\langle \delta \sigma ,\delta {\dot{\varepsilon }}\right\rangle _{X_{j}}&= \tau \left\langle \frac{{\dot{\varepsilon }}_{1}}{\left| {\dot{\varepsilon }}_{1}\right| }-\frac{{\dot{\varepsilon }}_{2}}{\left| {\dot{\varepsilon }}_{2}\right| },{\dot{\varepsilon }}_{1}-{\dot{\varepsilon }}_{2}\right\rangle _{X_{j}} \\&\ge \tau \int _{X_{j}}^{}(\left| \dot{\varepsilon }_{1}\right| +\left| {\dot{\varepsilon }}_{2}\right| -\left| {\dot{\varepsilon }}_{1}\right| -\left| {\dot{\varepsilon }}_{2}\right| ) = 0. \end{aligned}$$

• 2. If ${\dot{\varepsilon }}_{k} = 0$ for both $k$, then $\delta {\dot{\varepsilon }} = 0$.

• 3. If $\tau {\dot{\varepsilon }}_{1} =\left| {\dot{\varepsilon }}_{1}\right| \sigma _{1}$ and ${\dot{\varepsilon }}_{2} = 0$, then

  $$\begin{aligned} \left\langle \delta \sigma ,\delta \dot{\varepsilon }\right\rangle _{X_{j}} =\left\langle \tau \frac{\dot{\varepsilon }_{1}}{\left| {\dot{\varepsilon }}_{1}\right| }-\sigma _{2},{\dot{\varepsilon }}_{1}\right\rangle _{X_{j}}\ge \int _{X_{j}}^{}(\tau -\left| \sigma _{2}\right| )\left| {\dot{\varepsilon }}_{1}\right| \ge 0. \end{aligned}$$

Hence, $A_{j,\parallel }$ is monotone. Finally, we note that $I+A_{j,\parallel }$ is surjective by the following arguments:

• 1. If $\left| {\dot{\varepsilon }}\right| <\tau$, then $(\sigma ,{\dot{\varepsilon }})\in I+A_{j,\parallel }$ for $\sigma ={\dot{\varepsilon }}$.

• 2. If $\left| {\dot{\varepsilon }}\right| \ge \tau$, then $(\sigma ,{\dot{\varepsilon }})\in I+A_{j,\parallel }$ for $\sigma = \tau \frac{{\dot{\varepsilon }}}{\left| {\dot{\varepsilon }}\right| }$.

In turn, Theorem 1.6 from [37] ensures that $A_{j,\parallel }$ is maximal monotone. $\square$

For the (perpendicular) contact law, we use the relation for rough surfaces described in Example 4.2. Setting $C_{j}^{1} = 0$, $C_{j}^{3}\ge 0$, and $C_{j}^{4}\ge 1$, this law is given by

$$\begin{aligned} A_{j,\perp } :=\left\{ (\sigma ,{\dot{\varepsilon }})\in \iota _{j,\perp }\mathfrak {G\times }\iota _{j,\perp }{\mathfrak {G}}\mid \sigma = -C_{j}^{3}(-{\dot{\varepsilon }})_{+}^{C_{j}^{4}}\right\} \end{aligned}$$

(5.19)

Lemma 5.9

The contact law $A_{j,\perp }$ of (5.19) satisfies 5.6.

Proof

It is easy to verify that $(0,0)\in A_{j,\perp }$. Boundedness of the inverse relation follows from the fact that $\left| \sigma \right| \le \left| \dot{\varepsilon }\right| ^{C_{j}^{4}}$. Secondly, it is easy to see that the function $f(x) = -(-x)_{+}^{C_{j}^{4}}$ is monotone. Since $A_{j,\perp }$ is the inverse of the graph of $f$, it is a monotone relation as well. Finally, it is clear that no monotone extension of $A_{j,\perp }$ exists, and hence, it is maximal. $\square$

We finalize this subsection by stating the well-posedness, which is a direct result of Theorem 5.3.

Theorem 5.4

Problem (5.15), in which the relations are given by Example 5.1, (5.18) and (5.19), is well-posed in the space $L_{\nu }^{2}({\mathbb {R}},{\hat{U}})$.

5.6.2 Bounded, c-maximal monotone contact relations

On the other hand, our base model (5.10) contains assumption 5.2 which requires that the relation itself (instead of its inverse) is bounded. This is not the case for Tresca friction introduced in (5.18). Thus, in order to obtain a model that satisfies 5.2, we perform two regularizations. First, we introduce a maximal strain rate $c_{\infty }>0$ such that the regularized friction relation becomes

$$\begin{aligned} A_{j,\parallel }^{reg} :=\left\{ (\sigma ,\dot{\varepsilon })\in \iota _{j,\parallel }\mathfrak {G\times }\iota _{j,\parallel }{\mathfrak {G}}\left| \begin{array}{c} \tau {\dot{\varepsilon }} =\left| {\dot{\varepsilon }}\right| \sigma \text { or }\left| \sigma \right| {\dot{\varepsilon }} = c_{\infty }1_{\left\{ \left| \sigma \right| >\tau \right\} }\sigma \\ \left| {\dot{\varepsilon }}\right| \le c_{\infty } \\ \end{array}\right. \text {a.e. on } X_{j}\right\} \end{aligned}$$

(5.20)

with $1_{\left\{ \left| \sigma \right| >\tau \right\} }$ the indicator function of the set $\left\{ \left| \sigma \right| >\tau \right\}$. We emphasize that $c_{\infty }$ can be chosen sufficiently large to ensure that this bound is not reached in physical applications. Secondly, we add a constant $\check{c}$ to the regularized law in order to ensure $c$-maximal monotonicity.

Lemma 5.10

The friction law $(A_{j,\parallel }^{reg}+\check{c})$ satisfies 5.2 for any $\check{c}>0$. That is, it is bounded, $\check{c}$-maximal monotone, and $(0,0)\in A_{j,\parallel }^{reg}+\check{c}$.

Proof

The fact that $(0,0)\in A_{j,\parallel }^{reg}+\check{c}$ is clear. Moreover, the boundedness of the post-set is guaranteed by the inequality $\left| {\dot{\varepsilon }}\right| \le c_{\infty }.$ It remains to show that $A_{j,\parallel }^{reg}$ is maximal monotone. Let $(\sigma _{1},\dot{\varepsilon }_{1}),(\sigma _{2},{\dot{\varepsilon }}_{2})\in A_{j,\parallel }^{reg}$. We distinguish three cases:

• 1. If $\tau {\dot{\varepsilon }}_{k} =\left| {\dot{\varepsilon }}_{k}\right| \sigma _{k}$ for both $k$, then we use the same arguments as the first case in Lemma 5.8.

• 2. If $\left| \sigma _{k}\right| {\dot{\varepsilon }}_{k} = c_{\infty }1_{\left\{ \left| \sigma _{k}\right| >\tau \right\} }\sigma _{k}$ for both $k$, then

  $$\begin{aligned} \begin{aligned} \left\langle \delta \sigma ,\delta {\dot{\varepsilon }}\right\rangle _{X_{j}}&= c_{\infty }\left\langle 1_{\left\{ \left| \sigma _{1}\right|>\tau \right\} }\frac{\sigma _{1}}{\left| \sigma _{1}\right| }-1_{\left\{ \left| \sigma _{2}\right|>\tau \right\} }\frac{\sigma _{2}}{\left| \sigma _{2}\right| },\sigma _{1}-\sigma _{2}\right\rangle _{X_{j}} \\&\ge c_{\infty }\int _{X_{j}}^{}(1_{\left\{ \left| \sigma _{1}\right|>\tau \right\} }(\left| \sigma _{1}\right| -\left| \sigma _{2}\right| )+1_{\left\{ \left| \sigma _{2}\right| >\tau \right\} }(\left| \sigma _{2}\right| -\left| \sigma _{1}\right| ))\ge 0 \\ \end{aligned} \end{aligned}$$

• 3. If $\tau {\dot{\varepsilon }}_{1} =\left| {\dot{\varepsilon }}_{1}\right| \sigma _{1}$ and $\left| \sigma _{2}\right| {\dot{\varepsilon }}_{2} = c_{\infty }1_{\left\{ \left| \sigma _{2}\right| >\tau \right\} }\sigma _{2}$, then $\left| \sigma _{1}\right| = \tau$ and so we obtain

  $$\begin{aligned} \begin{aligned} \left\langle \delta \sigma ,\delta {\dot{\varepsilon }}\right\rangle _{X_{j}}&=\left\langle \sigma _{1}-\sigma _{2},\left| {\dot{\varepsilon }}_{1}\right| \frac{\sigma _{1}}{\left| \sigma _{1}\right| } -c_{\infty }1_{\left\{ \left| \sigma _{2}\right|>\tau \right\} }\frac{\sigma _{2}}{\left| \sigma _{2}\right| }\right\rangle _{X_{j}} \\&\ge \int _{X_{j}}^{}(\left| {\dot{\varepsilon }}_{1}\right| (\left| \sigma _{1}\right| -\left| \sigma _{2}\right| )+c_{\infty }1_{\left\{ \left| \sigma _{2}\right|>\tau \right\} }(\left| \sigma _{2}\right| -\left| \sigma _{1}\right| )) \\&= \int _{X_{j}}^{}(\left| {\dot{\varepsilon }}_{1}\right| (\tau -\left| \sigma _{2}\right| )+c_{\infty }1_{\left\{ \left| \sigma _{2}\right| >\tau \right\} }(\left| \sigma _{2}\right| -\tau )) \\ \end{aligned} \end{aligned}$$

  The case $\left| \sigma _{2}\right| \le \tau$ now follows directly. For the other case, we have

  $$\begin{aligned} \left| {\dot{\varepsilon }}_{1}\right| (\tau -\left| \sigma _{2}\right| )+c_{\infty }(\left| \sigma _{2}\right| -\tau ) =(c_{\infty }-\left| \dot{\varepsilon }_{1}\right| )(\left| \sigma _{2}\right| -\tau )\ge 0. \end{aligned}$$

Hence, $A_{j,\parallel }^{reg}$is monotone. Finally, maximality can now be shown by noting that $1+A_{j,\parallel }^{reg}$ is surjective and invoking Theorem 1.6 from [37]. $\square$

For the contact law, we encounter the same issue: The relation (5.19) is not a bounded relation because the post-set of $\sigma = 0$ contains all ${\dot{\varepsilon }}>0$. We therefore use the same regularization as in the friction law and define a $c_{\infty }>0$ that is set outside of the relevant limits of the physical model. This leads us to the following relation:

$$\begin{aligned} A_{j,\perp }^{reg} :=\left\{ (\sigma ,{\dot{\varepsilon }})\in \iota _{j,\perp }\mathfrak {G\times }\iota _{j,\perp }{\mathfrak {G}}\left| \begin{array}{c} (\sigma +C_{j}^{3}(-{\dot{\varepsilon }})_{+}^{C_{j}^{4}})(\varepsilon -c_{\infty }) = 0 \\ {\dot{\varepsilon }} \le c_{\infty } \\ \sigma \ge -C_{j}^{3}(-{\dot{\varepsilon }})_{+}^{C_{j}^{4}} \\ \end{array}\right. \text { a.e. on } X_{j}\right\} . \end{aligned}$$

(5.21)

Lemma 5.11

The contact law $(A_{j,\perp }^{reg}+{\hat{c}})$ satisfies 5.2 for any ${\hat{c}}>0$.

Proof

Boundedness and the fact that $(0,0)\in A_{j,\perp }^{reg}+\check{c}$ are easy to verify. We show monotonicity of $A_{j,\perp }^{reg}$ by considering three cases:

• 1. If ${\dot{\varepsilon }}_{k} = c_{\infty }$ for both $k = 1,2$ then we have $\delta {\dot{\varepsilon }} = 0$.

• 2. If $\sigma _{k} = -C_{j}^{3}(-{\dot{\varepsilon }}_{k})_{+}^{C_{j}^{4}}$ for both $k = 1,2$, then we use the fact that $f(x) = -(-x)_{+}^{c}$ is a monotone function and thus

  $$\begin{aligned} \left\langle \delta \sigma ,\delta \dot{\varepsilon }\right\rangle _{X_{j}} = C_{j}^{3}\left\langle -(-\dot{\varepsilon }_{1})_{+}^{C_{j}^{4}}+(-\dot{\varepsilon }_{2})_{+}^{C_{j}^{4}},{\dot{\varepsilon }}_{1}-\dot{\varepsilon }_{2}\right\rangle _{X_{j}}\ge 0. \end{aligned}$$

• 3. If $\sigma _{1} = -C_{j}^{3}(-{\dot{\varepsilon }}_{1})_{+}^{C_{j}^{4}}$ and ${\dot{\varepsilon }}_{2} = c_{\infty }$ then $\sigma _{1}\le 0$ and $\sigma _{2}\ge 0$ and so

  $$\begin{aligned} \left\langle \delta \sigma ,\delta \dot{\varepsilon }\right\rangle _{X_{j}} =\left\langle \sigma _{1}-\sigma _{2},{\dot{\varepsilon }}_{1}-c_{\infty }\right\rangle _{X_{j}}\ge 0. \end{aligned}$$

This shows that $A_{j,\perp }^{reg}$ is monotone. No monotone extension of this law exists and hence it is maximal. In turn, $(A_{j,\perp }^{reg}+{\hat{c}})$ is ${\hat{c}}$-maximal monotone. $\square$

Remark 5.5

Applying the bound $c_{\infty }$ is equivalent to applying a cut-off operator on the strain rate similar to [38, 39]. We note that, in practice, one can always check a posteriori whether the solution has attained the bound at any moment in time. If not, then the solution is independent of this bound, as desired. If the solution attains the bound, however, then the value of $c_{\infty }$ can be increased. If no bounded $c_{\infty }$ exists, then the physical model assumptions will at some point be violated, since this would indicate an arbitrarily large bulk velocity.

The well-posedness of the resulting problem is now a direct consequence of Theorem 5.2, which we present formally in the following theorem.

Theorem 5.5

Problem (5.10), in which the relations are given by Example 5.1 and those analyzed in Lemmas 5.10 and 5.11, is well-posed in the space $L_{\nu }^{2}({\mathbb {R}},U)$.

6 Summary and final remarks

In this manuscript, we have provided a general mixed-dimensional finite strain model for poromechanics in the presence of fracture (Sect. 4.4), its simplification in the case of linearized strain (Sect. 4.5), and a well-posedness theory allowing in the setting of linearized stain, still allowing for a generality in terms of the constitutive laws (Sect. 5.4). These developments are to the best of our knowledge all new. Our finite strain theory is rotationally invariant, and our mixed-dimensional model have several well-established models as its simplifications (Sect. 4.6).

As presented, the model allows for a range of physical phenomena, some of which may be desirable to neglect in various concrete applications. One such example arises in friction, where Tresca friction does not conform to a parameterization of the full model, but instead is in a different class of maximally monotone relations where there is no positive lower bound. This barely violates assumption 5.2, and is thus a degenerate limit of our model framework (using our general theory, as stated in Theorem 5.2, leads to a continuity constant that is unbounded). On the other hand, Tresca friction (and other models satisfying assumption 5.6) can nevertheless be allowed by a small perturbation of the spaces considered, as illustrated in Sect. 5.5.

An example of a different kind of degeneracy, which is not analyzed herein, is related to the presence of mechanical strain and stress terms at solid surfaces. This allows for modeling of surface effects, as may appear due to mineral processes in the subsurface, or through surface coating in industrial applications. On the other hand, such surface effects may be desirable to neglect for “clean" fractures. We consider this also as a degenerate limit of our model, but in this case it is assumption 5.3 that is violated (and possibly also 5.1).

Another important point which is not covered by our well-posedness analysis is the dependency of fracture permeability on aperture. In terms of the model structure, this enters in the sense of the permeability depending on the mixed-dimensional strain, as stated in (4.4). Such dependencies have recently been analyzed in the fixed-dimensional case [40], and we are optimistic that their approach can be extended to the mixed-dimensional setting.

We have previously shown how numerical methods can be derived for the simpler problem of flow in porous media (see references in Sect. 4.6.2). Development of numerical methods for the current problem is ongoing, and we look forward to reporting on this in future work.

Appendix A: Evolutionary equations and monotonicity

This work employs the theoretical framework from [10] in order to present the model and its analysis. This setting is more general than the conventional approach in which mappings in Sobolev spaces are used, and therewith provides us three key advantages.

• First, the domains of differential operators are derived from the operator, rather than vice versa. This relieves the need for characterizing the solution space and the theory is formulated in $L^{2}$-type spaces instead.

• Second, it allows us to include constitutive laws in our model that are not, strictly speaking, mappings between Sobolev spaces. In fact, we can consider the larger class of maximal monotone relations.

• Third, the theory of evolutionary equations naturally incorporates the time derivatives in a continuous setting. This provides an existence result for the entire time domain, rather than requiring arguments based on discrete time stepping.

In order to provide accessible reference, we recall the key concepts from this framework in this appendix. Let $H$ be a Hilbert space endowed with inner product $\left\langle \cdot , \cdot \right\rangle$.

Definition A.1

A binary relation between $H$ and $H$ is a subset $A\subseteq H\times H$.

This set-theoretical perspective on binary relations allows us to speak of closed relations (if the set $A$ is closed) and of the closure of a relation, which we denote by ${\overline{A}}$.

Binary relations have an algebraic structure in the sense that for $A,B\subseteq H\times H$ and $\lambda \in {\mathbb {R}}$, we have

$$\begin{aligned} A+\lambda B =\left\{ (x,y)\in H\times H\left| \begin{array}{c} \exists y_{A}\in H \text { with } (x,y_{A})\in A, \\ \exists y_{B}\in H \text { with } (x,y_{B})\in B, \\ \end{array} \right. \text { such that } y = y_{A}+\lambda y_{B}\right\} . \end{aligned}$$

(A1)

Example A.1

Given a mapping $f:H\rightarrow H$, then its graph $\textrm{Gr}(f)$ given by all pairs $(x,f(x))$ with $x\in H$ is a binary relation. The algebraic structure of binary relations naturally generalizes that of mappings in the sense that

$$\begin{aligned} \textrm{Gr}(f)+\lambda \textrm{Gr}(g) = \textrm{Gr}(f+\lambda g). \end{aligned}$$

We do not distinguish between a mapping and its graph and instead reuse the notation $f$, i.e., write $(x,y)\in f\subseteq H\times H$ when referring to the corresponding binary relation.

Example A.2

For a constant $c\in {\mathbb {R}}$, we have $(x,y)\in c$ if $y = cx$. The binary relation $c\subseteq H\times H$ therefore corresponds to the multiplication by $c$. Importantly, this multiplicative structure implies that $(x,\partial _t y) \in \partial _t c$.

Example A.3

We adopt a matrix/vector notation for binary relations between tuples of variables. Thus, let $H = H_{1}\times H_{2}$ and let $u :=\begin{bmatrix} u_{1} \\ u_{2} \end{bmatrix} \in H$ and $v :=\begin{bmatrix} v_{1} \\ v_{2} \end{bmatrix} \in H$. Moreover, let $A_{ij}\subseteq H_{j}\times H_{i}$ for $i,j\in \left\{ 1,2\right\}$. We denote

$$\begin{aligned} (u,v) =\left( \left[ \begin{array}{c} u_{1} \\ u_{2} \\ \end{array}\right] ,\left[ \begin{array}{c} v_{1} \\ v_{2} \\ \end{array}\right] \right) \in \begin{bmatrix} A_{11} &{}\quad A_{12} \\ A_{21} &{}\quad A_{22} \\ \end{bmatrix} \end{aligned}$$

(A2)

if and only if $v_{i} = \sum _{j}^{}v_{ij}$ with $(u_{j},v_{ij})\in A_{ij}$ for each $i$ and $j$.

Definition A.2

The domain and range of a binary relation $A\subseteq H\times H$ are denoted by $\textrm{dom}(A)$ and $\text {ran}(A)$, respectively, and are given by

$$\begin{aligned} \begin{aligned} \textrm{dom}(A)&:=\left\{ x\in H\mid \exists y\in H \text { s.t. }(x,y)\in A\right\} , \\ \text {ran}(A)&:=\left\{ y\in H\mid \exists x\in H \text { s.t. }(x,y)\in A\right\} . \\ \end{aligned} \end{aligned}$$

(A3)

As is common in a functional analysis setting, $\text {ran}(A)$ may be a proper subset of $H$, i.e., the range of an operator is allowed to be smaller than its codomain. In this setting, we typically also have that $\textrm{dom}(A)$ is a proper subset of $H$. Thus, it is important to remember that both the domain and range of $A$ can be proper subsets of $H$, despite $A$ being defined as $A\subseteq H\times H$.

Definition A.3

A binary relation $A\subseteq H\times H$ is monotone if for all $(x_{1},y_{1}),(x_{2},y_{2})\in A$ it holds that

$$\begin{aligned} \text {Re}(\left\langle x_{1}-x_{2},y_{1}-y_{2}\right\rangle )\ge 0. \end{aligned}$$

Definition A.4

A binary relation $A\subseteq H\times H$ is maximal monotone if it is monotone and if for all $B\subseteq H\times H$ with $B$ monotone and $A\subseteq B$, it follows that $A = B$.

Definition A.5

A binary relation $A\subseteq H\times H$ is $c$-maximal monotone for some $c>0$ if $A-c$ is maximal monotone.

The following example illustrates the generality that these definitions allow for.

Example A.4

Let $H\mathbb { :=R}$ and let the binary relation $A\mathbb {\subset R\times R}$ be given by

$$\begin{aligned} A :=\left\{ (x,y)\in \mathbb {R\times R}\left| \begin{array}{c} y = 0 \text { if } \left| x\right| <1 \\ xy\ge 0 \text { if } \left| x\right| = 1 \\ \end{array}\right. \right\} . \end{aligned}$$

Graphically, $A$ is given by the horizontal line segment $(x,0)$ for $-1\le x\le 1$ and the vertical half-lines $(\text {sgn}(y), y)$. In turn, $\textrm{dom}(A) =\left[ -1,1\right]$, which we emphasize is a proper subset of ${\mathbb {R}}$, and $\text {ran}(A) = {\mathbb {R}}$. The fact that $A$ is monotone can be verified by a straightforward computation:

• $\left| x_{1}\right| ,\left| x_{2}\right| <1.$ Then $y_{1} = y_{2} = 0$ and so $\left\langle x_{1}-x_{2},y_{1}-y_{2}\right\rangle = 0.$

• $\left| x_{1}\right| <1 =\left| x_{2}\right| .$Then $y_{1} = 0$ and $\left\langle x_{1}-x_{2},y_{1}-y_{2}\right\rangle =(x_{2}-x_{1})y_{2}\ge 0$ since the two terms have the same sign.

• $\left| x_{1}\right| =\left| x_{2}\right| = 1$. Then, $(x_{1}-x_{2})$ is either zero or has the same sign as $(y_{1}-y_{2})$, giving the result.

Moreover, $A$ is maximal monotone since there exists no monotone extension $B\mathbb {\subseteq R\times R}$ with $A\subseteq B$ except for $B = A$.

Definition A.6

The adjoint of a binary relation $A\subseteq H_{1}\times H_{2}$ is given by

$$\begin{aligned} A' :=\left\{ (u,v)\in H_{2}\times H_{1}\mid \left\langle u,y\right\rangle _{H_{2}} =\left\langle v,x\right\rangle _{H_{1}}, \forall (x,y)\in A\right\} . \end{aligned}$$

Definition A.7

The inverse of a binary relation $A\subseteq H\times H$ is given by

$$\begin{aligned} A^{-1} :=\left\{ (x,y)\in H\times H\mid (y,x)\in A\right\} . \end{aligned}$$

Importantly, this definition implies that any binary relation has an inverse relation. Moreover, it is easy to see that the inverse of a maximal monotone relation is itself maximal monotone as well (by the symmetry of the inner product).

Example A.5

Continuing with example A.4, we have

$$\begin{aligned} A^{-1} =\left\{ (x,y)\in \mathbb {R\times R}\left| \begin{array}{c} y = \text {sgn}(x) \quad \text { if } x\ne 0 \\ y\in \left[ -1,1\right] \quad \text { if } x = 0 \\ \end{array}\right. \right\} . \end{aligned}$$

As is apparent here, $A^{-1}$ can be described as a set-valued function. However, for the reasons mentioned at the beginning of the section, we herein prefer the use of binary relations.

Even though the inverse of a relation is always well-defined, we often require more properties of the inverse, in particular that it corresponds to a Lipschitz continuous mapping. This can be obtained by Minty’s theorem, which lies at the heart of the analysis of evolutionary equations.

Theorem A.1

(Minty’s theorem) Let $A\subseteq H\times H$ be a $c$-maximal monotone relation for some $c>0$. Then the inverse relation $A^{-1}$ defines a Lipschitz continuous mapping with domain $\textrm{dom}(A^{-1}) = H$ and a Lipschitz constant bounded by $\frac{1}{c}$.

Proof

See Theorem 1.1 of [10]. $\square$

Example A.6

Continuing with example A.5, let $B = A+\epsilon$ for some $\epsilon >0.$ We then compute

$$\begin{aligned} B^{-1} =\left\{ (x,y)\in \mathbb {R\times R}\left| \begin{array}{c} y = \text {sgn}(x) \quad \text { if }\left| x\right| \ge \epsilon \\ y =\frac{1}{\epsilon }x \quad \text { if }\left| x\right| <\epsilon \\ \end{array}\right. \right\} . \end{aligned}$$

Note that, by construction, $B-c = A$ is maximal monotone with $c = \epsilon >0$. In turn, Minty’s theorem implies that $B^{-1}$ is Lipschitz continuous with possible Lipschitz constant $\frac{1}{\epsilon }$ and this can easily be verified. On the other hand, we see that $A^{-1}$ is not Lipschitz continuous at the origin. In turn, the converse of Minty’s theorem implies that there is no $c>0$ such that $A$ is $c$-maximal monotone.

Minty’s theorem forms a powerful tool and it is therefore important to determine whether a binary relation is maximal monotone. It becomes a natural question to ask whether the sum of two maximal monotone relations is maximal monotone as well. To answer this, we first define bounded relations.

Definition A.8

A relation $A\subseteq H\times H$ is called bounded if for all bounded sets $M\subseteq H$, the post-set given by $\left\{ y\in H\mid \exists x\in M \text { s.t. } (x,y)\in A\right\}$ is bounded as well.

We remark that if $A$ corresponds to a mapping, the post-set is typically referred to as the image of $M$ under $A$. Using the boundedness property, we have the following result concerning the sum of maximal monotone relations.

Lemma 7.1

Let $A,B\subset H\times H$ be maximal monotone, let $A$ be bounded and let $\textrm{dom}(A)\cap \textrm{dom}B)\ne \varnothing$. Then, the sum $A+B$ is maximal monotone.

Proof

See Proposition 1.22 of [37]. $\square$

On the other hand, we will frequently encounter unbounded operators as well. An important class of these is given by densely defined, unbounded linear operators.

Definition A.9

An operator $A\subseteq H\times H$ is densely defined if $\textrm{dom}(A)$ is dense in $H$.

Differential operators are a primary example of unbounded operators in our work, defined as mappings between $L^{2}$-type spaces. We emphasize this using the following example

Example A.7

Let $\Omega \subset {\mathbb {R}}^{n}$ be a Lipschitz domain with $n\in \left\{ 2,3\right\}$. Let ${\tilde{\nabla }}_{0}$ be the gradient operator defined on $\mathring{C}^{\infty }(\Omega )$, i.e., on the space of infinitely differentiable functions with vanishing trace on $\partial \Omega$:

$$\begin{aligned} {\tilde{\nabla }}_{0}:\mathring{C}^{\infty }(\Omega )\subset L^{2}(\Omega ,{\mathbb {R}})\rightarrow L^{2}(\Omega ,{\mathbb {R}}^{n}). \end{aligned}$$

Now, let the divergence $(\nabla \cdot )$ and its dual be given by

$$\begin{aligned} (\nabla \cdot ) :=-({\tilde{\nabla }}_{0})' \text { and } \nabla _{0} :=-(\nabla \cdot )'. \end{aligned}$$

Let $H$ be the pair of spaces

$$\begin{aligned} H :=L^{2}(\Omega , {\mathbb {R}})\times L^{2}(\Omega ,{\mathbb {R}}^{n}). \end{aligned}$$

Finally, let the relation $A\subseteq H\times H$ be given by

$$\begin{aligned} A :=\begin{bmatrix} 0 &{} \quad \nabla \cdot \\ \nabla _{0} &{}\quad 0 \\ \end{bmatrix}. \end{aligned}$$

We emphasize that the domain of the operator $A$ follows from its definition. In fact, $\textrm{dom}(A)$ is given by all functions $(f, g)\in H$ with $f\in \textrm{dom}(\nabla _{0})\subseteq L^{2}(\Omega , {\mathbb {R}})$ and $g\in \textrm{dom}(\nabla \cdot )\subseteq L^{2}(\Omega ,{\mathbb {R}}^{n})$. For these functions, it follows by definition that $\nabla _{0}f\in L^{2}({\mathbb {R}}^{n})$ and $\nabla \cdot g\in L^{2}(\Omega , {\mathbb {R}})$. In turn, we have the following characterization in more conventional notation (see Fig. 8)

$$\begin{aligned} \textrm{dom}(A) = \textrm{dom}(\nabla _{0})\times \textrm{dom}(\nabla \cdot ) = \mathring{H}^{1}(\Omega )\times H(\nabla \cdot ,\Omega ). \end{aligned}$$

Finally, we note that $\textrm{dom}(A)$ is dense in $H$ and thus $A$ is a densely defined, unbounded linear operator.

Remark A.1

Analogously, we may start with $({\tilde{\nabla }}_{0} \cdot )$ as the divergence acting on vector-valued, infinitely differentiable functions with vanishing normal trace. By taking the appropriate adjoints, this leads us to the operators $\nabla$ and $(\nabla _{0} \cdot )$ with

$$\begin{aligned} \textrm{dom}(\nabla )&= H^{1}(\Omega ),&\textrm{dom}(\nabla _{0} \cdot )&= \mathring{H}(\nabla \cdot ,\Omega ). \end{aligned}$$

(A4)

Here, $\mathring{H}(\nabla \cdot ,\Omega )$ denotes the subspace of $H(\nabla \cdot ,\Omega )$ consisting of functions with zero normal trace on $\partial \Omega$.

Fig. 8

The gradient and divergence mappings defined as densely defined, unbounded linear operators. Note that both the domain and range of $\nabla _{0}$ are a proper subset of $L^{2}$. The fact that the range of $\nabla \cdot$ equals $L^{2}$ is classical

Example A.8

Continuing with example A.7, we note that $A$ is linear and skew self-adjoint, and therefore maximal monotone. In particular, for $[f_{1}, g_{1}]^T, [f_{2}, g_{2}]^T \in \textrm{dom}(A)$, we have

$$\begin{aligned} \begin{aligned} \left\langle \left[ \begin{array}{c} f_{1}-f_{2} \\ g_{1}-g_{2} \\ \end{array}\right] ,A\left[ \begin{array}{c} f_{1} \\ g_{1} \\ \end{array}\right] -A\left[ \begin{array}{c} f_{1} \\ g_{1} \\ \end{array}\right] \right\rangle&=\left\langle f_{1}-f_{2},\nabla \cdot (g_{1}-g_{2})\right\rangle +\left\langle \nabla _{0}(f_{1}-f_{2}),g_{1}-g_{2}\right\rangle \\&=\left\langle f_{1}-f_{2},\nabla \cdot (g_{1}-g_{2})\right\rangle -\left\langle f_{1}-f_{2},\nabla \cdot (g_{1}-g_{2})\right\rangle \\&= 0. \\ \end{aligned} \end{aligned}$$

(A5)

Here, we have used angled brackets for the inner products of both the product space $H$ and its components.

In order to define an evolutionary equation, we next consider the spatiotemporal setting. Following [10], we introduce an exponentially weighted Bochner-type function space.

Definition A.10

For $\nu > 0$, let

$$\begin{aligned} L_{\nu }^{2}({\mathbb {R}},H) :=\left\{ f\mathbb {:R\rightarrow }H\mid \int _{{\mathbb {R}}}^{}\left| e^{-\nu t}f(t)\right| _{H}^{2}\textrm{d}t<\infty \right\} \end{aligned}$$

Definition A.11

Given a binary relation $A\subseteq H\times H$, its temporal extension $A_{\nu }\subseteq L_{\nu }^{2}({\mathbb {R}},H)\times L_{\nu }^{2}({\mathbb {R}},H)$ is given by

$$\begin{aligned} A_{\nu } :=\left\{ u,v\in L_{\nu }^{2}({\mathbb {R}},H)\mid (u(t),v(t))\in A, \text {for almost every } t\in {\mathbb {R}}\right\} \end{aligned}$$

Definition A.12

An operator $A:\textrm{dom}(A)\subseteq L_{\nu }^{2}({\mathbb {R}},H)\rightarrow L_{\nu }^{2}({\mathbb {R}},H)$ is time-translation invariant if for each $(u,v)\in A$ and $h\in {\mathbb {R}}$, we have $(u( \cdot +h),v( \cdot +h))\in A.$

Definition A.13

A closed mapping $A:\textrm{dom}(A)\subseteq L_{\nu }^{2}({\mathbb {R}},H)\rightarrow L_{\nu }^{2}({\mathbb {R}},H)$ is called causal if for all $h\in {\mathbb {R}}$ and $f,g\in \textrm{dom}(A)$ with $f\left. \right| _{t\le h} = g\left. \right| _{t\le h}$, it follows that $(Af)\left. \right| _{t\le h} =(Ag)\left. \right| _{t\le h}$.

Finally, we present the time derivative on the weighted space $L_{\nu }^{2}({\mathbb {R}},H)$. Its definition is motivated by the following short calculation for differentiable $f\mathbb {:R\rightarrow R}$.

$$\begin{aligned} e^{\nu t}\partial _{t}(e^{-\nu t}f) = e^{\nu t}(-\nu e^{-\nu t}f+e^{-\nu t}\partial _{t}f) =(-v+\partial _{t})f. \end{aligned}$$

In turn, we have $e^{\nu t}(\partial _{t}+\nu )e^{-\nu t}f = \partial _{t}f$, which suggests the following definition.

Definition A.14

Given $\nu >0$, let $\partial _{0,\nu }:\textrm{dom}(\partial _{0,\nu })\subseteq L_{\nu }^{2}({\mathbb {R}},H)\rightarrow L_{\nu }^{2}({\mathbb {R}},H)$ be given by

$$\begin{aligned} \partial _{0,\nu } :=e^{\nu t}(\partial _{t}+\nu )e^{-\nu t}. \end{aligned}$$