Update aem_model.md

docs/theory/aemmodel/aem_model.md

@@ -12,51 +12,51 @@ In order to explain the AEM approach and in order compare/benchmark the 3 dimens
 
 ## Two-Dimensional Solution
 
-The steady state pressure conditions at distance r from the well radius is given by Thiem’s solution for an infinite aquifer:
-
-$$p(r) = \frac{\mu q}{2 \pi k H} ( \ln\left(\frac{r}{r_w}\right) + S)$$
-
-The pressure differences can be superposed for N wells. When choosing a very large reference radius , we can write the superposing effects of pressure for a well i with additional surrounding wells j as:
-
-$$P_i = \frac{\mu}{2 \pi k H} \left( \ln\left( \frac{r_e}{r_{w_i}} \right) q_i + \sum_{\substack{j=1 \\ j \ne i}}^{N} \ln\left( \frac{r_e}{L_j} \right) q_j \right)$$
-
-$r_wi$ is the radius of the well completion, and  $L_j$ is distance of well j to the well i. As the sum of flow rates is zero, this can be simplified to:
-
-
+The steady state pressure conditions at distance r from the well radius is given by Thiem’s solution for an infinite aquifer. In the well, a correction for the skin is applied:
 
 $$
-P_i = \frac{\mu}{2 \pi k H} \left( \sum_{j=1}^{N} \ln(r_e) q_j + \ln\left(\frac{1}{r_{w_i}}\right) q_i + \sum_{\substack{j=1 \\\\ j \ne i}}^{N} \ln\left(\frac{1}{L_j}\right) q_j \right)
+p(r) = 
+\left\\{ 
+\begin{matrix}
+\frac{\mu q}{2 \pi k H} \ln\frac{r_e}{r}                     & \text{(in reservoir)}\\
+\frac{\mu q}{2 \pi k H} \left( \ln\frac{r_e}{r_w} + S\right) & \text{(in well)}
+\end{matrix}
+\right.
 $$
 
-=
-
+The pressure differences can be superposed for N wells. When choosing a very large reference radius, we can write the superposing effects of pressure for a well i with additional surrounding wells j as:
 
 $$
-P_i = \frac{\mu}{2 \pi k H} \left( \ln\left(\frac{1}{r_{w_i}}\right) q_i + \sum_{\substack{j=1 \\ j \ne i}}^{N} \ln\left(\frac{1}{L_j}\right) q_j \right)
+P_i = \frac{\mu}{2 \pi k H} \left( q_i\left[\ln\frac{r_e}{r_{wi}} + S_i\right]  + \sum_{\substack{j=1 \\ j \ne i}}^{N} q_j\ln\frac{r_e}{L_j}  \right)
 $$
 
+$r_{wi}$ is the radius of the well completion, and  $L_j$ is distance of well j to the well i. As the sum of flow rates is zero, this can be simplified to:
 
 
-Consequently, we can write for the system of equation to solve:
+
+$$
+P_i = \frac{\mu}{2 \pi k H} \left( \sum_{j=1}^{N} q_j\ln r_e  + q_i\left[\ln\frac{1}{r_{wi}} + S_i\right]  + \sum_{\substack{j=1 \\ j \ne i}}^{N} q_j\ln\frac{1}{L_j}  \right)=
+\frac{\mu}{2 \pi k H} \left( q_i\left[\ln\frac{1}{r_{wi}} + S_i\right] + \sum_{\substack{j=1 \\ j \ne i}}^{N} q_j\ln\frac{1}{L_j}  \right)
+$$
 
 
+When the injectors have a pressure $P_i=P_{inj}$ and the producers are known to have the same (but yet unknown) pressure $P_P$, we have the following system of equations to solve:
 
 $$
 \begin{bmatrix}
-I_{i=1}(r_{w_i}) & I_{i=1}(L_1) & \cdots & I_{i=1}(L_N) & D_1 \\
-I_{i=2}(L_1) & I_{i=2}(r_{w_i}) & \cdots & I_{i=2}(L_N) & D_2 \\
+I_{11} & I_{12} & \cdots & I_{1N} & D_1    \\
+I_{22} & I_{22} & \cdots & I_{2N} & D_2    \\
 \vdots & \vdots & \ddots & \vdots & \vdots \\
-I_N(L_1) & I_N(r_{w_N}) & \cdots & I_N(L_N) & D_N \\
-1 & 1 & \cdots & 1 & 0
+I_{N1} & I_{N2} & \cdots & I_{NN} & D_N    \\
+1      & 1      & \cdots & 1      & 0
 \end{bmatrix}
 \begin{bmatrix}
 q_1 \\
 q_2 \\
 \vdots \\
 q_N \\
-P_p
-\end{bmatrix}
-=
+P_P
+\end{bmatrix}=
 \begin{bmatrix}
 P_1 \\
 P_2 \\
@@ -66,16 +66,21 @@ P_N \\
 \end{bmatrix}
 $$
 
-	(eq. 4)
+(eq. 4)
 
 with
 
 
 $$
-I_i(r) = A \mu_i \left( \ln\left( \frac{1}{r} \right) + S_i \right), \quad A = \frac{1}{2 \pi k H}
+I_{ij} = \left\\{
+\begin{matrix}
+\frac{\mu_i}{2 \pi k H} \left[ \ln\frac{1}{r_{wi}} + S_i \right] & (i=j) \\
+\frac{\mu_i}{2 \pi k H}  \ln\frac{1}{r_{ij}} & (i\ne j)
+\end{matrix}
+\right.
 $$
 
-where mu_i is viscosity at well i and is the well radius or the distance between the wells i and j
+where $\mu_i$ is viscosity at well $i$; $r_{wi}$ is the well radius and $r_{ij}$ the distance between the wells i and j
 
 $D_i$: 0 when the well is injector, -1 when producer
 
@@ -88,46 +93,44 @@ dependent of the sign of .
 
 $P_p$: resulting production pressure to maintain the steady state mass balance.
 
-For a doublet configuration, the system of equations reduces to the well-known steady state solution (e.g. Van Wees et al., 2012):
+For a doublet configuration with equal viscosities, equal wellbore radius and equal skin, the system of equations reduces to the well-known steady state solution (e.g. Van Wees et al., 2012):
 
-$$\Delta p = \frac{\mu q}{ \pi k H} ( \ln\left( \frac{L}{r_w} \right) + S)$$
+$$\Delta p = \frac{\mu q}{ \pi k H} \left[ \ln\frac{L}{r_w} + S\right]$$
 
 Please note that the viscosity is constant for each row of the equations(cf. eq. 3). In order to include the effects of the viscosity contrasts, the viscosity for each row is set to the viscosity around the well under consideration.
 
 ## Three-Dimensional Solution
 
-In three dimensions the system of equations is setup for N well line segments instead of N wells. 
+In three dimensions the system of equations is set up for $N$ well line segments instead of $N$ wells. 
 The structure of the system of equations remains the same, 
 and the system of equations and solution approach (excluding the last row and column) 
 is explained in detail in the annex of Egberts et al. (2013). 
-The M submatrix is a summation over line segment contribution of injectors and producers (k = 0) 
-and and multiple images below the base and above the top of the reservoir, marked by the indices k = -1 and k = 1. 
-These image contributions enforce no flow boundary conditions at the top and base of the reservoir. 
-The i and j denote pairs of control point location for each of the line segments (chosen at the middle) and is at least $r_w$ away from the well bore segment i.
+The $M$ submatrix is a summation over line segment contribution of injectors and producers 
+and multiple images below the base and above the top of the reservoir, marked by the index $k$. 
+These image contributions enforce no-flow boundary conditions at the top and base of the reservoir. 
+The $i$ and $j$ denote pairs of control point location for each of the line segments (chosen at the middle) and is at least $r_w$ away from the well bore segment $i$.
+The expression for the pressure field resulting from a finite well segment is different from the 2D expression for fully penetrating wells, and the skin of a well is accounted for by implementing an effective wellbore radius $r_w'=r_we^{-S}$. The expressions read:
 
-$$M_{ij} = \sum_{k=-1}^{1} \frac{\mu_i}{2 \pi k H} (\ln\left( \frac{1}{r_{ijk}} +S_i) \right)$$
+$$M_{ij} = \frac{\mu_i}{2 \pi k H} \sum_{k=-N_{im}}^{N_{im}} \frac{1}{L_i}\ln\left( \frac{d_{ijk}+e_{ijk}+2L_i}{d_{ijk}+e_{ijk}-2L_i}\right)$$
 
-Such that the full system of equations becomes:
+$L_i$ is the length of the well segment $i$; $d_{ijk}$ and $e_{ijk}$ are the distance of the point of evaluation (at the effective well radius of the segment source, or at one of the other well segments) to the centre of the source or one of its images. The full system of equations becomes:
 
 
 $$
 \begin{bmatrix}
-M_{ij} & \cdots & D_1 \\\\
-\vdots & \ddots & \vdots \\\\
-M_{ij} & \cdots & D_N \\\\
-1 & \cdots & 0
+M_{11} & \cdots & M_{1N} & D_1 \\\\
+\vdots & \ddots & \vdots & \vdots \\\\
+M_{N1} & \cdots & M_{NN} & D_N \\\\
+1 & \cdots & 1 & 0
 \end{bmatrix}
 \begin{bmatrix}
 q_1 \\\\
-q_2 \\\\
 \vdots \\\\
 q_N \\\\
 P_p
-\end{bmatrix}
-=
+\end{bmatrix}=
 \begin{bmatrix}
 P_1 \\\\
-P_2 \\\\
 \vdots \\\\
 P_N \\\\
 0
@@ -136,9 +139,8 @@ $$
 
 
 The method of Egberts et al. (2013) allows to incorporate an anisotropic permeability field.
-In addition additional images are included to enforce spatial constrains in flow (i.e. no flow barriers). 
-The AEM also allows to be extended to incorporate frictional losses in the right hand side based on the solution for .
-This requires an iterative solution and needs to model the exact well segment pressure. Please note that the solution is in 3D is based a single reference  well pressure  and  relative to the ambient pressure, and neglects potential vertical variation in pressure in the well bore compared to the ambient reservoir pressure. This is considered a valid  assumption for geothermal reservoirs.
+It also allows to be extended to incorporate frictional losses in the solution; this has not been implemented.
+The current solution in 3D is based a single reference well pressure relative to the ambient pressure, and neglects potential vertical variation in pressure in the wellbore compared to the ambient reservoir pressure. This is considered a valid  assumption for geothermal reservoirs.
 
 ## Benchmark